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