1.
up
2.
This page in PDF

HW 12 Mathematics 503, computer part, July 26, 2007

Nasser M. Abbasi

June 15, 2014
1.
This is the project description handout. It describes the problem to solve PDF
2.
Mathematica implementation
(a)
notebook
(b)
HTML
(c)
PDF
3.
Maple implementation
(a)
maple_FEM_solution_basic.mw
(b)
HTML
(c)
PDF
4.
Ada implementation
(a)
finite_elements.adb.txt
(b)
Test_Matrix.adb.txt
5.
small animation movie in swf format FEMdirect.swf

1 Derivation for the Ax=b

This is a suplement to the report for the computer project for Math 503. This includes the symbolic derivation of the A matrix and the b vector for the problem of Ax = b which is generated from the FEM formulation for this project. I also include a very short Mathematica program which implements the FEM solution.

For x = [0,L]  where L is the length, we define the shape functions (called tent function in this case) as shown below

PIC

The shape function is defined by          (x   )
ϕi(x) = ψ  h − i where

       (
       |
       ||{ 1 − z    0< z < 1
ψ (z) =
       |||
       (   0      z > 1
(1)

And ψ (z)= ψ (− z)  as shown in this diagram

PIC

Now the derivative of  ′
ϕi (x)  is given by

       (
       ||||   1   (i− 1)h < x≤ i h
       |||{   h
 ′          1
ϕi (x)= || − h  i h< x < (i+1 )h
       |||
       ||(   0      otherwise

Now we write the weak form in terms of the above shape function (which is our admissible direction). From part 1 we had

    ∫
     L ′    ′
I =  0 y(x)ϕ (x)+ q y (x)ϕ (x)− f ϕ (x) dx= 0

And Let

pict

Hence, now we pick one admissible direction at a time, and need to satisfy the above integral for each of these. Hence we write

    ∫  (         )         (         )
      L  N    ′     ′        N
Ij = 0   ∑i=1ciϕi (x) ϕj(x) +q   ∑i=1ciϕi(x) ϕj (x)− f ϕj(x) dx = 0    j = 1,2,⋅⋅⋅N

But due to sphere on influence of the ϕj(x)  extending to only xj−1⋅⋅⋅xj+1   the above becomes

        (           )         (          )
   ∫ xj+1   j+1                    j+1
Ij =        ∑  ciϕ′i (x) ϕ′j(x)+ q    ∑  ciϕi(x) ϕj(x)− f ϕj(x) dx= 0     j = 1,2,⋅⋅⋅N
     xj−1  i=j−1                  i=j−1

Hence we obtain N equations which we solve for the N coefficients cj

Now to evaluate Ij  we write

pict

Now we will show the above for j = 1  which will be sufficient to build the A matrix due to symmetry.

For j = 1
    ∫   (        )          (         )
      2h  2   ′      ′        2
I1 = 0   ∑i=1ciϕi (x) ϕ1(x)+ q  ∑i=1ciϕi(x) ϕ1(x)− f ϕ1(x) dx

Hence breaking the interval into 2 parts we obtain

pict

Hence

pict

Now set up a little table to do the above integral.

|-------|----|-----|-----------------|-------------------|
|       |    |     |                 |                   |
|Range  | ϕ′1 | ϕ′2  |       ϕ1        |        ϕ2         |
|-------|----|-----|-----------------|-------------------|
| [0,h] | 1  |N ∕A | ψ (− x+ 1)→  x  |       N∕A         |
|-------|-h--|-----|------h-------h--|-------------------|
|       | −1 |  1  |  (x    )      x |  (  x   )   x     |
--[h,2h]---h-----h---ψ--h-− 1-→-2-−-h--ψ--−-h +-2-→-h −-1-|

The above table was build by noting that for ϕj, it will have the equation   (x   )
ψ  h − i when x is under the left leg of tent. And it will have the equation   (  x   )
ψ  − h + i when x is under the right leg of the tent. This is because for x < 0  ,  the argument to ψ ()  is negative and so we flip the argument as per the definition for ψ shown in the top of this report.

Hence we obtain for the integral in (2)

pict

so the above becomes integral becomes

pict

Hence

pict

Which becomes

pict

or

pict

Therefore

                           (    )     (   )
I = c1+ c  qh+ c1 − c2 +qc   1-h + qc   1h  − fh
1   h    1 3    h    h    1  3       2  6

Hence

      (        )     (         )
        2-  2-          1-  1-
I1 = c1 h + 3qh  + c2 − h + 6qh  − fh = 0

Multiply by h we obtain

|------------------------------------------|
| I = c (2+ 2h2q)+ c (− 1+ 1h2q)− fh2 = 0  |
|  1   1    3       2      6               |
-------------------------------------------
(2)

Hence we now can set up the Ax= b system using only the above equation by taking advantage that A will be tridiagonal and there is symmetry along the diagonal.

⌊                                                                         ⌋ ⌊  ⌋       ⌊ ⌋
   (       )  (         )
||   2+ 23h2q    − 1+  16h2q        0            0           ⋅⋅⋅          0    || || c1||       ||1||
|                                                                         | |  |       | |
|| (− 1 + 1h2q) (2 + 2h2q)   (− 1+ 1h2q)       0           ⋅⋅⋅          0    || || c2||       ||1||
||       6           3            6                                        || ||  ||       || ||
||             (      12 )   (   2 2 )   (     1 2 )                       || ||  || = f h2|| ||
||      0       − 1+  6h q     2+ 3h q     − 1 + 6h q      ⋅⋅⋅          0    || || c3||       ||1||
||                                                                     .   || || .||       ||.||
|      0           0            0            0           ...           ..   | | ..|       |..|
||                                                                         || ||  ||       || ||
⌈                                                    (     1 2)  (    2 2)⌉ ⌈  ⌉       ⌈ ⌉
       0           0            0           ⋅⋅⋅      − 1+  6h q    2 + 3hq    cN         1

The following is the FEM program to implement the above, with few plots showing how close it gets to the real solution as N increases.

PIC

I also written a small Manipulate program to simulate the above. Here it is

PIC