Solution of a One-Dimension Heat Equation Using Higher-Order Finite Difference Methods and Their Stability

M. Emran Ali¹, Wahida Zaman Loskor², Samia Taher², Farjana Bilkis²

¹Department of Textile Engineering, Northern University Bangladesh, Dhaka, Bangladesh.
²Department of Science and Humanities, Bangladesh Army International University of Science & Technology, Cumilla, Bangladesh.

DOI: 10.4236/jamp.2022.103060                

Abstract

One-dimensional heat equation was solved for different higher-order finite difference schemes, namely, forward time and fourth-order centered space explicit method, backward time and fourth-order centered space implicit method, and fourth-order implicit Crank-Nicolson finite difference method. Higher-order schemes have complexity in computing values at the neighboring points to the boundaries. It is required there a specification of the values of field variables at some points exterior to the domain. The complexity was incorporated using Hicks approximation. The convergence and stability analysis was also computed for those higher-order finite difference explicit and implicit methods in case of solving a one dimensional heat equation. The obtained numerical results were compared with exact solutions. It is found that backward time and fourth-order centered space implicit scheme along with Hicks approximation performed well over the other mentioned higher-order approaches.

Keywords

Heat Equation, Boundary Condition, Higher-Order Finite Difference Methods, Hicks Approximation

1. Introduction

Partial differential equations (PDEs) are found to arise in illustrating numerous physical problems, such as, heat transfer, fluid flow, solid mechanics, economics, biological process, etc. [1]. The PDEs can often be classified into three categories, namely, parabolic type equation, hyperbolic type equation and elliptic type equation. Generally, they are related, respectively, to diffusion, advection and steady-states of either hyperbolic or parabolic problems. Heat equation a parabolic type PDE enable to explain numerous physical problems raised in scientific, engineering, economic sectors that is also known as diffusion equation [2]. The equation generally represents temperature distribution in a fixed region over elapsing of time [3]. It is often solved using different analytical and numerical tools [1] - [6]. For example, Gorguis and Chan [3] solved (1 + 1) dimensional heat equation analytically using Adomian decomposition method and compared obtained solution with that came through separation of variable method; Paul and Ali [4] and Ahmad and Yaacob [6] used method of lines in coordination with higher-order finite difference approximation to solve the equation; Muhiddin and Sulaiman [2] solved the equation using fourth-order Crank-Nicolson (CN4) finite difference method (FDM) and fourth-order standard implicit FDM (BTCS) and made a comparison of obtained results; fourth-order iterative alternating decomposition explicit method of Mitchell and Fairweather was exercised in the study due to Mansor et al. [4]. As we cannot find the analytical solution of most of the PDEs, efficient and faster numerical techniques are highly appreciable in research community for solving those. It is of interest to note here that higher-order finite difference approximation methods increase the accuracy of model results and it also processes local truncation errors (LTEs) of higher order [7] [8]. But use of higher-order FDM has complexity in generating values at the neighboring points to the boundaries of the domain (red points in Figure 1) wherein it is required to know the values at the points outside the boundaries (green or violet points in Figure 1). Muhiddin and Sulaiman [2] solve the complexity by considering a regular second-order FDM scheme for the points immediate to the boundaries. But we believe that this kind of consideration may weaken the beauty the use of higher-order FDM. Paul and Ali [7] incorporated the complexity by considering weighted average technique in storm surge simulation. Hicks and Wei [9] suggested an approximation technique (hereafter Hicks approximation) in incorporating those points outside the boundary. In this study, we have used Hicks approximation along with different explicit and implicit FDMs. The LTEs and stability for mentioned higher-order FDMs have also been studied here. The comparisons among the computed numerical results are also made with the exact solution and found them to perform well. It is important to mention here that as higher-order FDMs, we use forward time and fourth-order centered space (FTCS4), backward time and fourth-order centered space (BTCS4), and fourth-order Crank Nicolson (CN4) FDMs were used.

Figure 1. The idealized gridding whereas water blue represents boundary points.

2. Materials and Methods

2.1. Problem Statements

As a model equation, one-dimensional heat equation available in Ahmad and Yaacob [6], with auxiliary conditions, can be written in form

$v_t=c{v}_{xx}$ ; $0<x<1;t>0$, (1)

subject to initial and boundary conditions

$v\left(0<x<1,0\right)=70˚\text{C}$, (2)

$v\left(0,t\right)=50˚\text{C}$, (3)

$v\left(1,t\right)=20˚\text{C}$, (4)

where the state variable $v\left(x,t\right)$ stands for representing temperature distribution of rod of unit length at point x over time t, which is initially heated with 70˚C providing 50˚C and 20˚C boundary temperatures. Here c represents the diffusivity of the substance of the rod.

It is of interest noting here that Hicks approximation is applied only for homogeneous boundary (see [9] ). So, a convenient transformation $u\left(x,t\right)=v\left(x,t\right)+30x-50$ has been introduced in Equation (1) and it yielded

$u_t=c{u}_{xx}$ ; $0<x<1;t>0$, (5)

subject to transformed initial and boundary conditions

$u\left(0<x<1,0\right)=30x+20˚\text{C}$, (6)

$u\left(0,t\right)=0˚\text{C}$ (7)

$u\left(1,t\right)=0˚\text{C}$. (8)

Thus the boundary conditions reform into homogenous type boundary conditions.

2.2. Higher FDMs and Discretization of Model Equation

In this manuscript, we have discretized the model equation using three higher-order FDM techniques, namely, FTCS4, BTCS4, and CN4. Whereas all FDM schemes are of implicit type except the former one. Suppose $u_i^n=u\left(x_i,t_n\right)$ where $x_i=\left(i-1\right)\Delta x$ for $i=1,2,3,\cdots ,N$ are the discretized grid points in space and $t_n=n\Delta t,n=1,2,3,\cdots ,M$ discretized time. It is worth mentioning here that ∆x and ∆t stand for denoting step size in space and time, respectively. Thus the discretization of Equation (5) using FTCS4, BTCS4, and CN4 can respectively put to the following forms,

$\frac{u_i^{n+1}-u_i^n}{\Delta t}=c\frac{-u_{i-2}^n+16u_{i-1}^n-30u_i^n+16u_{i+1}^n-u_{i+2}^n}{12\Delta x^2}$ (9)

$\frac{u_i^{n+1}-u_i^n}{\Delta t}=c\frac{-u_{i-2}^{n+1}+16u_{i-1}^{n+1}-30u_i^{n+1}+16u_{i+1}^{n+1}-u_{i+2}^{n+1}}{12\Delta x^2}$ (10)

$\begin{array}{c}\frac{u_i^{n+1}-u_i^n}{\Delta t}=c[\frac{-u_{i-2}^n+16u_{i-1}^n-30u_i^n+16u_{i+1}^n-u_{i+2}^n}{24\Delta x^2}\\ +\frac{-u_{i-2}^{n+1}+16u_{i-1}^{n+1}-30u_i^{n+1}+16u_{i+1}^{n+1}-u_{i+2}^{n+1}}{24\Delta x^2}]\end{array}$ (11)

Thus the initial and boundary conditions symbolized by Equations (6)-(8) can be reformed, respectively, as

$u_i^1=30x_i+20˚\text{C}$, $u_1^n=0˚\text{C}$, and $u_N^n=0˚\text{C}$. (12)

3. Local Truncation Error for Heat Equation

The local truncation error (LTE) at (n, i)-th grid point represents the deviation of approximated scheme result from exact result at the grid point [10]. Smith et al. [11] defines LTE as follows:

Suppose $G_i^n\left(f\right)=0$ represents a difference equation approximating a PDE at (n, i)-th grid point. If we replace f by F, where F is the exact solution of the PDE, the value $G_i^n\left(F\right)$ represents the LTE,

$LT{E}_i^n=G_i^n\left(F\right)$, (13)

where $LT{E}_i^n$ represents LTE at (n, i)-th grid point. The right hand side of the expression can easily be expressed as a power of ∆t and ∆x using Taylor series expansion. The LTE has been made for FTCS4, BTCS4, and CR4 as follows:

For the FTCS4 explicit scheme (Equation (9)), the local truncation error can be calculated as

$LT{E}^{FTCS4}=\frac{u_i^{n+1}-u_i^n}{\Delta t}-c\frac{-u_{i-2}^n+16u_{i-1}^n-30u_i^n+16u_{i+1}^n-u_{i+2}^n}{12\Delta x^2}$. (14)

Considering u is an exact solution and using Taylor’s series expansion on the right side of Equation (14), it can be written as [12]

$LT{E}^{FTCS4}=\frac{\Delta t}{2}\frac{{\partial }^{2}u}{\partial t^2}+O\left(\Delta t^2\right)-\frac{2c\Delta x^4}{60}\frac{{\partial }^{6}u}{\partial x^6}+O\left(\Delta x^5\right)$,

or,

$LT{E}^{FTCS4}=O\left(\Delta t\right)+O\left(\Delta x^4\right)$. (15)

Similarly for BTCS4, we can get

$LT{E}^{BTCS4}=\frac{\Delta t}{2}\frac{{\partial }^{2}u}{\partial t^2}+O\left(\Delta t^2\right)-\frac{2c\Delta x^4}{60}\frac{{\partial }^{6}u}{\partial x^6}+O\left(\Delta x^5\right)$,

or,

$LT{E}^{BTCS4}=O\left(\Delta t\right)+O\left(\Delta x^4\right)$ (16)

Thus for CN4 finite difference scheme

$LT{E}^{CN4}=O\left(\Delta t^2\right)+O\left(\Delta x^4\right)$. (17)

The LTE may help to understand how farter the error will converse with proper choice of step size.

4. Stability Analysis for Heat Equation

In contrast, stability obtained from a numerical scheme is the necessity in a numerical computation. It is mentioned in the study due to Tavella and Randall [13] that a numerical scheme is said to be stable if the difference between the numerical solution and the exact solution remains bounded as the number to time steps tend to infinity.

Suppose error values emanated for a numerical scheme ${\epsilon }_i^n$ can be defined as ${\epsilon }_i^n=u_i^n-E_i^n$, where $u_i^n$ and $E_i^n$ are the numerical and exact solution, respectively, at any (i, n-th) grid point [14].

From the scheme in Equation (9), we can write

${\epsilon }_i^{n+1}={\epsilon }_i^n+r\left(-{\epsilon }_{i-2}^n+16{\epsilon }_{i-1}^n-30{\epsilon }_i^n+16{\epsilon }_{i+1}^n-{\epsilon }_{i+2}^n\right)$, (18)

where $r=\frac{c\Delta t}{12\Delta x^2}$.

Expanding ${\epsilon }_i^n$ using finite Fourier series, one can define it as (see [14] )

${\epsilon }_i^n={\text{e}}^{bn\Delta t}{\text{e}}^{Iki\Delta x}$, $I=\sqrt{-1}$ (19)

where ${\text{e}}^{bn\Delta t}$ and k are the wave amplitude and the wave number, respectively, of error propagation curve, and b is a constant.

Imposing Equation (19) in Equation (18), it can put to the following form

${\text{e}}^{b\Delta t}=1+r\left(-2\cos \left(2k\Delta x\right)+32\cos \left(k\Delta x\right)-30\right)$. (20)

A scheme is said to be stable if the error growth factor $G\left(k\right)=\left|\frac{{\epsilon }_i^{n+1}}{{\epsilon }_i^n}\right|=\left|{\text{e}}^{b\Delta t}\right|\le 1$ [14].

In worst case, suppose that $k\Delta x=\pi$, or $G\left(k\right)=1-64r$, implies $r\le \frac{1}{32}$. Thus the Von Neumann stability criteria for FTCS4 explicit finite difference scheme can be defined as $\frac{c\Delta t}{\Delta x^2}\le \frac{3}{8}$.

Similarly, BTCS4 implicit scheme (Equation (10)), the error growth factor

$G\left(k\right)=\left|\frac{1}{1+r\left(30+2\cos \left(2k\Delta x\right)-32\cos \left(k\Delta x\right)\right)}\right|$ (21)

Considering the maximum value of $\cos \left(k\Delta x\right)=1$, it is found that $G\left(k\right)\le 1$. Hence the system is unconditionally stable.

For the CN4 finite difference scheme, we can find the growth factor

$G\left(k\right)=\left|\frac{1+s\left(-2\cos \left(2k\Delta x\right)+32\cos \left(k\Delta x\right)-30\right)}{1+s\left(30+2\cos \left(2k\Delta x\right)-32\cos \left(k\Delta x\right)\right)}\right|$, (22)

where $s=\frac{c\Delta t}{24\Delta x^2}$. The right side of the Equation (22) is less than 1 for every

value of $\cos \left(k\Delta x\right)$. Thus the implicit BTCS4 and CN4 schemes are unconditionally stable.

5. Integration Procedure

It is to be noted here that Hicks approximation is applicable only for homogeneous boundary conditions [9]. So a transformation was imposed in Equation (1) and it turned into Equation (5) subject to initial conditions (Equation (6)) and homogeneous boundary conditions (Equations (7) and (8)). This equation was solved firstly solved using difference equation (Equation (9)) for getting numerical solution for FTCS4 scheme. Use of FTCS4 scheme has complexities in incorporating results at red grid points of Figure 1. It is required to know the values at the green points of Figure 1 to compute results at red points (Figure 1). This complexity was incorporated in different ways such as Muhiddin and Sulaiman [1] incorporated that using lower-order scheme, Paul and Ali [4] suggested maintaining boundary temperature at all grid points outside the boundary, Paul and Ali [7] used weighted average technique, and following Hicks and Wei [9], we have incorporated it as $u_{N-1}=-u_{N+1}$, where N is boundary points. In this study, we have incorporated Hicks approximation technique.

Similarly, for BTCS4 FDM (Equation (10)), it is required to approximate the values at violet points (Figure 1) for computing results at red points (Figure 1). Finally, when the values at red grid points (Figure 1) are computed using CN4 FDM in Equation (11), both green and violet grid points are required to approximate using Hicks approximation technique. It is important to note here that for all the mentioned methods, we have developed Matlab routines whereas same spatial step ( $\Delta x=0.10$ ) and tome step ( $\Delta t=0.0051$ ) size was maintained.

6. Result and Discussion

It is to be noted here that we have solved Equation (5) subject to initial and boundary conditions characterized by Equations (6)-(8) using FTCS4, BTCS4 and CN4 finite difference methods shown in Equations (9)-(11), respectively. At the first a Matlab routine was developed for the difference equation symbolized by Equation (9) with the help of boundary conditions represented by Equation (12) whereas the Hicks approximation $u_{N-1}=-u_{N+1}$ was used for generating the values of the field variables at the neighboring points and it yielded the solution of the solution of Equation (5). Taking the inverse transformation $v\left(x,t\right)=u\left(x,t\right)-30x+50$, we get the required solution of Equation (1) using FTCS4 scheme. It is worth to note here that similar techniques was implemented for other mentioned FDM techniques. The graphical representation of temperature profiles obtained from the mentioned numerical FDMs are presented, respectively, in Figures 2-4. The obtained numerical results are also compared with the analytic result (available in [6] ) in Figure 5 for $t=0.3$ second. One can perceive from the figure that the BTCS4 shows better results over the other numerical techniques (FTCS4 and CN4). It may have two reasons. In the first, it is an implicit scheme so it will be performed better over FTCS4 explicit scheme. Secondly, during numerical computation, it is required to use Hicks approximation in a single time (only in violet points) whereas for CN4 scheme, Hicks

Figure 2. Obtained temperature profile $u\left(x,t\right)$ (in ˚C) distribution from FTCS4 method at different points of domain with elapsing time t second.

Figure 3. Obtained temperature profile $u\left(x,t\right)$ (in ˚C) distribution from BTCS4 method at different points of domain with elapsing time t second.

approximation is needed for double points (both green and violet in Figure 1). Though, order of LTE for CN4 scheme is higher than the other two scheme, BTCS4 shows better performance over the CN4 scheme. It is worth noting here that we have tested the mentioned stability conditions during execution of numerical algorithms. It was found that the BTCS4 and CN4 were free from any restrictions on spatial or time step size and BTCS4 explicit scheme was separately run several values of r obtained for changing values of step size (both spatial and time) and the thermal diffusivity c.

For better understanding, we have also computed the absolute relative errors at different spatial grid points with $t=0.3$ second for mentions FDMs from the analytic results and presented in Table 1. It can be seen from the table that the

Figure 4. Obtained temperature profile $u\left(x,t\right)$ (in ˚C) distribution from CN4 method at different points of domain with elapsing time t second.

Figure 5. Comparison of the obtained results from different results with observed data at t = 0.3 second.

Table 1. Comparison of relative errors from analytic solution for FTCS4, BTCS4, and CN4 FDMs.

absolute relative errors at different points of the domain computed for the FTCS4 scheme are lower towards the left boundary and gradually raising up along the right boundary whereas CN4 shows those higher along the central points of the domain of interest. In contrast, the relative errors emanated from BTCS4 scheme show lower relative than other two schemes. Thus on the basis of the comparison among the errors made in Table 1, we can say the BTCS4 finite difference along with Hicks approximation could be an alternative approach in solving parabolic equations like heat equation. Thus, it could be an efficient alternative approach in such numerical computation.

7. Conclusion

In this study, we have solved a one-dimensional heat equation using three different FDMs and compared their computed results. Their stability and LTE are also studied and found that having lower LTE of BTCS4 than CN4, BTCS4 found to show better performance over CN4 and FTCS4 in case of Hicks approximation. Even at some grid points, FTCS4 shows a relatively lower relative error than CN4. Whereas FTCS4 explicit scheme is conditionally stable. This study may help the researchers in identifying more accurate and stable schemes and hence numerical solution of the parabolic equation and thus, the approach can be a better alternative in such computation.

Conflicts of Interest

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

References

[1]	Abdulkadir, Y. (2015) Comparison of Finite Difference Schemes for the Wave Equation Based on Dispersion. Journal of Applied Mathematics and Physics, 3, 1544-1562. https://doi.org/10.4236/jamp.2015.311179
[2]	Muhiddin, F.A. and Sulaiman, J. (2018) Performance Analysis of Fourth-Order Crank-Nicolson Scheme for Solving Diffusion Equations. AIP Publishing. https://doi.org/10.1063/1.5041630
[3]	Gorguis, A. and Chan, W.K.B.C. (2008) Heat Equation and Its Comparative Solutions. Computers & Mathematics with Applications, 55, 2973-2980. https://doi.org/10.1016/j.camwa.2007.11.028
[4]	Mansor, A., Zulkifle, A.K., Alias, N., Hasan, M.K. and Boyce, M.J.N. (2013) The Higher Accuracy Fourth-Order IADE Algorithm. Journal of Applied Mathematics, 2013, Article ID: 236548. https://doi.org/10.1155/2013/236548
[5]	Paul, G.C. and Ali, M.E. (2019) Solution of One-Dimensional Heat Equation: An Alternative Approach. Menemui Matematik (Discovering Mathematics), 41, 22-33.
[6]	Ahmad, R.R. and Yaacob, N. (2013) Arithmetic-Mean Runge-Kutta Method and Method of Lines for Solving Mildly Stiff Differential Equations. Menemui Matematik (Discovering Mathematics), 35, 21-29.
[7]	Paul, G.C. and Ali, M.E. (2019) Numerical Storm Surge Model with Higher Order Finite Difference Method of Lines for the Coast of Bangladesh. Acta Oceanologica Sinica, 38, 100-116. https://doi.org/10.1007/s13131-019-1385-7
[8]	Chapra, S.C. and Canale, R.P. (2015) Numerical Methods for Engineers. Seventh Edition, McGrawHill Education, New York, USA.
[9]	Hicks, J.S. and Wei, J. (1967) Numerical Solution of Parabolic Partial Differential Equations with Two Point Boundary Conditions by Use of the Method of Lines. Journal of the ACM, 14, 549-562. https://doi.org/10.1145/321406.321417
[10]	Bouwer, A. (2009) The Du Fort and Frankel Finite Difference Scheme Applied to and Adapted for a Class of Finance problems. Ph.D. Dissertation, University of Pretoria, Pretoria.
[11]	Smith, G.D., Smith, G.D. and Smith, G.D.S. (1985) Numerical Solution of Partial Differential Equations: Finite Difference Methods. Oxford University Press, New York, USA.
[12]	Li, N., Steiner, J. and Tang, S. (1994) Convergence and Stability Analysis of an Explicit Finite Difference Method for 2-Dimensional Reaction-Diffusion Equations. The ANZIAM Journal, 36, 234-241. https://doi.org/10.1017/S0334270000010377
[13]	Tavella, D. and Randall, C. (2000) Pricing Financial Instruments: The Finite Difference Method. John Wiley & Sons, Inc., Cananda.
[14]	Abdullah, A. (2006) Simplified von Neumann Stability Analysis of Wave Equation Numerical Solution. ARPN Journal of Engineering and Applied Sciences, 14, 4164-4168.