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

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 deﬁne the shape functions (called tent function in this case) as shown below

The shape function is deﬁned 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

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

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 inﬂuence 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$ coeﬃcients $cj$

Now to evaluate $Ij$ we write

Now we will show the above for $j = 1$ which will be suﬃcient 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

Hence

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 ﬂip the argument as per the deﬁnition for $ψ$ shown in the top of this report.

Hence we obtain for the integral in (2)

so the above becomes integral becomes

Hence

Which becomes

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.

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