sci.math.num-analysis SUNA45, The condition of 2-D diffusion ====================================== In SUNA42 for the integrand belonging to one arbitrary F.E. diffusion triangle the following discretization was found: (df/dx)^2 + (df/dy)^2 ] . dx.dy = 1/2.Det * 1/Det^2 * (y23.f1 + y31.f2 + y12.f3)^2 + (x23.f1 + x31.f2 + x12.f3)^2 Where Det = x21.y31 - x31.y21 = 2 * area triangle . Therefore the variational integral for one triangle is in fact equivalent with a Least Squares Finite Element Method for two equations with three unknowns. The accompanying finite element matrix, therefore, is expected to be singular: | y32.y32 + x32.x32 y32.y13 + x32.x13 y32.y21 + x32.x21 | | y13.y13 + x13.x13 y13.y21 + x13.x21 | / Det | symmetrical y21.y21 + x21.x21 | We can check this with help of MAPLE. Type " maple < problem.in " (UNIX), where "problem.in" is quoted without permission (-: > with(linalg); > y32 := y3 - y2; > y21 := y2 - y1; > y13 := y1 - y3; > x32 := x3 - x2; > x21 := x2 - x1; > x13 := x1 - x3; > E11 := y32*y32 + x32*x32; > E22 := y13*y13 + x13*x13; > E33 := y21*y21 + x21*x21; > E12 := y32*y13 + x32*x13; > E13 := y32*y21 + x32*x21; > E23 := y13*y21 + x13*x21; > E := array([[ E11 , E12 , E13 ] , > [ E12 , E22 , E23 ] , > [ E13 , E23 , E33 ]]); > D := det(E); # The answer of MAPLE will be: D := 0 > quit; Eigenvalues of the F.E. matrix may be determined according to: Det ( | y23.y23 + x23.x23 - L y23.y31 + x23.x31 y23.y12 + x23.x12 | | y31.y31 + x31.x31 - L y31.y12 + x31.x12 | ) = 0 | symmetrical y21.y21 + x21.x21 - L | Or: ( y32.y32 + x32.x32 - L ).( y13.y13 + x13.x13 - L ).( y21.y21 + x21.x21 - L ) + 2 * ( y32.y13 + x32.x13 ).( y13.y21 + x13.x21 ).( y32.y21 + x32.x21 ) - ( y32.y32 + x32.x32 - L ).( y13.y21 + x13.x21 )^2 - ( y13.y13 + x13.x13 - L ).( y32.y21 + x32.x21 )^2 - ( y21.y21 + x21.x21 - L ).( y32.y13 + x32.x13 )^2 = 0 Collect terms belonging to L^power : L^3 : - 1 L^2 : (y32.y32 + x32.x32) + (y13.y13 + x13.x13) + (y21.y21 + x21.x21) Define a , b , c as lengths of the three triangle sides, then for the case L^2 : a^2 + b^2 + c^2 . L^1 : - (y32.y32 + x32.x32).(y13.y13 + x13.x13) - ... + (y32.y13 + x32.x13)^2 + ... --------------------------------------- - y32.y32.y13.y13 - y32.y32.x13.x13 - x32.x32.y13.y13 - x32.x32.x13.x13 + y32.y13.y32.y13 + 2 * y32.y13.x32.x13 + x32.x13.x32.x13 ----------------------------------------------------------------------- - (x13.y32 - x32.y13)^2 = - Det^2 Same for the other three terms, giving a total for L^1 : - 3 . Det^2 | x1 y1 1 | Remark. The following holds too: Det = Det ( | x2 y2 1 | ) . | x3 y3 1 | L^0 : equals the determinant of the element-matrix itself, which was proven to be zero: 0 . The eigenvalues are thus found by solving the algabraic equation: - L^3 + (a^2 + b^2 + c^2).L^2 - 3.Det^2.L + 0 = 0 Therefore one eigenvalue is L1 = 0 and the other two are found by solving: L^2 - (a^2 + b^2 + c^2).L + 3.Det^2 = 0 Make dimensionless by calculating L/Det instead of L . And define: K = (a^2 + b^2 + c^2) / Det . Then: L^2 - K.L + 3 = 0 --> L = K/2 +/- sqrt[ (K/2)^2 - 3 ] The eigenvalues must be all (real valued) positive or zero, due to properties of positive (semi)definite symmetric matrices. Consequently, the minimum value of L is found for: K = 2.sqrt(3) ; In general: K >= 2.sqrt(3) . L2 = L3 = K/2 . -------------- It remains to be seen for which kind of triangles the above minimum is reached. /|\ The quantities K and L/Det are invariant for scaling b / | \ a (and for translation and for rotation). Therefore adopt a / |h \ triangle with side-length c = 1 , and height h . Then: / | \ -------|------- a^2 + b^2 + c^2 = x^2 + h^2 + (1 - x)^2 + h^2 + 1 1-x c x = 2.(x^2 + h^2 - x + 1) Det = h . Thus: K/2 = (x^2 + h^2 - x + 1)/h An extremum is reached for dK/dx = 0 and dK/dh = 0 (where "d" = partial). dK/dx = 0 --> 2.x - 1 = 0 --> x = 1/2 : a and b equal. dK/dh = 0 --> 2.h.h - (x^2 + h^2 - x + 1) = 0 --> h^2 - x^2 + x - 1 = 0 --> h^2 = x^2 - x + 1 = 1/4 - 1/2 + 1 = 3/4 --> h = sqrt(3)/2 --> a = b = c = 1 : equilateral triangle ! ------------- --> K = 3 / ( sqrt(3)/2 ) = 2.sqrt(3) , as was found before. This completes the proof that two _equal_ eigenvalues <> 0 are found if the FE matrix belongs to an equilateral triangle. Left as an exercise is another proof of the above, which uses the well known formula: Det/2 = area = sqrt[s.(s-a).(s-b).(s-c)] where s = (a+b+c)/2 . OK. So far so good for the eigenvalues of _one_ element matrix. What happens if these element matrices are assembled into a global matrix, for a whole mesh of diffusion triangles? It may be suspected, but not proved, that the _condition_ number (see "CNN.sum") of such a matrix is _best_ for a mesh consisting of all kind of equilateral triangles. Because the condition number is the largest divided by the smallest eigenvalue for positive (semi) definite symmetric matrices. And for a single element this C.N. is minimal for an equilateral triangle. Therefore we think ... - * Han de Bruijn; Applications&Graphics | "A little bit of Physics * No * TUD Computing Centre; P.O. Box 354 | would be NO idleness in * Oil * 2600 AJ Delft; The Netherlands. | Mathematics" (HdB). * for * E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood