ODE Finite Element Method

Download as pdf or txt
Download as pdf or txt
You are on page 1of 36

Examples solving an ODE using finite elements method.

Mathematica and Matlab implementation and animation


Nasser M. Abbasi

Sept 28,2006 Compiled on September 9, 2023 at 3:03pm

Contents
1 Introduction 1

2 Weighted Residual method. Global trial functions. 4


2.1 First example. First order ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Second example. 4th order ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3 Finite element method 12


3.1 Example one. First order ODE, linear interpolation . . . . . . . . . . . . . . . . . . 12
3.2 Example Two. 2nd order ODE, Boundary value problem. Linear interpolation.
Symmetric weak form. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

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

equations can be represented as an operator L


dy
L := (y, x) → a + by(x) − f (x) (1)
dx
If we know the exact solution, call it ye (x) then we write

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

2 Weighted Residual method. Global trial functions.


The best way to learn how to use WRM is by working over and programming some examples.
We analyze the solution in terms of errors and the effect of changing N on the result.

2.1 First example. First order ODE


Given the following ODE
dy
− y(x) = 0
dx
defined over 0 ≤ x ≤ 1 with the boundary condition y(0) = 1, we wish to solve this numerically
using the WRM. This ODE has an exact solution of y = ex .
The solution using WRM will always follow these steps.
STEP 1
Assume a solution that is valid over the domain 0 ≤ x ≤ 1 to be a series solution of trial (basis)
functions. We start by selecting a trial functions. The assumed solution takes the form of
N
X
ỹ(x) = ỹ0 + qj φj (x)
j=1

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

Hence from (1), we find the residual to be


   
N N
d X X
R(x) = 1+ qj xj  − 1 + qj x j 
dx
j=1 j=1
 
XN XN
= qj jxj−1 − 1 + qj x j 
j=1 j=1
N
X
qj jxj−1 − xj

= −1 + (2)
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 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

vi (x) = xi−1 (4)

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
  

Substituting the above in (3) gives


x=1
Z
xi−1 R(x) dx = 0 i = 1· · · N
x=0
x=1
Z
xi−1 −1 + q 1 − x1 + q2 2x − x2 + q3 3x2 − x3 dx = 0
  
i = 1· · · 3
x=0
6

The above generates N = 3 equations to solve. They are


x=1
Z
−1 + q1 1 − x1 + q2 2x − x2 + q3 3x2 − x3 dx = 0
  
i=1
x=0
x=1
Z
x −1 + q1 1 − x1 + q2 2x − x2 + q3 3x2 − x3 dx = 0
  
i=2
x=0
x=1
Z
x2 −1 + q1 1 − x1 + q2 2x − x2 + q3 3x2 − x3 dx = 0
  
i=3
x=0

Performing the integration above, gives the following 3 equations


     
1 2 3
−1 − q1 − − q2 − − q3 − =0
2 3 4
     
1 1 5 11
− − q1 − − q2 − − q3 − =0
2 6 12 20
     
1 1 3 13
− − q1 − − q2 − − q3 − =0
3 12 10 30

Which can be written in matrix form as


   
1 2 3
1
 
2 3 4  q1
1 5 11  q 
1
6 12 20   2 =
2

q3 1
 
1 3 13
12 10 30 3

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

Figure 2: Code used for N = 5


8

Figure 3: Result for N = 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

Figure 4: Error using N = 5

2.2 Second example. 4th order ODE


Now we will use a more complex example and repeat the above steps. We now want to numerically
solve the following
d4 y
+ y(x) = 1
dx4
defined over 0 ≤ x ≤ 1 with the boundary conditions y(0) = 0, y(1) = 0, y 0 (0) , y 0 (1) = 0.
This problem is taken from Professor S.N.Atluri text book ’Methods of computer modeling in
engineering and the sciences’ Volume 1, page 47-50. Professor Atluri used a trigonometric functions
for the trial function
N
X
ỹ = qi sin ((2i − 1) πx)
i=1

Which already satisfies the boundary conditions. For the test function, the same function as above
is used hence the test (weight) function is

vj (x) = sin ((2j − 1) πx)

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

Selecting the trial solution.


