Efficient Finite Difference Methods for the Numerical Analysis of One-Dimensional Heat Equation

Md. Shahadat Hossain Mojumder¹, Md. Nazmul Haque^1,2, Md. Joni Alam^1*

¹Department of Mathematics, Comilla University, Kotbari, Bangladesh.
²Department of Mathematical Sciences, Kent State University, Kent, OH, USA.

DOI: 10.4236/jamp.2023.1110204                

Abstract

In this paper, we investigate and analyze one-dimensional heat equation with appropriate initial and boundary condition using finite difference method. Finite difference method is a well-known numerical technique for obtaining the approximate solutions of an initial boundary value problem. We develop Forward Time Centered Space (FTCS) and Crank-Nicolson (CN) finite difference schemes for one-dimensional heat equation using the Taylor series. Later, we use these schemes to solve our governing equation. The stability criterion is discussed, and the stability conditions for both schemes are verified. We exhibit the results and then compare the results between the exact and approximate solutions. Finally, we estimate error between the exact and approximate solutions for a specific numerical problem to present the convergence of the numerical schemes, and demonstrate the resulting error in graphical representation.

Keywords

Explicit Scheme, Implicit Scheme, C-N Scheme, CFL Condition

1. Introduction

The heat equation is a parabolic second-order linear partial differential equation that describes the distribution of heat (or temperature variation) in a given region over time. It demonstrates the fact of how heat flows through an object. For example, if a heated body is placed in a box of cold water, the body’s temperature will decrease. Eventually, the temperature in the box will equalize at a time if no other heat source is applied there.

To solve the one-dimensional heat equation, many effective methods can be used, such as separation of variables [1] [2] [3] [4] , Fourier series methods [5] [6] , domain decomposition method [7] [8] [9] , finite difference method [4] [10] [11] [12] [13] , finite element method [14] [15] , finite volume method [16] [17] and many more. In this paper, we use the finite difference method to solve and analyze the one-dimensional time-dependent heat equation. The finite difference method (FDM) is a powerful numerical technique for obtaining the approximate solution to the initial boundary value problems. It is easy to solve PDEs using this method [18] . Taylor series approach is used to find the solutions of a partial differential equation in this method [19] . The technique of this method is that all the derivatives in a partial differential equation are replaced by their finite difference approximations, and the PDEs are converted to a set of linear algebraic equations. Any iterative procedure can be used to solve this system of linear equations. Then, the solution of this system of linear equations is referred to as the solution of the PDE [1] . The finite difference method has been extensively used by researchers over the past years due to the fact that most of the calculations of this method can be carried out on the computer, and hence the solution is easy to obtain [20] [21] .

The aim of this work is to investigate a one-dimensional heat equation by using explicit Forward Time Centered Space (FTCS) and implicit Crank-Nicolson (CN) finite difference schemes. The rest of this paper is organized as follows. We present the one-dimensional heat equation with Neumann boundary conditions in the next Section. In Section 3, a brief introduction of thermal diffusivity is reported. The derivatives of the temperature distribution $u\left(x\mathrm{,}t\right)$ are approximated by finite differences using the Taylor series in Section 4. The explicit and Crank-Nicolson schemes are developed for the one-dimensional heat equation in Section 5. The stability criterion is discussed, and the stability conditions for both schemes are verified in Section 6. We have exhibited the results and finally compared the results between the exact and the approximate solutions. The algorithm is executed by standard computer programming languages in Sections 6 and 7. We estimate error between the exact and the approximate solutions for a specific numerical problem and show it graphically in Section 7. Finally, a short conclusion is reported in Section 8.

2. Model Equations

In our model problem, we consider the one-dimensional time-dependent heat conduction equation of the form [3] [11] [12] [21] ,

$\frac{\partial u}{\partial t}-{\alpha }^2\frac{{\partial }^{2}u}{\partial x^2}=\mathrm{0,}\left(x\mathrm{,}t\right)\in \Omega \times \Gamma$ (1)

where, $\Omega =\left[\mathrm{0,}L\right]$ and $\Gamma =\left[\mathrm{0,}T\right]$ , with initial condition

$u\left(x\mathrm{,0}\right)=u_0\left(x\right)\text{for}0\le x\le L\mathrm{,}$ (2)

and boundary conditions

$u\left(\mathrm{0,}t\right)=f\left(x\right)\text{for}0\le t\le T\mathrm{,}$ (3)

$u\left(L\mathrm{,}t\right)=g\left(x\right)\text{for}0\le t\le T\mathrm{,}$ (4)

where, $u=u\left(x\mathrm{,}t\right)$ is the temperature of the substance, α² is called the thermal diffusivity of the substance. Equations (1)-(4) together are called initial boundary-value problem.

3. Thermal Diffusivity (α²)

Thermal diffusivity is a significant physical quantity that describes the heat transfer rate through a medium. Different materials have different thermal diffusivity. It tells us how fast heat transfer occurs within a material from the hotter to the colder side. Thus, thermal diffusivity measures how quickly a material reacts to temperature changes.

It is defined as the ratio between the thermal conductivity and the product of the density and the specific heat capacity of the body with the pressure held constant, mathematically,

${\alpha }^2=\frac{K_0}{\rho c_p}\mathrm{,}$ (5)

where,

K₀ = thermal conductivity of the substance;

ρ = mass per unit volume of the substance;

c_p = substance’s specific heat capacity.

4. Finite Difference Discretization

In order to describe the finite difference techniques, let us consider the rectangular domain $\Omega \times \Gamma$ into a finite number of nodes $\left(x_i\mathrm{,}{t}_j\right)$ , where Ω is spatial domain and Γ is temporal domain.

Now, we develop the approximations of finite difference formulas according to [22] .

Let I be an open interval and $u\mathrm{<mo>\; :\; </mo>}I\to ℝ$ be differentiable $\left(n+1\right)$ times. Now, we expand u in Taylor series about a point x, evaluate at $x+h$ and at $x-h$ , respectively as follows,

$u\left(x+h\right)=u\left(x\right)+h{u}^{\prime }\left(x\right)+\frac{h^2}{\mathrm{2\; <mo>\; !\; </mo>}}{u^{\prime }}^{\prime }\left(x\right)+\frac{h^3}{\mathrm{3\; <mo>\; !\; </mo>}}{{u^{\prime }}^{\prime }}^{\prime }\left(x\right)+\cdots$ (6)

$u\left(x-h\right)=u\left(x\right)-h{u}^{\prime }\left(x\right)+\frac{h^2}{\mathrm{2\; <mo>\; !\; </mo>}}{u^{\prime }}^{\prime }\left(x\right)-\frac{h^3}{\mathrm{3\; <mo>\; !\; </mo>}}{{u^{\prime }}^{\prime }}^{\prime }\left(x\right)+\cdots$ (7)

Rearranging (6), we get

$u^{\prime }\left(x\right)\approx \frac{u\left(x+h\right)-u\left(x\right)}{h}+O\left(h\right)\mathrm{,}$ (8)

which is forward difference formula for $u^{\prime }\left(x\right)$ .

Similarly, from (7), we get the backward difference formula as

$u^{\prime }\left(x\right)\approx \frac{u\left(x-h\right)-u\left(x\right)}{h}+O\left(h\right)\mathrm{,}$ (9)

By adding (6) and (7), we have

$\frac{u\left(x+h\right)-2u\left(x\right)+u\left(x-h\right)}{h^2}={u^{\prime }}^{\prime }\left(x\right)+\frac{h^2}{12}\left(u^{iv}\left(x\right)\right)+\cdots$

This gives us the finite difference formula for the second derivative with second order accuracy,

${u^{\prime }}^{\prime }\left(x\right)=\frac{u\left(x+h\right)-2u\left(x\right)+u\left(x-h\right)}{h^2}+O\left(h^2\right)$ (10)

Now, we can define $\frac{\partial u}{\partial t}$ as,

$\frac{\partial u}{\partial t}=\underset{k\to \infty }{lim}\frac{u\left(x\mathrm{,}t+k\right)-u\left(x\mathrm{,}t\right)}{k}\mathrm{.}$

Thus, we get,

$\frac{\partial u}{\partial t}\approx \frac{u\left(x\mathrm{,}t+k\right)-u\left(x\mathrm{,}t\right)}{k}\mathrm{,}$ (11)

which is called a forward difference approximation for $\frac{\partial u}{\partial t}$ .

In order to find an approximation to the second derivative, $\frac{{\partial }^{2}u}{\partial x^2}$ ; we start with the forward difference

$\frac{\partial u}{\partial x}\approx \frac{u\left(x+h\mathrm{,}t\right)-u\left(x\mathrm{,}t\right)}{h}\mathrm{.}$ (12)

We need to approximate the terms in the numerator. It is customary to use a backward difference approximation. This is given by letting $h\to -h$ in the forward difference form,

$\frac{\partial u}{\partial x}\approx \frac{u\left(x\mathrm{,}t\right)-u\left(x-h\mathrm{,}t\right)}{k}$ (13)

applying (13) to $\frac{\partial u}{\partial x}$ and evaluated at $x=x$ and $x=x+h$ we have,

$\frac{\partial u}{\partial x}\left(x\mathrm{,}t\right)\approx \frac{u\left(x\mathrm{,}t\right)-u\left(x-h\mathrm{,}t\right)}{h}\mathrm{,}$ (14)

and

$\frac{\partial u}{\partial x}\left(x+h\mathrm{,}t\right)\approx \frac{u\left(x+h\mathrm{,}t\right)-u\left(x\mathrm{,}t\right)}{h}\mathrm{.}$ (15)

5. Schemes for the Heat Equation

In order to develop the schemes for solving the one-dimensional heat Equation (1), both the time and space derivatives are required to replace by their finite differences. To perform this task, we need to specify both the time and spatial locations of the u values in the finite difference formulas [1] . Therefore, we need to introduce subscript i and j to designate the spatial and the time steps of the discrete solution, respectively. Hence, we can rewrite Equation (11) under these notations as

$\frac{\partial u}{\partial t}=\frac{u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}}{k}\mathrm{,}$ (16)

and Equation (10) as

$\frac{{\partial }^{2}u}{\partial x^2}=\frac{u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}}{h^2}\mathrm{.}$ (17)

For Crank-Nicolson scheme, the derivation of Equation (10) as

$\frac{{\partial }^{2}u}{\partial x^2}=\frac{u_{i+\mathrm{1,}j+1}-2u_{i\mathrm{,}j+1}+u_{i-\mathrm{1,}j+1}}{h^2}\mathrm{.}$ (18)

The second order partial derivative $\frac{{\partial }^{2}u}{\partial x^2}$ can also be expressed by using the central difference at $x_i$ and $t_j$ and for any constant $\lambda$ ( $0\le \lambda \le 1$ ),

${\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{i\mathrm{,}j}=\left(1-\lambda \right){\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{i\mathrm{,}j}+\lambda {\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{i\mathrm{,}j+1}$ (19)

${\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{i\mathrm{,}j}=\left(1-\lambda \right)\frac{u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}}{h^2}+\lambda \frac{u_{i+\mathrm{1,}j+1}-2u_{i\mathrm{,}j+1}+u_{i-\mathrm{1,}j+1}}{h^2}\mathrm{,}$ (20)

where, the notation $u_{i\mathrm{,}j}=u\left(x_i\mathrm{,}{t}_j\right)$ is used to represent the numerical solution on a finite rod (substance) $x\in \Omega$ and $t\in \Gamma$ .

5.1. Forward Time Centered Space (FTCS) Explicit Scheme

In this section, we consider the Forward in Time and Central in Space (FTCS) scheme, where we replace the time derivative in (1) by the forward difference scheme (16) and the space derivative in (1) by the central difference scheme (17). This yields,

$\frac{u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}}{k}-{\alpha }^2\frac{u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}}{h^2}=0$

$\text{or}\mathrm{,}\frac{u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}}{k}={\alpha }^2\frac{u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}}{h^2}$

$\text{or}\mathrm{,}{u}_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}=\frac{{\alpha }^{2}k}{h^2}\left(u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}\right)\mathrm{.}$ (21)

Now we use the Courant-Friedrichs-Lewy or CFL condition in Equation (21).

The CFL condition for one dimension heat equation is $r=\frac{{\alpha }^{2}k}{h^2}$ . By using this condition in Equation (21), we get

$u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}=r\left(u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}\right)$

$\text{or}\mathrm{,}{u}_{i\mathrm{,}j+1}=r{u}_{i+\mathrm{1,}j}-2r{u}_{i\mathrm{,}j}+r{u}_{i-\mathrm{1,}j}+u_{i\mathrm{,}j}$

$\text{or}\mathrm{,}{u}_{i\mathrm{,}j+1}=r{u}_{i-\mathrm{1,}j}+u_{i\mathrm{,}j}-2r{u}_{i\mathrm{,}j}+r{u}_{i+\mathrm{1,}j}$

$\therefore u_{i\mathrm{,}j+1}=r{u}_{i-\mathrm{1,}j}+\left(1-2r\right)u_{i\mathrm{,}j}+r{u}_{i+\mathrm{1,}j}\mathrm{.}$ (22)

Equation (22) is called the Forward Time, Centered Space or FTCS approximation to the heat equation. The FTCS solution can be extracted from the grid of Figure 1.

Figure 1. Forward time centered space (FTCS) explicit scheme.

We exhibit the nature of the FTCS scheme (22) in Figure 2 and Figure 3 for the one-dimensional heat Equation (1) by the standard computer programming languages with boundary conditions $u\left(\mathrm{0,}t\right)=u\left(\mathrm{1,}t\right)=0$ , and initial condition $u_0\left(x\right)=u\left(x\mathrm{,0}\right)=\sin \left(n\pi x\right)$ , $n=\mathrm{1,2}$ .

5.2. Crank-Nicolson Scheme

To develop the Crank-Nicolson scheme [23] for our model Equation (1), we put $\lambda =\frac{1}{2}$ in Equation (20), for each $i\mathrm{,}j=\mathrm{0,1,2,3},\cdots$, [24] , we have,

${\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{i\mathrm{,}j}=\frac{1}{2h^2}\left(u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}+u_{i+\mathrm{1,}j+1}-2u_{i\mathrm{,}j+1}+u_{i-\mathrm{1,}j+1}\right)\mathrm{,}$ (23)

Now, using (23) and (16) in Equation (1), we get

$\frac{u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}}{k}-\frac{{\alpha }^2}{2h^2}\left(u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}+u_{i+\mathrm{1,}j+1}-2u_{i\mathrm{,}j+1}+u_{i-\mathrm{1,}j+1}\right)=0$

$\Rightarrow u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}-\frac{k{\alpha }^2}{2h^2}\left(u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}+u_{i+\mathrm{1,}j+1}-2u_{i\mathrm{,}j+1}+u_{i-\mathrm{1,}j+1}\right)=0.$ (24)

For the simplicity of calculation, we use the CFL condition, $r=\frac{k{\alpha }^2}{h^2}$ , thus Equation (24) becomes,

$u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}-\frac{r}{2}\left(u_{i+\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i-\mathrm{1,}j}+u_{i+\mathrm{1,}j+1}-2u_{i\mathrm{,}j+1}+u_{i-\mathrm{1,}j+1}\right)=0$

$\Rightarrow u_{i\mathrm{,}j+1}=\frac{r{u}_{i+\mathrm{1,}j}+2\left(1+r\right)u_{i\mathrm{,}j}+r{u}_{i-\mathrm{1,}j}+r{u}_{i+\mathrm{1,}j+1}+r{u}_{i-\mathrm{1,}j+1}}{2\left(1+r\right)}\mathrm{.}$ (25)

Equation (25) is the Crank-Nicolson scheme for the one dimensional heat equation.

Setting $r=1$ , i.e., $k=\frac{h^2}{{\alpha }^2}$ then Equation (25) reduces to

$u_{i\mathrm{,}j+1}=\frac{1}{4}\left(u_{i+\mathrm{1,}j}+u_{i-\mathrm{1,}j}+4u_{i\mathrm{,}j}+u_{i+\mathrm{1,}j+1}+u_{i-\mathrm{1,}j+1}\right)\mathrm{.}$ (26)

(a)(b)

Figure 2. Temperature distribution and mesh plot for the explicit scheme in (a) and (b) are presented, respectively, for $r=0.5$ and spatial domain $x=0$ to $x=\pi$ . (a) Temperature distribution; (b) Temperature distribution in mesh.

(a)(b)

Figure 3. Temperature distribution and mesh plot for the explicit scheme in (a) and (b) are presented, respectively, for $r=0.5$ and spatial domain $x=0$ to $x=2\pi$ . (a) Temperature distribution; (b) Temperature distribution in mesh.

The CN solution (26) can be made sense from the grid of Figure 4. Moreover, We show the temperature distribution of the one-dimensional heat Equation (1) by the standard computer programming languages using the CN scheme (26) under the boundary conditions $u\left(0,t\right)=u\left(1,t\right)=0$ , and the initial conditions $u_0\left(x\right)=u\left(x,0\right)=\sin \left(n\pi x\right)$ , $n=1,2$ in Figure 5 and Figure 6, respectively.

6. Stability Analysis

Stability is one of the most important properties of a system. It is a measure of the performance of a system which refers to the ability of the system to perform satisfactorily. It tells us the fact of how good the performance of the system is, whether it is stable or unstable. In this section, we discuss the Von-Neumann stability analysis for the one-dimensional heat equation for the explicit and Crank-Nicolson Schemes. The stability analysis was based on the citation [20] [24] [25] [26] .

From the definition of the error in the numerical approximation, we get

${\xi }_{i\mathrm{,}j}=u_{i\mathrm{,}j}-U_{i\mathrm{,}j}\mathrm{,}$ (27)

where, $U_{i\mathrm{,}j}$ is the exact solution at grid point $\left(i\mathrm{,}j\right)$ and $u_{i\mathrm{,}j}$ is the approximate solution.

For linear differential equations with periodic boundary condition, the spatial variation of error may be expanded in a finite Fourier series, in the interval Ω as

$\xi \left(x\right)={\displaystyle \underset{m=1}{\overset{M}{\sum }}}{A}_m{\text{e}}^{i{k}_{m}x},$ (28)

where the wave number $k_m=\frac{\pi m}{L},m=1,2,\cdots ,M$ and $M=\frac{L}{h}$ , $i=\sqrt{-1}$ and x is the independent space variable, ${\text{e}}^{i{k}_{m}x}$ is the complex exponential.

The time dependence of the error is included by assuming the amplitude of error $A_m$ is a function of time. Since, the error tends to grow or decay exponentially with time, it is reasonable to assume that the amplitude varies exponentially with time; hence,

Figure 4. Crank-Nicolson implicit scheme.

(a)(b)

Figure 5. Temperature distribution and mesh plot for the Crank-Nicolson scheme in (a) and (b) are presented, respectively, for $r=1.0$ and spatial domain $x=0$ to $x=\pi$ . (a) Temperature distribution; (b) Temperature distribution in mesh.

(a)(b)

Figure 6. Temperature distribution and mesh plot for the Crank-Nicolson scheme in (a) and (b) are presented, respectively, for $r=1.0$ and spatial domain $x=0$ to $x=2\pi$ . (a) Temperature distribution; (b) Temperature distribution in mesh.

$\xi \left(x,t\right)={\xi }_{i,j}={\displaystyle \underset{m=1}{\overset{M}{\sum }}}{\text{e}}^{at}{\text{e}}^{i{k}_{m}x}.$ (29)

Since the difference equation for error is linear (the behavior of each term of the series is the same as the series itself), it is enough to consider the growth of error of a typical term:

${\xi }_m\left(x,t\right)={\xi }_{i,j}={\text{e}}^{at}{\text{e}}^{i{k}_{m}x}.$ (30)

The stability characteristics can be studied using just this form for the error with no loss in generality.

We now define the amplification factor as,

$G=\frac{{\xi }_{i\mathrm{,}j+1}}{{\xi }_{i\mathrm{,}j}}\mathrm{.}$ (31)

Then, the necessary and sufficient condition for the error to remain bounded, maintaining numerical stability is,

$\left|G\right|\le 1\text{or}-1\le G\le 1.$ (32)

6.1. Stability Criterion for FTCS Scheme

Now, we are going to show that Equation (32) holds for our FTCS numerical scheme (22). From Equation (30), we get by [20] [26] ,

${\xi }_{i,j}={\text{e}}^{at}{\text{e}}^{i{k}_{m}x}$ (33)

${\xi }_{i,j+1}={\text{e}}^{a\left(t+k\right)}{\text{e}}^{i{k}_{m}x}$ (34)

${\xi }_{i+1,j}={\text{e}}^{at}{\text{e}}^{i{k}_m\left(x+h\right)}$ (35)

${\xi }_{i-1,j}={\text{e}}^{at}{\text{e}}^{i{k}_m\left(x-h\right)}.$ (36)

Since, both the numerical solution $u_{i\mathrm{,}j}$ and the exact solution $U_{i\mathrm{,}j}$ satisfy Equation (22), therefore, the error ${\xi }_{i\mathrm{,}j}$ also follows the discretized PDE Equation (22). Thus we have,

${\xi }_{i\mathrm{,}j+1}={\xi }_{i\mathrm{,}j}+r\left({\xi }_{i+\mathrm{1,}j}-2{\xi }_{i\mathrm{,}j}+{\xi }_{i-\mathrm{1,}j}\right)\mathrm{.}$ (37)

Using Equations (33)-(36) in Equation (37), we get

${\text{e}}^{a\left(t+k\right)}{\text{e}}^{ik}{m}^x={\text{e}}^{at}{\text{e}}^{ik}{m}^x+r\left({\text{e}}^{at}{\text{e}}^{ik}{m}^{\left(x+h\right)}-2{\text{e}}^{at}{\text{e}}^{ik}{m}^x+{\text{e}}^{at}{\text{e}}^{ik}{m}^{\left(x-h\right)}\right)\mathrm{.}$

Divide by ${\text{e}}^{at}{\text{e}}^{ik}{m}^x$ to yield after simplification,

${\text{e}}^{ak}=1+r\left({\text{e}}^{ik}{m}^h+{\text{e}}^{-ik}{m}^h-2\right)$ (38)

We know that,

${\text{e}}^{i\theta }=\cos \theta +i\sin \theta \mathrm{,}$ (39)

and

${\text{e}}^{-i\theta }=\cos \theta -i\sin \theta \mathrm{.}$ (40)

Subtracting (39) from (40) we get,

${\text{e}}^{i\theta }-{\text{e}}^{-i\theta }=2i\sin \theta$

$\Rightarrow \sin \theta =\frac{{\text{e}}^{i\theta }-{\text{e}}^{-i\theta }}{2i}\mathrm{.}$ (41)

Putting $\theta =\frac{k_{m}h}{2}$ in (41), we get

$\sin \left(\frac{k_{m}h}{2}\right)=\frac{{\text{e}}^{\frac{i{k}_{m}h}{2}}-{\text{e}}^{-\frac{i{k}_{m}h}{2}}}{2i}\mathrm{,}$

squaring both sides, we get

${\sin }^2\left(\frac{k_{m}h}{2}\right)=-\left[\frac{{\text{e}}^{i{k}_{m}h}+{\text{e}}^{i{k}_{m}h}-2}{4}\right]$ (42)

From (38) and (42), we get

${\text{e}}^{ak}=1-4r{\sin }^2\left(\frac{k_{m}h}{2}\right)$ (43)

However, by Equations (33)-(36) in Equation (31), the last yield

$G={\text{e}}^{ak}$ (44)

Thus, by Equations (43), (44) plugged into (32), the condition for stability is given by

$\left|G\right|=\left|1-4r{\sin }^2\left(\frac{k_{m}h}{2}\right)\right|\le 1$ (45)

And from (45), we can write,

$1-4r{\sin }^2\left(\frac{k_{m}h}{2}\right)\ge -1$ (46)

And

$1-4r{\sin }^2\left(\frac{k_{m}h}{2}\right)\le -1$ (47)

For the above conditions, Equations (46) and (47), to hold for all m (and therefore all ${\sin }^2\left(\frac{k_{m}h}{2}\right)$ ), we need to consider (46) and (47) separately.

Equation (46) yields

$4r{\sin }^2\left(\frac{k_{m}h}{2}\right)\le 2$ (48)

Note that the term $4r{\sin }^2\left(\frac{k_{m}h}{2}\right)$ is always positive, i.e., always in the range [0, 1]. The worst case is when ${\sin }^2\left(\frac{k_{m}h}{2}\right)=1$ , therefore, Equation (48) yields

$r\le \frac{1}{2}$

$\Rightarrow r=\frac{{\alpha }^{2}k}{h^2}\le \frac{1}{2}$ (49)

Again Equation (47) yields

$4r{\sin }^2\left(\frac{k_{m}h}{2}\right)\ge 0$ (50)

Since ${\sin }^2\left(\frac{k_{m}h}{2}\right)$ range is [0, 1] Equation (50) yields

$r\ge 0$ (51)

which always holds. Combining (49) and (51), we get $0\le r\le \frac{1}{2}$ . Hence, the CFL condition $r=\frac{{\alpha }^{2}k}{h^2}$ plays an important role to the stability analysis [27] .

Figure 7 and Figure 8 show that the explicit scheme is stable for $r=0.4$ and $r=0.5$ respectively. But, the explicit scheme is unstable when the value of r exceeds 0.5. The instability of the explicit scheme is shown in Figure 9 for $r=0.6$ . The spatial domain of the substance is considered as $x=0$ to $x=\pi$ in Figures 7(a)-9(a) and the spatial domain $x=0$ to $x=2\pi$ is used in Figures 7(b)-9(b), respectively.

We observe that our explicit numerical scheme satisfy the stability property. Hence, we can conclude that our scheme is numerically efficient.

6.2. Stability Criterion for Crank-Nicolson Scheme

By the Crank-Nicolson scheme, the discretized form of Equation (1)

$u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}-\frac{r}{2}\left[u_{i-\mathrm{1,}j+1}-2u_{i\mathrm{,}j+1}+u_{i+\mathrm{1,}j+1}\right]-\frac{r}{2}\left[u_{i-\mathrm{1,}j}-2u_{i\mathrm{,}j}+u_{i+\mathrm{1,}j}\right]=0.$

Using similar arguments as in the FTCS method, the error ${\xi }_{k\mathrm{,}m}$ also follows the discretized ODE equation of the Crank-Nicolson scheme

${\xi }_{i\mathrm{,}j+1}-{\xi }_{i\mathrm{,}j}-\frac{r}{2}\left[{\xi }_{i-\mathrm{1,}j+1}-2{\xi }_{i\mathrm{,}j+1}+{\xi }_{i+\mathrm{1,}j+1}\right]-\frac{r}{2}\left[{\xi }_{i-\mathrm{1,}j}-2{\xi }_{i\mathrm{,}j}+{\xi }_{i+\mathrm{1,}j}\right]=0$ (52)

Divide (52) by ${\xi }_{i\mathrm{,}j}$ ,

$\frac{{\xi }_{i\mathrm{,}j+1}}{{\xi }_{i\mathrm{,}j}}-\frac{{\xi }_{i\mathrm{,}j}}{{\xi }_{i\mathrm{,}j}}-\frac{r}{2}\left[\frac{{\xi }_{i-\mathrm{1,}j+1}}{{\xi }_{i\mathrm{,}j}}-\frac{2{\xi }_{i\mathrm{,}j+1}}{{\xi }_{i\mathrm{,}j}}+\frac{{\xi }_{i+\mathrm{1,}j+1}}{{\xi }_{i\mathrm{,}j}}\right]-\frac{r}{2}\left[\frac{{\xi }_{i-\mathrm{1,}j}}{{\xi }_{i\mathrm{,}j}}-\frac{2{\xi }_{i\mathrm{,}j+}}{{\xi }_{i\mathrm{,}j}}+\frac{{\xi }_{i+\mathrm{1,}j}}{{\xi }_{i\mathrm{,}j}}\right]=0$ (53)

Simplifying each term in Equation (53), we get

$\frac{{\xi }_{i\mathrm{,}j+1}}{{\xi }_{i\mathrm{,}j}}=G=\frac{{\text{e}}^{a\left(t+k\right)}}{{\text{e}}^{at}}$ (54)

$\frac{{\xi }_{i-\mathrm{1,}j+1}}{{\xi }_{i\mathrm{,}j}}=G{\text{e}}^{-i{k}_{m}h}$ (55)

$\frac{{\xi }_{i+\mathrm{1,}j+1}}{{\xi }_{i\mathrm{,}j}}=G{\text{e}}^{i{k}_{m}h}$ (56)

$\frac{{\xi }_{i-\mathrm{1,}j}}{{\xi }_{i\mathrm{,}j}}={\text{e}}^{-i{k}_{m}h}$ (57)

$\frac{{\xi }_{i+\mathrm{1,}j}}{{\xi }_{i\mathrm{,}j}}={\text{e}}^{i{k}_{m}h}$ (58)

(a)(b)

Figure 7. Stability for the explicit numerical scheme for the spatial domain $x=0$ to $x=\pi$ in (a) and spatial domain $x=0$ to $x=2\pi$ in (b), respectively, using $r=0.4$ .

(a)(b)

Figure 8. Stability for the explicit numerical scheme for the spatial domain $x=0$ to $x=\pi$ in (a) and spatial domain $x=0$ to $x=2\pi$ in (b), respectively, using $r=0.5$ .

(a)(b)

Figure 9. Instability for the explicit numerical scheme for the spatial domain $x=0$ to $x=\pi$ in (a) and spatial domain $x=0$ to $x=2\pi$ in (b), respectively, using $r=0.6$ .

Plugging Equations (54)-(58) into Equation (53) yields

$G-1=\frac{r}{2}G\left({\text{e}}^{-i{k}_{m}h}-2+{\text{e}}^{i{k}_{m}h}\right)+\frac{r}{2}\left({\text{e}}^{-i{k}_{m}h}-2+{\text{e}}^{i{k}_{m}h}\right)$ (59)

Using a trigonometric identity involving complex exponentials, such as

$2\cos \theta ={\text{e}}^{i\theta }+{\text{e}}^{-i\theta }$ (60)

Using (60) into (59) yields

$G-1=\frac{r}{2}G\left(2\cos \left(k_{m}h\right)-2\right)+\frac{r}{2}\left(2\cos \left(k_{m}h\right)-2\right)\mathrm{.}$

Simplifying, we get

$G=\frac{1-r\left(1-\cos \left(k_{m}h\right)\right)}{1-r\left(1-\cos \left(k_{m}h\right)\right)}$ (61)

Remembering that the cosine range is [−1, 1], the value we should explore for cosine values are $\cos \left(k_{m}h\right)=-\mathrm{1,0,1}$ , hence we arrive at three facts.

Fact-01: If $\cos \left(k_{m}h\right)=-1$ , then from (61), we get

$G=\frac{1-2r}{1+2r}\mathrm{.}$ (62)

Since $r>0$ , then $G\le 1$ , always.

Fact-02: If $\cos \left(k_{m}h\right)=0$ , then from (61), we get

$G=\frac{1-r}{1+r}\mathrm{.}$ (63)

Since $r>0$ , then $G\le 1$ , always.

Fact-03: If $\cos \left(k_{m}h\right)=1$ , then from (61), we get

$G=1.$ (64)

Then $G\le 1$ , always. Therefore, we can simply say, the Crank-Nicolson scheme is unconditionally stable.

Figures 10-12 show that the Crank-Nicolson scheme is always stable. The mesh ratio $r=0.5$ , $r=1.0$ and $r=1.5$ are used in Figures 10-12, respectively. The spatial domain is considered as $x=0$ to $x=\pi$ in Figures 10(a)-12(a) and the spatial domain $x=0$ to $x=2\pi$ is used in Figures 10(b)-12(b), respectively.

The above figures illustrate that the Crank-Nicolson numerical scheme satisfies the stability conditions and thus we can conclude that the scheme is numerically efficient.

7. Error Analysis and Results Discussion

In this section, we focus on an error analysis between the exact solution and the approximate solution. To derive the relative error between the exact solution and the approximate solution for FTCS scheme, we compute the relative error in L₁-norm defined by

$\text{error}=\frac{\Vert u_e-u_a\Vert }{\Vert u_e\Vert }\mathrm{,}$

(a)(b)

Figure 10. Stability for Crank-Nicolson scheme for the spatial domain $x=0$ to $x=\pi$ in (a) and spatial domain $x=0$ to $x=2\pi$ in (b), respectively, using $r=0.5$ .

(a)(b)

Figure 11. Stability for Crank-Nicolson scheme for the spatial domain $x=0$ to $x=\pi$ in (a) and spatial domain $x=0$ to $x=2\pi$ in (b), respectively, using $r=1.0$ .

(a)(b)

Figure 12. Stability for Crank-Nicolson scheme for the spatial domain $x=0$ to $x=\pi$ in (a) and spatial domain $x=0$ to $x=2\pi$ in (b), respectively, using $r=1.5$ .

where, $u_e$ and $u_a$ are the exact solution and the approximate solution computed by finite difference scheme, respectively.

The exact solution of heat Equation (1) with boundary conditions $u\left(\mathrm{0,}t\right)=u\left(L\mathrm{,}t\right)=0$ , and initial condition $u_0\left(x\right)=u\left(x,0\right)=\sin \left(\frac{\pi x}{L}\right)$ is [1] ,

(a)(b)

Figure 13. Comparison between the exact solution and the explicit numerical solution (a) and L₁-norm of the error (b). (a) Comparison of the exact and the approximate solutions; (b) Error.

$u\left(x\mathrm{,}t\right)=\sin \left(\frac{\pi x}{L}\right){\text{e}}^{-\frac{{\alpha }^2{\pi }^{2}t}{L^2}}\mathrm{.}$ (65)

In Figure 13, we compare the exact and the approximate solution of the one-dimensional heat equation. In this case, we use FTCS explicit scheme to approximate the numerical solution. The L₁-norm of the error is plotted in Figure 13(b). From the figure, it is clear that our numerical scheme converges and efficiently accurate.

8. Conclusion

In this work, two finite difference schemes for the one-dimensional heat equation are presented. In the first step, the explicit and Crank-Nicolson schemes have been developed. Later, the stability criterion is set, and the stability conditions are verified for both the schemes, which shows the validity of the numerical schemes. By using standard computer programming codes, we examine the stability criteria. We find that the explicit scheme is conditionally stable with

$0\le r\le \frac{1}{2}$ , andthe Crank-Nicolson scheme is unconditionally stable. Finally, we

