sci.math.num-analysis

SUNA52, Fast 1-D FEM solving
============================

Disclaimer: this article is just for fun.
It has no theoretical nor practical significance whatsoever.

Consider the (one-dimensional) 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
in-between 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 element-matrix:

       | +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
* E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood