The Higher Accuracy Fourth-Order IADE Algorithm

Abstract

This study develops the novel fourth-order iterative alternating decomposition explicit (IADE) method of Mitchell and Fairweather (IADEMF4) algorithm for the solution of the one-dimensional linear heat equation with Dirichlet boundary conditions. The higher-order finite difference scheme is developed by representing the spatial derivative in the heat equation with the fourth-order finite difference Crank-Nicolson approximation. This leads to the formation of pentadiagonal matrices in the systems of linear equations. The algorithm also employs the higher accuracy of the Mitchell and Fairweather variant. Despite the scheme’s higher computational complexity, experimental results show that it is not only capable of enhancing the accuracy of the original corresponding method of second-order (IADEMF2), but its solutions are also in very much agreement with the exact solutions. Besides, it is unconditionally stable and has proven to be convergent. The IADEMF4 is also found to be more accurate, more efficient, and has better rate of convergence than the benchmarked fourth-order classical iterative methods, namely, the Jacobi (JAC4), the Gauss-Seidel (GS4), and the successive over-relaxation (SOR4) methods.

1. Introduction

Numerical methods with accuracy of the order, 𝑂(Δℎ)^𝑛, 𝑛>2, are referred to as higher-order methods (ℎ = mesh size). Recent developments seem to desire methods of higher-order for achieving higher accuracy numerical solutions for problems involving partial differential equations. When higher-order methods are evaluated, factors such as rate of convergence, stability, and boundary conditions have to be considered too. There are mainly two categories of finite difference higher-order schemes.

The first category is the noncompact stencils that utilize any grid point surrounding the grid points of interest where the difference schemes are implemented. The method involves larger matrix bandwidth, but in many cases, the increase is not large [1]. There will be a probable increase in execution time as more grid points are used. However, increasing the number of grid points provides advantages in terms of enhancement in accuracy and improvement in resolution [2]. To enhance accuracy, the approach should consider proper treatment of boundary conditions of comparable accuracy.

The second category is the compact stencils that use Lax-Wendroff’s idea [3] proposed by MacKinnon and Carey [4]. The method uses smaller number of stencils, making it computationally efficient and highly accurate. However, it requires inversion of matrix to obtain spatial derivative at each point. Also, the boundary stencil has a large effect on the stability and accuracy of the scheme [5]. Usually, compact schemes of order higher than four require formation of auxiliary equations due to boundary conditions. This results in large bandwidth matrices which are not symmetric [6]. Furthermore, the schemes can get fairly complicated for more complex equations.

Some of the current findings on higher-order methods include the work by Sulaiman et al. [7] who suggested a fourth-order quarter sweep modified successive over-relaxation (QSMSOR) iterative method for solving a one-dimensional parabolic equation. It is found to be superior in terms of rate of convergence and execution time as compared to other SOR methods. Jha [8] formulated the six-order accurate quarter sweep alternating group explicit (QSAGE) iterative finite difference method for solving nonlinear singular two-point boundary value problems. The method can be implemented in parallel and is found to be superior to the corresponding full sweep alternating group explicit (AGE) and SOR methods. Jin et al. [9] proposed an AGE iteration method of the fourth-order accuracy by integrating the grouping explicit method with numerical boundary conditions. The method was used to solve initial-boundary value problem of convection equations. Fu and Tan [10] showed that an unconditionally stable split-step FDTD method with higher-order spatial accuracy is more accurate than the lower-order methods. The dispersion error of the proposed method is comparable with the higher-order ADI-FDTD method.

Higher-order methods are also studied by Mohebbi and Deghan [11] who applied a compact finite difference approximation of fourth-order and the cubic C¹ spline collocation method to some one-dimensional heat and advection-diffusion equations. The scheme has fourth-order accuracy in both space and time and unconditionally stable. Liao [12] proposed an efficient and high-order accuracy of the fourth-order compact finite difference method to solve one-dimensional Burgers’ equation. Tian and Ge [13] studied a stable fourth-order compact ADI method for solving two-dimensional unsteady convection diffusion problems. The method is temporally second-order and spatially fourth-order accurate, which requires only a regular five-point 2D stencil similar to that in the standard second-order methods. Chun [14] solved some nonlinear equations by applying the fourth-order iterative methods containing the King’s fourth-order family. It was observed that the proposed method has at least equal performance compared to other methods of the same order. Zhu et al. [15] presented a high-order parallel finite difference algorithm based on the domain decomposition strategy. The study used the classical explicit scheme to calculate the interface values between subdomains and the fourth-order compact schemes for the interior values. The method has high accuracy and is convergent and stable. Gao and Xie [16] devised a fourth-order alternating direction implicit compact finite difference schemes for the solution of two-dimensional Schrödinger equations. The method is highly competitive as compared to other existing methods, and it achieves the expected convergence rate.

This study develops a higher-order finite difference algorithm, with noncompact stencils, that is capable of delivering highly accurate solutions for a one-dimensional heat equation with Dirichlet boundary conditions. The study focuses on modifying the unconditionally stable and convergent second-order iterative alternating decomposition explicit (IADE) method of Mitchell and Fairweather (IADEMF2). The IADEMF2, which was originally proposed by Sahimi et al. [17], employs the fractional splitting of the Mitchell and Fairweather (MF) variant that has an accuracy of the order 𝑂((Δ𝑡)²+(Δ𝑥)⁴). The scheme executes a two-stage process involving the solution of sets of tridiagonal equations along lines parallel to the first and second time steps, respectively. It is found to be unconditionally stable, convergent, and more accurate than the classical AGE class of methods, namely, the AGE method based on the Peaceman and Rachford variant (AGE-PR), whose spatial accuracy is of second order, and the AGE method based on the Douglas and Rachford variant (AGE-DR), whose temporal accuracy is only of first order. The detailed derivation of the IADEMF2 can be obtained from [17].

In this paper, a fourth-order Crank-Nicolson (CN) difference approximation is applied to the spatial derivative in the heat equation, and the MF variant is employed, leading to the formation of the fourth-order IADEMF (IADEMF4) numerical algorithm. The convergence of the IADEMF4 is analyzed and proven. Numerical experiments verify the potential of the IADEMF4 in enhancing the accuracy of the IADEMF2. The results of the proposed higher-order scheme are also compared with the benchmarked fourth-order classical iterative methods, such as the fourth-order Gauss-Seidel (GS4), the fourth-order Jacobi (JAC4), and the fourth-order successive over-relaxation (SOR4) methods.

This paper is organized as follows. Section 2 discusses the implementation and stability of the fourth-order CN approximation on the heat equation. In Section 3, the IADEMF4 algorithm is formulated. Section 4 analyses the convergence of the IADEMF4. Section 5 provides the equations for the benchmarked fourth-order classical iterative methods. The computational complexity of the methods considered in this paper is given in Section 6, while the pseudocode for the IADEMF4 sequential algorithm is presented in Section 7. The experiments conducted are given in Section 8. Sections 9 and 10 provide some discussion and conclusion based on the obtained numerical results.

2. Fourth-Order Crank-Nicolson Approximation to the Spatial Derivative in the Heat Equation

