140 
 

 American Scientific Research Journal for Engineering, Technology,  and Sciences  (ASRJETS) 
ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 

© Global Society of Scientific Research and Researchers  
http://asrjetsjournal.org/  

 

Application of Some Finite Difference Schemes for Solving 

One Dimensional Diffusion Equation 

Tsegaye Simona, Purnachandra Rao Koyab* 

a,bSchool of Mathematical and Statistical Sciences, Hawassa University, Hawassa, Ethiopia 
aEmail: tsegaye_simon@yahoo.com 

bEmail: drkpraocecc@yahoo.co.in 

 

 

Abstract 

In this paper the numerical solutions of one dimensional diffusion equation using some finite difference methods 

have been considered. For that purpose three examples of the diffusion equation together with different 

boundary conditions are examined. The finite difference methods applied on each example are (i) forward time 

centered space (ii) backward time centered space and (iii) Crank – Nicolson. In each case, we have studied 

stability of finite difference method and also obtained numerical result. The performance of each scheme is 

evaluated in accordance with both the accuracy of the solution and programming efforts. The implementation 

and behavior of the schemes have been compared and the results are illustrated pictorially. It is found in case of 

the test examples studied here that the Crank – Nicolson scheme gives better approximations than the two other 

schemes. 

Keywords: Crank – Nicolson; Diffusion equation; Forward time centered space; Backward time centered space; 

Stability.  

1. Introduction 

One dimensional diffusion equation plays an important role in modeling numerous physical phenomena. The 

application of such diffusion equation includes a wide range of areas such as physical, biological and financial 

sciences. One of the most common applications is propagation of heat in the environment, where 𝑢𝑢(𝑥𝑥, 𝑡𝑡) 

represents the temperature of some substance at point 𝑥𝑥 and time 𝑡𝑡.  

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

* Corresponding author.  

http://asrjetsjournal.org/


American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

141 
 

The diffusion or heat equation is applied when attempting to describe the density fluctuations in a material that 

undergoes diffusion [8]. In the diffusion equation there appear derivatives with respect to time and space 

coordinates.  

Several finite difference techniques have been employed to solve these equations so as to fit their physical 

nature [1 – 2, 5, 7, 13,  14]. The process of solving requires specification of a suitable boundary conditions viz., 

(i) Dirichlet (ii) Neumann or gradient (iii) Mixed or Robin (iv) Combination of Dirichlet and Neumann. 

Boundary conditions are applied on the spatial coordinate at 𝑥𝑥 = 0 and 𝑥𝑥 = 𝐿𝐿 and initial conditions are applied 

on the temporal coordinate when  𝑡𝑡 = 0. In the process of finding numerical solutions the continuous partial 

differential equations are replaced by their discrete approximations. In the present context the word discrete 

indicates that the numerical solution is known only at a finite number of points in the physical domain. The 

number of those discrete points can vary and can be fixed by the user of the numerical method. However, 

increment in the number of discrete points increases not only the resolution but also the accuracy of the 

numerical solution. 

The process of discretization of a diffusion equation leads to a set of algebraic equations. These algebraic 

equations are evaluated so as to obtain values for the unknown quantities of the discretization. In turn, the values 

of unknowns provide an approximate solution to the original diffusion equation. 

We now divide the  𝑥𝑥𝑡𝑡   domain into a mesh. The coordinate axes are divided into steps of uniform lengths ∆𝑥𝑥 

and ∆𝑡𝑡 along 𝑥𝑥 and 𝑡𝑡 axis respectively. When horizontal and vertical lines are drawn along the step nodes the 

resulting image would resemble as a net or mesh. The two key parameters of the mesh are ∆𝑥𝑥 and   ∆𝑡𝑡. The 

former denotes the local distance between adjacent points in space while the latter denotes local distance 

between adjacent time steps. The intersection points of this mesh are called nodes. The discrete solutions are 

computed at the mesh nodes.  

 The core idea of the finite difference scheme is to replace continuous derivatives with difference formulae. 

Difference formulae provide discrete values of the function at nodes of the mesh. A variety of finite difference 

schemes are popular and widely used. Use of different combinations of mesh points in the difference formulas 

results in different schemes. In the limit of the mesh step spacing tending to zero i.e., as  ∆𝑥𝑥 → 0  and  ∆𝑡𝑡 → 0 , 

the numerical solution obtained by any useful scheme will approach the true solution of the differential 

equation. However, the rate at which the numerical solution converges to the true solution varies with the 

scheme. In addition, there are some practically useful schemes which may fail to yield a solution for bad 

combinations of ∆𝑥𝑥 and  ∆𝑡𝑡. Hence, the selection of combinations of ∆𝑥𝑥 and ∆𝑡𝑡 plays an important role [3]. 

Further, considerable attention is required so that the physical interpretations of solutions of diffusion equation 

remain meaningful. That is, we have to select such a combination of ∆𝑥𝑥 and ∆𝑡𝑡 that provides physically 

meaningful approximate solution for the diffusion equation. 

2. Governing Equation and Finite Difference Schemes 

 Here we now introduce and discuss governing equation describing one dimensional diffusion. Also we present 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

