ODE Finite Element Method
ODE Finite Element Method
ODE Finite Element Method
Contents
1 Introduction 1
4 References 31
This report is a basic review of weighted residual methods and FEM. Some basic differential
equations are used to illustrate the method. Mathematica and Matlab code written to solve
numerically a first and second order ODE using FEM.
1 Introduction
This is a basic review of Finite Elements Methods from Mathematical point of view with examples
of how it can be used to numerically solve first and second order ODE’s. Currently I show how
to use FEM to solve first and second order ODE. I am also working on a detailed derivation and
implementation using FEM to solve the 2D Poisson’s equation but this work is not yet completed.
FEM is a numerical method for solving differential equations (ordinary or partial). It can also be
used to solve non-linear differential equations but I have not yet studied how this is done. FEM
is a more versatile numerical method than the finite difference methods for solving differential
equations as it supports more easily different types of geometry and boundary conditions, in
addition the solution of the differential equation found using FEM can be used at any point in the
domain and not just on the grid points as the case is with finite difference methods. On the other
hand FEM is more mathematically complex method, and its implementation is not as straight
forward as with finite difference methods.
Considering only ordinary differential equations with constant coefficients over the x domain (real
1
2
line).
dy
a + b y(x) = f (x)
dx
dy
a + b y(x) − f (x) = 0
dx
defined over 0 ≤ x ≤ 1 with the boundary condition y(0) = y0 . In the above, only when y(x) is
the exact solution, call it ye (x), do we have the above identity to be true.
In other words, only when y = ye we can write that a dy
dx + b ye (x) − f (x) = 0. Such a differential
e
dye
L(ye , x) → a + bye (x) − f (x)
dx
→0
FEM is based on the weighted residual methods (WRM) where we assume that the solution of an
differential equation is the sum of weighted basis functions represented by the symbol φi in here.
This is in essence is similar to Fourier series, where we represent a function as a weighted sums
of series made up of the basis functions which happened to be in that case the sin and cosine
functions.
So the first step in solving the differential equation is to assume that the solution, called ỹ(x) ,
can be written as
N
X
ỹ(x) = qj φj (x)
j=1
Where qj are unknown coefficients (the weights) to be determined. Hence the main computational
part of FEM will be focused on determining these coefficients.
When we substitute this assumed solution in the original ODE such as shown in the above example,
equation (1) now becomes
dỹ
L(ỹ, x) = a + bỹ(x) − f (x)
dx
= R(x)
Where R(x) is called the differential equation residual, which is a function over x and in general
will not be zero due to the approximate nature of our assumed solution. Our goal is to to determine
the coefficients qj which will make R(x) the minimum over the domain of the solution.
The optimal case is for R(x) to be zero over the domain. One method to be able to achieve this is
by forcing R(x) to meet the following requirement
Z
R(x) vi (x) dx = 0
Ω
3
for all possible sets of function vi (x) which are also defined over the same domain. The functions
vi (x) are linearly independent from each others. If we can make R(x) satisfy the above for each
one of these functions, then this implies that R(x) is zero. And the solution ỹ(x) will be as close as
possible to the exact solution. We will find out in FEM that the more elements we use, the closer
to the exact solution we get. This property of convergence when it comes to FEM is important,
but not analyzed here.
Each one of these functions vi (x) is called a test function (or a weight function), hence the name
of this method.
In the Galerkin method of FEM, the test functions are chosen from the same class of functions as
the trial functions as will be illustrated below.
By making R(x) satisfy the above integral equation (called the weak form of the original differential
equation) for N number of test functions, where N is the number of the unknown coefficients qi ,
then we obtain a set of N algebraic equations, which we can solve for qi .
The above is the basic outline of all methods based on the weighted residual methods. The choice
of the trial basis functions, and the choice of the test functions, determine the method used.
Different numerical schemes use different types of trial and test functions.
In the above, the assumed solution ỹ(x) is made up of a series of trial functions (the basis). This
solution is assumed to be valid over the whole domain. This is called a global trial function. In
methods such as Finite Elements and Finite volume, the domain itself is discretized, and the
assumed solution is made up of a series of solutions, each of which is defined over each element
resulting from the discretization process.
In addition, in FEM, the unknown coefficients, called qi above, have a physical meaning, they are
taken as the solution values at each node. The trial functions themselves are generated by using
polynomial interpolation between the nodal values. The polynomial can be linear, quadratic or
cubic polynomial or higher order.
Lagrangian interpolation method is normally used for this step. The order of the polynomial is
determined by the number of unknowns at the nodes. For example, if our goal is to determine the
axial displacement at each node, then we have 2 unknowns, one at each end of the element. Hence
a linear interpolation will be sufficient in this case, since a linear polynomial a0 + a1 x contain 2
unknowns, the a1 and a0 . If in addition to the axial displacement, we wish to also solve for the
rotation at each end of the element, hence we have a total of 4 unknowns, 2 at each end of the
element, which are the displacement and the rotation.
Hence in this case the minimum interpolating polynomial needed will be a cubic polynomial
a0 + a1 x + a2 x2 + a3 x3 . In the examples below, we assume that we are only solving for axial
displacement, hence a linear polynomial will be sufficient.
At first, we will work with global trial functions to illustrate how to use weighted residual method.
4
Where φj (x) is the trial function which we have to choose, and qj are the N unknown coefficients
to be determined subject to a condition which will be shown below. ỹ0 is the assumed solution
which needs to be valid only at the boundary conditions. Hence in this example, since we are given
that the solution must be 1 at the initial condition x = 0, then ỹ0 = 1 will satisfy this boundary
condition. Hence our trial solution is
N
X
ỹ = 1 + qj φj (x)
j=1
STEP 2
Now we decide on what trial function φ(x) to use. For this example, we can select the trial functions
to be polynomials in x or trigonometric functions. Let us choose a polynomial φj (x) = xj , hence
our assumed solution becomes
N
X
ỹ = 1 + qj x j (1)
j=1
Now we need to determine the coefficients qj , and then our solution will be complete. This is done
in the following step.
STEP 3
Substituting the above assumed solution back into the original ODE, we obtain the residual R(x)
dỹ
− ỹ(x) = R(x)
dx
R(x) is the ODE residual. This is the error which will result when the assumed solution is used
in place of the exact solution.
5
Our goal now is to reduce this residual to minimum. The way we achieve this is by requiring that
the residual satisfies the following integral equation
Z1
vi (x) R(x) dx = 0 i = 1· · · N (3)
0
The above is a set of N equations. The integration is carried over the whole domain, and vi (x) is
a weight (test) function, which we have to also select. Depending on the numerical scheme used,
the test function will assume different forms.
For the Galerkin method, we select the test function to be from the same family of functions as
the trial (basis) functions. Hence in this example, let us select the test function to the following
polynomial
STEP 4
We now choose a value for N and solve the set of equations generated from (3). Let us pick N = 3,
hence R(x) becomes
3
X
qj jxj−1 − xj
R(x) = −1 +
j=1
= −1 + q1 1 − x1 + q2 2x − x2 + q3 3x2 − x3
The solution is
72 30 20
q1 = , q2 = , q3 =
71 71 71
Hence our assumed series solution is now complete, using the above coefficients, and from equation
(1) we write
N
X
ỹ = 1 + qj xj
j=1
ỹ = 1 + q1 x + q2 x2 + q3 x3
Hence
72 30 20
ỹ(x) = 1 + x + x2 + x3
71 71 71
Let use compare the above solution to the exact solution y = ex by comparing the values of the
solution over a number of points. This is done using the following small Mathematica code
7
Figure 1: Using N = 3
To make this more useful, we can examine how the error changes as N changes. The following
Mathematica code determines the solution and calculates the same table as above for N = 1· · · 5
This code below plots the absolute error as N changes. Notice that the number of peaks in the
error plot is also N which is the polynomial order (the trial solution) used to approximate the
exact solution, which is to be expected.
9
Which already satisfies the boundary conditions. For the test function, the same function as above
is used hence the test (weight) function is
The book above then reduces the residual equation to a symmetric form by doing integration by
parts before solving it for the coefficients. In here, we will use the unsymmetrical weak form and
compare the results with those shown in the above textbook. We now start again with the same
steps as we did in the above example.
step 1
10
step 2
Selecting trial basis function φj (x). As mentioned above, we select φj (x) = sin ((2j − 1) πx) ,
hence the trial solution is
N
X
ỹ = qj sin ((2j − 1) πx)
j=1
step 3
Substituting the above assumed solution into the original ODE, gives the differential equation
residual R(x)
d4 ỹ
+ ỹ(x) − 1 = R(x)
dx4
4 N N
d X X
qj φj (x) + qj φj (x) − 1 = R(x)
dx4
j=1 j=1
Notice the requirement above that the trial basis functions must be 4 times differentiable, which
is the case here. From above we obtain
N
X
R(x) = 1 + ((2j − 1) π)4 qj sin ((2j − 1) πx) − 1
j=1
Our goal now is to reduce this residual to minimum. The way we achieve this is by requiring that
the residual satisfies the following weak form integral equation
Z
vi (x) R(x) dx = 0 i = 1· · · N (3)
Ω
The above is a set of N equations. The integration is carried over the whole domain, and vi (x) is
a weight (test) function, which we have to also select. As mentioned above, in this problem we
select the test function to be
vi (x) = sin ((2i − 1) πx) (4)
step 4
Deciding on a value for N and solving the set of equations generated from (3). Let us pick N = 3,
hence R(x) becomes
X3
R(x) = 1 + ((2j − 1) π)4 (qj sin (2j − 1) πx) − 1
j=1
1 + π 4 (q1 sin πx) + 1 + (3π)4 (q2 sin (3πx)) + 1 + (5π)4 (q3 sin 5πx) − 1
=
11
Z1 n o
1 + π 4 (q1 sin πx) + 1 + (3π)4 (q2 sin 3πx) + 1 + (5π)4 (q3 sin 5πx) − 1 dx = 0
(sin πx) i=
0
Z1 n o
1 + π 4 (q1 sin πx) + 1 + (3π)4 (q2 sin (3πx)) + 1 + (5π)4 (q3 sin 5πx) − 1 dx = 0
(sin 3πx) i=
0
Z1 n o
1 + π 4 (q1 sin πx) + 1 + (3π)4 (q2 sin 3πx) + 1 + (5π)4 (q3 sin 5πx) − 1 dx = 0
(sin 5πx) i=
0
Carrying the integration above and simplifying and solving for qi gives the numerical solution.
This below is a Mathematica code which solves this problem for different N values, and compares
the error as N changes. The error shown is the percentage error in the solution (approximate
compared to exact) for up to N = 10. The result below agrees well with the result in Professor’s
Atluri textbook.
12
Figure 5: Error for different N using FEM for second order ODE
The function φj (x) is called the j th basis function. qj φj (x) is called a trial function.
The function φj (x) is made up of functions called the shape functions Nk as they are normally
called in structural mechanics books.
qj are the unknown coefficients Z
which are determined by solving N set of equations generated by
setting N integrals of the form R(x) vi (x) dx to zero. Where R(x) is the differential equation
Ω
residual. In all what follows N is the taken as the number of nodes.
13
In FEM, we also carry the same basic process as was described above, the differences are the
following:
Now we divide the domain itself into a number of elements. Earlier we did not do this.
Next, the φj (x) function is found by assuming the solution to be an interpolation between the
nodes of the element. The solution values at the nodes are the qi and are of course unknown except
at the boundaries as given by the problem.
We start by deciding on what interpolation between the nodes to use. We will use polynomial
interpolation here. Then qj φj (x) will become the interpolation function.
In addition, the coefficients qj represent the solution at the node j. These are the unknowns, which
we will solve for by solving the weak form integral equation as many times as there are unknowns
to solve for.
By solving for the nodal values, we can then use the interpolating function again to find the
solution at any point between the nodes.
This diagram illustrates the above, using the first example given above to solve a differential
equation du
dx − u(x) = 0 with the given boundary condition of u(0) = 1 and defined over 0 ≤ x ≤ 1
This is the given These are the unknowns we are solving for.
solution value at
the boundary
y
q3
q1 x
q2 q4
Using linear interpolation, the solution u(x) when x is located inside first element, is found from
u(x) = q1 + slope (x − x0 )
q2 −q1
But the slope of the linear interpolating line over the first element is x2 −x1 , hence the above
becomes
q2 − q1
u(x) = q1 + (x − x1 )
x2 − x1
(x2 − x) (x − x1 )
= q1 + q2
(x2 − x1 ) (x2 − x1 )
14
The above is the linear interpolating polynomial. We could also have used the formula of Lagrangian
interpolation to arrive at the same result.
The above is the approximate solution which is valid over the first element only. Using superscript
to indicate the element number, and assuming we have equal division between nodes of length
say h (i.e. element length is h) then we write
(x2 − x) (x − x1 )
u1 (x) = q1 + q2
h h
Again, the above is valid for x1 ≤ x ≤ x2 . We now do the same for the second element
(x3 − x) (x − x2 )
u2 (x) = q2 + q3
h h
The above is valid for x2 ≤ x ≤ x3 . And finally for the 3rd element
(x4 − x) (x − x3 )
u3 (x) = q3 + q4
h h
The above is valid for x2 ≤ x ≤ x3 . This is now illustrated in the following diagram
1
u2
u
u3
y
q2 x
q1
q2
q3
Since our goal is to express the global approximate solution u(x) as a series sum of basis functions
each multiplied by qj , we now rewrite each of the uj to allow this, as follows
The above is valid for x1 ≤ x ≤ x2 . Notice the use of the following notation: Since each element will
have defined on it two shape functions, N1 (x) and N2 (x), one per node, then we use a superscript
to indicate the element number. Hence for element 1, we will write its two shapes functions as
N11 (x) and N21 (x).
We now do the same for the second element
N12 (x) N22 (x)
z }| { z }| {
(x 3 −x) (x−x2 )
u2 = q2 + q3
h h
= q2 N12 (x) + q3 N22 (x)
The above is valid for x2 ≤ x ≤ x3 . And finally for the 3rd element
N13 (x) N23 (x)
z }| { z }| {
(x 4 −x) (x−x3 )
u3 = q3 + q4
h h
= q3 N13 + q4 N23
u(x) = u1 + u2 + u3
= q1 N11 + q2 N21 + q2 N12 + q3 N22 + q3 N13 + q4 N23
φ1 = N11
x2 − x
=
h
And the shape function for node 2 is
φ2 = N21 + N12
(x − x1 ) (x3 − x)
= +
h h
The shape function for node 3 is
φ3 = N22 + N13
(x − x2 ) (x4 − x)
= +
h h
The shape function for the last node is
φ4 = N23
x − x3
=
h
We see that the shape function for any internal node is
x − xj−1 xj+1 − x
φj = +
h h
16
u 2 q 2N 21x q 3N 22x
u 1 q 1N 11x q 2N 12x u 3 q 3N 31 q 4N 32
y
q3
q1 x
q2
q4
u1 u2 uj u N1
u j1
q3 qN
q1
q2 qj
q j1 q N1
element N-1 X=1
X=0 element 1 element 2 element j-2 element j
Next the residual R(x) is found by substituting the above global solution into the original differ-
ential equation as we did before.
Let The differential equation to solve be
du
a + bu(x) = f (x)
dx
Defined over 0 ≤ x ≤ 1 with the initial condition u(0) = u0 .
N is the number of nodes, therefore the residual is
N N
dX X
R(x) = a qi φi + b qi φi − f (x)
dx
i=1 i=1
N
X
qi aφ0i + bφi − f (x)
R(x) =
i=1
N equations are obtained from the above, which are solved for the qi .
The derivatives of the shape functions are found first. Assuming in this example that the domain
is divided into 3 elements results in
i φi φ0i
(x2 −x) −1
1 h h
n o
(x−x1 ) (x3 −x) 1 −1
2 h , h h, h
(x−x2 ) 1
3 h h
(x−x1 )
Since φ1 is not zero only over x1 ≤ x ≤ x2 , and given that φ1 = x2h−x , φ01 = −1h , φ2 = h due
to the range of integration limits, and that φ02 = h1 , then the above can simplifies to
x=x
Z 2
φ1 q1 aφ01 + bφ1 + q2 aφ02 + bφ2 − f (x) dx = 0
x=x1
x=x
Z 2
(x2 − x) −a (x2 − x) 1 (x − x1 )
q1 +b + q2 a + b − f (x) dx = 0
h h h h h
x=x1
The above equation gives the first row in the global stiffness matrix for any first order linear ODE
of the form a du
dx + bu(x) = f (x). The above shows that numerical integration is only needed to be
x=x
R 2
performed on the term (x2 − x) f (x).
x=x1
Next the last equation is found, which will be the last row of the stiffness matrix. For the last
node only j = N the following results
x=x
Z N "N # !
X
0
φN qi aφi + bφi − f (x) dx = 0
x=xN −1 i=1
x=xN −1
Which simplifies to
x=x
Z N
(3a + b(xN −1 − xN )) (xN −1 − xN )2 (−3a + 2b(xN −1 − xN )) (xN −1 − xN )2 1
N −1 −q N − (x − xN −1 ) f (x) d
6h2 6h2 h
x=xN −1
Hence the last row of the stiffness matrix can be determined directly except for the term under
the integral which needs to be evaluated using integration.
Now the equation that represents any internal node is found. This will be any row in the global
stiffness matrix between the first and the last row.
For any j other than 1 or N the following results
x=x
Z j x=x
Z j+1
N N
" # ! " # !
X X
qi aφ0i + bφi qi aφ0i + bφi − f (x) dx = 0
φj − f (x) dx + φj
x=xj−1 i=1 x=xj i=1
Where the integral was broken into two parts to handle the domain of influence of the shape
functions.
x=x
Z j x=x
Z j+1
qj−1 aφ0j−1 qj aφ0j φj qj aφ0j + bφj + qj+1 aφ0j+1 + bφj+1 − f (
φj + bφj−1 + + bφj − f (x) dx+
x=xj−1 x=xj
20
For
x − xj−1 0 1 xj − x 0 −1
xj−1 ≤ x ≤ xj , φj = , φj = , φj−1 = , φj−1 =
h h h h
xj+1 − x 0 −1 x − xj 0 1
xj ≤ x ≤ xj+1 , φj = , φj = , φj+1 = , φj+1 =
h h h h
Hence the weak form integral can be written as
x=x
Z j
x − xj−1 −a xj − x a x − xj−1
qj−1 +b + qj +b − f (x) dx
h h h h h
x=xj−1
x=x
Z j+1
xj+1 − x −a xj+1 − x a x − xj
+ qj +b + qj+1 +b − f (x) dx = 0
h h h h h
x=xj
Which simplifies to
x=x
(3a − 2b(xj−1 − xj )) (xj−1 − xj )2 Z j
!
2
(−3a − b(xj−1 − xj )) (xj−1 − xj ) 1
qj−1 +qj − (x − xj−1 ) f (x) dx
6h2 6h2 h
x=xj−1
! ! x=x
Z j+1
(−3a − 2b(xj − xj+1 )) (xj − xj+1 )2 (3a − b(xj − xj+1 )) (xj − xj+1 )2 1
qj +qj+1 − (xj+1 − x) f (x) dx
6h2 6h2 h
x=xj
For equal distance between elements, xj − xj+1 = −h, xj−1 − xj = −h the above simplifies to
x=x
Z j
(3a + 2bh) h2
(−3a + bh) h2 (−3a + 2bh) h2
1
qj−1 + qj − (x − xj−1 ) f (x) dx + qj
6h2 6h2 h 6h2
x=xj−1
x=x
Z j+1
(3a + bh) h2
1
+ qj+1 − (xj+1 − x) f (x) dx = 0
6h2 h
x=xj
Combining gives
x=x
Z j x=x
Z j+1
−3a + bh 2bh 3a + bh 1
qj−1 +qj +qj+1 − (x − xj−1 ) f (x) dx + (xj+1 − x) f (x) dx = 0
6 3 6 h
x=xj−1 x=xj
The above gives the expression for any row in the stiffness matrix other than the first and the last
21
The above shows that the stiffness matrix can be build quickly without the use of any numerical
integration in this example. Integration is needed to evaluate the force (or load) vector. Depending
on the forcing function, this can be simple or difficult to perform.
Once the load vector is calculated, the unknowns qi are solved for. But before doing that, q1 is
first replaced by the initial condition given in the problem
x=x
R 2
1
h (x2 − x) f (x
x=x1
x=x
R 2 x=x
R 3
1
h (x − x1 ) f (x) dx + (x
2bh−3a bh+3a
0 0 0 ... 0 x=x1 x=x2
6 6
u0
−3a+bh 2bh 3a+bh
0 0 ... 0
x=x x=x
q2 R 3 R 4
6 3 6 1
h (x − x2 ) f (x) dx + (x
0 −3a+bh 2bh 3a+bh
0 ... 0 q 3
x=x2 x=x3
6 3 6
.
..
0 0 −3a+bh 2bh 3a+bh
... 0 ..
6 3 6 .
=
...
.. .. .. .. .. .. .. x=x
R j x=x
R j+1
. . . . . . . 1
. h (x − xj−1 ) f (x) dx + (
−3a+bh 2bh 3a+bh
0 ..
x=xj−1 x=xj
0 0 0 6 3 6
q ..
−3a+bh 2bh 3a+bh N −1
0 0 0 0 6 3 6
.
q
N
bh−3a 3a+2bh x=x
0 0 0 0 0 6 6
1 RN −1 x=x
RN
h (x − xN −2 ) f (x) dx +
x=xN −2 x=xN −
x=x
1
RN
h (x − xN −1 ) f
x=xN −1
22
Now the first row is removed (and remembering to multiply u0 by the first entry in the second
row) gives
x=x
R 2
1
h (x − x1 ) f (x) dx + x
x=x1
x=x x=x
−3a+bh 2bh 3a+bh
1
R 3 R 4
u0 0 0 . . . 0
(x − x ) f (x) dx +
q2
2
6 3 6 h
x=x2 x=x3
−3a+bh 2bh 3a+bh q
0 6 3 6 0 . . . 0 3
..
..
0 0 −3a+bh 2bh 3a+bh
. . . 0 . .
6 3 6
. 1 x=x R j x=xRj
.. .. .. .. .. .. .. . =
. . . . . . . .
h (x − x j−1 ) f (x) dx +
x=xj−1 x=x
−3a+bh 2bh 3a+bh .
0 0 0 0 .
6 3 6 .
..
0 0 0 0 −3a+bh 2bh 3a+bh
q N −1
.
6 3 6
x=x x=
bh−3a 3a+2bh qN N −1
0 0 0 0 0 1 R R
6 6 h (x − x N −2 ) f (x) dx +
x=xN −2 x=x
x=x
1
RN
h (x − xN −1
x=xN −1
The first column is now removed after moving the first entry in the first row to the RHS to become
part of the load vector
x=x
R 2 x=x
R 3
1 (x − x 1 ) f (x) dx + (x3 − x) f (x) dx
h
x=x1 x=x2
x=x x=x
2bh 3a+bh
1
R 3 R 4
3 6 0 0 . . . 0
q
h (x − x 2 ) f (x) dx + (x4 − x)
2
−3a+bh
x=x2 x=x3
2bh 3a+bh
6 3 6 0 ... 0 q3
..
..
0 −3a+bh 2bh 3a+bh
. . . 0
. .
6 3 6
x=x
R j x=xR j+1
... =
. .. .. .. .. .. 1
..
. . . . . h x=x
(x − x j−1 ) f (x) dx + (xj+1 −
x=xj
j−1
−3a+bh 2bh 3a+bh
0 ...
0 0 6 3 6
..
0 0 0 −3a+bh 2bh 3a+bh
q
.
N −1
6 3 6
x=x x=x
bh−3a 3a+2bh qN RN −1 RN
0 0 0 0 6 6
1
h (x − xN −2 ) f (x) dx + (xN −
x=xN −2 x=xN −1
x=x
1
RN
h (x − xN −1 ) f (x)
x=xN −1
Now the system above is solved for q2 , q3 , · · · , qN . Once the qi are found, the solution is
N umber
XN odes
u(x) = qi φ i
i=1
The Appendix shows a Mathematica code which solve the above general first order ODE. 4 first
order ODE solved for illustration and the solution is compared to the exact solution. Animation
23
was made to help illustrate this. the RMS error was calculated. The animations can be access by
clicking on the links below (shows only in the HTML version of this report).
3.2 Example Two. 2nd order ODE, Boundary value problem. Linear
interpolation. Symmetric weak form.
Now the FEM formulation is given for a boundary value, second order ODE with constant
coefficients of the form
d2 u du
a 2 +b + cu(x) = f (x)
dx dx
with the initial conditions u(x0 ) = u0 (called essential Dirichlet condition), and the Neumann
boundary condition du dx = β at L, for x0 ≤ x ≤ L
As before, the solution is assumed to be of the form
N
X
u= qi φ i
i=1
ZL
φj R(x) dx = 0 j = 1· · · N
x0
Where R(x) is the ODE residual obtained by substituting the assumed solution into the original
ODE. Hence the above becomes (For simplicity, j = 1· · · N is not written each time as it is
assumed to be the case)
ZL
d2 u
du
φj a 2 +b + cu(x) − f (x) dx = 0
dx dx
x0
ZL ZL
d2 u
du
φj a 2 dx + φj b + cu(x) − f (x) dx = 0
dx dx
x0 x0
RL 2
Applying integration by parts on φj a ddxu2 dx. Since
x0
d2 u
d du du
a φj = aφj 2 + aφ0j
dx dx dx dx
ZL ZL ZL
d2 u
d du du
a φj dx = aφj 2 dx + aφ0j dx
dx dx dx dx
x0 x0 x0
L ZL ZL
d2 u
du du
aφj = aφj 2 dx + aφ0j dx
dx x0 dx dx
x0 x0
24
Hence
ZL ZL
d2 u du L
du
aφj 2 dx = aφj − aφ0j dx
dx dx x0 dx
x0 x0
L ZL ZL
du du du
aφj − aφ0j dx + φj b + cu(x) − f (x) dx = 0 (A2)
dx x0 dx dx
x0 x0
Hence the symmetric weak form equation (A2) can now be simplified more giving
(
aq1
) ZL ZL
j=1
du du
h
− aφ0j dx + φj b + cu(x) − f (x) dx = 0 j = 1· · · N (A3)
aβ j=N dx dx
x0 x0
The trial function obtain by linear interpolation are used. These are shown in this table
i φi φ0i
(x2 −x) −1
1 h h
n o
(x−xi−1 ) (xi+1 −x) 1 −1
1<i<N h , h h, h
(x−xN −1 ) 1
N h h
25
The global stiffness matrix is now constructed. For the first equation (which corresponds to the
first row in the global stiffness matrix) the result is
ZL ZL
aq1 du du
− aφ01 dx + φ1 b + cu(x) − f (x) dx = 0
h dx dx
x0 x0
N
P
Since u = qi φi , and since φ1 domain of influence is only from x1 to x2 then above becomes
i=1
Zx2 Zx2
aq1 0 d d
− aφ1 (q1 φ1 + q2 φ2 ) dx + φ1 b (q1 φ1 + q2 φ2 ) + c(q1 φ1 + q2 φ2 ) − f (x) dx = 0
h dx dx
x1 x1
The above gives the first row in the stiffness matrix. Since it is assumed that each element will
have the same length, hence x1 − x2 = −h, and the above becomes
Zx2
bh2 ch3 bh2 ch3
q1 q2 1
− 2 −ah + ax1 − + − ax2 + 2 −ax1 + + + ax2 − (x2 − x) f (x) = 0
h 2 3 h 2 6 h
x1
Zx2
a ax1 b ch ax2 ax1 b ch ax2 1
−q1 − + 2 − + − 2 + q2 − 2 + + + 2 − (x2 − x) f (x) = 0
h h 2 3 h h 2 6 h h
x1
For the last equation, which will be the last row in the global stiffness matrix
ZxN ZxN
du du
aβ − aφ0N dx + φN b + cu(x) − f (x) dx = 0
dx dx
xN −1 xN −1
N
P
Since u = qi φi , and since φN domain of influence is only from xN −1 to xN The above becomes
i=1
ZxN ZxN
d d
aβ− aφ0N (qN −1 φN −1 + qN φN ) dx+ φN b (qN −1 φN −1 + qN φN ) + c(qN −1 φN −1 + qN φN ) − f (x) dx
dx dx
xN −1 xN −1
26
ZxN ZxN
a −1 1 (x − xN −1 ) −1 1 (xN − x) (x − xN
aβ− qN −1 + qN dx+ b qN −1 + qN + c qN −1 + qN
h h h h h h h h
xN −1 xN −1
The above represents the last row in the global stiffness matrix. Since it is assumed that each
element will have the same length, hence xN −1 − xN = −h, and the above becomes
ZxN
bh2 ch3 bh2 ch3
qN −1 qN 1
aβ − 2 −axN −1 − + + axN + 2 axN −1 + + − axN − (x − xN −1 ) f (x) = 0
h 2 6 h 2 3 h
xN −1
ZxN
axN −1 b ch axN axN −1 b ch axN 1
aβ − qN −1 − 2
− + + 2 + qN 2
+ + − 2 − (x − xN −1 ) f (x) = 0
h 2 6 h h 2 3 h h
xN −1
The expression for any row in between the first and the last rows is now found. For a general node
j the result is
xZj+1 xZj+1
0 du du
− aφj dx + φj b + cu(x) − f (x) dx = 0
dx dx
xj−1 xj−1
Breaking the integral into halves to make it easier to write the trial functions over the domain of
influence gives
Zxj xZj+1 Zxj xZj+1
0 du 0 du du du
− aφj dx + aφj dx+ φj b + cu(x) − f (x) dx+ φj b + cu(x) − f (x) dx = 0
dx dx dx dx
xj−1 xj xj−1 xj
Zxj Zxj
du du
− aφ0j dx + φj b + cu(x) − f (x) dx
dx dx
xj−1 xj−1
0 0 xj −x 0
Over this range, u = qj−1 φj−1 + qj φj hence dudx = qj−1 φj−1 + qj φj where φj−1 = h , φj−1 =
−1 x−xj−1
h , φj = h , φ0j = h1 hence the above becomes
Zxj Zxj
1 −1 1 x − xj−1 −1 1 xj − x x − xj−1
− a qj−1 + qj dx + b qj−1 + qj + c qj−1 + qj − f (x) dx
h h h h h h h h
xj−1 xj−1
(A4)
27
0 0 xj+1 −x
Over this range, u = qj φj + qj+1 φj+1 hence du
dx = qj φj + qj+1 φj+1 where φj = h , φ0j =
−1 x−xj 0 1
h , φj+1 = h , φj+1 = h hence the above becomes
xZj+1 xZj+1
−a −qj qj+1 xj+1 − x −qj 1 xj+1 − x x − xj
− + dx + b + qj+1 + c qj + qj+1 − f (x) dx
h h h h h h h h
xj xj
Since we are assuming each element will have the same length, hence xj−1 − xj = −h, xj − xj+1 =
−h and the above becomes
axj−1 b ch axj axj−1 2ch axj+1
qj−1 − 2 − + + 2 + qj + −
h 2 6 h h2 3 h2
Zxj xZj+1
−axj b ch axj+1 x − xj−1 xj+1 − x
+ qj+1 + + + − f (x) dx − f (x) dx = 0
h2 2 6 h2 h h
xj−1 xj
The above gives any row in the stiffness matrix other than the first and the last row. Now we can
28
ax1 ax2
− ha + b ch
− ax b ch ax2
h2
− 2 + 3 − h2 h2
1
+ 2 − 6 + h2 0 0
−ax2
− ax b ch ax2 ax1
+ 2ch ax3 b ch ax3
h2
1
− 2 + 6 + h2 h2 3 − h2 h2
+ 2 + 6 + h2
0
axj−1 b ch axj axj−1 axj+1 −axj axj+1
+ 2ch + 2b + ch
0 − h2
− 2 + 6 + h2 h2 3 − h2 h2 6 + h2
..
0 0 . ···
ax axN axN −1 axN
− hN2−1 − 2b + ch b ch
0 0 6 + h2 h2
+ 2 + 3 − h2
Zx2
1
h (x2 − x) f (x)
x1
.
..
..
.
..
.
Zxj xZj+1
=
h1 (x − xj−1 ) f (x) dx + h1 (xj+1 − x) f (x) dx
xj−1 xj
..
.
..
.
ZxN
1
(x − x ) f (x) − aβ
h N −1
xN −1
We must first replace q1 by the initial condition given in the problem, and we must remove the
29
Zx2
1
h (x2 − x) f (x)
x1
..
.
..
.
..
.
Zxj xZj+1
=
1 1
(x − xj−1 ) f (x) dx + h (xj+1 − x) f (x) dx
h
xj−1 xj
..
.
..
.
ZxN
1
(x − x ) f (x) − aβ
h N −1
xN −1
30
Multiply the first element in the second row by u0 and removing the first row gives
−ax2
− ax b ch ax2 ax1
+ 2ch ax3
+ 2b + ch ax3
h2
1
− 2 + 6 + h2
u0 h2 3 − h2 h2 6 + h2 0
−ax3
− ax − 2b + ch ax3 ax2
+ 2ch ax4 b ch ax4
0 6 + h2 3 − h2 + 2 + 6 + h2
2
h2 h2 h2
.. ..
0 0 . .
ax b ch axN
axN −1 b ch axN
0 0 − hN2−1 − 2 + 6 + h2 h2
+ 2 + 3 − h2
q
Zx2 Zx3
1
h (x − x2 ) f (x) dx + h1 (x3 − x) dx
x1 x2
..
.
..
.
Zxj xZj+1
1 (x − xj−1 ) f (x) dx + h1 (xj+1 − x) f (x) dx
= h
xj−1 xj
..
.
..
.
x
Z N
1
h (x − xN −1 ) f (x) − aβ
xN −1
Now moving the first element of the first row above to the RHS and removing the first column
31
gives
q2
ax1 2ch ax3
−ax2 b ch ax3
..
h2
+ 3 − h2 h2
+ 2 + 6 + h2
0 0
.
..
−axj axj+1
− ax b ch ax3 ax2 2ch ax4 b ch
0 h2
2
− 2 + 6 + h2 h2
+ 3 − h2 h2
+ 2 + 6 + h2
.
qj
.. .. ..
. . ··· . ..
.
axN −1 b ch axN
ax N −1 b ch axN
− h2 − 2 + 6 + h2 + 2 + 3 − h2
0 0 h2 qN −1
qN
Zx2 Zx3
1 (x − x ) f (x) dx + 1 (x − x) dx − − ax1 − b + ch + ax2 u
h 2 h 3 h2 2 6 h2 0
x1 x2
..
.
..
.
Zx j x
Z j+1
1 1
=
h (x − xj−1 ) f (x) dx + h (xj+1 − x) f (x) dx
xj−1 xj
..
.
..
.
x
Z N
1
h (x − xN −1 ) f (x) − aβ
xN −1
Now we solve for the vector [q2 , q3 , · · · , qN ]T and this completes our solution.
A Matlab implementation is below which solves any second order ODE. This code was used to
generate an animation. This animation can be accessed by clicking on the link below.
4 References
1. Methods of computer modeling in engineering and the sciences. Volume 1. By Professor
Satya N. Atluri. Tech Science Press.
2. Class lecture notes. MAE 207. Computational methods. UCI. Spring 2006. Instructor: Pro-
fessor SN Atluri.
3. Computational techniques for fluid dynamics, Volume I. C.A.J.Fletcher. Springer-Verlag
Matlab driver file This file sets up the call to do FEM and the animation
function nma_fem_driver
%script to do a simulation of FEM solution of a second order ODE
%by Nasser Abbasi
%This script calls the m file called nma_fem.m
32
%
%this solves the ODE
%
% a u''(t) + b u'(t) + c u(t) = f(t)
%
% with t over the range t0 to len
% and with initial conditions u(0)=u0
% and with u'(len)=beta
%
%To use this to solve an ODE, edit the values below for the variables
% a,b,c,f,len,u0,x0
a=1;
b=1;
c=1;
x0=0; % starting domain point
u0=-10;
len=10;
beta=1;
ymax=4; ymin=-2; % for plotting only. To make the plot looks better
%This below solves the equation exactly to compare the FEM solution
%against.
sol= dsolve('a*D2y+b*Dy+c*y=cos(t)*sin(t)','y(0)=-10','Dy(10)=1');
sol=subs(sol);
directory='matlab_ANIMATION/';
extension='.png';
MAX_NODES = 100;
for n = 3:MAX_NODES
ezplot(sol,[0:len,-ymin:ymax]);
hold on;
image = getframe(gcf);
p=frame2im(image);
file_name=[directory num2str(n) extension] ;
imwrite(p,file_name,'png');
hold off;
% pause(.5)
33
end
end
%%%%%%%%%%%%%%%%%%%%%%%%
%
% edit this below to change the forcing function.
%%%%%%%%%%%%%%%%%%%%%%%%
function v=f(x)
v=sin(x).*cos(x);
end
Matlab function to solve ODE using FEM This file does FEM solution on general second order
ODE
function rmserror=nma_fem(a,b,c,f,x0,u0,beta,len,exact_sol,nNodes)
% Solves numerically a second order ODE with constant coeff.
% using the symmetric Garlekin Finite Elements Methods.
%
% see the file nma_fem_driver.m on how to call this function.
%
% Solve a u''(t) + b u'(t) + c u(t) = f(t)
%
% with t over the range t0 to len
% and with initial conditions u(0)=u0
% and with u'(len)=beta
xc=linspace(x0,len,nNodes);
A = build_stiffness_matrix(nNodes,xc,a,b,c);
load = build_load_vector(a,u0,nNodes,xc,f,beta);
% Now remove the first row and column from the stiffness matrix
A(2,1) = A(2,1)*u0;
load(2) = load(2)-A(2,1);
A = A(2:end,2:end);
load = load(2:end);
34
q = A\load;
q = [u0;q];
plot(xc,y,'ro');
hold on;
line(xc,y);
rmserror = rmserror/NPOINTS;
rmserror = sqrt(rmserror);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = trial(x,x0,len,xc,nNodes,q)
if x<x0||x>len
error('in Trial. x outside range');
end
v = 0;
for i = 1:nNodes
[s,d] = phi(i,x,nNodes,x0,len,xc); %notice ignore d here.
v = v+ q(i)*s;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [v,d]=phi(i,x,nNodes,x0,len,xc)
if(i<1 || i>nNodes)
error('node number outside range');
end
if(x<x0 || x>len)
error('x outside range');
end
35
h = xc(2)-xc(1);
if i==1
if(x>xc(2))
v = 0;
d = 0;
else
v = (xc(2)-x)/h;
d = -1/h;
end
return;
end
if i==nNodes
if(x<xc(nNodes-1))
v = 0;
d = 0;
else
v = (x-xc(nNodes-1))/h;
d = 1/h;
end
return;
end
if(x>xc(i+1) || x<xc(i-1) )
v = 0;
d = 0;
else
if(x<=xc(i))
v = (x-xc(i-1))/h;
d = 1/h;
else
v = (xc(i+1)-x)/h;
d = -1/h;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEE my report for description of this function.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=build_stiffness_matrix(nNodes,xc,a,b,c)
A = zeros(nNodes,nNodes);
h = xc(2)-xc(1);
for i = 1:nNodes
if i==1
A(1,1) = -a/h+ a*xc(1)/h^2 - b/h +c*h/3 - a*xc(2)/h^2;
A(1,2) = -a*xc(1)/h^2 +b/2 -c*h/6 + a*xc(2)/h^2;
else
if i==nNodes
A(nNodes,nNodes-1) = -a*xc(nNodes-1)/h^2 -b/2 +c*h/6 +a*xc(nNodes)/h^2;
A(nNodes,nNodes) = a*xc(nNodes-1)/h^2 +b/2 +c*h/3 -a*xc(nNodes)/h^2;
36
else
A(i,i-1) = -a*xc(i-1)/h^2 -b/2 + c*h/6 + a*xc(i)/h^2;
A(i,i) = a*xc(i-1)/h^2 + 2*c*h/3 -a*xc(i+1)/h^2;
A(i,i+1) = -a*xc(i)/h^2 + b/2 + c*h/6 + a*xc(i+1)/h^2;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Called to integrate the 'force' function
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v=integrand(x,f,xj)
v = (x-xj).*f(x);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SEE my report for more description of this
% This build the load vector.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function load = build_load_vector(a,u0,nNodes,xc,f,beta)
load = zeros(nNodes,1);
h = xc(2)-xc(1);