discuss the error between the exact and the approximate solution of the one-dimensional heat equation, and validate the efficiency of our numerical schemes.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Recktenwald, G.W. (2004) Finite Difference Approximations to the Heat Equation. Portland State University, Portland.
[2]	Raishinghania, M.D. (2000) Advanced Differential Equations. S. Chand & Company, New Delhi.
[3]	Loskor, W.Z. and Sarkar, R. (2022) A Numerical Solution of Heat Equation for Several Thermal Diffusivity Using Finite Difference Scheme with stability Conditions. Journal of Applied Mathematics and Physics, 10, 449-465. https://doi.org/10.4236/jamp.2022.102034
[4]	Kabir, M.H., Afroz, A., Hossain, M.M. and Gani, M.O. (2010) A Finite Difference Method for One Dimensional Heat Equation. Journal of Scientific and Engineering Research, 1, 54-59.
[5]	Bernatz, R. (2010) Fourier Series and Numerical Methods for Partial Differential Equations. John Wiley & Sons, Inc., Hoboken. https://doi.org/10.1002/9780470651384
[6]	Makhtoumi, M. (2010) Numerical Solutions of Heat Diffusion Equation over One Dimensional Rod Region. International Journal of Science and Advanced Technology, 7, 10-13.
[7]	Gorguis, A. and Chan, W.K.B. (2008) Heat Equation and Its Comparative Solutions. Computers & Mathematics with Applications, 55, 2973-2980. https://doi.org/10.1016/j.camwa.2007.11.028
[8]	Adomian, G. (1998) A Review of the Decomposition Method in Applied Mathematics. Journal of Mathematical Analysis and Applications, 135, 501-544. https://doi.org/10.1016/0022-247X(88)90170-9
[9]	Adomian, G. (1994) On the Solution of Partial Differential Equations with Specified Boundary Conditions. Journal of Mathematical Analysis and Applications, 140, 569-581. https://doi.org/10.1016/0022-247X(89)90084-X
[10]	Maturi, D.A., Alsulami, N.M. and Alaidarous, E.S. (2020) Finite Difference Approximation for Solving Transient Heat Conduction Equation of Copper. Advances in Pure Mathematics, 10, 350-358. https://doi.org/10.4236/apm.2020.105021
[11]	Mebrate, B. (2015) Numerical Solution of a One Dimensional Heat Equation with Dirichlet Boundary Conditions. American Journal of Applied Mathematics, 3, 305-311. https://doi.org/10.11648/j.ajam.20150306.20
[12]	Dhawan, S. and Kumar, S. (2009) A Numerical Solution of One Dimensional Heat Equation Using Cubic B-Spline Basis Functions. International Journal of Research and Reviews in Applied Science, 1, 71-77.
[13]	Mu’lla, M.A.M., Gaweash, A.M.A. and Bakur, H.Y.I. (2022) Numerical Solution of Parabolic in Partial Differential Equations (PDEs) in One and Two Space Variable. Journal of Applied Mathematics and Physics, 10, 311-321. https://doi.org/10.4236/jamp.2022.102024
[14]	Zhou, C., Zou, Y., Chai, S., Zhang, Q. and Zhu, H. (2018) Weak Galerkin Mixed Finite Element Method for Heat Equation. Applied Numerical Mathematics, 123, 180-199. https://doi.org/10.1016/j.apnum.2017.08.009
[15]	Naveros, I., Ghiaus, C., Ordoñez, J. and Ruíz, D.P. (2016) Thermal Networks from the Heat Equation by Using the Finite Element Method. WIT Transactions on Engineering Sciences, 106, 1743-3533. https://doi.org/10.2495/HT160041
[16]	Saptaningtyas, F.Y. and Setyarsi, A.D. (2018) Finite Volume Method with Explicit Scheme Technique for Solving Heat Equation. Journal of Physics: Conference Series, 1097, Article ID: 012089. https://doi.org/10.1088/1742-6596/1097/1/012089
[17]	Asish, M. (2016) Finite-Volume Solution of Diffusion Equation and Application to Model Problems. International Journal of Research in Engineering and Technology, 2, 46-55.
[18]	Causon, D.M. and Mingham, C.G. (2010) Introductory Finite Difference Methods for PDEs.
[19]	Mahardika, D.P. and Haryani, F.F. (2020) Numerical Analysis of One Dimensional Heat Transfer on Varying Metal. Journal of Physics: Conference Series, 1511, Article ID: 012049. https://doi.org/10.1088/1742-6596/1511/1/012049
[20]	Hasnain, S., Saqib, M., and Mashat, D.S. (2017) Numerical Study of One Dimensional Fishers KPP Equation with Finite Difference Schemes. American Journal of Computational Mathematics, 7, 70-83. https://doi.org/10.4236/ajcm.2017.71006
[21]	Burden, R.L. and Faires, J.D. (2012) Numerical Analysis. Ninth Edition, Brooks/Cole, Cengage Learning, Boston, MA, USA.
[22]	Herman, R.L. (2015) Introduction to Partial Differential Equations.
[23]	Crank, J. and Nicolson, P. (1996) A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat-Conduction Type. Advances in Computational Mathematics, 6, 207-226. https://doi.org/10.1007/BF02127704
[24]	Fadugba, S.E., Edogbanya, O.H. and Zelibe, S.C. (2013) Crank Nicolson Method for Solving Parabolic Partial Differential Equations. International Journal of Applied Mathematics and Modeling, 1, 8-23.
[25]	Chung, T.J. (2002) Computational Fluid Dynamics. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511606205
[26]	Gerald, C.F. and Wheatley, P.O. (1915) Applied Numerical Analysis. Seventh Edition, Pearson, London.
[27]	Alam, M.J., Loskor, W.Z. and Mohiuddin, M. (2021) Time Dependent Wave Propagation Modeling Using Finite Difference Scheme of 2D Wave Equation Based on Absorbing and Reflecting Boundaries. Journal of Applied Mathematics and Physics, 9, 2334-2344. https://doi.org/10.4236/jamp.2021.99148