sci.math.numanalysis
SUNA52, Fast 1D FEM solving
============================
Disclaimer: this article is just for fun.
It has no theoretical nor practical significance whatsoever.
Consider the (onedimensional) ordinary second order differential equation with
constant coefficients:
a.x'' + b.x' + c.x = d ; x(0) = A ; x(L) = B
The equation will be solved numerically by employing the Finite Element Method.
Due to the constant coefficients, a uniform mesh with identical elements may be
employed. Such gives rise to one and only one Element matrix and vector:
 e11 e12   x1  =  r1 
 e21 e22   x2   r2 
Assembly is performed for two identical elements:
 e11 e12 0   x1   r1 
 e21 e22+e11 e12   x2  =  r1+r2 
 0 e21 e22   x3   r2 
Solve for x2 : x2 = ( r1+r2  e21.x1  e12.x3 ) / ( e22+e11 )
Substitute in: e11.x1 + e12.x2 =
e11.x1 + e12.( r1+r2  e21.x1  e12.x3 ) / ( e22+e11 ) = r1
Giving:
[ e11  e12.e21/(e22+e11) ].x1 + [  e12.e12/(e22+e11) ].x3 =
= r1  e12.(r1+r2)/(e22+e11)
Substitute in: e21.x2 + e22.x3 =
e21.( r1+r2  e21.x1  e12.x3 ) / ( e22+e11 ) + e22.x3 = r2
Giving:
[  e21.e21/(e22+e11) ].x1 + [ e22  e21.e12/(e22+e11) ].x3 =
= r2  e21.(r1+r2)/(e22+e11)
A fast algorithm is found: Symmetric case e21 = e12 :
^^^^
tra := 1/(e11+e22) tra := 1/(e11+e22)
r1 := r1  e12.(r1+r2).tra add :=  e12.(r1+r2).tra
r2 := r2  e21.(r1+r2).tra r1 := r1 + add
sto :=  e21.e12.tra r2 := r2 + add
e11 := e11 + sto sto :=  e12.e12.tra
e22 := e22 + sto e11 := e11 + sto
e12 :=  e12.e12.tra e22 := e22 + sto
e21 :=  e21.e21.tra e12 := sto
The above calculations are carried out repeatdly. All results e(i,j) and b(i)
must be saved for later usage. At each step, the element is enlarged twice.
For a 8 = 2^3 elements mesh, the process is visualized (ASCII) as follows:
o  o  o  o  o  o  o  o  o : step 0
o  o  o  o  o : step 1
o  o  o : step 2
o  o : step 3
At step 3, the resulting equations shouls be solved, at hand of the boundary
conditions present. After that, a backsubstitution process will evaluate the
inbetween values. Here the results e(i,j) and b(i) previously stored must
be used again, together with:
x2 = ( r1+r2  e21.x1  e12.x3 ) / ( e22+e11 ) : see above.
The algorithm has time performance log(N) , if N is the number of elements.
That's what we mean by "fast".
As an example, the diffusion problem (d/dx)^2.T = 0 may be discretized at an
uniform mesh by the following elementmatrix:
 +1 1 
 1 +1 
Hence the algorithm for the symmetric case e21 = e12 :
tra := 1/2
sto :=  1/2 ==> 1/2 .  +1 1 
e11 := 1  1/2 = 1/2  1 +1  : step 1
e22 := 1  1/2 = 1/2
e12 :=  1/2
tra := 1
sto :=  1/4 ==> 1/4 .  +1 1 
e11 := 1/2  1/4 = 1/4  1 +1  : step 2
e22 := 1/2  1/4 = 1/4
e12 :=  1/4
tra := 2
sto :=  1/8 ==> 1/8 .  +1 1 
e11 := 1/4  1/8 = 1/8  1 +1  : step 3
e22 := 1/4  1/8 = 1/8
e12 :=  1/8
Guess what, after M steps the boundaries are reached with: 1/2^M .  +1 1 
 1 +1  .
Backsubtitution happens to be always by: x2 = (x1 + x3)/2 .
The net result is a linear function between the two boundary values, as has to
be expected for this (quite trivial) problem.
But, of course it's better to solve a second order differential equation with
constant coefficients, as required by the above method, by: ANALYTICAL means.
Typical, if a problem is easy, then there exist so MANY methods for solving it
in an efficient way. I wonder if there are more serious applications which may
employ a technique like the one above, but which cannot be solved analytically.

* 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