142 
 

the finite difference methods viz., forward time centered space or FTCS, backward time centered space or BTCS 

and Crank – Nicolson schemes. 

2.1 Governing Equation 

The more general diffusion equation is a partial differential equation and it describes the density fluctuations in 

the material undergoing diffusion. The equation can be expressed as: 

𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

=  ∇. (𝐷𝐷∇𝑢𝑢)                                                                               (1) 

Here in (1), 𝑢𝑢 ≡  𝑢𝑢(𝒙𝒙, 𝑡𝑡)  denotes the density of the diffusing material at location 𝒙𝒙 = (𝑥𝑥, 𝑦𝑦, 𝑧𝑧) and at time  𝑡𝑡. 

Also 𝐷𝐷 ≡ 𝐷𝐷(𝑢𝑢(𝒙𝒙, 𝑡𝑡), 𝑥𝑥) denotes the collective diffusion coefficient for the density 𝑢𝑢 at location 𝐱𝐱.  

Now for simplicity let the diffusion coefficient be independent of both density and location i.e., 𝐷𝐷 is a constant. 

Thus equation (1) reduces to a linear form as 

𝜕𝜕𝜕𝜕 
𝜕𝜕𝜕𝜕

=  𝐷𝐷∇2𝑢𝑢 .                                                                              (2) 

  Equation (2) is a diffusion or heat equation and describes the distribution of material or heat in a given region 

over time with constant diffusion coefficient. In the present study we further simplify (2) and consider one – 

dimensional diffusion with constant diffusion coefficient  𝐷𝐷. This consideration simplifies (2) to give 

𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕

 =   𝐷𝐷 𝜕𝜕2𝜕𝜕
𝜕𝜕𝑥𝑥2

                                                                                        (3) 

Here in (3),  𝑢𝑢 is the density of the diffusing material at location 𝑥𝑥 and time  𝑡𝑡. 𝐷𝐷 is the diffusivity coefficient in 

the 𝑥𝑥 – direction. We impose appropriate initial and boundary conditions on (3) in order to evaluate a numerical 

or approximate solution of it. Several combinations of boundary conditions are possible. We consider three 

distinct cases and apply them in a finite domain as follows: 

Case 1: Dirichlet type of boundary conditions 

𝑢𝑢(𝑥𝑥, 0) = 𝑓𝑓(𝑥𝑥)                                                                                           (4) 

                                   𝑢𝑢(0, 𝑡𝑡) =  𝑢𝑢0      

                                   𝑢𝑢(𝐿𝐿, 𝑡𝑡) =  𝑢𝑢𝐿𝐿  

Case 2: Neumann type of boundary conditions 

𝑢𝑢(𝑥𝑥, 0) = 𝑓𝑓(𝑥𝑥)                                                                                           (5) 

                                   𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥

(0, 𝑡𝑡)  =  𝑢𝑢0   



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

143 
 

                                   𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥

(𝐿𝐿, 𝑡𝑡)  =  𝑢𝑢𝐿𝐿 

Case 3: One end Neumann and other end Dirichlet type of boundary conditions: 

𝑢𝑢(𝑥𝑥, 0) = 𝑓𝑓(𝑥𝑥)                                                                                           (6) 

                                   𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥

(0, 𝑡𝑡)  =  0   

                                   𝑢𝑢(𝐿𝐿, 𝑡𝑡)  =  𝑢𝑢𝐿𝐿 

Here in (4) to (6), the quantities 𝑢𝑢0 and 𝑢𝑢𝐿𝐿 represent constant densities of the diffusing material respectively at 

 𝑥𝑥 = 0  and  𝑥𝑥 = 𝐿𝐿. The function 𝑓𝑓(𝑥𝑥) is an arbitrary one of its argument. 

2.2 Forward Time Centered Space Explicit Scheme 

The equation (3) does not always exhibit an analytical solution and even if it exhibits finding is not easy. Hence, 

finite difference schemes are applied for finding approximate solutions. We now approximate equation (3) by 

applying forward difference on time derivative and central difference on space derivative. Thus equation (3) 

takes the form as 

𝜕𝜕𝑖𝑖
𝑛𝑛+1− 𝜕𝜕𝑖𝑖

𝑛𝑛

∆𝜕𝜕
  =  𝐷𝐷 �𝜕𝜕𝑖𝑖−1

𝑛𝑛 −2 𝜕𝜕𝑖𝑖
𝑛𝑛+ 𝜕𝜕𝑖𝑖+1

𝑛𝑛  
∆𝑥𝑥2

� + 𝑂𝑂(∆𝑡𝑡) + 𝑂𝑂(∆𝑥𝑥2).                                           (7) 

The temporal error 𝑂𝑂(∆𝑡𝑡) and spatial error  𝑂𝑂(∆𝑥𝑥2) have different orders. Their values are very small and their 

influence on the solution is negligible. The big O notation expresses the rate at which the truncation error goes 

to zero. Hence drop the truncation error terms from equation (7) and after a rearrangement it leads to 