Consider the following one-dimensional parabolic heat equation (1) that has been suitably assumed to be nondimensionalised. It models the flow of heat in a homogeneous unchanging medium of finite extent in the absence of heat source 𝜕𝑈=𝜕𝜕𝑡²𝑈𝜕𝑥² (1)

subject to given initial and Dirichlet boundary conditions 𝑈()(𝑥)𝑈()(𝑡)𝑈()(𝑡)𝑥,0=𝑓,0≤𝑥≤1,0,𝑡=𝑔,0<𝑡≤𝑇,1,𝑡=ℎ,0<𝑡≤𝑇. (2)

For the problem in (1), the finite difference approach discretizes the time-space domain by placing a rectangular grid over the domain, with grid spacing of Δ𝑡and Δ𝑥 in the 𝑡- and 𝑥-directions, respectively. The grid consists of the set of lines parallel to the𝑡-axis given by 𝑥_𝑖=𝑖Δ𝑥, 𝑖=0,1,…,𝑚,𝑚+1 and a set of lines parallel to the 𝑥-axis given by 𝑡_𝑘=𝑘Δ𝑡, 𝑘=0,1,…,𝑛,𝑛+1. For simplicity, the grid spacing is taken to be uniform, so that Δ𝑥=1/(𝑚+1), and Δ𝑡=𝑇/(𝑛+1). At a grid-point 𝑃(𝑥_𝑖,𝑡_𝑘) in the solution domain, the dependent variable 𝑈(𝑥,𝑡)which represents the nondimensional temperature at time 𝑡 and at position 𝑥 is approximated by 𝑢^𝑘_𝑖.

At the grid-point𝑃(𝑥_𝑖,𝑡_𝑘+1/2), the IADEMF4 scheme replaces the spatial derivative in the heat equation with a higher-order, particularly the fourth-order Crank-Nicolson (CN) difference approximation [18]. This is shown in the expression given in (3), with the central difference operators defined as 𝛿²_𝑥𝑢^𝑘_𝑖=𝑢^𝑘_𝑖−1−2𝑢^𝑘_𝑖+𝑢^𝑘_𝑖+1 and 𝛿²_𝑥𝑢_𝑖^𝑘+1=𝑢^𝑘+1_𝑖−1−2𝑢_𝑖^𝑘+1+𝑢^𝑘+1_𝑖+1. The approach gives the fourth-order scheme a spatial truncation error of the order 𝑂(Δ𝑥)⁴. It enhances the accuracy of its second-order counterpart, which has a larger error of the order 𝑂(Δ𝑥)²1𝑢Δ𝑡^𝑖𝑘+1−𝑢^𝑘_𝑖=12(Δ𝑥)²𝛿²_𝑥−1𝛿12^4𝑥𝑢_𝑖^𝑘+1+𝑢^𝑘_𝑖. (3)

To determine the stability of (3), the von Neumann stability analysis can be applied. Let 𝜆=Δ𝑡/(Δ𝑥)², and since 𝛿⁴_𝑥𝑢^𝑘_𝑖=𝛿²_𝑥(𝛿²_𝑥𝑢^𝑘_𝑖), the discretization of (3) becomes 𝑢_𝑖^𝑘+1−𝑢^𝑘_𝑖=𝜆2−1𝑢12^{𝑘+1𝑖−2}+43𝑢^𝑘+1_𝑖−1−52𝑢_𝑖^𝑘+1+43𝑢^𝑘+1_𝑖+1−1𝑢12^𝑘+1𝑖+2−1𝑢12^𝑘𝑖−2+43𝑢^𝑘_𝑖−1−52𝑢^𝑘_𝑖+43𝑢^𝑘_𝑖+1−1𝑢12^𝑘𝑖+2. (4)

The discretized equation in (4) is assumed to have a solution in the form of a Fourier harmonic function; that is, 𝑢^𝑘_𝑖=𝜌^𝑘𝑒^{𝑠𝛽𝑖Δ𝑥}, where 𝜌 is referred to as the amplification factor, 𝛽 is an arbitrary constant, and 𝑠=√−1. The amplification factor represents the time dependence of the solution. If the Fourier function is substituted into (4) and then solve for 𝜌, the result will be ()𝜌=1−𝜆/6cos𝛽−4²−91+(𝜆/6)cos𝛽−4²−9,with𝛽=𝛽Δ𝑥. (5)

Since |𝜌|≤1, then the fourth-order CN approximation is unconditionally stable for any choice of 𝛽, Δ𝑡, and Δ𝑥.

3. The Formulation of the IADEMF4

The IADEMF4 is firstly developed based on the execution of the unconditionally stable fourth-order CN approximation (3) on the heat equation.

Equation (4) can also be expressed as in (6), using the definitions of the constants given in (7) 𝑎𝑢^𝑘+1_𝑖−2+𝑏𝑢^𝑘+1_𝑖−1+𝑐𝑢_𝑖^𝑘+1+𝑑𝑢^𝑘+1_𝑖+1+𝑒𝑢^𝑘+1_𝑖+2=−𝑎𝑢^𝑘_𝑖−2−𝑏𝑢^𝑘_𝑖−1+𝑐𝑢̂^𝑘_𝑖−𝑑𝑢^𝑘_𝑖+1−𝑒𝑢^𝑘_𝑖+2,𝜆𝑖=2,3,…,𝑚−1,𝑎=24,𝑏=−2𝜆3,𝑐=4+5𝜆4,𝑑=−2𝜆3𝜆,𝑒=,̂24𝑐=4−5𝜆4. (6)(7) The approximation in (6) can be displayed in a matrix form such as 𝐴𝐮=𝐟 (8), where 𝐴 is a sparse pentadiagonal coefficient matrix and 𝐮=(𝑢₂,𝑢₃,…,𝑢_𝑚−2,𝑢_𝑚−1)^𝑇 is the column vector containing the unknown values of 𝑢 at the time level 𝑘+1. The column vector 𝐟=(𝑓₂,𝑓₃,…,𝑓_𝑚−2,𝑓_𝑚−1)^𝑇 consists of boundary values and known 𝑢 values at the previous time level 𝑘. The definitions for every entry in 𝐟 are given in (9) 𝐴𝐮=𝐟𝑐𝑑𝑒𝑏𝑐𝑑𝑒𝐎𝑎𝑏𝑐𝑑𝑒⋱⋱⋱⋱𝑎𝑏𝑐𝑑𝑒𝐎𝑎𝑏𝑐𝑑𝑎𝑏𝑐_{(𝑚−2)×(𝑚−2)}𝑢₂𝑢₃⋮𝑢_𝑚−2𝑢_{𝑚−1𝑘+1}=𝑓₂𝑓₃⋮𝑓_𝑚−2𝑓_𝑚−1,𝑓₂=−𝑏𝑢^𝑘₁+𝑢₁^𝑘+1+𝑐𝑢̂^𝑘₂−𝑑𝑢^𝑘₃−𝑒𝑢^𝑘₄,𝑓₃=−𝑎𝑢^𝑘₁+𝑢₁^𝑘+1−𝑏𝑢^𝑘₂+𝑐𝑢̂^𝑘₃−𝑑𝑢^𝑘₄−𝑒𝑢^𝑘₅,𝑓_𝑖=−𝑎𝑢^𝑘_𝑖−2−𝑏𝑢^𝑘_𝑖−1+𝑐𝑢̂^𝑘_𝑖−𝑑𝑢^𝑘_𝑖+1−𝑒𝑢^𝑘_𝑖+2𝑓for𝑖=4,5,…,𝑚−3,_𝑚−2=−𝑎𝑢^𝑘_𝑚−4−𝑏𝑢^𝑘_𝑚−3+𝑐𝑢̂^𝑘_𝑚−2−𝑑𝑢^𝑘_𝑚−1−𝑒𝑢^𝑘_𝑚+𝑢_𝑚^𝑘+1,𝑓_𝑚−1=−𝑎𝑢^𝑘_𝑚−3−𝑏𝑢^𝑘_𝑚−2+𝑐𝑢̂^𝑘_𝑚−1−𝑑𝑢^𝑘_𝑚+𝑢_𝑚^𝑘+1. (8)(9)

The evaluations of 𝑓₂, 𝑓₃, 𝑓_𝑚−2, and 𝑓_𝑚−1 require the values of𝑢 at the boundaries 𝑖=1 and 𝑖=𝑚. However, these values cannot be obtained numerically because their computations involve nodes at 𝑖=−1  and 𝑖=𝑚+2, which are exterior to the considered solution domain. If the exact solutions of 𝑢 are available at 𝑖=1 and 𝑖=𝑚, then it is appropriate to consider them as the required boundary values. Otherwise, the boundary conditions have to be formulated, bearing in mind that they should be of comparable accuracy [1].

The IADEMF4 scheme secondly employs the higher-order accuracy formula of MF [19]. The variant, whose accuracy is of the order 𝑂((Δ𝑡)²+(Δ𝑥)⁴), is as given in (10) and (11) 𝑟𝐼+𝐺₁𝐮^(𝑝+1/2)=𝑟𝐼−𝑔𝐺₂𝐮^(𝑝)+𝐟,𝑟𝐼+𝐺₂𝐮^(𝑝+1)=𝑟𝐼−𝑔𝐺₁𝐮^(𝑝+1/2)+𝑔𝐟, (10)(11) where 𝑟, 𝑝, and 𝐼 represent an acceleration parameter, the iteration index, and an identity matrix, respectively. 𝐺₁ and 𝐺₂ are two constituent matrices. The vectors 𝐮^(𝑝+1) and 𝐮^(𝑝+1/2) represent the required solution at the iteration level (𝑝+1) and at some intermediate level (𝑝+1/2), respectively. The relation of 𝑔 and 𝑟 is given by 𝑔=(6+𝑟)/6, 𝑟>0.

Substitute the following expression obtained from (10) into (11) 𝐮^(𝑝+1/2)=𝑟𝐼+𝐺₁⁻¹𝑟𝐼−𝑔𝐺₂𝐮^(𝑝)+𝐟, (12) and then simplify to obtain 𝑟𝐼+𝐺₂𝐮^(𝑝+1)=𝑟𝐼−𝑔𝐺₁𝑟𝐼+𝐺₁⁻¹𝑟𝐼−𝑔𝐺₂𝐮^(𝑝),+𝑟𝐼−𝑔𝐺₁𝑟𝐼−𝑔𝐺₁⁻¹+𝑔𝐟. (13)As𝑝in (13) becomes sufficiently large, the temperature solution reaches a steady state that is, as 𝑝→∞, then 𝐮^(𝑝+1)→𝐮 and 𝐮^(𝑝)→𝐮. Simplify the equation, and then multiply by (𝑟𝐼+𝐺₁). After some algebraic manipulations, the following expression is obtained: 𝑟𝐺₁+𝐺₂1+𝑔+1−𝑔²𝐺₁𝐺₂𝐮=𝑟1+𝑔𝐟. (14)Multiply (14) by 1/𝑟(1+𝑔), and use the definition of 𝑔 to finally obtain the following form: 𝐺₁+𝐺₂−16𝐺₁𝐺₂𝐮=𝐟. (15)

The comparison between the expression 𝐴𝐮=𝐟 in (8) and the form given in (15) suggests that the coefficient matrix𝐴for the IADEMF4 can be decomposed into 𝐴=𝐺₁+𝐺₂−16𝐺₁𝐺₂. (16)The IADEMF4 requires the constituent matrices 𝐺₁ and 𝐺₂ in (16) to be in the form of lower and upper tridiagonal matrices, respectively, in order to retain the pentadiagonal structure of 𝐴. Thus, 𝐺₁=1𝑙₁𝑚1𝑂₁𝑙₂⋱𝑚₂⋱⋱⋱𝑙_𝑚−41𝑂𝑚_𝑚−4𝑙_𝑚−31_{(𝑚−2)×(𝑚−2)},𝐺₂=𝑒₁̂𝑢₁𝑣₁̂𝑒₂̂𝑢₂𝑣₂̂𝑂𝑂𝑒̂⋱⋱⋱_𝑚−4𝑢_𝑚−4𝑣̂_𝑚−4𝑒̂_𝑚−3𝑢_𝑚−3𝑒̂_{𝑚−2(𝑚−2)×(𝑚−2)}. (17)

By substituting 𝐺₁ and 𝐺₂ in (17) into the formula for the decomposition of 𝐴, the entries of the resultant matrix are compared with the entries of 𝐴 in (8), yielding the following definitions: 𝑒₁̂=6(𝑐−1)5,𝑢₁=6𝑑5,𝑙₁=6𝑏6−𝑒₁̂,𝑒₂̂=6(𝑐−1)+𝑙₁𝑢₁5,𝑣_𝑖̂=6𝑒5for𝑖=1,2,…,𝑚−4. (18) And for 𝑖=2,3,…,𝑚−3, 𝑢_𝑖=6𝑑+𝑙_𝑖−1𝑣̂_𝑖−15,𝑚_𝑖−1=6𝑎6−𝑒̂_𝑖−1,𝑙_𝑖=6𝑏+𝑚_𝑖−1𝑢_𝑖−16−𝑒_𝑖̂,𝑒̂_𝑖+1=6(𝑐−1)+𝑙_𝑖𝑢_𝑖+𝑚_𝑖−1𝑣̂_𝑖−15. (19)

Since 𝐺₁ and 𝐺₂ are three banded matrices, then (𝑟𝐼+𝐺₁) and (𝑟𝐼+𝐺₂) can be inverted easily. The equations in (10) and (11) are rearranged as in (20) and (21), respectively, 𝐮^(𝑝+1/2)=𝑟𝐼+𝐺₁⁻¹𝑟𝐼−𝑔𝐺₂𝐮^(𝑝)+𝑟𝐼+𝐺₁⁻¹𝐮𝐟,_(𝑝+𝟏)=𝑟𝐼+𝐺₂⁻¹𝑟𝐼−𝑔𝐺₁𝐮^(𝑝+1/2)+𝑔𝑟𝐼+𝐺₂⁻¹𝐟. (20)(21)The above two equations are computed and simplified, leading to the computational formulae at each of the half iteration levels as given in (22) and (23). (i)At the (𝑝+1/2) iteration level, 𝑢_𝑖^(𝑝+1/2)=1𝑅𝐸_𝑖−1𝑢_𝑖^(𝑝)+𝑊_𝑖−1𝑢^(𝑝)_𝑖+1+𝑉_𝑖−1𝑢^(𝑝)_𝑖+2−𝑚_𝑖−3𝑢^(𝑝+1/2)_𝑖−2−𝑙_𝑖−2𝑢^(𝑝+1/2)_𝑖−1+𝑓_𝑖,𝑖=2,3,…,𝑚−2,𝑚−1 (22)(ii)At the (𝑝+1) iteration level, 𝑢_𝑖^(𝑝+1)=1𝑍_𝑖−1𝑆_𝑖−3𝑢^(𝑝+1/2)_𝑖−2+𝑄_𝑖−2𝑢^(𝑝+1/2)_𝑖−1+𝑃𝑢_𝑖^(𝑝+1/2)−𝑢_𝑖−1𝑢^(𝑝+1)_𝑖+1−𝑣̂_𝑖−1𝑢^(𝑝+1)_𝑖+2+𝑔𝑓_𝑖,𝑖=𝑚−1,𝑚−2,…,3,2 (23)

with 𝑚₋₁=𝑚₀=𝑙₀=𝑉_𝑚−2=𝑉_𝑚−3=𝑊_𝑚−2=𝑢_𝑚−2=𝑣̂_𝑚−2=𝑣̂_𝑚−3=𝑄₀=𝑆₋₁=𝑆₀𝐸=0,𝑅=1+𝑟,𝑃=𝑟−𝑔,_𝑖=𝑟−𝑔𝑒_𝑖̂,𝑍_𝑖=𝑟+𝑒_𝑖̂𝑊,𝑖=1,2,…,𝑚−2,_𝑖=−𝑔𝑢_𝑖,𝑄_𝑖=−𝑔𝑙_𝑖𝑉,𝑖=1,2,…,𝑚−3,_𝑖=−𝑔𝑣_𝑖̂,𝑆_𝑖=−𝑔𝑚_𝑖,𝑖=1,2,…,𝑚−4. (24)

The IADEMF4 algorithm is regarded as a two-stage process involving two iteration levels, (𝑝+1/2) and (𝑝+1). It is completed explicitly by using the required equations at the two levels in alternate sweeps along all the grid points in the interval (0,1) until convergence is reached. In (22), the calculation to determine the unknown 𝑢_𝑖^(𝑝+1/2) begins at the left boundary and then moves to the right. In a similar manner, 𝑢_𝑖^(𝑝+1) in (23) is calculated by proceeding from the right boundary towards the left (Figure 1).

Figure 1 

The two-stage IADEMF4 algorithm. The directions of the sweeps at the (𝑝+1/2) and (𝑝+1) iteration levels.

The computational molecules are depicted in Figures 2 and 3.

Figure 2 

Computational molecule of the IADEMF4 at the (𝑝+1/2) iteration level.

Figure 3 

Computational molecule of the IADEMF4 at the (𝑝+1) iteration level.

At each level of iteration, the computational molecules involve two known grid points at the new level and another three known ones at the old level. Clearly, the method is explicit.

4. Convergence Analysis of the IADEMF4

This section proves the convergence of the IADEMF4. Since 𝜆>1 will not guarantee an accurate approximation for 𝜕𝑈/𝜕𝑡 [18], the values of 𝜆 that are considered appropriate for the proof are 0<𝜆≤1.

From the definitions in ((18) and (19)), the following results are obtained: 𝑒₁̂=6(𝑐−1)5=3𝜆2,implyingthat0<𝑒₁̂≤32,𝑢₁=6𝑑5<0,𝑙₁=6𝑏6−𝑒₁̂<0,𝑣_𝑖̂=6𝑒5𝑢>0for𝑖=1,2,…,𝑚−4,₂=6𝑑+𝑙₁𝑣₁̂5=𝑢₁+𝑙₁𝑣₁̂5<𝑢₁,since𝑙₁𝑣₁̂𝑒̂<0,₂=6(𝑐−1)+𝑙₁𝑢₁5=𝑒₁̂+𝑙₁𝑢₁5>𝑒₁̂,since𝑙₁𝑢₁>0. (25) Simple computation shows that 𝑙₁𝑢₁/5<1. Clearly, 𝑒₂̂<6.

Lemma 1. 𝑒_𝑖̂ in ((18) and (19)) is such that, 0<𝑒_𝑖̂<6, for all 𝑖=1,2,3,…,𝑚−2.

Proof. The results in (25) show that 6>𝑒₂̂>𝑒₁̂>0. Assume it is true that 6>𝑒_𝑘̂>𝑒̂_𝑘−1>0 for all 𝑘=3,4,5,…,𝑚−3. Since the assumption implies that 6>𝑒̂_𝑘−1>𝑒̂_𝑘−2>0, then 𝑚_𝑘−1=6𝑎/(6−𝑒̂_𝑘−1)>𝑚_𝑘−2=6𝑎/(6−𝑒̂_𝑘−2)>0. Therefore, 𝑚_𝑘−1𝑣̂_𝑘−1>𝑚_𝑘−2𝑣̂_𝑘−2.
It has also been shown that 𝑢₂<𝑢₁<0. Assume it is true that 𝑢_𝑗<𝑢_𝑗−1<0for all 𝑗=3,4,…,𝑚−4. The assumption implies that 𝑢_𝑗−1<𝑢_𝑗−2<0. Thus, 𝑚_𝑗−1𝑢_𝑗−1<𝑚_𝑗−2𝑢_𝑗−2<0𝑢_𝑗+1=6𝑑+𝑙_𝑗𝑣_𝑗̂5=𝑢₁+𝑙_𝑗𝑣_𝑗̂5=𝑢₁+6𝑏+𝑚_𝑗−1𝑢_𝑗−1𝑣_𝑗̂56−𝑒_𝑗̂<𝑢₁+6𝑏+𝑚_𝑗−2𝑢_𝑗−2𝑣̂_𝑗−156−𝑒̂_𝑗−1=𝑢₁+𝑙_𝑗−1𝑣̂_𝑗−15=𝑢_𝑗<0. (26)
Since it is true that, for 𝑗+1, 𝑢_𝑗+1<𝑢_𝑗<0, then, by induction, the assumption 𝑢_𝑗<𝑢_𝑗−1<0 is true for all 𝑗. An equivalent to the preceding statement would be 𝑢_𝑘<𝑢_𝑘−1<0 for all 𝑘=3,4,5,…,𝑚−3. It follows that 𝑙_𝑘𝑢_𝑘=6𝑏+𝑚_𝑘−1𝑢_𝑘−1𝑢_𝑘6−𝑒_𝑘̂>6𝑏+𝑚_𝑘−2𝑢_𝑘−2𝑢_𝑘−16−𝑒̂_𝑘−1=𝑙_𝑘−1𝑢_𝑘−1,𝑒̂_𝑘+1=6(𝑐−1)+𝑙_𝑘𝑢_𝑘+𝑚_𝑘−1𝑣̂_𝑘−15=𝑒₁̂+𝑙_𝑘𝑢_𝑘+𝑚_𝑘−1𝑣̂_𝑘−15>𝑒₁̂+𝑙_𝑘−1𝑢_𝑘−1+𝑚_𝑘−2𝑣̂_𝑘−25=𝑒_𝑘̂>0. (27)
Since it is true that, for 𝑘+1, 𝑒̂_𝑘+1>𝑒_𝑘̂>0, then, by induction, 𝑒_𝑘̂>𝑒̂_𝑘−1>0 is also true for all 𝑘.
Suppose there is a 𝑘 such that 𝑒_𝑘̂>6. Then, 𝑒̂_𝑘+1=𝑒₁̂+𝑙_𝑘𝑢_𝑘+𝑚_𝑘−1𝑣̂_𝑘−15=𝑒₁̂+6𝑏+𝑢_𝑘56−𝑒_𝑘̂+𝑚_𝑘−15𝑢_𝑘−1𝑢_𝑘6−𝑒_𝑘̂+𝑣̂_𝑘−1<𝑒₁̂. (28)
This is a contradiction since 𝑒̂_𝑘+1>𝑒₁̂. This verifies that 6>𝑒_𝑘̂>𝑒̂_𝑘−1>0, 𝑘=3,4,5,…,𝑚−3. Let 𝑖=2,3,4,…,𝑚−2, and then 6>𝑒_𝑖̂>𝑒̂_𝑖−1>0. Therefore, 0<𝑒_𝑖̂<6,𝑖=1,2,3,…,𝑚−2. Lemma 1 is proved.