N
X
ỹ(x) = ỹ0 + qj φj (x)
j=1
ỹ0 = 0 as this will satisfy the boundary conditions. Hence the trial solution is
N
X
ỹ(x) = qj φj (x)
j=1

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

Hence (3) becomes


Z
vi (x) R(x) dx = 0 i

Z1
(sin (2i − 1) πx) R(x) dx = 0 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 (2i − 1) πx) i
0

The above generates N equations to solve for the coefficients qi

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

3 Finite element method


3.1 Example one. First order ODE, linear interpolation
Let us first summarize what we have done so far. Given a differential equation defined over domain
N
X
Ω, we assume its solution to be of the form ũ(x) = qj φj (x).
j=1

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

element 1 element 2 element 3


X=1
node 1 node 2 node 3 node 4
x1 x2 x3 x4
Here we assume linear interpolation between nodal
values of the solution

Figure 6: Finite elements with linear interpolation

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

element 1 element 2 element 3


X=1
node 1 node 2 node 3 node 4
x1 x2 x3 x4
Here we assume linear interpolation between nodal
values of the solution

Figure 7: Finite elements with linear interpolation second diagram

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

N11 (x) N21 (x)


z }| { z }| {
(x 2 −x) (x−x1 )
u1 (x) = q1 + q2
h h
= q1 N11 (x) + q2 N21 (x)
15

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

The above is valid for x2 ≤ x ≤ x3 . The global trial function is

u(x) = u1 + u2 + u3
= q1 N11 + q2 N21 + q2 N12 + q3 N22 + q3 N13 + q4 N23
  

= q1 N11 + q2 N21 + N12 + q3 N22 + N13 + q4 N23


  

The shape function for node 1 is

φ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

The approximate solution is therefore


N umber
XN odes
u(x) = qi φ i
i=1
This completes the first part, which was to express the global approximate solution as sum of
basis functions, each multiplied by an unknowns q coefficients.
The diagram below illustrates the above.

u 2  q 2N 21x  q 3N 22x

u 1  q 1N 11x  q 2N 12x u 3  q 3N 31  q 4N 32
y
q3
q1 x
q2
q4

element 1 element 2 element 3


X=1
node 1 node 2 node 3 node 4
x1 x2 x3 x4

Figure 8: Showing solutions at each element

The diagram below illustrates the numbering used.

u1 u2 uj u N1
u j1
q3 qN
q1
q2 qj
q j1 q N1
  element N-1 X=1
X=0 element 1 element 2 element j-2 element j

node 1 node 2 node 3 node j-1 node j node (N-1) node N


x1 x2 x3 x j1 xj x N1 xN
Numbering used for N number of nodes and N-1 number of elements

Figure 9: Numbering used for interpolation and elements


17

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

The test functions are


vj = φj j = 1· · · N
The weak form of the differential equation is
x=x
Z N
φj (x) R(x) dx = 0 j = 1· · · N
x=x1
x=x
Z N " N
# !
X
qi aφ0i + bφi − f (t) dx = 0

φj (x) j = 1· · · N
x=x1 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

In general, if there are N nodes, then for the first node


(x2 − x)
φ1 =
h
−1
φ01 =
h
And for the last node
(x − xN −1 )
φN =
h
1
φ0N =
h
18

And for any node in the middle


(x − xi−1 )
φi =
h
1
φ0i =
h
for xi−1 ≤ x ≤ xi and
(xi+1 − x)
φi =
h
−1
φ0i =
h
for xi ≤ x ≤ xi+1 .
Looking back at the weak form integral above, it is evaluated as follows
x=x
Z N "N # !
X
0

φj qi aφi + bφi − f (x) dx = 0 j = 1· · · N
x=x1 i=1

For the first node only, j = 1, the following results


x=x
Z N "N # !
X
qi aφ0i + bφi − f (x) dx = 0

φ1
x=x1 i=1

(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 simplifies further to


! x=x
Z 2
(3a + 2b(x1 − x2 )) (x1 − x2 )2 (−3a + b(x1 − x2 )) (x1 − x2 )2
−q1 − q2 − φ1 f (x) = 0
6h2 6h2
x=x1
! x=x
Z 2
(3a + 2b(x1 − x2 )) (x1 − x2 )2 (−3a + b(x1 − x2 )) (x1 − x2 )2 1
−q1 − q2 − (x2 − x) f (x) = 0
6h2 6h 2 h
x=x1

Since x2 − x1 = h and x1 − x2 = −h the above reduces to


x=x
Z 2
(3a − 2bh) h2 (−3a − bh) h2
 
1
−q1 − q2 − (x2 − x) f (x) = 0
6h2 6h2 h
x=x1
    x=x
Z 2
2bh − 3a bh + 3a 1
q1 + q2 − (x2 − x) f (x) = 0
6 6 h
x=x1
19

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

Since φN domain of influence is xN −1 ≤ x ≤ xN , the above simplifies to


x=x
Z N
φN qN −1 aφ0N −1 + bφN −1 + qN aφ0N + bφN − f (x) dx = 0
   

x=xN −1

(xN −x) −1 x−xN −1


Since φN −1 = h , φ0N −1 = h , φN = h , φ0N = 1
h the above becomes
x=x
Z N      
x − xN −1 −a (xN − x) 1 x − xN −1
qN −1 +b + qN a +b − f (x) dx = 0
h h h h h
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

Letting xN −1 − xN = −h in the above becomes


  x=x
Z N
bh − 3a (3a + 2bh) 1
qN −1 + qN − (x − xN −1 ) f (x) = 0
6 6 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

row. Hence the global stiffness matrix is



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
 q1 
 −3a+bh 2bh 3a+bh
0 0 ... 0   q2  
  
 x=x
R 3 x=x
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 (x − xj−1 ) f (x) dx + (
 . 
 h
−3a+bh 2bh 3a+bh
0   ..  
 x=xj−1 x=xj
 0 0 0 6 3 6


 
..
 

0 0 0 0 −3a+bh 2bh 3a+bh  q
 N −1 .
  
6 3 6
 
qN
  
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

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

And the unknowns qi are found by solving N algebraic equations

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

then integrating both sides of the above gives

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

Substituting the above in equation A1 gives

 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

Considering the term


 L
du du du
aφj = aφj − aφj
dx x0 dx x=L dx x=x0

First considering the term aφj du du


dx x=L . Since dx = β at x = L then this term becomes aβ[φj ]x=L .
th
But [φj ]x=L is non zero only for the N shape function evaluated at x = L. Since linear inter-
x−xN −1
polation is used, the N th shape function is φN = h which have the value of 1 at x = xN .
Hence
du
aφj = aβ whenj = N, otherwise,0
dx x=L
Now considering the term aφj dudx x=x0 , since at x = x0 all shape functions will be zero except for
φ1 which has the value of 1 at x = x0 . Hence this simplifies to a du
dx x=x0
N N
du d P
qi φ0i . But at x = x0 only φ1 is defined and its has the slope of
P
Recalling that dx = dx qi φ i =
i=1 i=1
−1
h , hence
du −aq1
a = when j = 1, otherwise 0
dx x=x0 h
However, q1 = u0 since that is by definition the initial condition. Finally one obtains
( )
aq1
du L j
 
h = 1
aφj =
dx x0 aβ j=N

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

But φ1 = (x2h−x) and φ01 = −1


h and over the domain from x1 to x2 , φ2 = (x−x1 )
h , φ02 = h1 , hence the
above becomes
Zx2   Zx2      
aq1 −a −q1 q2 (x2 − x) −q1 q2 (x2 − x) (x − x1 )
− + dx+ b + + c q1 + q2 − f (x) dx = 0
h h h h h h h h h
x1 x1

The above simplifies to


! Zx2
b(x1 − x2 )2 c(x1 − x2 )3 b(x1 − x2 )2 c(x1 − x2 )3
 
aq1 q1 q2 1
− 2 ax1 − − − ax2 + 2 −ax1 + − + ax2 − (
h h 2 3 h 2 6 h
x1
! ! Zx2
q1 b(x1 − x2 )2 c(x1 − x2 )3 q2 b(x1 − x2 )2 c(x1 − x2 )3 1
− 2 −ah + ax1 − − − ax2 + 2 −ax1 + − + ax2 − (
h 2 3 h 2 6 h
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

(x−x ) (xN −x)


But φN = N −1
h and φ0N = 1
h and over the domain from xN −1 to xN , φN −1 = h , φ0N −1 =
−1
h hence the above becomes

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 simplifies to


!
qN −1 b(xN −1 − xN )2 c(xN −1 − xN )3 qN b(xN −1 − xN )2 c(xN −1 − x
aβ− 2 −axN −1 − − + axN + 2 +axN −1 + −
h 2 6 h 2 3

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

Considering the first domain xj−1 ≤ x ≤ xj gives

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

Now considering the second domain xj ≤ x ≤ xj+1 gives


xZj+1 xZj+1  
du du
− aφ0j dx + φj b + cu(x) − f (x) dx
dx dx
xj xj

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

Combine A4 and A5 and simplifying gives


!
qj−1 b(xj−1 − xj )2 c(xj−1 − xj )3 qj b(xj−1 − xj )2 c(xj−1 − xj )3 b(xj −
−axj−1 − − + ax j + ax j−1 + − −
h2 2 6 h2 2 3 2
! Zxj xZj+1
qj+1 b(xj − xj+1 )2 c(xj − xj+1 )3 x − xj−1 xj+1 − x
+ 2 −axj + − + axj+1 − f (x) dx− dx = 0
h 2 6 h h
xj−1 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

write the global stiffness matrix as

 
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

first row after that as follows



 
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
  
1
− + + 3 − h2 + + + 0

 h2 2 6 h2 h2 h2 2 6 h2

 
 .. 
. 
0 0 ···

 
 
 .. 
 0 0 . ··· 
 
ax b ch axN
 axN −1 b ch axN
 
0 0 − hN2−1 − 2 + 6 + h2 h2
+ 2 + 3 − h2 

 
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

close all; clear all;

% equation 1. u''+u'+u=sin(x)*cos(x). u(0)=1, u'(len)=beta, len=10


% maxy=4 miny=-1
% a=1, b=1, c=1, x0=0, u0=1,len=10

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;

r = nma_fem(a,b,c,@f,x0,u0,beta,len,sol,n); %SOLVES by FEM

title(sprintf('Solving y''''+y''+y(x)=sin(x)*cos(x). y(0)=1, y''(10)=1. by FEM. \nN=%d\nRMSER


drawnow;

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

%by Nasser Abbasi. Sept 26,2006.

xc=linspace(x0,len,nNodes);

% This plots the shape functions.


% x=linspace(x0,len,1000);
% y=linspace(x0,len,1000);
% for i=1:nNodes
% for j=1:length(y)
% [v,d]=phi(i,x(j),nNodes,x0,len,xc);
% y(j)=v;
% end
% plot(x,y);
% hold on;
% end

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

% SOLVE for unknowns

q = A\load;
q = [u0;q];

% Plot the solution


y = zeros(length(xc),1);
for i=1:length(xc)
y(i)=trial(xc(i),x0,len,xc,nNodes,q);
end

plot(xc,y,'ro');
hold on;
line(xc,y);

% Calculate RMSerror. Use 50 points. Should be more than enough.


NPOINTS = 50;
x = linspace(x0,len,NPOINTS);
rmserror = 0;
for i = 1:length(x)
y = trial(x(i),x0,len,xc,nNodes,q);
t=x(i); % USED for subs below. do not remove.
yexact = real(double(subs(exact_sol)));
rmserror = rmserror+(y-yexact)^2;
end

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);

load(1) = -(1/h)* quad(@(x)integrand(x,f,xc(2)),xc(1),xc(2)) ;


load(nNodes) = (1/h)* quad(@(x)integrand(x,f,xc(nNodes-1)),xc(nNodes-1),xc(nNodes)) - a*beta;
for i = 2:nNodes-1
load(i) = quad(@(x)integrand(x,f,xc(i-1)),xc(i-1),xc(i)) - ...
quad(@(x)integrand(x,f,xc(i+1)),xc(i),xc(i+1));
load(i) = (1/h)*load(i) ;
end
end
 

You might also like