𝑢𝑢𝑖𝑖𝑛𝑛+1 = 𝑟𝑟𝑢𝑢𝑖𝑖−1𝑛𝑛 + (1 − 2𝑟𝑟)𝑢𝑢𝑖𝑖𝑛𝑛 + 𝑟𝑟𝑢𝑢𝑖𝑖+1𝑛𝑛                                                          (8) 

In (8), we have used the notation  𝑟𝑟 = 𝐷𝐷(∆𝑡𝑡 ∆𝑥𝑥2⁄ ). Furthermore, equation (8) can be expressed in terms of 

matrix multiplication as 

𝑢𝑢𝑛𝑛+1 = (𝐼𝐼 + 𝑟𝑟𝑟𝑟) 𝑢𝑢𝑛𝑛                                                                           (9) 

In (9), we have used the following matrix notations 

𝑢𝑢𝑛𝑛+1 =  [𝑢𝑢1𝑛𝑛+1,  𝑢𝑢2𝑛𝑛+1, . . . ,  𝑢𝑢𝑁𝑁−1𝑛𝑛+1  ]𝑇𝑇 

𝑢𝑢𝑛𝑛 =  [𝑢𝑢1𝑛𝑛 + 𝑟𝑟𝑢𝑢0𝑛𝑛 ,  𝑢𝑢2𝑛𝑛 , . . . ,  𝑢𝑢𝑁𝑁−1𝑛𝑛 +   𝑢𝑢𝑁𝑁𝑛𝑛 ]𝑇𝑇  

𝐼𝐼 =  �

1   0    0  .  .  .    0
0    1    0  .  .  .    0
.      .     .    .  .  .   .

0    0    0  .  .  .    1 

� 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

144 
 

𝑟𝑟 =  

⎣
⎢
⎢
⎢
⎡
−2    1    0  .  .  .  0    0    0
1   − 2    1  .  .  .  0    0    0
.          .      .  .  .  .   .     .     .

0    0    0  .  .  .  1   − 2    1
0    0    0  .  .  .  0    1   − 2⎦

⎥
⎥
⎥
⎤
.                                          (10) 

In (9), both  𝑢𝑢𝑛𝑛 and 𝑢𝑢𝑛𝑛+1 are column vectors of dimension  𝑁𝑁 − 1. The superscript 𝑇𝑇 denotes transpose of the 

matrix.  The  𝑁𝑁 − 1 dimensional square matrices  𝐼𝐼 and 𝑟𝑟 respectively denote an identity and a tridiagonal 

matrices. 

To find the stability condition of equation (8), we substitute a trial solution or a Fourier mode 

𝑢𝑢𝑖𝑖𝑛𝑛 =  𝜆𝜆𝑛𝑛𝑒𝑒[𝑛𝑛𝑖𝑖𝑛𝑛 𝑝𝑝⁄ ]                                                                         (11) 

to get  

𝜆𝜆 = 𝑟𝑟�𝑒𝑒𝑛𝑛𝑛𝑛 𝑝𝑝⁄ + 𝑒𝑒   −𝑛𝑛𝑛𝑛 𝑝𝑝⁄  � + 1 − 2𝑟𝑟 = 1 − 2𝑟𝑟[1 − cos (𝜋𝜋 𝑝𝑝⁄ )] 

Here 𝑝𝑝 is any non-zero integer. Since we must have |𝜆𝜆|  ≤ 1 for non – divergence, the stability condition turns 

out to be  𝑟𝑟 = [𝐷𝐷∆𝑡𝑡 ∆𝑥𝑥2⁄ ] ≤ 1 2⁄  . That is, the equation (8) is stable as long as the spatial interval ∆𝑥𝑥 satisfies 

the condition ∆𝑥𝑥 ≤ √2𝐷𝐷 ∆𝑡𝑡 for any given time interval   ∆𝑡𝑡. Otherwise (8) will not be stable [13]. 

Further, the truncation error  𝑇𝑇𝑖𝑖,𝑛𝑛  for the equation (8) can be derived and expressed as [14] 

𝑇𝑇𝑖𝑖,𝑛𝑛 =   �
𝜕𝜕𝑢𝑢
𝜕𝜕𝑡𝑡

 −  𝐷𝐷
𝜕𝜕2𝑢𝑢
𝜕𝜕𝑥𝑥2

�
𝑖𝑖,𝑛𝑛

+  
1
2
∆𝑡𝑡 �

𝜕𝜕2𝑢𝑢
𝜕𝜕𝑡𝑡2

�
𝑖𝑖,𝑛𝑛

+ 
1

12
(∆𝑥𝑥)2 �

𝜕𝜕4𝑢𝑢
𝜕𝜕𝑥𝑥4

�
𝑖𝑖,𝑛𝑛

+ 𝑂𝑂((∆𝑡𝑡)2) + 𝑂𝑂((∆𝑥𝑥)4). 

It can be observed that the truncation error  𝑇𝑇𝑖𝑖,𝑛𝑛  goes to zero as both temporal and spatial intervals go to zero. 

That is 𝑇𝑇𝑖𝑖,𝑛𝑛 → 0 as ∆𝑡𝑡 → 0 and   ∆𝑥𝑥 → 0.  It shows that the FTCS scheme is consistent with partial differential 