Lemma 2. If 𝑟>0 and 𝑔=(6+𝑟)/6, then max|(𝑟−𝑔𝑒_𝑖̂)/(𝑟+𝑒_𝑖̂)|<1,i=1,2,3,…,𝑚−2.

Proof. Let max|(𝑟−𝑔𝑒_𝑖̂)/(𝑟+𝑒_𝑖̂)|=|(𝑟−𝑔𝑒_𝑗̂)/(𝑟+𝑒_𝑗̂)| for some𝑗.
Assume |(𝑟−𝑔𝑒_𝑗̂)/(𝑟+𝑒_𝑗̂)|≥1. Since 𝑟>0 and 𝑒_𝑗̂>0, then 𝑟+𝑒_𝑗̂>0.
If (𝑟−𝑔𝑒_𝑗̂)/(𝑟+𝑒_𝑗̂)≥1, then −((6+𝑟)/6)𝑒_𝑗̂≥𝑒_𝑗̂, which implies that 𝑟≤−12. This contradicts the fact that 𝑟>0.
If (𝑟−𝑔𝑒_𝑗̂)/(𝑟+𝑒_𝑗̂)≤−1, then 2𝑟−((6+𝑟)/6)𝑒_𝑗̂≤−𝑒_𝑗̂, which implies that 𝑒_𝑗̂≥12. This contradicts Lemma 1.
So, the assumption that |(𝑟−𝑔𝑒_𝑗̂)/(𝑟+𝑒_𝑗̂)|≥1 is false. Hence, max|(𝑟−𝑔𝑒_𝑖̂)/(𝑟+𝑒_𝑖̂)|=|(𝑟−𝑔𝑒_𝑗̂)/(𝑟+𝑒_𝑗̂)|<1.

Proposition 3. ‖(𝑟𝐼−𝑔𝐺₁)(𝑟𝐼+𝐺₁)⁻¹‖₂<1.

Proof. Let 𝐹=(𝑟𝐼−𝑔𝐺₁)(𝑟𝐼+𝐺₁)⁻¹; that is, 𝐹=𝑟−𝑔𝑑⋱𝑟−𝑔𝑑𝑂⋱⋱⋱⋮⋮⋱⋱⋯⋯⋯⋯𝑟−𝑔𝑑. (29)𝐹 is a lower triangular matrix, with all the diagonal entries (whence the eigenvalues of 𝐹) equal to (𝑟−𝑔)/𝑑, where 𝑑=𝑟+1. Denote all the eigenvalues of 𝐹 by 𝜆_𝐹. If 𝜌[𝐹] is defined as the spectral radius of 𝐹, then 𝜌[𝐹]=max𝜆_𝐹=𝜆_𝐹=𝑟−𝑔≤‖𝐹‖𝑟+1². (30)But, by definition of 2-norm, ‖𝐹‖₂=max_{‖𝐱‖2≠0}‖𝐹𝐱‖₂‖𝐱‖₂=max_{‖𝐲‖2=1}‖𝐹𝐲‖₂‖𝐲‖₂=max_{‖𝐲‖2=1}‖𝐹𝐲‖₂≤max_{‖𝐰‖2=1}‖𝐹𝐰‖₂=max_{‖𝐰‖2=1}𝜆_𝐹𝐰₂. (31)Since all the eigenvalues of 𝐹 are equal, then ‖𝐹‖₂≤𝜆_𝐹max_{‖𝐰‖2=1}‖𝐰‖₂=𝜆_𝐹=𝜌[𝐹]. (32)Thus, from (30) and (32), ‖𝐹‖₂=|𝜆_𝐹|=|(𝑟−𝑔)/(𝑟+1)|. By Lemma 2 with 𝑒_𝑗 replaced by 1, |(𝑟−𝑔)/(𝑟+1)|<1 is obtained, leading to the result of Proposition 3, which is ‖𝐹‖₂=‖(𝑟𝐼−𝑔𝐺₁)(𝑟𝐼+𝐺₁)⁻¹‖₂<1.

Proposition 4. ‖(𝑟𝐼−𝑔𝐺₂)(𝑟𝐼+𝐺₂)⁻¹‖₂<1.

Proof. Let𝐾=(𝑟𝐼−𝑔𝐺₂)(𝑟𝐼+𝐺₂)⁻¹; that is, 𝑒̂𝐾=𝑟−𝑔₁𝑧₁𝑒̂⋱⋯⋯⋯𝑟−𝑔₂𝑧₂𝑒̂⋱⋯⋮⋱⋮⋮𝑂⋱⋱⋮𝑟−𝑔_𝑚𝑧_𝑚. (33)𝐾 is an upper triangular matrix, with all the diagonal entries equal to (𝑟−𝑔𝑒_𝑖̂)/𝑧_𝑖,𝑖=1,2,…,𝑚−2, where 𝑧_𝑖=𝑟+𝑒_𝑖̂. Since all the eigenvalues of 𝐾 are distinct, then 𝐾 is similar to a diagonal matrix 𝐷.
By the Schur triangularization theorem [20], there is an orthogonal matrix 𝑂 such that 𝑂^𝑇𝐾𝑂=𝐷. The diagonal entries of 𝐷 are the eigenvalues of 𝐾‖𝐾‖₂=𝑂^𝑇𝐾𝑂₂=‖𝐷‖₂=maximumeigenvalueof𝐷^𝑇𝐷=maximumeigenvalueof𝐷²=𝜌𝐷²=𝜌𝐾²=𝜌(𝐾). (34) By Lemma 2, ‖𝐾‖₂(𝐾)𝑒̂=𝜌=max𝑟−𝑔_𝑖𝑟+𝑒_𝑖̂=𝑟−𝑔𝑒_𝑗̂𝑟+𝑒_𝑗̂<1, (35) for some 𝑗. Hence, this proves Proposition 4.

From (20) and (21), 𝐮^(𝑝+1)=𝑟𝐼+𝐺₂⁻¹𝑟𝐼−𝑔𝐺₁×𝑟𝐼+𝐺₁⁻¹𝑟𝐼−𝑔𝐺₂𝐮^(𝑝)+𝑟𝐼+𝐺₁⁻¹𝐟+𝑔𝑟𝐼+𝐺₂⁻¹𝐟. (36)

Let𝑀(𝑟)=𝑟𝐼+𝐺₂⁻¹𝑟𝐼−𝑔𝐺₁𝑟𝐼+𝐺₁⁻¹𝑟𝐼−𝑔𝐺₂, (37) And let 𝑞(𝑟)=𝑟𝐼+𝐺₂⁻¹𝑟𝐼−𝑔𝐺₁𝑟𝐼+𝐺₁⁻¹+𝑔𝑟𝐼+𝐺₂⁻¹𝐟. (38)

Then, 𝐮^(𝑝+1)=𝑀(𝑟)𝐮^(𝑝)+𝑞(𝑟).

Theorem 5. The IADEMF4 is convergent if 𝜌[𝑀(𝑟)]<1, for 𝑟>0.

Proof. Define 𝑀(𝑟)=(𝑟𝐼+𝐺₂)𝑀(𝑟)(𝑟𝐼+𝐺₂)⁻¹; then 𝑀(𝑟)=𝑟𝐼−𝑔𝐺₁𝑟𝐼+𝐺₁⁻¹𝑟𝐼−𝑔𝐺₂𝑟𝐼+𝐺₂⁻¹. (39)Thus, by similarity, 𝑀(𝑟) and 𝑀(𝑟) have the same set of eigenvalues. Therefore, 𝜌[𝑀(𝑟𝑀(𝑟≤𝑀(𝑟)]=𝜌))₂=𝑟𝐼−𝑔𝐺₁𝑟𝐼+𝐺₁⁻¹𝑟𝐼−𝑔𝐺₂𝑟𝐼+𝐺₂⁻¹₂≤𝑟𝐼−𝑔𝐺₁𝑟𝐼+𝐺₁⁻¹₂𝑟𝐼−𝑔𝐺₂𝑟𝐼+𝐺₂⁻¹₂<1. (40)The last inequality is due to Propositions 3 and 4. The proven Theorem 5 assures the convergence of the IADEMF4.

5. The Fourth-Order Classical Iterative Methods

The system of linear equations in (8) may also be solved by using the classical iterative methods of the fourth order. They include the JAC4, the GS4, and the SOR4 methods. These iterative methods are capable of exploiting the sparse structure of the pentadiagonal matrix.

The JAC4 algorithm can be represented by 𝑢_𝑖^(𝑝+1)=𝑓_𝑖−𝑎𝑢^(𝑝)_𝑖−2−𝑏𝑢^(𝑝)_𝑖−1−𝑑𝑢^(𝑝)_𝑖+1−𝑒𝑢^(𝑝)_𝑖+2𝑐,𝑖=2,3,…,𝑚−1. (41)

The approximation of 𝑢_𝑖^(𝑝+1) at the (𝑝+1)th iteration level is computed using the relevant values in the𝑝th iteration.

The GS4 method uses the most recent values of 𝑢^(𝑝+1)_𝑖−2 and 𝑢^(𝑝+1)_𝑖−1 to update the approximation value of 𝑢_𝑖^(𝑝+1). The algorithm for the GS4 is as expressed in (42) 𝑢_𝑖^(𝑝+1)=𝑓_𝑖−𝑎𝑢^(𝑝+1)_𝑖−2−𝑏𝑢^(𝑝+1)_𝑖−1−𝑑𝑢^(𝑝)_𝑖+1−𝑒𝑢^(𝑝)_𝑖+2𝑐,𝑖=2,3,…,𝑚−1. (42)

The SOR4 iterative method accelerates the convergence rate of the GS4. If𝜔is a relaxation parameter, then for any 𝜔≠0, (42) can be rewritten as 𝑢_𝑖^(𝑝+1)=(1−𝜔)𝑢_𝑖^(𝑝)+𝜔𝑓_𝑖−𝑎𝑢^(𝑝+1)_𝑖−2−𝑏𝑢^(𝑝+1)_𝑖−1−𝑑𝑢^(𝑝)_𝑖+1−𝑒𝑢^(𝑝)_𝑖+2𝑐,𝑖=2,3,…,𝑚−1. (43)

The SOR4 algorithm reduces to the GS4 if 𝜔=1. Except for very special cases, it is difficult to obtain the analytic expression for 𝜔. According to Young [21], the precise determination of the optimal 𝜔 in the SOR is only known for a small class of matrices. The value of 𝜔 generally lies in the range of 1<𝜔<2.

The computational molecules for the JAC4 and the GS4 are as illustrated in Figures 4 and 5, respectively. The SOR4 has the same form as the GS4.

Figure 4 

Computational molecule of the JAC4.

Figure 5 

Computational molecule of the GS4/SOR4.

6. Computational Complexity

The cost of implementing an algorithm can be assessed by examining its computational complexity. The number of arithmetic operations such as additions (including subtractions) and multiplications (including divisions) that is needed to perform by the algorithm can be straightforwardly counted. Table 1 gives the number of sequential arithmetic operations per iteration that is required to evaluate the algorithms.

Table 1 

Sequential arithmetic operations per iteration (𝑚: problem size).

For the higher-order schemes, trade-off between accuracy and speed usually happens. The computational complexity has an effect on the efficiency of a particular scheme. Thus, this factor will be taken into account when discussing the results of the numerical experiments conducted in this study.

7. The Pseudocode of the IADEMF4 Sequential Algorithm

Algorithm 1 illustrates the pseudocode of the IADEMF4 sequential algorithm. This algorithm can also be generalized and applied to other methods discussed in this paper.

begin   determine the parameters 𝑚,𝑟,Δ𝑡,Δ𝑥,𝜆,𝜀,𝑇   determine initial conditions 𝑢𝑖(𝑝)   determine exact solutions 𝑈𝑖 at time level 𝑇   while (time level <𝑇)      determine boundary conditions      for 𝑖=2 to 𝑖=𝑚−1         compute 𝑓𝑖 (refer to (9))      end-for         set iteration = 0         while (convergence conditions are not satisfied)          for 𝑖=2 to 𝑖=𝑚−1         compute 𝑢𝑖(𝑝+1/2) (refer to (22))          end-for          for 𝑖=2 to 𝑖=𝑚−1         compute 𝑢𝑖(𝑝+1) (refer to (23))          end-for          test for convergence:          compute 𝑒𝑖←𝑢𝑖(𝑝+1)−𝑢𝑖(𝑝)          if (max𝑒𝑖<𝜀)          then 𝑢𝑖(𝑝)←𝑢𝑖(𝑝+1)          add 1 to iteration (if necessary)         end-while   end-while   for 𝑖=2 to 𝑖=𝑚−1      determine average absolute error, root mean      square error and maximum error   end-for end

Algorithm 1 

The IADEMF4 sequential algorithm.

8. Numerical Experiments

Two experiments were conducted to test the sequential numerical performance of the proposed IADEMF4 method against those of the benchmarked classical iterative methods, namely, the JAC4, the GS4, and the SOR4. Comparison is also made with the corresponding IADE method of the second order.

Experiment 1. This problem was taken from Saulev [22] 𝜕𝑈=𝜕𝜕𝑡²𝑈𝜕𝑥²,0≤𝑥≤1 (44) subject to the initial condition 𝑈(𝑥,0)=4𝑥(1−𝑥),0≤𝑥≤1 (45) and the boundary conditions 𝑈(0,𝑡)=𝑈(1,𝑡)=0,𝑡≥0. (46) The exact solution to the given problem is given by 𝑈(𝑥,𝑡)=32𝜋^3∞_(2)𝑘=1,1𝑘³𝑒^{−𝜋2𝑘2𝑡}sin(𝑘𝜋𝑥). (47)

Experiment 2. This problem was taken from Johnson and Reiss [23], 𝜕𝑈=𝜕𝜕𝑡²𝑈𝜕𝑥²,0≤𝑥≤1 (48) subject to the initial condition 𝑈(𝑥,0)=sin(𝜋𝑥)(1+6cos(𝜋𝑥)),0≤𝑥≤1 (49) and the boundary conditions 𝑈(0,𝑡)=𝑈(1,𝑡)=0,𝑡≥0. (50) The exact solution to the given problem is given by 𝑈(𝑥,𝑡)=sin(𝜋𝑥)𝑒^−𝜋2𝑡+3sin(2𝜋𝑥)𝑒^{−4𝜋2𝑡}. (51)

9. Results and Discussion

In each experiment, the exact solutions at 𝑖=1 and 𝑖=𝑚 were taken as boundary values for the IADEMF4 and the other fourth-order classical iterative methods. The convergence criterion used in the testing of each method was taken as ‖𝑢^(𝑝+1)−𝑢^(𝑝)‖_∞≤𝜀, where𝜀is the convergence tolerance. The selections of the optimum𝑟 or 𝜔were determined by experiments. As for the IADEMF2, the CN scheme with 𝜃=1/2 was employed.

Figures 6, 7, 8, and 9 visualize the behavior of the one-dimensional parabolic heat solutions for both experiments. By using 𝑚=10 and a tolerance requirement of 𝜀=10⁻⁴, two different mesh sizes, 𝜆=0.5 and 𝜆=1.0, were considered in each experiment. The exact solutions are compared with the IADEMF4 and the IADEMF2 numerical solutions. The optimum value of𝑟chosen for each method, as well as their corresponding outcome of the number of iterations, 𝑛, is stated in the legend of each figure. Every figure reveals that the IADEMF4 is more accurate than the IADEMF2 and the former converges with fewer numbers of iterations in comparison to the latter. The numerical solutions of the IADEMF4 seem to be in very good agreement with the exact solutions. For example, in Figure 7, at 𝑥=0.5, the difference between the exact solution and the numerical solution using the IADEMF2 and using the IADEMF4 is about 3.6% and 0%, respectively. These results imply that the accuracy and convergence rate of the second-order IADE method are enhanced by the implementation of the fourth-order CN approximation that leads to the formation of the corresponding fourth-order IADE scheme.

Figure 6 

Numerical and exact solutions for Experiment 1. 𝜆=0.5, Δ𝑥=0.1, Δ𝑡=0.005, 𝑡=0.25, and 𝜀=10−4.

Figure 7 

Numerical and exact solutions for Experiment 1. 𝜆=1.0, Δ𝑥=0.1, Δ𝑡=0.01, 𝑡=0.5, and 𝜀=10−4.

Figure 8 

Numerical and exact solutions for Experiment 2. 𝜆=0.5, Δ𝑥=0.1, Δ𝑡=0.005, 𝑡=0.25, and 𝜀=10−4.

Figure 9 

Numerical and exact solutions for Experiment 2. 𝜆=1.0, Δ𝑥=0.1, Δ𝑡=0.01, 𝑡=0.5, and 𝜀=10−4.

Figure 10 displays the graph of log (RMSE) versus log(Δ𝑥) for decreasing values of Δ𝑥 implemented on Experiment 2. By considering 𝜀=10⁻⁴ and fixing Δ𝑡=0.0001, the value of Δ𝑥 was initially taken as Δ𝑥=0.125. It was then successively halved into Δ𝑥/2, Δ𝑥/4, and Δ𝑥/8.

Figure 10 

Log(RMSE) versus Log(Δ𝑥) Experiment 2. Δ𝑡=0.0001, 𝑡=0.002, and 𝜀=10−4.

The figure shows that amongst the tested methods, the root mean square error (RMSE) of the IADEMF4 and the IADEMF2 decreases linearly as the values of Δ𝑥 decreases. The slope of the IADEMF4 is approximately equal to 4, which corresponds to its fourth-order spatial accuracy. The IADEMF2 with a second-order spatial accuracy has a slope that is approximately equal to 2. It is clear that the IADEMF4 is always more accurate than the IADEMF2, for the different considered mesh sizes. The graphs of the SOR4 (𝜔=1.05), the GS4, and the JAC4 show that their accuracies of fourth order tend to lack as Δ𝑥 decreases, largely due to the effect of increasing round-off errors as the value of 𝑚 increases.

Tables 2, 3, 4, 5, 6, 7, 8, and 9 provide numerical results in terms of the average absolute error (AAE), root mean square error (RMSE), maximum error (ME), number of iterations (𝑛), and execution time (ET) measured in seconds (s). The results are obtained from both experiments for two different values of 𝑚 and mesh size, 𝜆. It is generally observed that when 𝑚=700, 𝜆=0.5, and 𝜀=10⁻⁶, the IADEMF4 has the least average absolute error, root mean square error, and maximum error in comparison with the other methods under consideration (Tables 2–5). When the size of 𝑚 was ten times bigger (𝑚=7000) and a more stringent tolerance criterion was set (𝜀=10⁻¹⁰), the accuracy of the IADEMF4 clearly outperforms the other methods for a mesh size of 𝜆=0.5 (Tables 6 and 8). However, the difference in the errors amongst the tested methods is not so obvious for the case of 𝜆=1.0 (Tables 7 and 9). The achievement of the fourth-order IADE method in Experiments 1 and 2 can be clearly seen from the results in Tables 2 and 8, respectively, where it has caused a huge 94% reduction of RMSE from its corresponding second-order IADE method.

Table 2 

Experiment 1. 𝑚=700, 𝜆=0.5, Δ𝑥=1.43×10⁻³, Δ𝑡=1.02×10⁻⁶, 𝑡=5.09×10⁻⁵, and 𝜀=10⁻⁶.

Table 3 

Experiment 1. 𝑚=700, 𝜆=1.0, Δ𝑥=1.43×10⁻³, Δ𝑡=2.04×10⁻⁶, 𝑡=5.09×10⁻⁵, and 𝜀=10⁻⁶.

Table 4 

Experiment 2. 𝑚=700, 𝜆=0.5, Δ𝑥=1.43×10⁻³, Δ𝑡=1.02×10⁻⁶, 𝑡=5.09×10⁻⁵, and 𝜀=10⁻⁶.

Table 5 

Experiment 2. 𝑚=700, 𝜆=1.0, Δ𝑥=1.43×10⁻³, Δ𝑡=2.04×10⁻⁶, 𝑡=5.09×10⁻⁵, and 𝜀=10⁻⁶.

Table 6 

Experiment 1. 𝑚=7000, 𝜆=0.5, Δ𝑥=1.43×10⁻⁴, Δ𝑡=1.02×10⁻⁸, 𝑡=5.10×10⁻⁷, and 𝜀=10⁻¹⁰.

Table 7 

Experiment 1. 𝑚=7000, 𝜆=1.0, Δ𝑥=1.43×10⁻⁴, Δ𝑡=2.04×10⁻⁸, 𝑡=5.1×10⁻⁷, and 𝜀=10⁻¹⁰.

Table 8 

Experiment 2. 𝑚=7000, 𝜆=0.5, Δ𝑥=1.43×10⁻⁴, Δ𝑡=1.02×10⁻⁸, 𝑡=5.10×10⁻⁷, and 𝜀=10⁻¹⁰.

Table 9 

Experiment 2. 𝑚=7000, 𝜆=1.0, Δ𝑥=1.43×10⁻⁴, Δ𝑡=2.04×10⁻⁸, 𝑡=5.1×10⁻⁷, and 𝜀=10⁻¹⁰.

The IADEMF4 has the advantage of featuring higher accuracy due to the execution of the fourth-order CN approximation coupled with the fourth-order accurate MF variant. The IADEMF2 is only derived from the second-order CN approximation, but its combination with its corresponding fourth-order MF variant managed to produce errors which are better than the GS4, JAC4, and the SOR4. Even though the classical iterative methods are also derived from the fourth-order accurate CN type approximation, they are lacking in accuracy due to the round-off errors that have been accumulated from the time the execution starts till it ends.

With regards to the rate of convergence, the results from each table demonstrate that the number of iterations produced by the IADEMF2 and the fourth-order classical iterative methods is greater than or at least equal to that of the IADEMF4. The convergence rate of the latter surpasses the others in the case of 𝑚=7000 and 𝜀=10⁻¹⁰ (Tables 6–8). Even though the operational count of the IADEMF4 is relatively quite large (Table 1), due to its higher level of accuracy, its increasing number of correct digits at each iteration causes it to converge at a faster rate. The IADEMF2 has the advantage of having less mathematical operations; thus, it is competitive in terms of convergence rate but at the expense of accuracy. The application of the fourth-order CN approximation on the heat equation proves that the three classical iterative methods also converge, with JAC4 appearing to be the slowest and the least accurate amongst all.

It is found that the convergence rate of the fourth-order and second-order methods improves with the application of larger mesh size, that is, from 𝜆=0.5 to 𝜆=1.0. The coarser meshes cause a reduction in the computational operations, thus giving better rate of convergence. However, a more accurate solution is obtained by using finer mesh. For example, Table 6 shows that, for 𝜆=0.5, the IADEMF4 has RMSE=4.33𝑒−10 and 𝑛=173, whereas Table 7 shows that, for 𝜆=1.0, its RMSE=1.63𝑒−7 and 𝑛=112. In general, amongst all the tested methods, the IADEMF4 still maintains its greater accuracy characteristic, even with coarser meshes.

In terms of execution time, the results from every table display shorter execution time for the IADEMF2 in comparison to the IADEMF4. This is expected, since the IADEMF2 has lower computational complexity (Table 1). Despite the achievement in accuracy, the IADEMF4 has to perform more computational work than the IADEMF2 since the former has to utilize values of 𝑢 at two grid points on either side of the point (𝑖Δ𝑥,𝑘Δ𝑡) along the 𝑘th time level. Thus, if accuracy is desired, then the preferred sequential numerical algorithm would be the IADEMF4. On the other hand, if execution time matters, then the choice would be the IADEMF2.

Amongst the fourth-order methods, the IADEMF4 executes in the least amount of time. Even though its computational complexity is relatively large, it operates with the least number of iterations, thus enabling it to be the most efficient.

10. Conclusion

This paper proposes the development of the novel fourth-order IADEMF4 finite difference scheme. The higher-order scheme is developed by representing the spatial derivative in the heat equation with the fourth-order finite difference CN approximation. This leads to the formation of pentadiagonal matrices in the systems of linear equations with larger computational stencils. The algorithm also employs the higher accuracy of the Mitchell and Fairweather variant. Despite the fourth-order IADE scheme’s higher computational complexity, this technique is proved to be valuable because it enhances the accuracy of its second-order counterpart, namely, the IADEMF2. In addition, the higher accuracy of IADEMF4 is verified as a convergent and unconditionally stable scheme and is superior in terms of rate of convergence. It has also been proven to execute more efficiently in comparison to the other benchmarked fourth-order classical iterative methods, such as the GS4, the SOR4, and the JAC4. The increasing number of correct digits at each iteration serves as an advantage for the IADEMF4, thus yielding faster rate of convergence with higher level of accuracy.

In conclusion, the proposed IADEMF4 scheme affords users many advantages with respect to higher accuracy, stability, and rate of convergence, and it serves as an alternative, efficient technique for the solution of a one-dimensional heat equation with Dirichlet boundary conditions.

The proposed fourth-order scheme can be modified and adapted to more general multidimensional linear and nonlinear parabolic, elliptic, and hyperbolic partial differential equations, using different types of boundary conditions. Besides, it can also be considered for applications in problems that require higher-order accuracy with high resolution, such as problems in nanocomputing that require the solution of very large sparse systems of equation.

The explicit and high accuracy feature of the IADEMF4 can be exploited in two- and three-dimensional heat problems, by applying the scheme as pentadiagonal solvers for the system of linear equations arising in the sweeps of a higher-order alternating direction implicit (ADI) scheme. The lower-order ADI scheme was initially proposed by Peaceman and Rachford [24]. The higher space dimensions are expected not to cause unsatisfactory performance for the ADI-IADEMF4, as long as the proposed scheme is stable, convergent, and efficient as the process of its implementation continues to advance from one-time step to another. Pathirana et al. [25] implemented the ADI in developing a two-dimensional model for incorporating flood damage in urban drainage planning. The proposed model is found to be stable, numerically accurate, and computationally efficient. Mirzavand et al. [26] proposed the ADI-FDTD method for the physical modeling of high-frequency semiconductor devices. The approach is able to reduce significantly the full-wave simulation time.

It is possible to parallelize the IADEMF4 algorithm since the calculation of the new iteration only depends on the known values from the last iteration (Figures 2 and 3). Future work is to exploit the explicit computational properties and efficiency of the IADEMF4 for parallelization and execution on distributed parallel computing systems. The idea is to speed up the execution time without compromising its accuracy, especially on problems involving very large linear systems of equations.