NUMERICAL ANALYSIS

University of Babylon/ College of Engineering/ Mechanical Engineering Dep.

Lecturer : Dr. Rafel Hekmat Class : 3rd B.Sc

Solution of Differential Equation by Finite

Difference Method
What is a Differential Equation?
An equation that consists of derivatives is called a differential equation. Differential
equations have applications in all areas of science and engineering. Mathematical formulation
of most of the physical and engineering problems leads to differential equations.

Types of Differential Equations

a- Ordinary Differential Equations (ODE)
b- Partial Differential Equations (PDE)

An ordinary differential equation is that in which all the derivatives are with respect to a single
independent variable. Examples of ordinary differential equations include

Ordinary differential equations are classified in terms of order and degree. Order of an ordinary
differential equation is the same as the highest derivative and the degree of an ordinary
differential equation is the power of highest derivative.

Thus the differential equation,

is of order 3 and degree 1, whereas the differential equation

is of order 1 and degree 2.

1
PDEs are differential equations containing two or more independent variables, so the
derivatives are partial derivatives. For example ,given a function of two variables
(x, y) f(x, y), the partial derivatives with respect to x and y are 𝜕𝑓 ⁄𝜕𝑥 and 𝜕𝑓 ⁄𝜕𝑦

PDE example I: Laplace equation

The Laplace equation is a second ordered PDE appearing for example in Fluid Mechanics
(potential flows) and Electromagnetics (electromagnetic field in a charge free region):

where  represents, for example, the velocity potential or the temperature.

PDE example II: Heat equation

The heat equation is a second order PDE describing how temperature T diffuses through a
medium of thermal diffusivity α:

PDE example III: Wave equation

The wave equation describes how the sound pressure p propagates at speed c0 through a
medium (at rest), and plays a role in acoustics, fluid mechanics, and quantum mechanics:

Classification of Partial Differential Equations

For analyzing the equations for fluid flow problems, it is convenient to consider the case of a
second-order differential equation given in the general form as:

In the coefficients A, B, C, D, E and F are either constants or functions of only (x, y) (do not
contain φ or its derivatives), it is said to be a linear equation, otherwise it is a non-linear
equation. If G = 0, the aforesaid equation is homogeneous, otherwise it is non-homogeneous.

2
Again for the above mentioned equation

if B2 − 4AC = 0, the equation is parabolic

if B2 − 4AC < 0, the equation is elliptic
if B2 − 4AC > 0, the equation is hyperbolic

Example:

Classify the following linear second order partial differential equations (PDEs) with solution
u(x,y)in the xy-plane.

3
Example: Consider the one-dimensional damped wave equation 9uxx= utt+ 6ut

Initial and Boundary Conditions

Most PDEs have an infinite number of admissible solutions. Thus the PDE alone is not
sufficient to get a unique solution. Usually some boundary conditions and initial conditions are
required.
For the heat equation the simplest boundary conditions are fixed temperatures at both ends:

Where l is the length of the rod and h1(t) the temperature at the first end and h2(t) the
temperature at the second end.
The initial conditions specify an arbitrary initial temperature distribution inside the rod:

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.

Steps of finite difference solution:

 Divide the solution region into a grid of nodes,
 Approximate the given differential equation by finite difference equivalent,

4
  Solve the differential equations subject to the boundary conditions and/or initial
conditions.

First Derivative Approximations

The derivative of a given function f(x) can be approximated as:

df ( x0 ) f ( x0  x)  f ( x0 )
(1) Forward Difference Formula
dx x
df ( x0 ) f ( x0 )  f ( x0  x)
(2) Backword Difference Formula
dx x
df ( x0 ) f ( x0  x)  f ( x0  x)
(3) Central Difference Formula
dx 2x

Second Derivative Approximation

2 f f ( x0  x)  2 f ( x0 )  f ( x0  x)

x 2
(x)2

Example:
Solve

2 f
 4 f  0 , 0  x  1.0
x 2
Subject to the boundary conditions:
f (0)  0
f (1)  1
Assume x=0.5

Solution: Discretization of the equation:

f (i  1)  2 f (i)  f (i  1)
 4 f (i)  0
(x) 2
f (2)  2 f (1)  f (0)
 4 f (1)  0
(0.5) 2
Since there is one unknown, one equation is enough to find this unknown.

5
 f (1)  f (0)  f (2)
f (1)  1
If x=0.25
f (2)  2 f (1)  f (0)
 4 f (1)  0
(0.25) 2
f (3)  2 f (2)  f (1)
 4 f (2)  0
(0.25) 2
f (4)  2 f (3)  f (2)
 4 f (3)  0
(0.25) 2
f(1), f(2) and f(3) are unknowns. Therefore, three equations are written to find three unknowns.

Example:
Consider the diffusion equation:
  2
k  2
t x
Where k is a constant. Discretized equation:

(i, j  1)  (i, j ) (i  1, j )  2(i, j )  (i  1, j )

k 
t (x)2
Where

x  ix, i  0,1, 2,..., n

t  jt , j  0,1, 2,..., n
 Forward difference for t and central difference for x is used.

t
r
k (x)2
 We can write:
(i, j  1)  r(i  1, j )  (1  2r )(i, j )  r(i 1, j ) `

 This equation can be calculated in terms of boundary and initial conditions.

 This explicit formula can be used to find (x,t+)explicitly in terms of (x,t)

Example:
The deflection in a simply supported beam with a uniform load q and a tensile axial load T is
given by

where
x:Location along the beam (in)
T:Tension applied (lbs)
E:Young’s modulus of elasticity of the beam (psi)
6
I:Second moment of area (in4)
q:Uniform loading intensity (lb/in)
L:Length of beam (in)

Given,

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.
Solution:
Substituting the given values,

𝜕2 𝑦
Approximating the derivative at node i by the central divided difference approximation
𝜕𝑥 2

We can rewrite the equation as:

Since Δx=25 , we have 4 nodes as given in Figure below

Writing the equation at each node, we get

7
Equations above are 4 simultaneous equations with 4 unknowns and can be written in matrix
form as

Solving the equations we get,

Example
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

The inner radius a= 5′′ , and the outer radius b= 8’’ , and the material of the pressure vessel is
ASTM A36 steel
Divide the radial thickness of the pressure vessel into 6 equidistant nodes, and find the radial
displacement profile

At node i along the radial thickness of the pressure vessel

8
Let us break the thickness, b-a , of the pressure vessel into n+1 nodes, that is is node and is
node . That r=a is nodes 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

We have been asked to do the calculations for n=5, that is a total of 6 nodes. This gives

9
Writing the above Equations in matrix form gives

After solving

Explicit versus implicit Finite Difference Schemes

During the lecture we solved the transient (time-dependent) heat equation in 1D

In explicit finite difference schemes, the temperature at time n+1 depends explicitly on the
temperature at time n. The explicit finite difference discretization of above equation is

This can rearranged in the following manner (with all quantities at time n+1 on the left-hand-
side and quantities at time n on the right-hand-side)

𝑛
Since we know𝑇𝑖+1 , 𝑇𝑖𝑛 and 𝑇𝑖−1
𝑛
, we can compute 𝑇𝑖𝑛 . The major advantage of explicit finite
difference methods is that they are relatively simple and computationally fast. However, the
main drawback is that stable solutions are obtained only when

𝜕2 𝑇
In implicit finite difference schemes, the spatial derivatives 2 are evaluated (at least partially)
𝜕𝑥
at the new time step. The simplest implicit discretization of heat equation in 1D is

This can be rearranged so that unknown terms are on the left and known terms are on the right

𝐾∆𝑡
Where 𝑠 = . The main advantage of implicit finite difference methods is that there are
∆𝑥 2

10
no restrictions on the time step, which is good news if we want to simulate geological processes
at high spatial resolution. Taking large time steps, however, may result in an inaccurate solution.

Solving an implicit finite difference scheme

We solve the transient heat equation on the domain –l/2 ≤ x ≤l/2 with the following boundary
conditions

with the initial condition

As usual, the first step is to discretize the spatial domain with nx finite difference points. The
implicit finite difference discretization of the temperature equation is

and the known right-hand-side vector rhs is

11