equation (3). 

2.3 Backward Time Centered Space Implicit Scheme 

To discretize equation (3), we substitute the backward difference approximation for the first partial derivative 

and central difference approximation for the second partial derivative, 

𝜕𝜕𝑖𝑖
𝑛𝑛− 𝜕𝜕𝑖𝑖

𝑛𝑛−1

∆𝜕𝜕
  =  𝐷𝐷 �𝜕𝜕𝑖𝑖−1

𝑛𝑛 −2 𝜕𝜕𝑖𝑖
𝑛𝑛+ 𝜕𝜕𝑖𝑖+1

𝑛𝑛  
∆𝑥𝑥2

� + 𝑂𝑂(∆𝑡𝑡) + 𝑂𝑂(∆𝑥𝑥2)                                            (12) 

To present the system of equations in a simple manner let us drop the truncation error terms from (12) and 

rearrange the resulting equation to get 

𝑢𝑢𝑖𝑖𝑛𝑛−1 = − 𝑟𝑟𝑢𝑢𝑖𝑖−1𝑛𝑛 + (1 + 2𝑟𝑟)   𝑢𝑢𝑖𝑖𝑛𝑛 − 𝑟𝑟𝑢𝑢𝑖𝑖+1𝑛𝑛                                                           (13) 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

145 
 

Here in (12) and (13), we denote  𝑖𝑖  is an index number taking the values 𝑖𝑖 = 1, 2, 3, . . ., 𝑁𝑁 − 1.  

If the Dirichlet type of boundary condition are given i.e., values of the end points  𝑢𝑢0  
𝑛𝑛  and 𝑢𝑢𝑁𝑁 

𝑛𝑛  are given, then 

(13) can be reduced into a compact form as 

𝑢𝑢𝑛𝑛−1 = (𝐼𝐼 − 𝑟𝑟𝑟𝑟)𝑢𝑢𝑛𝑛,                                                                       (14) 

Here in (14), both 𝑢𝑢𝑛𝑛−1  and  𝑢𝑢𝑛𝑛 represent (𝑁𝑁 − 1) dimensional column vectors. That is,  𝑢𝑢𝑛𝑛−1 =

 [𝑢𝑢1𝑛𝑛−1,  𝑢𝑢2𝑛𝑛−1, . . . ,  𝑢𝑢𝑁𝑁−1𝑛𝑛−1]𝑇𝑇 and  𝑢𝑢𝑛𝑛 =  [(𝑢𝑢1𝑛𝑛 −  𝑟𝑟𝑢𝑢0𝑛𝑛),  𝑢𝑢2𝑛𝑛, . . . ,  𝑢𝑢𝑁𝑁−2𝑛𝑛 , ( 𝑢𝑢𝑁𝑁−1𝑛𝑛 −   𝑢𝑢𝑁𝑁𝑛𝑛)]𝑇𝑇 

Now up on substitution of the Fourier mode (11) in (13) yields  1 𝜆𝜆⁄   =  −𝑟𝑟𝑒𝑒−(𝑛𝑛𝑛𝑛 𝑝𝑝⁄ ) +  (1 + 2𝑟𝑟) – 𝑟𝑟𝑒𝑒(𝑛𝑛𝑛𝑛 𝑝𝑝⁄ ) or 

equivalently  𝜆𝜆 =   {1 +  2𝑟𝑟[1 −  cos (𝜋𝜋 𝑝𝑝⁄ )]}−1 and thus |𝜆𝜆| ≤ 1. This indicates that the BTCS scheme is 

unconditionally stable. Further, the advantage of this scheme is it removes the stability limitation associated 

with the diffusion operator. The disadvantage is that the problem becomes more expensive to solve numerically 

[3]. 

The truncation error  𝑇𝑇𝑖𝑖,𝑛𝑛  for the BTCS solution of the diffusion equation is [14] 

𝑇𝑇𝑖𝑖,𝑛𝑛 =  �
𝜕𝜕𝑢𝑢
𝜕𝜕𝑡𝑡

 −  𝐷𝐷
𝜕𝜕2𝑢𝑢
𝜕𝜕𝑥𝑥2

�
𝑖𝑖,𝑛𝑛

+ 
1
2
∆𝑡𝑡 �

𝜕𝜕2𝑢𝑢
𝜕𝜕𝑡𝑡2

�
𝑖𝑖,𝑛𝑛
−  

1
12

(∆𝑥𝑥)2 �
𝜕𝜕4𝑢𝑢
𝜕𝜕𝑥𝑥4

�
𝑖𝑖,𝑛𝑛

+ 𝑂𝑂((∆𝑡𝑡)2) + 𝑂𝑂((∆𝑥𝑥)4). 

It can be observed that the truncation error  𝑇𝑇𝑖𝑖,𝑛𝑛  goes to zero as both temporal and spatial intervals go to zero. 

That is 𝑇𝑇𝑖𝑖,𝑛𝑛 → 0 as ∆𝑡𝑡 → 0 and   ∆𝑥𝑥 → 0.  It shows that the BTCS scheme is consistent with partial differential 

equation (3). 

2.4 Crank – Nicolson Scheme 

As we have already seen both the FTCS and BTCS schemes have temporal truncation errors of order  𝑂𝑂(∆𝑡𝑡). 

However, the Crank – Nicolson scheme has the error of order  𝑂𝑂(∆𝑡𝑡2). Also the Crank – Nicolson scheme is not 

significantly more difficult to implement than the BTCS schemes. Further, the Crank – Nicolson scheme has 

significant advantages whenever the time – accurate solutions play an important role. The Crank – Nicolson 

scheme, like BTCS, is also implicit and unconditionally stable.  

The Crank – Nicolson scheme approximates equation (3) using central differences of time intervals. The spatial 

derivatives are estimated by the average of their values at time steps 𝑛𝑛   and   𝑛𝑛 + 1  as 

𝜕𝜕𝑖𝑖
𝑛𝑛+1− 𝜕𝜕𝑖𝑖

𝑛𝑛

∆𝜕𝜕
=  𝐷𝐷

2
�𝜕𝜕𝑖𝑖+1

𝑛𝑛+1−2𝜕𝜕𝑖𝑖
𝑛𝑛+1+ 𝜕𝜕𝑖𝑖−1

𝑛𝑛+1

(∆𝑥𝑥)2
+ 𝜕𝜕𝑖𝑖+1

𝑛𝑛 −2𝜕𝜕𝑖𝑖
𝑛𝑛+ 𝜕𝜕𝑖𝑖−1

𝑛𝑛

(∆𝑥𝑥)2
� + 𝑂𝑂((∆𝑡𝑡)2) + 𝑂𝑂((∆𝑥𝑥)2)                            (15) 

The Crank – Nicolson scheme has a truncation error of temporal order  𝑂𝑂(∆𝑡𝑡2) and spatial order   𝑂𝑂(∆𝑥𝑥2). Drop 

the truncation error terms from (15) and rearranging the terms so that values of 𝑢𝑢 at time 𝑛𝑛 are on the left and 

values of  𝑢𝑢  at time  𝑛𝑛 + 1 are on the right gives 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

146 
 

𝑟𝑟𝑢𝑢𝑖𝑖−1𝑛𝑛  + (2 − 2𝑟𝑟)𝑢𝑢𝑖𝑖𝑛𝑛 + 𝑟𝑟𝑢𝑢𝑖𝑖+1𝑛𝑛  = −𝑟𝑟𝑢𝑢𝑖𝑖−1𝑛𝑛+1 + (2 + 2𝑟𝑟)𝑢𝑢𝑖𝑖𝑛𝑛+1  − 𝑟𝑟𝑢𝑢𝑖𝑖+1𝑛𝑛+1                                           (16) 

The left hand side of (16) contains three known values and the right hand contains three unknowns. This scheme 

generates a set of  (𝑁𝑁 − 1) linear equation and those have to be solved at each time level. Applying (16) for all 

the internal mesh points at 𝑖𝑖 = 1,2,3, . . . ,𝑁𝑁 − 1 and using the boundary conditions  𝑥𝑥 = 0 and 𝑥𝑥 = 𝐿𝐿 we obtain a 

tridiagonal set of linear algebraic equations. These equations have to be solved at each time level.  

The compact form of such tridiagonal set of linear algebraic equations can be written as 

𝑢𝑢𝑛𝑛+1 = (2𝐼𝐼 − 𝑟𝑟𝑟𝑟)−1(2𝐼𝐼 + 𝑟𝑟𝑟𝑟)𝑢𝑢𝑛𝑛                                                               (17) 

The system of equations (17) can be solved very efficiently.  Also its unconditional stability can be shown by 

substituting (11) into (16). Thus,  

                          2𝜆𝜆 �1 +  𝑟𝑟 �1 −  cos �𝑛𝑛
𝑝𝑝
��� = 2�1 −  𝑟𝑟 �1 − cos �𝑛𝑛

𝑝𝑝
���, 

                                                   𝜆𝜆 =  
1 −  𝑟𝑟�1 − cos �𝜋𝜋𝑝𝑝��

1 + 𝑟𝑟�1 −cos�𝜋𝜋𝑝𝑝��
 

|𝜆𝜆|  ≤ 1 

The truncation error  𝑇𝑇𝑖𝑖,𝑛𝑛  for the Crank – Nicolson solution of the diffusion equation is [14] 

𝑇𝑇𝑖𝑖,𝑛𝑛 =  �
𝜕𝜕𝑢𝑢
𝜕𝜕𝑡𝑡

 −  𝐷𝐷
𝜕𝜕2𝑢𝑢
𝜕𝜕𝑥𝑥2

�
𝑖𝑖,𝑛𝑛

+  
∆𝑡𝑡
2
𝜕𝜕
𝜕𝜕𝑡𝑡
�
𝜕𝜕𝑢𝑢
𝜕𝜕𝑡𝑡

 −  𝐷𝐷
𝜕𝜕2𝑢𝑢
𝜕𝜕𝑥𝑥2

�
𝑖𝑖,𝑛𝑛

+ 
(∆𝑡𝑡)2

