Department of Mechanical and Aerospace Engineering

MAE3456 – MEC3456

The Finite Element Method: one dimensional problems

Prof. Murray Rudman

Department of Mechanical and Aerospace Engineering,
Monash University

Lecture 29
Outline

•  Review of the weighted residual method

•  Problems with the weighted residual method

•  The finite element method

•  Dividing solution domain into sections, or elements
•  Piecewise trial functions
•  1 dimensional problem as an example

•  Application to an ODE
Review of MWR

•  Typical we can write a PDE or ODE as

∂2Q ∂2Q ∂2Q ∂2Q

2
+ 2 = f (x, y) L(Q) = 2 + 2 − f (x, y) = 0.0
∂x ∂y ∂x ∂y

•  The exact solution will have L(Q) = 0.

•  Recall in weighted residual method, we make an

approximation of solution by,
N
Q(x, y) = Q0 (x, y) + ∑ a jφ j (x, y)
j=1
•  Where φj(x, y) are known analytical function, and aj are
unknown coefficients.
Review of MWR
N
Q(x, y) = Q0 (x, y) + ∑ a jφ j (x, y)
j=1

•  Substituting the approximate solution into the PDE, we

obtain the residual R,

∂2Q ∂2Q
L (Q ) = 2 + 2 − f (x, y) = R
∂x ∂x

•  The residual will not be exactly zero for Q(x,y)

•  But given sufficient number of trial functions φj(x, y), the
residual will approach zero.

L(Q) = R → 0.0
Review of MWR
N
Q(x, y) = Q0 (x, y) + ∑ a jφ j (x, y)
j=1

•  In general, the unknown coefficients aj are determined

by requiring the integral of residual over solution domain
equal to zero, i.e. (for 2-D)

∫∫ W m (x, y) R(x, y)dx dy = 0

where, Wm is known as the weight function.

•  This is the ‘weak’ form of our original problem.

Review of MWR

∫∫ W m (x, y) R(x, y)dx dy = 0

•  The Galerkin method

•  Is a subclass of MWR
•  In this method, the weight function are selected as the
trial function, i.e.,
Wm (x, y) = φ m (x, y)

•  It can be shown that for certain trial functions, the Galerkin

method will approach the analytical answer as
m à ∞.
•  The Galerkin method is usually used with finite element
method.
Problems with MWR

•  Let us take the heat

equation as an example.
∂Q ∂2Q
−κ 2 = 0
∂t ∂x

•  If we used the weighted residual method

introduced previously, the solutions will be
approximated by using some trial basis function.
•  These trial functions are defined over the whole
solution domain.
Problems with MWR

•  In MWR, the trial functions are often defined over the

whole solution domain.

•  This leads to two issues/problems:

•  It seems unnecessary to include distant positions to
determine the solution at any given point.
•  If we use trial functions defined over the whole domain,
we end up with a full matrix to determine the unknown
coefficients aj (see the example in the previous lecture).
Finite Element Method

•  We can solve this problem by sub-dividing the domain

into several ‘finite elements’, and solving the problem
by MWR over each element.

el. 1 el. 2 el. 3 el. 4 el. 5 el. 6 el. 7 el. 8

•  This means that the solution is approximated over each

element, i.e., each trial basis function is defined only in
a single element, NOT over the entire domain.
Finite Element Method

•  In the (3-D) finite element method, we need to solve:

∫∫∫ W ( x, y, z) R dx dy dz =0
m

OR ∫∫∫ W ( ) ( x, y, z) R (
m el1 m el1 )
dxdydz +

∫∫∫ W ( ) ( x, y, z) R (
m el2 m el2 )
dxdydz +

+ ∫∫∫ Wm(elN ) ( x, y, z ) Rm(elN ) dxdydz = 0

el. 1 el. 2 el. 3 el. 4 el. 5 el. 6 el. 7 el. 8

Finite Element Method

•  We will consider a generic element

•  If we can determine the weak form of the equation on a
generic element, we can apply it to other elements

