Solution of ODE by FEM

Solution of ODE by FEM

Jerzy Pamin

e-mail:jpamin@L5.pk.edu.pl

Problem formulation

Problem to be solved

Solve boundary value problem

u ⁰⁰ (x) + 6x ² = 0 x ∈ (0, 1) , bcs: u(0) = 1 , u ⁰ (1) = − 1 2 using Galerkin formulation of FEM and 2 elements with linear interpolation.

Analytical solution

u ⁰⁰ (x) = −6x ² u ⁰ (x) = −2x ³ + C u(x) = − 1

2 x ⁴ + Cx + D u ^analit = − 1

2 x ⁴ + 3

2 x + 1

Weighted residual method

Global model via WRM

R = u ⁰⁰ (x) + 6x ² , Z 1

0

w(x)R(x)dx = 0 ∀w 6= 0 Z 1

0

wu ⁰⁰ dx + Z 1

0

w 6x ² dx = 0 ∀w

Weak formulation

− Z 1

0

w ⁰ u ⁰ dx + [wu ⁰ ] ¹ ₀ + Z 1

0

w 6x ² dx = 0 ∀w | · (−1) Z 1

0

w ⁰ u ⁰ dx − w(1)u ⁰ (1) + w(0)u ⁰ (0) − Z 1

0

w 6x ² dx = 0 , u(0) = 1

Note that u ⁰ (1) = − ¹ ₂ and u ⁰ (0) is unknown

FE discretization

2 elements with linear interpolation

0 0.5 1 x

1 2

i

x _i x _j

u _i u _j

3 2 j

e 1

x ¹ x ²

Topology e = 1 i = 1 j = 2 e = 2 i = 2 j = 3 Transformation x ^e ∈ (0, l ^e )

x = x ^e + d ^e d ¹ = 0, d ² = 0.5

Shape functions N i = 1 − ^x _l e ^e = 1 − 2x ^e N _j = ^x _{l e} ^e = 2x ^e N = [N _i , N _j ] q ^e =

 u i

u j



Bubnov-Galerkin approximation

u ≡ u ^e = N q ^e , w ≡ w ^e = N b = b ^T N ^T

Finite element equations

Integral equation for FE

Z l ^e 0

w ⁰ u ⁰ dx ^e − w(l ^e )u ⁰ (l ^e ) + w(0 ^e )u ⁰ (0 ^e ) − Z l ^e

0

w 6x ² dx ^e = 0 ∀w Z l ^e

0

w ⁰ u ⁰ dx ^e − w(l ^e )u ⁰ (l ^e ) + w(0 ^e )u ⁰ (0 ^e ) − Z l ^e

0

w 6(x ^e + d ^e ) ² dx ^e = 0 Substitute interpolation u = N q ^e , w = b ^T N ^T , invoke ∀b

Z l ^e 0

b ^T N ^0T N ⁰ q ^e dx ^e −b ^T N ^T (l ^e )u ⁰ (l ^e )+b ^T N ^T (0 ^e )u ⁰ (0 ^e )−

Z l ^e 0

b ^T N ^T 6(x ^e +d ^e ) ² dx ^e = 0 ∀b

Z l ^e 0

N ^0T N ⁰ q ^e dx ^e −N ^T (l ^e )u ⁰ (l ^e )+N ^T (0 ^e )u ⁰ (0 ^e )−

Z l ^e 0

N ^T 6(x ^e +d ^e ) ² dx ^e = 0

Notice that N ^T (l ^e ) =

 0 1



, N ^T (0 ^e ) =

 1 0



Z l ^e 0

N ^0T N ⁰ dx ^e q ^e −

 −u ⁰ (0 ^e ) u ⁰ (l ^e )



− Z l ^e

0

N ^T 6(x ^e + d ^e ) ² dx ^e = 0

Finite element equations

Finite element matrices

Z l ^e 0

N ^0T N ⁰ dx ^e q ^e −

 −u ⁰ (0 ^e ) u ⁰ (l ^e )



− Z l ^e

0

N ^T 6(x ^e + d ^e ) ² dx ^e = 0

K ^e = Z l ^e

0

N ^0T N ⁰ dx ^e , p ^e = Z l ^e

0

N ^T 6(x ^e +d ^e ) ² dx ^e , p ^e _b =

 −u ⁰ (0 ^e ) u ⁰ (l ^e )



Notice that N ⁰ = [−2 , 2]

Matrix equation for FE

K ^e q ^e − p ^e _b − p ^e = 0

K ^e q ^e = p ^e + p ^e _b

Numerical model at element level

Computations

Compute matrices for each element

K ¹ = K ² = Z 0.5

0

 −2 2



 −2 2  dx ^e =

 2 −2

−2 2



p ¹ = Z 0.5

0

 1 − 2x ¹ 2x ¹



6(x ¹ ) ² dx ¹ =

 0.0625 0.1875



p ² = Z 0.5

0

 1 − 2x ⁽²⁾ 2x ⁽²⁾



6(x ⁽²⁾ + 0.5) ² dx ⁽²⁾ =

 0.6875 1.0625



p ¹ _b =

 −u ⁰¹ (0 ¹ ) u ⁰¹ (l ¹ )



, p ² _b =

 −u ⁰² (0 ² ) u ⁰² (l ² )



Global set of equations

Assembly

Add element matrices to zeroed global arrays according to topology K = X

e

K ^e , q = X

e

q ^e , p = X

e

p ^e , p b = X

e

p ^e _b ,

Kq = p + p b

K =





2 −2 0

−2 2 + 2 −2

0 −2 2



 q =



 u ₁ u ₂ u ₃



 p =



 0.0625 0.8750 1.0625





p b =





−u ⁰¹ (0 ¹ ) u ⁰¹ (l ¹ ) − u ⁰² (0 ² )

u ⁰² (l ² )



 =





−u ⁰ (0) 0 u ⁰ (1)





Boundary conditions and solution

Set of 3 equations in 5 unknowns





2 −2 0

−2 4 −2

0 −2 2







 u 1

u 2

u 3



 =



 0.0625 0.8750 1.0625



 +





−u ⁰ (0) 0 u ⁰ (1)



 but we have boundary conditions u 1 = u(0) = 1 and u ⁰ (1) = −0.5!

Notice that until now the solution is independent of the boundary conditions.

Set of 3 equations in 3 unknowns





2 −2 0

−2 4 −2

0 −2 2







 1.0

u 2

u 3



 =



 0.0625 0.8750 1.0625



 +





−u ⁰ (0) 0

−0.5



 First solve equations 2 and 3, then equation 1 to obtain

u ₂ = 1.71875 , u ₃ = 2 , u ⁰ (0) = 1.5

Solution

Comparison of approximate and analytical solutions

Solutions Solution derivatives