6
�
𝜕𝜕3𝑢𝑢
𝜕𝜕𝑡𝑡3

�
𝑖𝑖,𝑛𝑛
−   

1
12

(∆𝑥𝑥)2 �
𝜕𝜕4𝑢𝑢
𝜕𝜕𝑥𝑥4

�
𝑖𝑖,𝑛𝑛

+ 𝑂𝑂((∆𝑡𝑡)3)

+ 𝑂𝑂((∆𝑥𝑥)3). 

It can be observed that the truncation error  𝑇𝑇𝑖𝑖,𝑛𝑛  goes to zero as both temporal and spatial intervals go to zero. 

That is 𝑇𝑇𝑖𝑖,𝑛𝑛 → 0 as ∆𝑡𝑡 → 0 and   ∆𝑥𝑥 → 0.  It shows that the Crank – Nicolson scheme is consistent with partial 

differential equation (3). 

3. Numerical Results  

Consider the special case of the diffusion equation (3) with  𝐷𝐷 = 1 to obtain  

𝜕𝜕𝑢𝑢 𝜕𝜕𝑡𝑡⁄ =  𝜕𝜕2𝑢𝑢 𝜕𝜕𝑥𝑥2⁄                                                                         (18) 

The numerical solutions of (18) together with appropriate initial and boundary conditions are drawn using finite 

difference schemes (i)  FTCS (ii) BTCS and (iii) Crank – Nicolson. The details are given in what follows from 

Sections 3.1 to 3.3.  In our computations 0.1 and 0.125 are assigned to the spatial interval  ∆𝑥𝑥 and the values 0.4 

and 0.625 for parameter  𝑟𝑟  are used.  



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

147 
 

3.1 Problem 1 

In this problem, consider (18) together with the following initial and boundary conditions  

𝑢𝑢(𝑥𝑥, 0) = sin(𝜋𝜋𝑥𝑥) ,      0 < 𝑥𝑥 < 1,                                                           (19) 

𝑢𝑢(0, 𝑡𝑡) = 𝑢𝑢(1, 𝑡𝑡) = 0,     𝑡𝑡 ≥ 0.                                                           (20) 

However the analytical solution of (18) with the initial boundary value problem (19) and (20) can be easily 

found using the method of separation of variables as  

𝑢𝑢(𝑥𝑥, 𝑡𝑡) =  𝑒𝑒−𝑛𝑛2𝜕𝜕 sin(𝜋𝜋𝑥𝑥).                                                                        (21) 

The numerical solutions and errors of (18) together with (19) and (20) are illustrated in the following plots.  

 

Figure 1: Numerical solutions and errors of (18) together with (19) and (20) obtained using the three different 

schemes when  𝑟𝑟 = 0.4 and  ∆𝑥𝑥 = 0.125. 

0 0.2 0.4 0.6 0.8 1

0
0.5

1
0

0.5

1

x

(c) BTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-0.02

0

0.02

x

(d) Error of BTCS Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
0

0.5

1

x

(e) Crank-Nicholson Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-5

0

5

x 10
-3

x

(f) Error of Crank-Nicolson Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
0

0.5

1

x

(a) FTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
0

0.005

0.01

x

(b) Error of FTCS Scheme

t

E
rr

o
r



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

148 
 

 

Figure 2: Numerical solutions and errors of (18) together with (19) and (20) obtained using the three different 

schemes when  𝑟𝑟 = 0.625 and  ∆𝑥𝑥 = 0.1. 

3.2 Problem 2 

In this problem, consider (18) together with the following initial and boundary conditions  

𝑢𝑢(𝑥𝑥, 0) = cos(𝜋𝜋𝑥𝑥) ,      0 ≤ 𝑥𝑥 ≤ 1,                                                    (22) 

𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥

(0, 𝑡𝑡) = 𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥

(1, 𝑡𝑡) = 0,    0 <  𝑡𝑡 ≤ 1.                                               (23) 

 

0 0.2 0.4 0.6 0.8 1

0
0.5

1
0

0.5

1

x

(c) BTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-0.02

0

0.02

x

(d) Error of BTCS Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
0

0.5

1

x

(e) Crank-Nicholson Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-5

0

5

x 10
-3

x

(f) Error of Crank-Nicolson Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-5

0

5

x 10
8

x

(a) FTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-5

0

5

x 10
8

x

(b) Error of FTCS Scheme

t

E
rr

o
r



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

149 
 

However, the analytical solution of (18) with the initial boundary value problem (22) and (23) can be easily 

found using the method of separation of variables as  

𝑢𝑢(𝑥𝑥, 𝑡𝑡) =  𝑒𝑒−𝑛𝑛2𝜕𝜕 cos(𝜋𝜋𝑥𝑥).                                (24) 

The numerical solutions and errors of (18) together with (22) and (23) are illustrated in the following plots: 

 

Figure 3: Numerical solutions and errors of (18) together with (22) and (23) obtained using the three different 

schemes when  𝑟𝑟 = 0.4 and  ∆𝑥𝑥 = 0.125. 

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

0

1

x

(a) FTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

0

1

x

(b) Error of FTCS Scheme

t

E
r
r
o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

0

1

x

(c) BTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

-0.5

0

x

(d) Error of BTCS Scheme

t

E
r
r
o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

0

1

x

(e) Crank-Nicolson Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

-0.5

0

x

(f) Error of Crank-Nicolson Scheme

t

E
r
r
o
r



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

150 
 

 

Figure 4: Numerical solutions and errors of (18) together with (22) and (23) obtained using the three different 

schemes when  𝑟𝑟 = 0.625 and  ∆𝑥𝑥 = 0.1. 

3.3  Problem 3 

In this problem, consider (18) together with the following initial and boundary conditions  

𝑢𝑢(𝑥𝑥, 0) = cos(𝑥𝑥) ,      0 < 𝑥𝑥 < 1                                                                        (25) 

𝜕𝜕𝜕𝜕
𝜕𝜕𝑥𝑥

(0, 𝑡𝑡) = 𝑢𝑢 �𝑛𝑛
2

, 𝑡𝑡� = 0,     𝑡𝑡 ≥ 0.                                                                                  (26) 

However the analytical solution of (18) with the initial boundary value problem (25) and (26) can be easily 

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-5

0

5

x 10
25

x

(a) FTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-5

0

5

x 10
25

x

(b) Error of FTCS Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

0

1

x

(c) BTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

-0.5

0

x

(d) Error of BTCS Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

0

1

x

(e) Crank-Nicolson Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.5

1
-1

-0.5

0

x

(f) Error of Crank-Nicolson Scheme

t

E
rr

o
r



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

151 
 

found using the method of separation of variables as  

𝑢𝑢(𝑥𝑥, 𝑡𝑡) =  𝑒𝑒−𝜕𝜕 cos(𝑥𝑥).                                                                              (27) 

The numerical solutions and errors of (18) together with (25) and (26) are illustrated in the following plots: 

 

Figure 5: Numerical solutions and errors of (18) together with (25) and (26) obtained using the three different 

schemes when  𝑟𝑟 = 0.4 and  ∆𝑥𝑥 = 0.125. 

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
0

0.5

1

x

(a) FTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
-5

0

5

x 10
-5

x

(b) Error of FTCS Scheme

t

E
r
r
o
r

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
0

0.5

1

x

(c) BTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
-5

0

5

x 10
-3

x

(d) Error of BTCS Scheme

t

E
r
r
o
r

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
0

0.5

1

x

(e) Crank-Nicolson Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
-2

0

2

x 10
-3

x

(f) Error of Crank-Nicolson Scheme

t

E
r
r
o
r



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

152 
 

 

Figure 6: Numerical solutions and errors of (18) together with (25) and (26) obtained using the three different 

schemes when  𝑟𝑟 = 0.625 and  ∆𝑥𝑥 = 0.1 

Performance of the three finite difference schemes can be evaluated up on comparing some factors viz. (i) 

stability of numerical schemes and (ii) error. Stability criteria and numerical solutions show that the three 

schemes work well and each produces reasonable results in case of all test examples. The numerical results are 

obtained using Matlab programming and the outcomes are illustrated in Figures 1 – 6.  

The selection of spatial interval ∆𝑥𝑥 = 0.125 and the time step   ∆𝑡𝑡 = 0.00625, leads to  𝑟𝑟 =  𝐷𝐷 [∆𝑡𝑡 (∆𝑥𝑥)2⁄ ] =

0.4 <  1 2⁄ . The stability analysis FTCS scheme guarantees that if  𝑟𝑟 ≤  1 2⁄   then the condition is stable. This 

stability criterion is well supported by the solutions provided by all the three schemes in case of all problems. 

The same fact has observed in Figures 1, 3 and 5.  

But if we decrease the spatial interval to ∆𝑥𝑥 = 0.1 for better resolution and considering the same time 

step   ∆𝑡𝑡 = 0.00625, leads to  𝑟𝑟 =  𝐷𝐷 [∆𝑡𝑡 (∆𝑥𝑥)2⁄ ] = 0.625 >  1 2⁄ . The stability analysis of FTCS scheme 

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
0

0.5

1

x

(a) FTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
0

2

4

x 10
-4

x

(b) Error of FTCS Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
0

0.5

1

x

(c) BTCS Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
-2

0

2

x 10
-3

x

(d) Error of BTCS Scheme

t

E
rr

o
r

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
0

0.5

1

x

(e) Crank-Nicolson Scheme

t

u

0 0.2 0.4 0.6 0.8 1

0
0.51

1.52
-1

0

1

x 10
-3

x

(f) Error of Crank-Nicolson Scheme

t

E
rr

o
r



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

153 
 

predicts that if  𝑟𝑟 >  1 2⁄   then the condition is unstable. Hence, the FTCS scheme provides unstable solution 

and as a result the error blows up as shown in Figures 2a – 2b and 4a – 4b. However, the BTCS and Crank – 

Nicolson schemes work quite well and provide stable solutions as shown in Figures 2c – 2f and 4c – 4f.  

With the same spatial interval ∆𝑥𝑥 = 0.1 and time step  ∆𝑡𝑡 = 0.00625, all the three schemes including FTCS 

produce stable solutions. Despite of the fact that 𝑟𝑟 = 0.625 doesn’t strictly satisfy the stability condition 𝑟𝑟 ≤

1 2⁄  as shown in Figures 6a – 6f. This implies that the stability of the FTCS scheme is not a necessary condition, 

but only sufficient one. Besides, if it converges then the accuracy of FTCS may be better than that of BTCS 

scheme, but generally not better than that of the Crank – Nicolson scheme. The observations regarding stability 

of the solutions made from the figures are summarized in a tabular form as below: 

Table 1: Stability of the solutions of one dimensional diffusion equation 

 

Problems 

Finite 

Difference 

Schemes 

Stability condition 

for 𝒓𝒓 ≤  𝟎𝟎.𝟓𝟓 

Stability condition 

for 𝒓𝒓 >  0.5 

 

Problem 1 

FTCS Stable (Fig. 1a) Unstable (Fig. 2a) 

BTCS Stable (Fig. 1c) Stable  (Fig. 2c) 

CN Stable (Fig. 1e) Stable (Fig. 2e) 

 

Problem 2 

FTCS Stable (Fig. 3a) Unstable (Fig. 4a) 

BTCS Stable (Fig. 3c) Stable (Fig. 4c) 

CN Stable (Fig. 3e) Stable (Fig. 4e) 

 

Problem 3 

FTCS Stable (Fig. 5a) Stable (Fig. 6a) 

BTCS Stable (Fig. 5c) Stable (Fig. 6c) 

CN Stable (Fig. 5e) Stable (Fig. 6e) 

 

4. Conclusions 

In this study, we have discussed application of numerical schemes on one dimensional diffusion equation.  It is 

observed from numerical computation that all the three schemes worked well according to the stability criteria 

and each scheme produced reasonable approximation for the density variable   𝑢𝑢(𝑥𝑥, 𝑡𝑡). The three schemes are 

compared based on the results of the three test problems. The comparison indicates that the approximate 

solution provided by Crank – Nicolson scheme is better than the other two schemes. Hence, Crank – Nicolson 

scheme is recommended for solving one – dimensional equation for a better approximation. 

References 

[1]  Bergara A. Finite – Difference Numerical methods of Partial Differential equation in Finance with 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 26, No  3, pp 140-154 

154 
 

Matlab, EPV/EHU, http://www.ehu.es/aitor  

[2]  Evans G., Blackdge J. and Yardley P. (2000). Numerical Methods for Partial Differential Equations, 

Springer – Verlag, Berlin, Heidelberg, New York. ISBN 3-540-76125-X. 

[3]  Gerald W. (2011). Finite – Difference Approximations to the Heat Equation, Portland State University, 

Portland, Oregon. 

[4]  H. Saberi Najafi (2008). Solving One – Dimensional Advection – Dispersion with Reaction Using 

Some Finite – Difference Methods, Applied Mathematical Sciences, Vol. 2, Pp. 2611 – 18. 

[5]  Karantonis A. (2001). Numerical Solution of Reaction – Diffusion Equations by the Finite Difference 

Method. 

[6]  Kværnø A. (2009). Numerical Mathematics, Lecture Notes, NTNU, University of Oslo. 

[7]   Owren B. (2011). Numerical Solution of Partial Differential Equations with Finite Difference 

Methods, Lecture Notes, NTNU, University of Oslo. 

[8]   Langtangen H. P. (2013). Finite difference methods for diffusion processes, University of Oslo. 

[9]  Shahraiyni Taheri H. and Ataie B. (2009). Comparison of Finite Difference Schemes for Water Flow in 

Unsaturated Soils, International Journal of Aerospace and Mechanical Engineering, Vol. 3(1), Pp.1 – 5.  

[10]  Sweilam N. H., Khader M. M. and Mahdy A. M. S. (2012). Crank – Nicolson Finite Difference 

Method for Solving Time – Fractional Diffusion Equation, Journal of Fractional Calculus and 

Applications, Vol. 2, No. 2, Pp. 1 – 9, http://www.fcaj.webs.com 

[11]  Thongmoon M. and McKibbin R.  (2006). A Comparison of Some Numerical Methods for the 

Advection – Diffusion Equation, Res. Lett. Inf. Math. Sci., Vol. 10, Pp 49 – 62. 

http://iims.massey.ac.nz/research/letters/ 

[12]  Tsegaye Simon (2013). Numerical Simulation of Diffusion – Reaction Equations: Application from 

River Pollution Model. Hawassa University, Hawassa, Ethiopia (Unpublished M. Sc. Thesis). 

[13]  Won Y., Wenwu C., Tae – Sang C. and John M. (2005). Applied Numerical Methods Using 

MATLAB, ISBN 0-471-69833-4 (cloth), John Wiley & Sons, Inc., Hoboken, New Jersey. 

[14]  Zana S. (2014). Numerical Solution of Diffusion Equation in One Dimension, Eastern Mediterranean 

University, Gazimağusa, North Cyprus. 

 

 

 

http://www.ehu.es/aitor
http://www.fcaj.webs.com/
http://iims.massey.ac.nz/research/letters/