el. 1 el. 2 el. 3 el. 4 el. 5 el. 6 el. 7 el. 8

-  Within each element we

have to choose ‘trial’
functions ϕ(xE) so that we
can write:

Q ( xE , t ) = ∑ a j (t ) φ j ( xE )
Finite Element Method – trial functions

Q ( xE , t ) = ∑ a j (t ) φ j ( xE )

•  Typically we use the Lagrange Q(x)

Q1
interpolation polynomials as the Q2
trial functions.
•  The simplest case is the linear
interpolation function. φ1(x)

•  Therefore, we can approximate

the solution in the element by:
φ2(x)
Q ( xE , t ) = ∑ Q jE (t ) φ j ( xE )
= Q1E (t ) φ1 (x) + Q2E (t ) φ2 (x)
Finite Element Method

Q ( xE , t ) = ∑ Q jE (t ) φ j ( xE )
= Q1E (t ) φ1 (x) + Q2E (t ) φ2 (x)
Q(x)
Q1
Q2
•  Because we have chosen the Lagrange
polynomials, the unknown coefficients
aj are simply Q1 and Q2 whose values
we are seeking at the two nodes. φ1(x)

•  Mathematically, the trial functions (i.e.

Lagrange polynomials) are written:
(E ) (E ) φ2(x)
x − x2 x − x1
φ1 = (E ) φ 2 = (E )
x1 − x 2( E ) x 2 − x1( E )
Department of Mechanical and Aerospace Engineering

MAE3456 – MEC3456

The Finite Element Method:

A one dimensional example

Prof. Murray Rudman

Department of Mechanical and Aerospace Engineering,
Monash University

Lecture 29
Application to a One-Dimensional Example

•  Example:
•  apply the FEM to solve the 1D heat equation

d 2T
= − f (x) With 0 ≤ x ≤ L
2
dx

€
Application to One-Dimensional Example

d 2T
2
= − f (x) over 0 ≤ x ≤ L
dx

•  Solution
•  The first step is to approximate the solution on a generic
element.
T ( xE ) = T1Eφ1 (x) + T2Eφ2 (x)

•  (E) refers to local “Element”

coordinates, not global coordinates

" (E) %
x − x
T ( x ) = T1E $ ( E ) 2 ( E ) '
# x1 − x2 &
" (E) %
x − x
+ T2E $ ( E ) 1 ( E ) '
# x2 − x1 &
Application to One-Dimensional Example

d 2T " x − x (E ) % " x − x (E ) %
= − f (x) T ( x ) = T1E $ (E ) 2 (E ) ' + T2E $ (E ) 1 (E ) '
dx 2 # x1 − x2 & # x2 − x1 &

•  The residual is defined as:

€
d 2T
R = L(T ) = 2 + f (x)
dx

•  By using MWR, the weak form of the original problem can

be expressed as:
x2 x2 " d 2T %
∫ x1
Wm (x) R(x)dx = ∫ x1
Wm (x)$ 2 + f (x)' dx = 0
# dx &
m = 1, 2
Application to One-Dimensional Example

d 2T x2 ! d 2T $
dx 2
= − f (x) ∫ x1
Wm (x)# 2 + f (x)& dx = 0
" dx %

•  The choice of weight function Wm(x) depends on methods.

€
•  The finite element method uses the Galerkin method, in
which Wm(x) is selected to be a trial basis function.

•  Therefore, we have:

# x − x2( E ) & # x − x1( E ) &

W1 (x) = φ1 (x) = % ( E ) (E) ( W2 (x) = φ2 (x) = % ( E ) (E) (
x
$ 1 − x 2 ' x
$ 2 − x1 '
Application to One-Dimensional Example

d 2T x2 ! d 2T $
dx 2
= − f (x) ∫ x1
Wm (x)# 2 + f (x)& dx = 0
" dx %

•  Substituting the weight function into the above weak form

€ of the equation, we have:
x2 ! d 2T $
∫ x1
φ m (x)# 2 + f (x)& dx = 0
" dx %

•  Since f(x) is known, we move it to the RHS,

x2 d 2T x2
∫ x1
φ m (x) 2 dx = − ∫ x φ m (x) f (x)dx
dx 1
Application to One-Dimensional Example

d 2T x2 d 2T x2

dx 2 = − f (x) ∫ x1
φ m (x) 2 dx = − ∫ x φ m (x) f (x)dx
dx 1

•  Here we manipulate the LHS by integrating by parts:

€ b b b
∫ udv = uv
a a
− ∫ vdu
a

•  Therefore the 2nd derivative term is written

€ x2 d 2T x2 d " dT %
∫ x1
φ m (x) 2 dx = ∫ x φ m (x) $ ' dx
dx 1 dx # dx &
x2
dT x2 dT dφ m
= φ m (x)
dx
− ∫ x1
dx dx
dx
x1
Application to One-Dimensional Example

d 2T x2 d 2T x2

dx 2 = − f (x) ∫ x1
φ m (x) 2 dx = − ∫ x φ m (x) f (x)dx
dx 1

•  For m=1, we have:

€ x2
x2 d 2T dT x2 dT dφ1
∫ x1 φ1 (x) dx 2 dx = φ1 (x) dx − ∫ x1
dx dx
dx
x1

•  The first term on the RHS immediately above is:

x2
dT dT (x2 ) dT (x1 )
φ1 (x) = φ1 (x2 ) − φ1 (x1 )
dx x1 dx dx
dT (x2 ) dT (x1 )
= (0) − (1)
dx dx
dT (x1 )
=−
dx
Application to One-Dimensional Example

•  Therefore the LHS of the equation is written

x2 d 2T dT (x1 ) x2 dT dφ1
∫ x1
φ1 (x) 2 dx = −
dx dx
− ∫ x1
dx dx
dx

•  Substituting the trial solution T and trial function φ1(x)

into the second term above (after some simple algebra)
x2 dT dφ1
∫ x1 dx dx dx
x2 d
) # x − x (E ) & # x − x ( E ) &,
d # x − x (E ) &
=∫ +T1E %% (E ) 2 (E ) (( + T2E %% (E ) 1 (E ) ((. %% (E ) 2 (E ) ((dx
x1
dx +* $ x1 − x2 ' $ x2 − x1 '.- dx $ x1 − x2 '
1
= (T1 − T2 )
x2 − x1
Application to One-Dimensional Example

x2 d 2T x2
∫ x1
φ m (x) 2 dx = − ∫ x f (x)φ m (x)dx m = 1, 2
dx 1

•  Therefore, we obtain the first equation, i.e., for m=1

1 dT (x1 ) x2
(T1 − T2 ) = − + ∫ f (x)φ1 (x)dx
x2 − x1 dx x1

•  If we do the mathematical derivation in a similar fashion,

we can obtain the equation for m=2

1 dT (x2 ) x2
( −T1 + T2 ) = + ∫ f (x)φ2 (x)dx
x2 − x1 dx x1
Application to One-Dimensional Example

1 dT (x1 ) x2
(T1 − T2 ) = − + ∫ f (x)φ1 (x)dx
x2 − x1 dx x1

1 dT (x2 ) x2
( −T1 + T2 ) = + ∫ f (x)φ2 (x)dx
x2 − x1 dx x1

•  These two equations can be expressed in a matrix form,

( dT (x1 ) + ( x2 +
− f (x)φ (x)dx
1 " 1 −1 % (T1 + . . dx . . ∫x1
. 1
.
$ ') , ) = +
, ) ,
x2 − x1 # −1 1 & *T2 - . dT (x2 ) . . x2 f (x)φ (x)dx .
.- * ∫x1
2
.* dx -
•  This matrix problem is the ‘weak’ governing equation for a
generic element.
Application to One-Dimensional Example

•  But to solve the problem, we need to consider the entire

domain

T3 T4
T2
T1 T5

x1 x2 x3 x4 x5

El. 1 El. 2 El. 3 El. 4

x1( El. 2 ) , T1( El. 2 ) x2( El. 2 ) , T2( El. 2 )

Application to One-Dimensional Example

•  We must sum up each elemental equation together to

determine the final answer.
•  Recall the total problem is

∫∫∫ W ( x, y, z) R dx dy dz =0
m

•  Which can be written on an element by element way

∫∫∫ W ( ) ( x, y, z) R (
m el1 m el1 )
dxdydz +

∫∫∫ W ( ) ( x, y, z) R (
m el2 m el2 )
dxdydz +

+ ∫∫∫ Wm(elN ) ( x, y, z ) Rm(elN ) dxdydz = 0
Application to One-Dimensional Example

•  We write the problem in matrix-

vector form " T1 %
$ '
•  Note, we will need to specify $ T2 '
a numbering system in 2D
and 3D problems A⋅$ T3 ' = [ f ]
$ '
$ T4 '
$ '
$# T5 '
&

•  Where A and f are built from the equations in each element

( dT (x1(E ) ) , ( x2 ,
− f (x)φ (x)dx
* dx * * ∫x1
*
(E )
1 " 1 −1 % (*T1 ,* * 1
*
$ ') =
- ) +
- ) -
x2(E ) − x1(E ) # −1 1 & +*T2(E ) .* * dT (x2(E ) ) * * x2 f (x)φ (x)dx *
*+ dx
∫
*. + x1
2
.

•  Note that T1(E) is T on the left side of element E and T2(E) is T on the
right side (not to be confused with T1, T2)
Application to One-Dimensional Example

•  To make the analysis clearer, we specify the geometry

and heat source as

f(x) = 10.0
0.0 10.0

•  Therefore, the governing equation for a generic element

is:
( dT (x (E ) ) ,
" 0.4 −0.4 %(*T (E ) ,* **− dx ** (12.5,
1

1
$ ') (E ) - = ) (E )
-+) -
# −0.4 0.4 &*+T *. * dT (x ) * +12.5.
2 2
*+ dx *.
Application to One-Dimensional Example

•  Look at the equation for the first element:

0.0
✔ 10.0

# 0.4 −0.4 0 0 0 &)T1 - )−dT(x1 ) /dx + 12.5-

%−0.4 0.4 0 0 0 (+T + +dT(x ) /dx + 12.5 +
% (++ 2 ++ ++ 2 ++
% 0 0 0 0 0 (*0 . = *0 .
% 0 0 0 0 0 (+0 + +0 +
% (+ + + +
%$ 0 0 0 0 0 ('+,0 +/ +,0 +/

( dT (x1(E ) ) ,
−
" 0.4 −0.4 % (*T1(E ) ,* ** dx ** (12.5 ,
$ ' ) (E ) - = ) -+ ) -
€ # −0.4 0.4 (E )
& *+T2 *. * dT (x2 ) * + 12.5 .
*+ dx *.
Application to One-Dimensional Example

•  Now consider the equation of the second element:

✔
0.0 10.0

" 0 0 0 0 0 %( 0 , ( 0 ,
$ '* * * *
$ 0 0.4 −0.4 0 0 '*T2 * *−dT(x2 ) / dx +12.5*
* * * *
$ 0 −0.4 0.4 0 0 ')T3 - = )dT(x3 ) / dx +12.5 -
$ 0 0 0 0 0 '* 0 * * 0 *
$ '* * * *
$# 0 0 0 0 0 '&*0 * *0 *.
+ . +
( dT (x1(E ) ) ,
−
" 0.4 −0.4 % (*T1(E ) ,* ** dx ** (12.5 ,
$ ' ) (E ) - = ) (E ) -+ ) -
# −0.4 0.4 & *+T2 *. * dT (x2 ) * + 12.5 .
*+ dx *.
Add equations for Node 1 and 2 together
" 0.4 −0.4 0 0 0 %(T1 , " 0 0 0 0 0 % (0 ,
$ ' * * $ '* *
$ −0.4 0.4 0 0 0 '** 2 ** $ 0 0.4 −0.4 0 0 '**T2 **
T
$ 0 0 0 0 0 ')0 - + $ 0 −0.4 0.4 0 0 ')T3 - =
$ '* * $ '* *
$ 0 0 0 0 0 0
'* * $ 0 0 0 0 0 ' *0 *
$ 0 0 0 0 0 ' *0 * $ 0 0 0 0 0 ' *0 *
# &+ . # &+ .
(−dT (x ) / dx +12.5, (0 ,
1
* * * *
* dT (x 2
) / dx +12.5 * * −dT (x 2
) / dx +12.5*
* * * *
) 0 +
- ) dT (x3
) / dx +12.5 -
*0 * * *
* * * 0 *
*0 * *+0 *.
+ .

# 0.4 −0.4 0 0 0 &)T1 - )−dT(x1 ) /dx + 12.5-

%−0.4 0.4 + 0.4 −0.4 0 0 (+T + +12.5 + 12.5 +
•  This gives %
% 0 −0.4
(++ 2 ++ ++ ++
0.4 0 0 (*T3 . = *dT(x 3 ) /dx + 12.5 .
% 0 0 0 0 0 (+0 + +0 +
% (+ + + +
%$ 0 0 0 0 0 ('+,0 +/ +,0 +/
Application to One-Dimensional Example

•  Equations for first 2 elements:

✔ ✔
0.0 10.0

# 0.4 −0.4 0 0 0 &)T1 - )−dT(x1 ) /dx + 12.5-

%−0.4 0.4 + 0.4 −0.4 0 0 (+T + +12.5 + 12.5 +
% (++ 2 ++ ++ ++
% 0 −0.4 0.4 0 0 (*T3 . = *dT(x 3 ) /dx + 12.5 .
% 0 0 0 0 0 (+0 + +0 +
% (+ + + +
%$ 0 0 0 0 0 ('+,0 +/ +,0 +/

( dT (x1(E ) ) ,
−
" 0.4 −0.4 % (*T1(E ) ,* ** dx ** (12.5 ,
$ ' ) (E ) - = ) -+ ) -
€ # −0.4 0.4 (E )
& *+T2 *. * dT (x2 ) * + 12.5 .
*+ dx *.
Application to One-Dimensional Example

•  Now consider the equation of the third element:

✔
0.0 10.0

" 0 0 0 0 0 %( 0 , ( 0 ,
$ '* * * *
$ 0 0 0 0 0 '* 0 * * 0 **
* * *
$ 0 0 0.4 −0.4 0 ')T3 - = )−dT(x3 ) / dx +12.5-
$ 0 0 −0.4 0.4 0 '*T4 * *dT(x4 ) / dx +12.5 *
$ '* * * *
$# 0 0 0 0 0 '&*0 * *0 *.
+ . +

( dT (x1(E ) ) ,
−
" 0.4 −0.4 % (*T1(E ) ,* ** dx ** (12.5 ,
$ ' ) (E ) - = ) (E ) -+ ) -
# −0.4 0.4 & *+T2 *. * dT (x2 ) * + 12.5 .
*+ dx *.
Application to One-Dimensional Example

•  Now we add the 3rd equation the first two:

✔
0.0 10.0

# 0.4 −0.4 0 0 0&)T1 - )−dT(x1) /dx + 12.5-

%−0.4 0.8 −0.4 0 0(+T2 + +25 +
% (++ ++ ++ ++
% 0 −0.4 0.4 + 0.4 −0.4 0(*T3 . = *12.5 + 12.5 .
% 0 0 −0.4 0.4 0(+T4 + +dT(x 4 ) /dx + 12.5 +
% (+ + + +
%$ 0 0 0 0 0('+,0 +/ +,0 +/

( dT (x1(E ) ) ,
−
" 0.4 −0.4 % (*T1(E ) ,* ** dx ** (12.5 ,
€ $ ' ) (E ) - = ) (E ) -+ ) -
# −0.4 0.4 & *+T2 *. * dT (x2 ) * + 12.5 .
*+ dx *.
Application to One-Dimensional Example

•  Similarly, add the equation of the 4th element to the

others:
✔
0.0 10.0

" 0.4 −0.4 0 0 0 %(T1 , (−dT(x1 ) / dx +12.5,

$ '* * * *
$ −0.4 0.8 −0.4 0 T
0 '* 2 * * 25 **
* * *
$ 0 −0.4 0.8 −0.4 0 ')T3 - = )25 -
$ 0 0 −0.4 0.4 + 0.4 −0.4 '*T4 * *12.5+12.5 *
$ '* * * *
$# 0 0 0 −0.4 0.4 '&*T * *dT(x ) / dx +12.5 *
+ 5. + 5 .

( dT (x1(E ) ) ,
−
" 0.4 −0.4 % (*T1(E ) ,* ** dx ** (12.5 ,
$ ' ) (E ) - = ) (E ) -+ ) -
# −0.4 0.4 & *+T2 *. * dT (x2 ) * + 12.5 .
*+ dx *.
Application to One-Dimensional Example

•  The final matrix problem looks like

" 0.4 −0.4 0 0 0 %(T1 , (−dT (x1 ) / dx +12.5,

$ '* * * *
$ −0.4 0.8 −0.4 0 T
0 '* 2 * * 25 **
* * *
$ 0 −0.4 0.8 −0.4 0 ')T3 - = )25 -
$ 0 0 −0.4 0.8 −0.4 '*T4 * *25 *
$ '* * * *
$# 0 0 0 −0.4 0.4 '&*T * * dT (x ) / dx +12.5 *
+ 5. + 5 .

•  We have almost everything we need, except there are

two derivative terms in the RHS vector

•  However, we haven’t yet mentioned BC’s

Application to One-Dimensional Example

•  Boundary condition
•  If we apply the boundary condition by specifying the heat
flux at the two ends of the solution domain as

dT (x1 ) dT (x5 )
= 66 = −34
dx 0.0 10.0 dx

•  Then we can replace the derivative terms below

" 0.4 −0.4 0 0 0 % (T1 , (−dT (x1 ) / dx + 12.5 ,

$ ' *T * *25 *
$ −0.4 0.8 −0.4 0 0 '* 2 * *
* * * **
$ 0 −0.4 0.8 −0.4 0 ' )T3 - = )25 -
$ 0 0 −0.4 0.8 −0.4 ' *T4 * *12.5 + 12.5 *
$ '* * * *
# 0 0 0 −0.4 0.4 & *+T5 *. *+dT (x5 ) / dx + 12.5 *.
Application to One-Dimensional Example

•  Then we have a tridiagonal matrix problem to solve to

obtain the final answer.
T −53.5
" 0.4 −0.4 0 0 0 %( 1 , ( ,
$ ' *T2 * *25 *
$ −0.4 0.8 −0.4 0 0 '* * * **
* * *
$ 0 −0.4 0.8 −0.4 0 ' )T3 - = )25 -
$ 0 0 −0.4 0.8 −0.4 ' *T4 * *12.5 + 12.5 *
$ '* * * *
# 0 0 0 −0.4 0.4 & *T * *−21.5 *.
+ 5. +

•  In comparison to the global MWR method we discussed

in last lecture, this can be solved much more efficiently.
Application to One-Dimensional Example

•  After solving the matrix problem, we can plot the

solution as:

T1 = 40
T2 = 173.75
T3 = 245
T4 = 253.75
T5 = 200
Summary

•  The implementation of finite element method can be

summarized as a step-by-step procedure

•  Dividing into simple shape elements

•  Approximate solution of PDE on

each element (Galerkin + local basis
function)
•  Sum up the weak form equation of
each element
•  Applying the boundary conditions

•  Solving the linear algebraic system

•  Displaying results graphically

Conclusion

•  We have considered a new method today.

•  You should be able to describe the each of the following

after this lecture
•  Weighted Residual method
•  Galerkin technique
•  The general idea of finite element method
•  The steps we used to develop the finite element method