sci.math.numanalysis
SUNA45, The condition of 2D 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 elementmatrix 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 sidelength c = 1 , and height h . Then:
/  \
 a^2 + b^2 + c^2 = x^2 + h^2 + (1  x)^2 + h^2 + 1
1x 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.(sa).(sb).(sc)] 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
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood