Subject: Numerical Methods (PGME-2nd sem.

-2021)
Chapter- 4
Application of Finite Difference Method to selected model of
Equations:
1. Parabolic types of Equations: Transient Heat Conduction
Equation
2. Hyperbolic Equation: 1D wave Equation
3. Equilibrium Equation: Laplace Equation

Class No.

Kalyan Kalita
[Pick the date]
Application of Finite Difference Method to selected model of Equations

Parabolic types of Equations:

Parabolic types of equations are associated with domains usually in the time
direction ( or time like direction). Computational (ie. Numerical Schemes) schemes for
such problems are march forward (ie. move forward) step by step,[ ie. After completing the
calculation in one step (say step no. n) calculation for next step (step no. n+1) starts] in the
direction of time(say along y axis) The stability consideration is important for such
numerical scheme.

The Numerical schemes applied for the solution of parabolic equations may be
categorised in to two types:

1. Explicit and
2. Implicit
Explicit Scheme:

An explicit scheme is one for which only one unknown appears in the finite difference
equation (FDE) in a manner which permits evaluation of the in terms of known quantities.
Example: Forward Time Central Space (FTCS) Scheme

Implicit Scheme:

An implicit scheme is one for which more than one unknown appears in the finite
difference equation (FDE) and therefore, need to solve a set of simultaneous equations.
Example: Crank-Nicholson Scheme.

Finite Difference Schemes for Transient Heat conduction Equation:

1. Forward Time Central Space (FTCS) Scheme:
Problem: Find the temperature distribution in a thin long (of unit length) insulated
bar of uniform cross-section.
The initial condition is u(x,0)=sin πx for 0 ≤ 𝑥 ≤ 1
And boundary conditions 𝑢(0, 𝑡) = 0 , 𝑢(1, 𝑡) = 0 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑡 ≥ 0

Problem Description:

The temperature distribution in a thin long insulated bar of uniform cross-section is

𝜕𝑢 𝜕2 𝑢
given as = 𝛼 where u be the temperature, x be the space coordinate, t be the time
𝜕𝑡 𝜕𝑥 2
and 𝛼 be the coefficient of thermal diffusion.
The bar (OA) is placed along x-axis between x = 0 and x = 1 (the bar length have been made
dimensionless with respect to maximum length of the bar).

The initial condition (i.e. the starting temperature at various location of the bar
at starting time t = 0 is give as u(x,0)=sin πx for 0 ≤ 𝑥 ≤ 1

The boundary conditions (i.e. at the end of the bar at any time is given as 𝑢(0, 𝑡) =
0 , 𝑢(1, 𝑡) = 0 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑡 ≥ 0
Type equation here.

𝑓𝑖𝑔. 1 ℎ𝑒𝑎𝑡 𝑐𝑜𝑛𝑑𝑢𝑐𝑡𝑖𝑜𝑛 𝑖𝑛 𝑎 𝑙𝑜𝑛𝑔 𝑏𝑎𝑟 𝑂𝐴

Step 1

Grid Generation/Mesh Generation: To solve the problem the domain 0 ≤ 𝑥 ≤ 1, 𝑡 > 0 is

subdivided into a network by drawing straight lines parallel to coordinate axis x and t, such that

𝑥 = 𝑖∆𝑥 𝑤ℎ𝑒𝑟𝑒 ∆𝑥 𝑏𝑒 𝑡ℎ𝑒 𝑚𝑒𝑠ℎ 𝑠𝑝𝑎𝑐𝑖𝑔 𝑎𝑙𝑜𝑛𝑔 𝑡ℎ𝑒 𝑏𝑎𝑟, & 𝑖 = 1,2,3 … … . 𝑁

𝑡 = 𝑛∆𝑡 𝑤ℎ𝑒𝑟𝑒 ∆𝑡 𝑏𝑒 𝑚𝑒𝑠ℎ 𝑠𝑝𝑎𝑐𝑖𝑔𝑠 𝑎𝑙𝑜𝑔 𝑡 𝑑𝑖𝑟𝑒𝑐𝑡𝑖𝑜𝑛 & 𝑛 = 1,2,3 … … …

In this problem let us divide the bar in to 4 equal division then

𝑙𝑒𝑛𝑔𝑡ℎ 𝑜𝑓 𝑡ℎ𝑒 𝑏𝑎𝑟 𝑥𝑚𝑎𝑥 1 1

∆𝑥 = = =𝑁=
𝑁𝑜𝑠.𝑜𝑓 𝑑𝑖𝑣𝑖𝑠𝑖𝑜𝑛 𝑚𝑎𝑑𝑒 𝑖𝑚𝑎𝑥 4
 𝐹𝑖𝑔. 2: 𝐺𝑟𝑖𝑑 𝑝𝑜𝑖𝑡𝑠 𝑓𝑜𝑟 𝑡ℎ𝑒 𝑝𝑟𝑜𝑏. 𝑜𝑓 ℎ𝑒𝑎𝑡 𝑐𝑜𝑑𝑢𝑐𝑡𝑖𝑜𝑛 𝑖𝑛 𝑎 𝑏𝑎𝑟

Notation used for Parabolic Problems: Say temperature at any point in (x,t) field can be written
𝑛
as u(x,t). Thus 𝑢(𝑥, 𝑡 ) = 𝑢(𝑖, 𝑛) = 𝑢𝑖

Step 2

Formation of Finite Difference Equation:

The governing equation (Partial differential Equation,

𝜕𝑢 𝜕2 𝑢
PDE) for the problem is given as = 𝛼 𝜕𝑥 2 where u be the temperature, x be the space
𝜕𝑡
coordinate, t be the time and 𝛼 be the coefficient of thermal diffusion (Considered as a
constant here) together with the Initial condition (IC) and boundary conditions (BC) as stated.
 For the formation of the formation of Finite Difference Equation (FDE) in this case we are using
𝜕𝑢
the FORWARD TIME CENTRAL SPACE (FTCS) SCHEME. In this scheme the TIME derivative ( )
𝜕𝑡
is represented by FORWARD DIFFERECE REPRESETATION and the SPACE derivative is
represented by CENTRAL DIFFERENCE representation. Thus (neglecting TE)

𝜕𝑢 𝑢𝑖𝑛+1 − 𝑢𝑖𝑛
=
𝜕𝑡 ∆𝑡
𝑛
𝜕2𝑢 𝑢𝑖+1 − 2𝑢𝑖𝑛 + 𝑢𝑖−1
𝑛
=
𝜕𝑥 2 ∆𝑥 2
Then the FDE can be written as

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛 𝑛
𝑢𝑖+1 − 2𝑢𝑖𝑛 + 𝑢𝑖−1
𝑛
= 𝛼 (Assuming 𝛼 = 1)
∆𝑡 ∆𝑥 2
𝜶 ∆𝒕
Or, 𝑢𝑖𝑛+1 = 𝑢𝑖𝑛 + ∆𝒙𝟐 ( 𝑢𝑖+1
𝑛
− 2𝑢𝑖𝑛 + 𝑢𝑖−1
𝑛
)
𝜶 ∆𝒕
Let us put =𝑟 in the above equation and will get
∆𝒙𝟐

𝑢𝑖𝑛+1 = 𝑢𝑖𝑛 + 𝑟 𝑢𝑖+1

𝑛
− 2𝑟 𝑢𝑖𝑛 + 𝑟 𝑢𝑖−1
𝑛

Or, 𝑢𝑖𝑛+1 = 𝑟 𝑢𝑖+1

𝑛
+ (1-2r) 𝑢𝑖𝑛 + 𝑟 𝑢𝑖−1
𝑛
(1)

n+1 time level

𝑢𝑖𝑛+1

n time level
𝑛
𝑢𝑖−1 𝑢𝑖𝑛 𝑛
𝑢𝑖+1

Fig.3 Computational Molecule for the FTCS Scheme

Calculation of Initial Conditions:

The IC is given as u(x,0)=sin πx for 0 ≤ 𝑥 ≤ 1 , at ∆𝑡 = 0, Thus for

𝑖 = 0, 𝑢(0.0) = 𝑢00 = sin 0 = 0
π
𝑖 = 1, 𝑢(1.0) = 𝑢10 = sin(π∆𝑥) = sin ( 4 ) = 0.707
π
𝑖 = 2, 𝑢(2.0) = 𝑢20 = sin(π(2∆𝑥)) = sin ( 2 ) = 1.0
3π
𝑖 = 3, 𝑢(3.0) = 𝑢30 = sin(π(3∆𝑥)) = sin ( ) = 0.707
4
𝑖 = 4, 𝑢(4.0) = 𝑢40 = sin(π(4∆𝑥)) = sin(π) = 0

Calculation Boundary Conditions:

The BCs are given as 𝑢(0, 𝑡) = 0 , 𝑢(1, 𝑡) = 0 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑡 ≥ 0. Thus for

the left end of the bar i.e at point O where 𝑖 = 0, 𝑢(0, 𝑡) = 𝑢0𝑛 =
0 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑣𝑎𝑙𝑢𝑒𝑠 𝑜𝑓 𝑡𝑖𝑚𝑒 𝑠𝑡𝑒𝑝

the right end of the bar i.e. at point A where 𝑖 = 𝑁 (ℎ𝑒𝑟𝑒 𝑁 = 4), 𝑢(4, 𝑡) = 𝑢4𝑛 =
0 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑣𝑎𝑙𝑢𝑒𝑠 𝑜𝑓 𝑡𝑖𝑚𝑒 𝑠𝑡𝑒𝑝

Right
Left Boundary Boundary

Time level n=3

Time level n=2

𝑢0𝑛+1 = 0 (BC) 𝑢4𝑛+1 = 0

Time level n+1=1
𝑢0𝑛 = 0 (BC)
i=0 i=1 i=2 i=3 𝑢4𝑛 = I=N
Initial Time 0 𝑢40 = 0
level 𝑢00 = 0 (IC) 𝑢10 = .707 𝑢20 = 1 𝑢30 = .707

Fig.4: IC, BC and time level

Solution of Governing Equation:

The governing equation is given as 𝑢𝑖𝑛+1 = 𝑟 𝑢𝑖+1

𝑛
+ (1-2r) 𝑢𝑖𝑛 + 𝑟 𝑢𝑖−1
𝑛

1
Let us assume 𝑟 =
2

-------------------------------------------------------------------------------------------------------------------------------

NOTE: Solution of this equation by FTCS scheme is sensitive to the value of the ratio r where

𝜶 ∆𝒕 1
𝑟= , involving the mesh lengths. It can be proved that if 0<𝑟 ≤ ,
∆𝒙𝟐 2

The FTCS scheme is both convergent and stable.

1
For 𝑟 > the scheme is neither convergent nor stable
2

------------------------------------------------------------------------------------------------------------
Now putting the value of r = ½ in the governing equation we get

1 1 1
𝑢𝑖𝑛+1 = 2 𝑢𝑖+1
𝑛
+ (1-2(2)) 𝑢𝑖𝑛 + (2) 𝑢𝑖−1
𝑛

1 1
Or, 𝑢𝑖𝑛+1 = 2 𝑢𝑖+1
𝑛 𝑛
+ 2 𝑢𝑖−1 (2)

Let us now solve the equation (2) for n+1=1 i.e. for n=0 for all ‘i’ values in that time level
1 1
Here 𝑢𝑖𝑛+1 = 𝑛
𝑢𝑖+1 + 𝑛
𝑢𝑖−1
2 2

For i =0 we get from the boundary value 𝑢10 = 0

1 1+0
for ‘i’ = 1 we get 𝑢11 = [𝑢20 + 𝑢00 ] = = 0.5
2 2

1 0.707+0.707
for ‘i’ = 2 we get 𝑢12 = [𝑢30 + 𝑢10 ] = = 0.707
2 2

1 0+1
for ‘i’ = 3 we get 𝑢13 = [𝑢40 + 𝑢20 ] = = 0.5
2 2

for i=4 we get from the boundary value 𝑢14 = 0

For Time level n=1

For i=0 we have 𝑢02 = 0 (BC)

1 .0.707+0
For i=1 we get 𝑢12 = [𝑢12 + 𝑢10 ] = = 0.354
2 2

1 0.5+0.5
For i=2 we get 𝑢22 = [𝑢13 + 𝑢11 ] = = 0.5
2 2

1 0+0.707
For i=3 we get 𝑢32 = [𝑢14 + 𝑢12 ] = = 0.354
2 2

For i=4 we have 𝑢42 = 0 (BC)

For Time level n=2

For i=0 we have 𝑢03 = 0 (BC)

1 0.5+0
For i=1 we get 𝑢13 = [𝑢22 + 𝑢02 ] = = 0.25
2 2

1 0.354+.354
For i=2 we get 𝑢23 = [𝑢32 + 𝑢12 ] = = 0.354
2 2

1 0+0.5
For i=3 we get 𝑢33 = [𝑢42 + 𝑢22 ] = = 0.25
2 2

For i=4 we have 𝑢43 = 0 (BC)

For Time level n=3

For i=0 we have 𝑢04 = 0 (BC)

1 0.354+0
For i=1 we get 𝑢14 = [𝑢23 + 𝑢03 ] = = 0.177
2 2

1 0.25+0.25
For i=2 we get 𝑢24 = [𝑢33 + 𝑢13 ] = = 0.25
2 2

1 0+0.354
For i=3 we get 𝑢34 = [𝑢43 + 𝑢23 ] = = 0.177
2 2

For i=4 we have 𝑢44 = 0 (BC)

Like this marching in the time direction at one stage say at n=P we get that u(i,p) = u(i,p-1) then we can
say that we have attained the steady state temperature distribution in the bar and then the iteration is
terminated, considering that solution is achieved.
 Right
Left Boundary Boudary
0 0.177 0.25 0.177
Time level n+1=4 0

0 0.25 0.354 0.25

Time level n+1=3
0

Time level n+1=2 0 0.354 0.5 0.354

0

Time level n+1=1 𝑢0𝑛+1 = 0 (BC) 0.5 0.707 0.5 0

𝑢0𝑛 = 0 (BC)
Time level i=0 i=1 I=2 I=3 𝑢4𝑛 =
Initial 0 𝑢4𝑛+1 = 0
𝑢00 = 0 (IC) 𝑢10 = .707 𝑢20 = 1 𝑢30 = .707 i=N
𝑢40 = 0
Fig.5: Computed solution for heat conduction Equation using FTCS Scheme
Problem 2.

𝜕𝑢 𝜕 2𝑢
Obtain the FDE for the heat conduction Equation 𝜕𝑡 = 𝛼 𝜕𝑥 2 for a bar of unit length using
FTCS scheme. where u be the temperature, x be the space coordinate, t be the time and 𝛼 be the
coefficient of thermal diffusion. Compute the solution for the first two time step where 𝛼 = 1, 𝑟 =
𝜶 ∆𝒕 1 1
= 2 and ∆𝑥 = 4.
∆𝒙𝟐
The initial condition is give as:
𝑢(𝑥, 0) = 𝑥 𝑓𝑜𝑟 0 ≤ 𝑥 ≤ 0.5 𝑎𝑛𝑑 𝑢(𝑥, 0) = (1 − 𝑥)𝑓𝑜𝑟 0.5 < 𝑥 ≤ 1
The Boundary Condition is given as:
At the ends of the bar at any time, is given as 𝑢(0, 𝑡) = 0 , 𝑢(1, 𝑡) = 0 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑡 ≥ 0.

Problem 3.
𝜕𝑢 𝜕 2𝑢
Obtain the FDE for the heat conduction Equation 𝜕𝑡 = 𝛼 𝜕𝑥 2 for a bar of unit length using
FTCS scheme. where u be the temperature, x be the space coordinate, t be the time and 𝛼 be the
coefficient of thermal diffusion. Compute the solution for the first two time step where 𝛼 = 1, 𝑟 =
𝜶 ∆𝒕 1 1
= and ∆𝑥 = 4.
∆𝒙𝟐 2
The initial condition is give as: u(x,0)=sin πx for 0 ≤ 𝑥 ≤ 1

The Boundary Condition is given as:

𝜕𝑢 𝜕𝑢
= 𝑢 𝑎𝑡 𝑥 = 0 𝑎𝑛𝑑 = −𝑢 𝑎𝑡 𝑥 = 1
𝜕𝑥 𝜕𝑥

Problem 4.
𝜕𝑢 𝜕 2𝑢
Obtain the FDE for the heat conduction Equation 𝜕𝑡 = 𝛼 𝜕𝑥 2 for a bar of unit length using FTCS
scheme. where u be the temperature, x be the space coordinate, t be the time and 𝛼 be the coefficient
𝜶 ∆𝒕 1
of thermal diffusion. Compute the solution for the first two time step where 𝛼 = 1, 𝑟 = ∆𝒙𝟐 = 2 and
1
∆𝑥 = .
4
The initial condition is give as:

𝑢(𝑥, 0) = 𝑥 𝑓𝑜𝑟 0 ≤ 𝑥 ≤ 0.5 𝑎𝑛𝑑 𝑢(𝑥, 0) = (1 − 𝑥)𝑓𝑜𝑟 0.5 < 𝑥 ≤ 1

The Boundary Condition is given as:

𝜕𝑢 𝜕𝑢
= 𝑢 𝑎𝑡 𝑥 = 0 𝑎𝑛𝑑 = −𝑢 𝑎𝑡 𝑥 = 1
𝜕𝑥 𝜕𝑥
Solution Prob 3:

Step 1: similar to prob.1

Step 2

Formation of Finite Difference Equation: Using FTCS scheme the FDE can be written as

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛 𝑛
𝑢𝑖+1 − 2𝑢𝑖𝑛 + 𝑢𝑖−1
𝑛

∆𝑡
= 𝛼 ∆𝑥 2
(here 𝛼 = 1, 𝑎𝑠 𝑔𝑖𝑣𝑒𝑛)

𝜶 ∆𝒕
Or, 𝑢𝑖𝑛+1 = 𝑢𝑖𝑛 + ∆𝒙𝟐 ( 𝑢𝑖+1
𝑛
− 2𝑢𝑖𝑛 + 𝑢𝑖−1
𝑛
)
𝜶 ∆𝒕
Let us put =𝑟 in the above equation and will get
∆𝒙𝟐

𝑢𝑖𝑛+1 = 𝑢𝑖𝑛 + 𝑟 𝑢𝑖+1

𝑛
− 2𝑟 𝑢𝑖𝑛 + 𝑟 𝑢𝑖−1
𝑛

Or, 𝑢𝑖𝑛+1 = 𝑟 𝑢𝑖+1

𝑛
+ (1-2r) 𝑢𝑖𝑛 + 𝑟 𝑢𝑖−1
𝑛
(1)
Again, we have r=1/2, Thus from equation 1 we get

1
𝑢𝑖𝑛+1 = 2 [𝑢𝑖+1
𝑛 𝑛
+ 𝑢𝑖−1 ] (2)

Calculation of Initial Conditions:

The IC is given as u(x,0)=sin πx for 0 ≤ 𝑥 ≤ 1 , at ∆𝑡 = 0, Thus for

𝑖 = 0, 𝑢(0,0) = 𝑢00 = sin 0 = 0
π
𝑖 = 1, 𝑢(1,0) = 𝑢10 = sin(π∆𝑥) = sin ( 4 ) = 0.707
π
𝑖 = 2, 𝑢(2,0) = 𝑢20 = sin(π(2∆𝑥)) = sin ( 2 ) = 1.0
3π
𝑖 = 3, 𝑢(3,0) = 𝑢30 = sin(π(3∆𝑥)) = sin ( 4 ) = 0.707
𝑖 = 4, 𝑢(4,0) = 𝑢40 = sin(π(4∆𝑥)) = sin(π) = 0
Boundary Conditions:
𝜕𝑢
At the left edge (at x=0) of the Bar BC is given as =u (3)
𝜕𝑥

Let us descritized the above equation using central differencing and then we get
𝑛 𝑛
𝑢𝑖+1 − 𝑢𝑖−1
= 𝑢𝑖𝑛
2∆𝑥
𝑢1𝑛 −𝑢−1
𝑛
Now for i=0 the above equation becomes = 𝑢0𝑛
2∆𝑥

𝑛 1
Or, 𝑢−1 = −2∆𝑥𝑢0𝑛 + 𝑢1𝑛 = − 2 𝑢0𝑛 + 𝑢1𝑛 (a) (as ∆𝑥 = 1/4)

Thus for this case we have to incorporate a false boundary at i = -1 and we assume that the PDE
( so also the FDE) is valid at the false boundary.
1
Again equation 2 is given as 𝑢𝑖𝑛+1 = 2 [𝑢𝑖+1
𝑛 𝑛
+ 𝑢𝑖−1 ],

1
For i=0 this equation gives 𝑢0𝑛+1 = 2 [𝑢1𝑛 + 𝑢−1
𝑛
] (b)

Now replacing the value of 𝑢𝑛

−1 from equation (b) using equation (a), we may get

1 𝑛 1 1
𝑢0𝑛+1 = [𝑢1 + (− 𝑢0𝑛 + 𝑢1𝑛 )] = 𝑢1𝑛 − 𝑢0𝑛
2 2 4
𝜕𝑢
At the right edge (at x=1) of the Bar BC is given as = -u (4). Descritizig eqn. 4 we get
𝜕𝑥

𝑛 𝑛
𝑢𝑖+1 −𝑢𝑖−1
= −𝑢𝑖𝑛
2∆𝑥

𝑢5𝑛 −𝑢3𝑛
For i=4 we get = −𝑢4𝑛 ,
2∆𝑥

1
or, 𝑢5𝑛 = −2∆𝑥 𝑢4𝑛 + 𝑢3𝑛 = − 2 𝑢4𝑛 + 𝑢3𝑛 (c)

for this case also we have to incorporate a false boundary at i = 5 and we assume that the PDE (
so also the FDE) is valid at the false boundary.
1 1
Again equation 2 is given as 𝑢𝑖𝑛+1 = 𝑛
[𝑢𝑖+1 𝑛
+ 𝑢𝑖−1 ], for i=4 we get 𝑢4𝑛+1 = [𝑢5𝑛 + 𝑢3𝑛 ] (d)
2 2
 1 1 1
From (d) and (c ) 𝑢4𝑛+1 = 2 [− 2 𝑢4𝑛 + 𝑢3𝑛 + 𝑢3𝑛 ] = 𝑢3𝑛 − 4 𝑢4𝑛

Solution of Governing Equations:

1
The governing equation is 𝑢𝑖𝑛+1 = 2 [𝑢𝑖+1
𝑛 𝑛
+ 𝑢𝑖−1 ] with IC and BC

Computation for n+1 = 0 i.e for n=0

1 1
For i=0 𝑢10 = 𝑢10 − 𝑢00 = 0.707 − 4 x0 = 0.707 (Boundary value)
4

1 1
For i=1 𝑢11 = 2 [𝑢20 + 𝑢00 ] = 2[1+0] = 0.5

1 1
For i=2 𝑢12 = 2 [𝑢30 + 𝑢10 ] = 2[0.707+0.707] = 0.707

1 1
For i=3 𝑢13 = 2 [𝑢40 + 𝑢20 ] = 2[0+1] = 0.5
1 1
For i=4 𝑢14 = 𝑢30 − 4 𝑢40 = 0.707 − 4 x 0 = 0.707 (Boundary value)

Computation for n = 1
1 1
For i=0 𝑢02 = 𝑢11 − 𝑢10 = 0.5 − 4 x0.707 = 0.323 (Boundary value)
4

1 1
For i=1 𝑢12 = 2 [𝑢12 + 𝑢10 ] = 2[0.707+0.707] = 0.707

1 1
For i=2 𝑢22 = 2 [𝑢13 + 𝑢11 ] = 2[0.5+0.5] = 0.5
1 1
For i=3 𝑢32 = 2 [𝑢14 + 𝑢12 ] = 2[0.707+0.707] = 0.707
1 1
For i=4 𝑢42 = 𝑢13 − 4 𝑢14 = 0.5 − 4 x 0.707 = 0.323 (Boundary value)

Computation for n = 2
1 1
For i=0 𝑢03 = 𝑢12 − 𝑢02 = 0.707 − 4 x0.323 = 0.626 (Boundary value)
4
1 1
For i=1 𝑢13 = 2 [𝑢22 + 𝑢02 ] = 2[0.5+0.323] = 0.412
1 1
For i=2 𝑢23 = 2 [𝑢32 + 𝑢12 ] = 2[0.707+0.707] = 0.707
1 1
For i=3 𝑢33 = [𝑢42 + 𝑢22 ] = [0.323+0.5] = 0.412
2 2
 1 1
For i=4 𝑢43 = 𝑢32 − 4 𝑢42 = 0.707 − 4 x 0.323 = 0.626 (Boundary value)

2=2

n=1
Time step

n=0

i=-1 i=0 i=1 i=2 i=3 i=4 i=5

Some other Explicit Schemes:

2. Leap – Frog Scheme:

In this scheme both the time derivative and space derivative of the 1-D heat
𝜕𝑢 𝜕 2𝑢
conduction equation 𝜕𝑡
= 𝛼 𝜕𝑥 2 are represented by Central Difference representation.
Thus the FDE ca be written as-

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛−1 𝑛
𝑢𝑖−1 − 2 𝑢𝑖𝑛 + 𝑢𝑖+1
𝑛
= 𝛼
2∆𝑡 ∆𝑥 2

This scheme is second order accurate I both time and space. However the scheme is
UCONDITIONALLY UNSTABLE for parabolic types of equation and hence is not used
for solution of Parabolic types of equation.

3. Dufort – Frankel Scheme:

Dufort & Frankel (1953) observed that in the Leap – Frog scheme, if the middle term
of Right hand side of the FDE i.e. 𝑢𝑖𝑛 term is approximated by the average value of
𝑢𝑖 at (n-1) and (n+1) time levels, the stability behavior of the scheme is improved
𝑢𝑖𝑛+1 +𝑢𝑖𝑛−1
significantly. i.e the term 𝑢𝑖𝑛 is replaced by
2

Thus the FDE is modified as

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛−1 𝑛
𝑢𝑖−1 − (𝑢𝑖𝑛+1 + 𝑢𝑖𝑛−1 ) + 𝑢𝑖+1
𝑛
= 𝛼
2∆𝑡 ∆𝑥 2
And this scheme is termed as Dufort – Frankel Scheme. The scheme uses previous
data of two time steps i.e n & n-1 steps while calculating values at present step (i.e
n+1 step).
The scheme is unconditionally stable, however under certain conditions it is not
consistent.
 IMPLICIT SCHEMES for Parabolic Equations:

Crank – Nicholson scheme: (for 1-D Heat Conduction Equation)

In this scheme the time derivative is represented by forward differences while the space
derivative is represented by the average central difference at the present time level nth
and next time level (n+1)th ( here, solution is sought at n+1 time level). Neglecting the
truncation errors (TE) the Crank-Nicholson discretisation of the 1-D heat conduction
𝜕𝑢 𝜕 2𝑢
equation 𝜕𝑡
= 𝛼 𝜕𝑥 2 is give as –

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛 𝑛+1

𝛼 𝑢𝑖−1 − 2 𝑢𝑖𝑛+1 + 𝑢𝑖+1
𝑛+1 𝑛
𝑢𝑖−1 − 2 𝑢𝑖𝑛 + 𝑢𝑖+1
𝑛
= [ + ]
∆𝑡 2 ∆𝑥 2 ∆𝑥 2

Simplifying we get
𝑛+1
−𝑟 𝑢𝑖+1 + 2(1 + 𝑟)𝑢𝑖𝑛+1 − 𝑟 𝑢𝑖−1
𝑛+1 𝑛
= 𝑟 𝑢𝑖+1 + 2(1 − 𝑟) 𝑢𝑖𝑛 + 𝑟𝑢𝑖−1
𝑛
(A1)
𝛼 ∆𝑡
𝑤ℎ𝑒𝑟𝑒 𝑟 = 2∆𝑥 2

Simplification:

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛 𝑛+1

𝛼 𝑢𝑖−1 − 2 𝑢𝑖𝑛+1 + 𝑢𝑖+1
𝑛+1 𝑛
𝑢𝑖−1 − 2 𝑢𝑖𝑛 + 𝑢𝑖+1
𝑛
= [ + ]
∆𝑡 2 ∆𝑥 2 ∆𝑥 2

𝛼 ∆𝑡
Or, 𝑢𝑖𝑛+1 − 𝑢𝑖𝑛 = 𝑛+1
[𝑢𝑖−1 − 2 𝑢𝑖𝑛+1 + 𝑢𝑖+1
𝑛+1 𝑛
+ 𝑢𝑖−1 − 2 𝑢𝑖𝑛 + 𝑢𝑖+1
𝑛
}
2∆𝑥 2

𝑟 𝑛+1 𝑟 𝑟 𝑟
or, − 𝑢𝑖+1 + 𝑟𝑢𝑖𝑛+1 + 𝑢𝑖𝑛+1 − 𝑛+1
𝑢𝑖−1 = 𝑢𝑖𝑛 + 𝑛
𝑢𝑖+1 − 𝑟 𝑢𝑖𝑛 + 2 𝑢𝑖−1
𝑛
2 2 2

multiplying both side by 2

𝑛+1
− 𝑟 𝑢𝑖+1 + 2(1 + 𝑟)𝑢𝑖𝑛+1 − 𝑟 𝑢𝑖−1
𝑛+1 𝑛
= 𝑟 𝑢𝑖+1 + 2(1 − 𝑟 ) 𝑢𝑖𝑛 + 𝑟 𝑢𝑖−1
𝑛
The internal mesh points are 𝑖 = 1,2,3, … … … . . (𝑁 − 1) and boundary points are i=0 and i=i
maximum = N. The initial values for all ‘i’ at starting time level are prescribed or given values.

Moreover, boundaries values at ‘i’ = 0 and at ‘i’ = N are prescribed values (as shown in the
previous problems).
In this scheme the equation (A1) represent a set of equation for various values of ‘i’ . The
quantities at time level (n+1) are unknown while those of at the nth level are known
𝑛+1
quantities. Here, each equation consists of three unknown namely, 𝑢𝑖+1 , 𝑢𝑖𝑛+1 𝑎𝑛𝑑 𝑢𝑖−1
𝑛+1

except the first equation for i=0 and the last equation for i=i max – 1 . These two equations
contain only two unknowns.
The system of equation so obtained are form a ‘Tri-diagonal System of Equation’ which ca be
solved using ‘Thomas Algorithm’.

Thomas Algorithm for Tri-diagonal System of Equation:

To illustrate the algorithm let consider the example of problem 1. The governing equation for
𝜕𝑢 𝜕 2𝑢
the heat conduction in the bar is given as 𝜕𝑡
= 𝛼 𝜕𝑥 2 . The FDE for the PDE using C-N scheme can
be obtained as
𝑛+1
−𝑟 𝑢𝑖+1 + 2(1 + 𝑟)𝑢𝑖𝑛+1 − 𝑟 𝑢𝑖−1
𝑛+1 𝑛
= 𝑟 𝑢𝑖+1 + 2(1 − 𝑟) 𝑢𝑖𝑛 + 𝑟𝑢𝑖−1
𝑛
(A1)

From this equation (A1) for different values of I we have:

For i=1 𝑎1 𝑢0 + 𝑏1 𝑢1 + 𝑐1 𝑢2 = 𝐾𝑛𝑜𝑤 𝑣𝑎𝑙𝑢𝑒𝑠 , here ′𝑢0 ′ be a boundary value, Hence we can write

for i=1 𝑏1 𝑢1 + 𝑐1 𝑢2 = 𝑑1

For i=2 𝑎2 𝑢1 + 𝑏2 𝑢2 + 𝑐2 𝑢3 = 𝑑2
For i=3 𝑎3 𝑢2 + 𝑏3 𝑢3 + 𝑐3 𝑢4 = 𝑑3
For any i 𝑎𝑖 𝑢𝑖−1 + 𝑏𝑖 𝑢𝑖 + 𝑐𝑖 𝑢𝑖+1 = 𝑑𝑖
For i= imax - 1 𝑎𝑀 𝑢𝑀−1 + 𝑏𝑀 𝑢𝑀 = 𝑑𝑀
In the above set of equation u’s are unknown. The coefficient 𝑎𝑖 , 𝑏𝑖 , 𝑐𝑖 and 𝑑𝑖 are prescribed
quantities.
To solve the set of equation let us start with first equation
𝑏1 𝑢1 + 𝑐1 𝑢2 = 𝑑1
𝑐1 𝑑1
𝑑𝑖𝑣𝑖𝑑𝑖𝑛𝑔 𝑏𝑜𝑡ℎ 𝑠𝑖𝑑𝑒 𝑏𝑦 𝑏1 , 𝑤𝑒 𝑔𝑒𝑡 𝑢1 + 𝑢2 =
𝑏1 𝑏1
𝑑1 𝑐1
𝑜𝑟, 𝑢1 = − 𝑢2
𝑏1 𝑏1
𝑑1 −𝑐1
𝑜𝑟, 𝑢1 = 𝑝1 + 𝑞1 𝑢2 where 𝑝1 = 𝑏1
𝑎𝑛𝑑 𝑞1 = 𝑏1
 Now substituting 𝑢1 value in the second equation 𝑎2 𝑢1 + 𝑏2 𝑢2 + 𝑐2 𝑢3 = 𝑑2 we get

𝑎2 (𝑝1 + 𝑞1 𝑢2 ) + 𝑏2 𝑢2 + 𝑐2 𝑢3 = 𝑑2
𝑑2 − 𝑝1 𝑎2 −𝑐2
𝑜𝑟, 𝑢2 = 𝑝2 + 𝑞2 𝑢3 where 𝑝2 = and 𝑞2 =
𝑏2 + 𝑞1 𝑎2 𝑏2 + 𝑞1 𝑎2

𝑑𝑖 − 𝑝𝑖−1 𝑎𝑖 −𝑐𝑖
Similarly, , 𝑢𝑖 = 𝑝𝑖 + 𝑞𝑖 𝑢𝑖+1 where 𝑝𝑖 = and 𝑞2 =
𝑏𝑖 + 𝑞𝑖−1 𝑎𝑖 𝑏𝑖 + 𝑞𝑖−1 𝑎𝑖
For i= M-1, M-2, M-3, ……….. 1

𝑑𝑀 − 𝑝𝑀−1 𝑎𝑀
Finally, we get, 𝑢𝑀 = 𝑝𝑀 where 𝑝𝑀 =
𝑏𝑀 + 𝑞𝑀−1 𝑎𝑀
Now, finding the last value of u i.e. 𝑢𝑀 we can find the other values of u.

Problem no.4:

𝜕𝑢 𝜕2 𝑢
Obtain the FDE for the heat conduction Equation 𝜕𝑡 = 𝛼 𝜕𝑥 2 for a bar of unit length using
Crank-Nicholson scheme. where u be the temperature, x be the space coordinate, t be the time
and 𝛼 be the coefficient of thermal diffusion. Compute the solution for the first two time step
𝜶 ∆𝒕 1
where 𝛼 = 1, 𝑟 = ∆𝒙𝟐 = 1 and ∆𝑥 = 4.
The initial condition is give as:
𝑢(𝑥, 0) = sin 𝜋 𝑥 𝑓𝑜𝑟 0 ≤ 𝑥 ≤ 1
The Boundary Condition is given as:
At the ends of the bar at any time, is given as 𝑢(0, 𝑡) = 0 , 𝑢(1, 𝑡) = 0 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑡 ≥ 0.

Solution:

Step1: Grid point Generation:

1 𝐿𝑒𝑛𝑔𝑡ℎ 𝑜𝑓 𝑡ℎ𝑒 𝑏𝑎𝑟
It is given that ∆𝑥 = 4 , Again we have ∆𝑥 = 𝑇𝑜𝑡𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝐷𝑖𝑣𝑖𝑠𝑖𝑜𝑛
𝐿𝑒𝑛𝑔𝑡ℎ 𝑜𝑓 𝑡ℎ𝑒 𝑏𝑎𝑟 1
Therefore, Total number of Division = = 1 = 4 and total grid point will be 5
∆𝑥
4
Viz. i= 0,1,2,3 and 4 .
 Right
Left Boundary Boundary

Time level n=3

Time level n=2

𝑢0𝑛+1 = 0 (BC) 𝑢4𝑛+1 = 0

Time level n+1=1

i=0 i=1 i=2 i=3 I=N

Initial Time 𝑢00 = 0 (IC) 𝑢10 = .707 𝑢20 = 1 𝑢30 = .707 𝑢40 = 0
level

Step 2: formation of FDE:

The governing equation is given as

𝜕𝑢 𝜕2𝑢
= 𝛼 2
𝜕𝑡 𝜕𝑥

The finite difference equation using Crank-Nicholson scheme will be

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛 𝑛+1

𝛼 𝑢𝑖−1 − 2 𝑢𝑖𝑛+1 + 𝑢𝑖+1
𝑛+1 𝑛
𝑢𝑖−1 − 2 𝑢𝑖𝑛 + 𝑢𝑖+1
𝑛
= [ + ]
∆𝑡 2 ∆𝑥 2 ∆𝑥 2

Simplifying we get
 𝑛+1
−𝑟 𝑢𝑖+1 + 2(1 + 𝑟)𝑢𝑖𝑛+1 − 𝑟 𝑢𝑖−1
𝑛+1 𝑛
= 𝑟 𝑢𝑖+1 + 2(1 − 𝑟) 𝑢𝑖𝑛 + 𝑟𝑢𝑖−1
𝑛
(A1)
𝛼 ∆𝑡 1
𝑤ℎ𝑒𝑟𝑒 𝑟 = , In this it is given that r=1 and ∆𝑥 = 4
2∆𝑥 2
Putting the value of r=1 in the equation (A1), we get

𝑛+1
− 𝑢𝑖+1 + 4𝑢𝑖𝑛+1 − 𝑢𝑖−1
𝑛+1 𝑛
= 𝑢𝑖+1 𝑛
+ 𝑢𝑖−1 (1)

Computation of Initial values and Boundary values:

The grid points are designated as i= 0,1,2,3 and 4.

Calculation of Initial values at various grid points:

The IC is given as u(x,0)=sin πx for 0 ≤ 𝑥 ≤ 1 , at ∆𝑡 = 0, Thus for

𝑖 = 0, 𝑢(0,0) = 𝑢00 = sin 0 = 0
π
𝑖 = 1, 𝑢(1,0) = 𝑢10 = sin(π∆𝑥) = sin ( 4 ) = 0.707
π
𝑖 = 2, 𝑢(2,0) = 𝑢20 = sin(π(2∆𝑥)) = sin ( 2 ) = 1.0
3π
𝑖 = 3, 𝑢(3,0) = 𝑢30 = sin(π(3∆𝑥)) = sin ( 4 ) = 0.707
𝑖 = 4, 𝑢(4,0) = 𝑢40 = sin(π(4∆𝑥)) = sin(π) = 0

The boundary value will be

at i=0 𝑢0𝑛+1 = 0 and at i=4 𝑢4𝑛+1 = 0 for all values of n

Formation of set of FDE for a certain time level ad its Solution by Thomas Algorithm:

Let us compute the values of for the first time step n+1 = 1 where initial time step is n=0
Putting n+1 = 1 and =0 in equation (1) we get
𝑛+1
− 𝑢𝑖+1 + 4𝑢𝑖𝑛+1 − 𝑢𝑖−1
𝑛+1 𝑛
= 𝑢𝑖+1 𝑛
+ 𝑢𝑖−1
For i=1 −𝑢12 + 4𝑢11 − 𝑢10 = 𝑢20 + 𝑢00
or, 4𝑢11 − 𝑢12 = 1 here 𝑢10 = 0 (𝐵𝐶), 𝑢20 = 1 (𝐼𝐶), and 𝑢00 = 0 (𝐼𝐶)
or, 4𝑢1 − 𝑢2 = 1 [for simplicity super script is omitted]

for i=2 −𝑢3 + 4𝑢2 − 𝑢1 = 𝑢3 + 𝑢1

or, −𝑢3 + 4𝑢2 − 𝑢1 = 0.707 + 0.707 here 𝑢3 = 0.707 and 𝑢1 = 0.707
or, −𝑢1 + 4𝑢2 − 𝑢3 = 1.14

for i=3 −𝑢4 + 4𝑢3 − 𝑢2 = 𝑢4 + 𝑢2

or, −𝑢2 + 4𝑢3 = 1 here 𝑢14 = 𝑢4 = 0, (BC), 𝑢20 = 1 (𝐼𝐶), 𝑢40 = 0 (𝐼𝐶)

Thus the set of equation obtained is

, 4𝑢1 − 𝑢2 = 1
−𝑢1 + 4𝑢2 − 𝑢3 = 1.14
−𝑢2 + 4𝑢3 = 1

Comparing the present equation with standard equation of Thomas Algorithm , we can have
𝑑1 −𝑐1 −1 1
𝑢1 = 𝑝1 + 𝑞1 𝑢2 where 𝑝1 = 𝑏1
= 14 and 𝑞1 = =− =4
𝑏1 4
𝑢2 = 𝑝2 + 𝑞2 𝑢3 where
1×(−1)
𝑑2 − 𝑝1 𝑎2 1.42− −𝑐2 −(−1)
4
𝑝2 = 𝑏2 + 𝑞1 𝑎2
= 1 = 0.445 and 𝑞2 = = 1 = 0.267
4+ ×(−1) 𝑏2 + 𝑞1 𝑎2 4+ ×(−1)
4 4

𝑑3 − 𝑝2 𝑎3 1−0.445×(−1)
𝑢3 = 𝑝3 where 𝑝3 = = 4+(0.267×(−1) = 0.387
𝑏3 + 𝑞2 𝑎3
Thus
𝑢3 = 𝑝3 = 0.387
𝑢2 = 𝑝2 + 𝑞2 𝑢3 = 0.445 + 0.267 × 0.387 = 0.548
1 1
𝑢1 = 𝑝1 + 𝑞1 𝑢2 = + × 0.548 = 0.387
4 4

Computation for time level n+1=2 (HW)