Finite Difference Method for Ordinary Differential
Equations
What is the finite difference method?
The finite difference method is used to solve ordinary differential equations that have
conditions imposed on the boundary rather than at the initial point. These problems are called
boundary-value problems. In this chapter, we solve second-order ordinary differential
equations of the form
d2y
= f ( x, y, y ' ), a x b , (1)
dx 2
with boundary conditions
y ( a ) = y a and y (b) = yb (2)
Many academics refer to boundary value problems as position-dependent and initial value
problems as time-dependent. That is not necessarily the case as illustrated by the following
examples.
The differential equation that governs the deflection y of a simply supported beam under
uniformly distributed load (Figure 1) is given by
d 2 y qx( L − x )
= (3)
dx 2 2 EI
where
x = location along the beam (in)
E = Young’s modulus of elasticity of the beam (psi)
I = second moment of area (in4)
q = uniform loading intensity (lb/in)
L = length of beam (in)
The conditions imposed to solve the differential equation are
y ( x = 0) = 0 (4)
y ( x = L) = 0
Clearly, these are boundary values and hence the problem is considered a boundary-value
problem.
Figure 1 Simply supported beam with uniform distributed load.
Now consider the case of a cantilevered beam with a uniformly distributed load (Figure 2).
The differential equation that governs the deflection y of the beam is given by
d 2 y q( L − x ) 2
= (5)
dx 2 2 EI
where
x = location along the beam (in)
E = Young’s modulus of elasticity of the beam (psi)
I = second moment of area (in4)
q = uniform loading intensity (lb/in)
L = length of beam (in)
The conditions imposed to solve the differential equation are
y ( x = 0) = 0 (6)
dy
( x = 0) = 0
dx
Clearly, these are initial values and hence the problem needs to be considered as an initial value
problem.
Figure 2 Cantilevered beam with a uniformly distributed load.
2
Example 1
The deflection y in a simply supported beam with a uniform load q and a tensile axial load
T is given by
d 2 y Ty qx( L − x)
− = (E1.1)
dx 2 EI 2 EI
where
x = location along the beam (in)
T = tension applied (lbs)
E = Young’s modulus of elasticity of the beam (psi)
I = second moment of area (in4)
q = uniform loading intensity (lb/in)
L = length of beam (in)
Figure 3 Simply supported beam for Example 1.
Given,
T = 7200 lbs, q = 5400 lbs/in, L = 75 in , E = 30 Msi , and I = 120 in 4 ,
a) Find the deflection of the beam at x = 50" . Use a step size of x = 25" and approximate
the derivatives by central divided difference approximation.
b) Find the relative true error in the calculation of y (50) .
Solution
a) Substituting the given values,
d2y 7200 y (5400) x(75 − x )
− =
dx 2
(30 10 )(120) 2(30 106 )(120)
6
d2y
2
− 2 10 −6 y = 7.5 10 −7 x(75 − x) (E1.2)
dx
d2y
Approximating the derivative at node i by the central divided difference approximation,
dx 2
Figure 4 Illustration of finite difference nodes using
central divided difference method.
d 2 y yi +1 − 2 yi + yi −1
(E1.3)
dx 2 ( x ) 2
We can rewrite the equation as
yi +1 − 2 yi + yi −1
− 2 10 −6 yi = 7.5 10 −7 xi (75 − xi ) (E1.4)
(x) 2
Since x = 25 , we have 4 nodes as given in Figure 3
Figure 5 Finite difference method from x = 0 to x = 75 with x = 25 .
The location of the 4 nodes then is
x0 = 0
x1 = x0 + x = 0 + 25 = 25
x 2 = x1 + x = 25 + 25 = 50
x3 = x 2 + x = 50 + 25 = 75
Writing the equation at each node, we get
Node 1: From the simply supported boundary condition at x = 0 , we obtain
y1 = 0 (E1.5)
Node 2: Rewriting equation (E1.4) for node 2 gives
y3 − 2 y 2 + y1
2
− 2 10 −6 y 2 = 7.5 10 −7 x 2 (75 − x 2 )
(25)
0.0016 y1 − 0.003202 y 2 + 0.0016 y 3 = 7.5 10 −7 (25)(75 − 25)
0.0016 y1 − 0.003202 y 2 + 0.0016 y 3 = 9.375 10 −4 (E1.6)
Node 3: Rewriting equation (E1.4) for node 3 gives
y 4 − 2 y3 + y 2
2
− 2 10 −6 y3 = 7.5 10 −7 x3 (75 − x3 )
(25)
0.0016 y 2 − 0.003202 y 3 + 0.0016 y 4 = 7.5 10 −7 (50)(75 − 50)
0.0016 y 2 − 0.003202 y 3 + 0.0016 y 4 = 9.375 10 −4 (E1.7)
Node 4: From the simply supported boundary condition at x = 75 , we obtain
y4 = 0 (E1.8)
4
Equations (E1.5-E1.8) are 4 simultaneous equations with 4 unknowns and can be written in
matrix form as
1 0 0 0 y1 0
0.0016 − 0.003202 y −4
0.0016 0 2 9.375 10
=
0 0.0016 − 0.003202 0.0016 y 3 9.375 10 − 4
0 0 0 1 y 4 0
The above equations have a coefficient matrix that is tridiagonal (we can use Thomas’
algorithm to solve the equations) and is also strictly diagonally dominant (convergence is
guaranteed if we use iterative methods such as the Gauss-Siedel method). Solving the
equations we get,
y1 0
y
2 = − 0.5852
y 3 − 0.5852
y 4 0
y (50) = y ( x 2 ) y 2 = −0.5852"
The exact solution of the ordinary differential equation is derived as follows. The
homogeneous part of the solution is given by solving the characteristic equation
m 2 − 2 10 −6 = 0
m = 0.0014142
Therefore,
y h = K 1e 0.0014142x + K 2 e −0.0014142x
The particular part of the solution is given by
y p = Ax 2 + Bx + C
Substituting the differential equation (E1.2) gives
d 2 yp
2
− 2 10 −6 y p = 7.5 10 −7 x(75 − x)
dx
d2
2
( Ax 2 + Bx + C ) − 2 10 −6 ( Ax 2 + Bx + C ) = 7.5 10 −7 x(75 − x)
dx
2 A − 2 10 −6 ( Ax 2 + Bx + C ) = 7.5 10 −7 x(75 − x)
− 2 10 −6 Ax 2 − 2 10 −6 Bx + (2 A − 2 10 −6 C ) = 5.625 10 −5 x − 7.5 10 −7 x 2
Equating terms gives
− 2 10 −6 A = −7.5 10 −7
− 2 10 −6 B = −5.625 10 −5
2 A − 2 10 −6 C = 0
Solving the above equation gives
A = 0.375
B = −28.125
C = 3.75 10 5
The particular solution then is
y p = 0.375x 2 − 28.125x + 3.75 105
The complete solution is then given by
y = 0.375x 2 − 28.125x + 3.75 10 5 + K 1e 0.0014142x + K 2 e −0.0014142x
Applying the following boundary conditions
y ( x = 0) = 0
y ( x = 75) = 0
we obtain the following system of equations
K 1 + K 2 = −3.75 10 5
1.1119K 1 + 0.89937K 2 = −3.75 10 5
These equations are represented in matrix form by
1 1 K1 − 3.75 10 5
1.1119 0.89937 K = 5
2 − 3.75 10
A number of different numerical methods may be utilized to solve this system of equations
such as the Gaussian elimination. Using any of these methods yields
K1 − 1.775656226 10 5
K = 5
2 − 1.974343774 10
Substituting these values back into the equation gives
y = 0.375x 2 − 28.125x + 3.75 10 5 − 1.775656266 10 5 e 0.0014142 x − 1.974343774 10 5 e −0.0014142 x
Unlike other examples in this chapter and in the book, the above expression for the deflection
of the beam is displayed with a larger number of significant digits. This is done to minimize
the round-off error because the above expression involves subtraction of large numbers that
are close to each other.
b) To calculate the relative true error, we must first calculate the value of the exact solution at
y = 50 .
y (50) = 0.375(50) 2 − 28.125(50) + 3.75 10 5 − 1.775656266 10 5 e 0.0014142(50)
− 1.974343774 10 5 e −0.0014142( 50)
y (50) = −0.5320
The true error is given by
E t = Exact Value – Approximate Value
Et = −0.5320 − (−0.5852)
E t = 0.05320
The relative true error is given by
True Error
t = 100%
True Value
0.05320
t = 100%
− 0.5320
t = −10%
6
Example 2
Take the case of a pressure vessel that is being tested in the laboratory to check its ability to
withstand pressure. For a thick pressure vessel of inner radius a and outer radius b , the
differential equation for the radial displacement u of a point along the thickness is given by
d 2 u 1 du u
+ − =0 (E2.3)
dr 2 r dr r 2
The inner radius a = 5 and the outer radius b = 8 , and the material of the pressure vessel is
ASTM A36 steel. The yield strength of this type of steel is 36 ksi. Two strain gages that are
bonded tangentially at the inner and the outer radius measure normal tangential strain as
t / r =a = 0.00077462
t / r =b = 0.00038462 (E2.4a,b)
at the maximum needed pressure. Since the radial displacement and tangential strain are related
simply by
u
t = , (E2.5)
r
then
u r =a = 0.00077462 5 = 0.0038731' '
u r =b = 0.00038462 8 = 0.0030769' '
The maximum normal stress in the pressure vessel is at the inner radius r = a and is given by
E u du
max =
2
+ (E2.7)
1 − r r =a dr r = a
where
E = Young’s modulus of steel (E= 30 Msi)
= Poisson’s ratio ( = 0.3)
The factor of safety, FS is given by
Yield strength of steel
FS = (E2.8)
max
a) Divide the radial thickness of the pressure vessel into 6 equidistant nodes, and find the
radial displacement profile
b) Find the maximum normal stress and factor of safety as given by equation (E2.8)
c) Find the exact value of the maximum normal stress as given by equation (E2.8) if it is
given that the exact expression for radial displacement is of the form
C
u = C1 r + 2 .
r
Calculate the relative true error.
Solution
b
i+1
i-1
0…… ……n
a i+1 b
i-1 i
Figure 4 Nodes along the radial direction.
a) The radial locations from r = a to r = b are divided into n equally spaced segments, and
hence resulting in n + 1 nodes. This will allow us to find the dependent variable u numerically
at these nodes.
At node i along the radial thickness of the pressure vessel,
d 2 u u i +1 − 2u i + u i −1
(E2.9)
dr 2 (r )2
du u i +1 − u i
(E2.10)
dr r
Such substitutions will convert the ordinary differential equation into a linear equation (but
with more than one unknown). By writing the resulting linear equation at different points at
which the ordinary differential equation is valid, we get simultaneous linear equations that can
be solved by using techniques such as Gaussian elimination, the Gauss-Siedel method, etc.
Substituting these approximations from Equations (E2.9) and (E2.10) in Equation (E2.3)
u i +1 − 2u i + u i −1 1 u i +1 − u i u i
+ − 2 =0 (E2.11)
(r )2 ri r ri
1 1
+ u + − 2 − 1 − 1 u i + 1 u i −1 = 0 (E2.12)
(r )2 r r i +1 (r )2 r r r 2 (r )2
i i i
Let us break the thickness, b − a , of the pressure vessel into n + 1 nodes, that is r = a is node
i = 0 and r = b is node i = n . That means we have n + 1 unknowns.
We can write the above equation for nodes 1,..., n − 1 . This will give us n − 1 equations. At the
edge nodes, i = 0 and i = n , we use the boundary conditions of
u0 = u r = a
8
u n = u r =b
This gives a total of n + 1 equations. So we have n + 1 unknowns and n + 1 linear equations.
These can be solved by any of the numerical methods used for solving simultaneous linear
equations.
We have been asked to do the calculations for n = 5, that is a total of 6 nodes. This
gives
b−a
r =
n
8−5
=
5
= 0.6 "
At node i = 0, r0 = a = 5" , u 0 = 0.0038731" (E2.13)
At node i = 1, r1 = r0 + r = 5 + 0.6 = 5.6" (E2.14)
1 2 1 1 1 1
u + − − − u + 2 + u = 0
0.6 2 0
0.6
2
(5.6)(0.6) (5.6) 0.6 (5.6)(0.6) 2
2 1
2.7778u 0 − 5.8851u1 + 3.0754u 2 = 0 (E2.15)
At node i = 2, r2 = r1 + r = 5.6 + 0.6 = 6.2"
1 2 1 1 1 1
u + − − − u + 2 + u = 0
0.6 2 1
0.6
2 2 2
(6.2)(0.6) 6.2 0.6 (6.2)(0.6) 3
2.7778u1 − 5.8504u 2 +3.0466u 3 = 0 (E2.16)
At node i = 3, r3 = r2 + r = 6.2 + 0.6 = 6.8"
1 2 1 1 1 1
u + − − − u + 2 +
2 3
u = 0
0.6 2 2
0.6 (6.8)(0.6) 6.8
2
0.6 (6.8)(0.6) 4
2.7778u 2 − 5.8223u 3 + 3.0229u 4 = 0 (E2.17)
At node i = 4, r4 = r3 + r = 6.8 + 0.6 = 7.4 ″
1 2 1 1 1 1
u + − − − u + 2 + u = 0
0.6 2 3 2 2 4
0.6 (7.4)(0.6) (7.4) 0.6 (7.4)(0.6) 5
2.7778u 3 − 5.7990u 4 + 3.0030u 5 = 0 (E2.18)
At node i = 5, r5 = r4 + r = 7.4 + 0.6 = 8 ″
u5 = u r =b = 0.0030769 ″ (E2.19)
Writing Equation (E2.13) to (E2.19) in matrix form gives
1 0 0 0 0 0 u 0 0.0038731
2.7778 − 5.8851 3.0754 0 0 0 u 0
1
0 2.7778 − 5.8504 3.0466 0 0 u 2 0
=
0 0 2.7778 − 5.8223 3.0229 0 u 3 0
0 0 0 2.7778 − 5.7990 3.0030 u 4 0
0 0 0 0 0 1 u 5 0.0030769
The above equations are a tri-diagonal system of equations and special algorithms such as
Thomas’ algorithm can be used to solve such a system of equations.
u 0 = 0.0038731″
u1 = 0.0036165″
u 2 = 0.0034222 ″
u 3 = 0.0032743 ″
u 4 = 0.0031618 ″
u 5 = 0.0030769 ″
b) To find the maximum stress, it is given by Equation (E2.7) as
E u du
max =
2
+
1 − r r =a dr r = a
E = 30 10 6 psi
= 0.3
u r =a = u 0 = 0.0038731″
du u1 − u 0
r =a
dr r
0.0036165 − 0.0038731
=
0.6
= −0.00042767
The maximum stress in the pressure vessel then is
30 106 0.0038731
max = 2
+ 0.3(− 0.00042767)
1 − 0.3 5
= 2.1307 10 psi
4
So the factor of safety FS from Equation (E2.8) is
36 103
FS = = 1.6896
2.1307 10 4
c) The differential equation has an exact solution and is given by the form
C
u = C1 r + 2 (E2.20)
r
where C1 and C 2 are found by using the boundary conditions at r = a and r = b .
C
u (r = a ) = u (r = 5) = 0.0038731 = C1 (5) + 2
5
C
u (r = b) = u (r = 8) = 0.0030769 = C1 (8) + 2
8
giving
C1 = 0.00013462
C2 = 0.016000
Thus
10
0.016000
u = 0.00013462r + (E2.21)
r
du 0.016000
= 0.00013462 − (E2.22)
dr r2
E u du
max = +
1 − 2 r r = a dr r = a
0.00013462(5) +
0.01600
30 10 6 5 0.016000
= + 0.3 0.0013462 −
1 − 0.3 2 5 52
= 2.0538 10 psi4
The true error is
Et = 2.0538 10 4 − 2.1307 10 4
= −7.6859 102
The absolute relative true error is
2.0538 10 4 − 2.1307 10 4
t = 100
2.0538 10 4
= 3.744 %
Example 3
The approximation in Example 2
du u i +1 − u i
dr r
is first order accurate, that is , the true error is of O(r ) .
The approximation
d 2 u u i +1 − 2u i + u i −1
(E3.1)
dr 2 (r )2
(
is second order accurate, that is , the true error is O (r )
2
)
Mixing these two approximations will result in the order of accuracy of O(r ) and
( )
O (r ) , that is O(r ) .
2
So it is better to approximate
du u i +1 − u i −1
(E3.2)
dr 2(r )
because this equation is second order accurate. Repeat Example 2 with the more accurate
approximations.
Solution
a) Repeating the problem with this approximation, at node i in the pressure vessel,
d 2 u u i +1 − 2u i + u i −1
(E3.3)
dr 2 (r ) 2
du u i +1 − u i −1
(E3.4)
dr 2r
Substituting Equations (E3.3) and (E3.4) in Equation (E2.3) gives
u i +1 − 2u i + u i −1 1 u i +1 − u i −1 u i
+ − 2 =0
(r )2 ri 2(r ) ri
1 1
− + u + − 2 − 1 u i + 1 + 1 u i +1 = 0 (E3.5)
2r (r ) (r )2 i −1 (r )2 r 2 (r )2 2r r
i i i
At node i = 0, r0 = a = 5 "
u 0 = 0.0038731" (E3.6)
At node i = 1, r1 = r0 + r = 5 + 0.6 = 5.6"
1 1 2 1 1 1
− u1 + 2 + u 2 = 0
2(5.6)(0.6) + (0.6)2 u 0 + − (0.6)2 − (5.6)2 2(5.6)(0.6)
0.6
2.6297u 0 − 5.5874u1 + 2.9266u 2 = 0 (E3.7)
At node i = 2, r2 = r1 + r = 5.6 + 0.6 = 6.2 "
1 1 2 1 1 1
− + u + −
2 1
− u + 2 +
2 2
u 3 = 0 (E3.8)
2(6.2)(0.6) 0.6 0.6
2
6.2 0.6 2(6.2)(0.6)
2.6434u1 − 5.5816u 2 +2.9122u3 = 0
At node i = 3, r3 = r2 + r = 6.2 + 0.6 = 6.8 "
1 1 2 1 1 1
− + u + −
2 2
− u + 2 +
2 3
u 4 = 0 (E3.9)
2(6.8)(0.6) 0.6 0.6
2
6.8 0.6 2(6.8)(0.6)
2.6552u 2 − 5.5772u 3 + 2.9003u 4 = 0
At node i = 4, r4 = r3 + r = 6.8 + 0.6 = 7.4 "
1 1 2 1 1 1
− + u 3 + − − u 4 + 2 + u 5 = 0 (E3.10)
2(7.4)(0.6) 0.6
2
0.6
2
(7.4)2
0.6 2(7.4)(0.6)
2.6651u3 − 5.5738u 4 + 2.8903u5 = 0
At node i = 5, r5 = r4 + r = 7.4 + 0.6 = 8 "
u 5 = u / r =b = 0.0030769 " (E3.11)
Writing Equations (E3.6) thru (E3.11) in matrix form gives
1 0 0 0 0 0 u 0 0.0038731
2.6297 − 5.5874 2.9266 0 0 0 u 0
1
0 2.6434 − 5.5816 2.9122 0 0 u 2 0
=
0 0 2.6552 − 5.5772 2.9003 0 u 3 0
0 0 0 2.6651 − 5.5738 2.8903 u 4 0
0 0 0 0 0 1 u 5 0.0030769
12
The above equations are a tri-diagonal system of equations and special algorithms such as
Thomas’ algorithm can be used to solve such equations.
u 0 = 0.0038731"
u1 = 0.0036115"
u 2 = 0.0034159 "
u 3 = 0.0032689 "
u 4 = 0.0031586 "
u 5 = 0.0030769 "
du − 3u 0 + 4u1 − u 2
b)
dr r =a 2(r )
− 3 0.0038731+ 4 0.0036115 − 0.0034159
=
2(0.6)
−4
= −4.925 10
30 10 6 0.0038731
+ 0.3(− 4.925 10 −4 )
max = 2
1 − 0.3 5
= 2.0666 10 psi
4
Therefore, the factor of safety FS is
36 10 3
FS =
2.0666 10 4
= 1.7420
c) The true error in calculating the maximum stress is
Et = 2.0538 10 4 − 2.0666 10 4
= −128 psi
The relative true error in calculating the maximum stress is
− 128
t = 100
2.0538 10 4
= 0.62323%
Table 1 Comparisons of radial displacements from two methods.
r u exact u1st order t u 2nd order t
5 0.0038731 0.0038731 0.0000 0.0038731 0.0000
5.6 0.0036110 0.0036165 1.5160 10 −1 0.0036115 1.4540 10 −2
6.2 0.0034152 0.0034222 2.0260 10 −1 0.0034159 1.8765 10 −2
6.8 0.0032683 0.0032743 1.8157 10 −1 0.0032689 1.6334 10 −2
7.4 0.0031583 0.0031618 1.0903 10 −1 0.0031586 9.5665 10 −3
8 0.0030769 0.0030769 0.0000 0.0030769 0.0000
ORDINARY DIFFERENTIAL EQUATIONS
Topic Finite Difference Methods of Solving Ordinary Differential Equations
Summary Textbook notes of Finite Difference Methods of solving ordinary
differential equations
Major General Engineering
Authors Autar Kaw, Cuong Nguyen, Luke Snyder
Date February 28, 2021
Web Site http://numericalmethods.eng.usf.edu
14
