Fourth-Order Deferred Correction Scheme for Solving Heat Conduction Problem

Abstract

A deferred correction method is utilized to increase the order of spatial accuracy of the Crank-Nicolson scheme for the numerical solution of the one-dimensional heat equation. The fourth-order methods proposed are the easier development and can be solved by using Thomas algorithms. The stability analysis and numerical experiments have been limited to one-dimensional heat-conducting problems with Dirichlet boundary conditions and initial data.

1. Introduction

The desired properties of finite difference schemes are stability, accuracy, and efficiency. These requirements are in conflict with each other. In many applications a high-order accuracy is required in the spatial discretization. To reach better stability, implicit approximation is desired. For a high-order method of traditional type (not a high-order compact (HOC)), the stencil becomes wider with increasing order of accuracy. For a standard centered discretization of order 𝑝, the stencil is 𝑝+1 points wide. This inflicts problems at the fictional boundaries, and using an implicit method results in the solution of an algebraic system of equations with large bandwidth. In light of conflict requirements of stability, accuracy, and computational efficiency, it is desired to develop schemes that have a wide range of stability and highorder of accuracy and lead to the solution of a system of linear equations with a tri-diagonal matrix, that is, the system of linear equations arising from a standard second-order discretization of heat equation.

The development of high-order compact (HOC) schemes [1–18] is one approach to overcome the antagonism between stability, accuracy, and computational cost. However, the HOC becomes complicated when applie to multidimensional problems or to non-Cartesian coordinate cases.

Another way of preserving a compact stencil at higher time level and reaching high-order spatial accuracy is the deferred correction approach [11]. A classical deferred correction procedure is developed in [19, 20].

In this paper we use the deferred correction technique to obtain fourth-order accurate schemes in space for the one-dimensional heat-conducting problem with Dirichlet boundary conditions. The linear system that needs to be solved at each time step is similar to the standard Crank-Nikolson method of second order which is solved by using Thomas algorithms. The fourth-order deferred (FOD) correction schemes are compared with the fourth-order semi-implicit (FOS) schemes and fourth-order compact (FOC) schemes for the Dirichlet boundary value problems.

A set of schemes are constructed for the one-dimensional heat-conducting problem with Dirichlet boundary conditions and initial data: 𝑢_𝑡=𝛽𝑢_𝑥𝑥()𝑢()+𝑓𝑥,𝑡,0<𝑥<𝑙,𝑡>0,𝑥,0=𝑢₀(𝑥)(),0<𝑥<𝑙,DirichletBC∶𝑢0,𝑡=𝛾₁(𝑡),𝑢(𝑙,𝑡)=𝛾₂(𝑡),𝑡>0, (1)(2)(3) where the diffusion coefficient 𝛽 is positive, 𝑢(𝑥,𝑡) represents the temperature at point (𝑥,𝑡), and 𝑓(𝑥,𝑡), 𝛾₁(𝑡), 𝛾₂(𝑡) are sufficiently smooth functions.

The rest of this paper is organized as follows. Section 2 presents an FOD scheme which we use to compare performance of proposed scheme with FOS and FOC schemes. Section 3 provides examples of comparisons. Although FOD schemes have a higher computational cost than FOS and FOC schemes, it is evident from these examples that the FOD schemes have the advantage of accuracy in the uniform norm, robustness, and the ability to be extended easily to the multidimensional case. We conclude the paper in Section 4.

2. The Fourth-Order Schemes

Let Δ𝑡 denote the temporal mesh size. For simplicity, we consider a uniform mesh consisting of 𝑁 points: 𝑥₁,𝑥₂,…,𝑥_𝑁 where 𝑥_𝑖=(𝑖−1)Δ𝑥 and the mesh size is Δ𝑥=𝑙/(𝑁−1). Below we use the notations 𝑢^𝑚_𝑖 and (𝑢_𝑥𝑥)^𝑚_𝑖 to represent the numerical approximations of 𝑢(𝑥_𝑖,𝑡^𝑚) and 𝑢_𝑥𝑥(𝑥_𝑖,𝑡^𝑚) where 𝑡^𝑚=𝑚Δ𝑡 and 𝑢^(𝑝) is the value of the 𝑝th derivative of the given function 𝑢.

2.1. Fourth-Order Semi-Implicit Scheme

The application to the well-known Crank-Nikolson scheme to (1) results in the following expression: 𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖=𝛽Δ𝑡2𝑢_𝑥𝑥𝑖^𝑛+1+𝑢_𝑥𝑥^𝑛_𝑖+𝑓_𝑖^𝑛+1/2, (4) where 𝑓_𝑖^𝑛+1/2=(𝑓_𝑖^𝑛+1+𝑓^𝑛_𝑖)/2,𝑖=2,…,𝑁−1. The Dirichlet boundary conditions 𝑢(0,𝑚Δ𝑡)=𝛾₁𝑡^𝑚=𝑢^𝑚₁,𝑢(𝑙,𝑚Δ𝑡)=𝛾₂𝑡^𝑚=𝑢^𝑚_𝑁 (5) are used to derive the following fourth-order approximation of second derivative terms: 𝑢_𝑥𝑥^𝑚₂=112Δ𝑥²𝑎₁𝑢^𝑚₁+𝑎₂𝑢^𝑚₂+𝑎₃𝑢^𝑚₃+𝑎₄𝑢^𝑚₄+𝑎₅𝑢^𝑚₅+𝑎₆𝑢^𝑚₆=𝑎₁12Δ𝑥²𝛾₁𝑡^𝑚+112Δ𝑥²×𝑎₂𝑢^𝑚₂+𝑎₃𝑢^𝑚₃+𝑎₄𝑢^𝑚₄+𝑎₅𝑢^𝑚₅+𝑎₆𝑢^𝑚₆,𝑢_𝑥𝑥^𝑚_𝑖=112Δ𝑥²×−𝑢^𝑚_𝑖−2+16𝑢^𝑚_𝑖−1−30𝑢^𝑚_𝑖+16𝑢^𝑚_𝑖+1−𝑢^𝑚_𝑖+2,𝑢𝑖=3,…,𝑁−2,_{𝑥𝑥𝑚𝑁−1}=112Δ𝑥²𝑎₁𝑢^𝑚_𝑁+𝑎₂𝑢^𝑚_𝑁−1+𝑎₃𝑢^𝑚_𝑁−2+𝑎₄𝑢^𝑚_𝑁−3+𝑎₅𝑢^𝑚_𝑁−4+𝑎₆𝑢^𝑚_𝑁−5=𝑎₁12Δ𝑥²𝛾₂𝑡^𝑚+112Δ𝑥²×𝑎₂𝑢^𝑚_𝑁−1+𝑎₃𝑢^𝑚_𝑁−2+𝑎₄𝑢^𝑚_𝑁−3+𝑎₅𝑢^𝑚_𝑁−4+𝑎₆𝑢^𝑚_𝑁−5, (6) where the coefficients can be found by matching the Taylor series expansion of left-hand-side terms up to order 𝑂(Δ𝑥⁴)𝑢⁽⁶⁾ which gives the following values of coefficients: 𝑎₁=10,𝑎₂=−15,𝑎₃𝑎=−4,₄=14,𝑎₅=−6,𝑎₆=1. (7) Schemes (6) can be combined and expressed in the following matrix form: 𝐮^𝑚_𝑥𝑥=1Δ𝑥²Λ_ℎ𝐮^𝑚+𝛾𝑡^𝑚, (8) where Λ_ℎ is the corresponding triangular and sparse (𝑁−2)×(𝑁−2) matrix, 𝐮^𝑚_𝑥𝑥=𝑢_𝑥𝑥^𝑚₂,𝑢_𝑥𝑥^𝑚₃,…,𝑢_𝑥𝑥^𝑚_𝑁−1^𝑇,𝐮^𝑚=𝑢^𝑚₂,𝑢^𝑚₃,…,𝑢^𝑚_𝑁−1^𝑇,𝛾𝑡^𝑚=𝛾₁𝑡^𝑚,0,…,0,𝛾₂𝑡^𝑚𝑇. (9) Substituting (6) into (4) gives us the following matrix form: 𝐸−𝛼Λ_ℎ𝐮^𝑛+1=𝐸+𝛼Λ_ℎ𝐮^𝑛+Δ𝑡𝛾𝑡^𝑛+1+𝛾𝑡^𝑛+Δ𝑡𝐟^𝑛+1/2, (10) where 𝛼=𝛽Δ𝑡/(2Δ𝑥²), 𝐟^𝑛+1/2=(𝑓₂^𝑛+1/2,𝑓₃^𝑛+1/2,…,𝑓^𝑛+1/2_𝑁−1)^𝑇, and 𝐸 denote the (𝑁−2)×(𝑁−2) identity matrix. The scheme (10) is FOSs for the heat-conducting problem with Dirichlet boundary condition. The order of approximation is 𝑂(Δ𝑡²,Δ𝑥⁴) in the uniform norm. The triangular and sparse (𝑁−2)×(𝑁−2) coefficient matrix in FOSs are time independent; hence, we have to store the inverse of the coefficient matrix 𝐸−𝛼Λ_ℎ before the time marching in the implementation for computational efficiency.

2.2. Fourth-Order Deferred Correction Schemes

A set of fourth-order deferred correction schemes is based on the well-known Crank-Nikolson type of scheme in the following form: 𝑢_𝑖^{𝑛+1,𝑠+1}−𝑢^𝑛_𝑖=𝛽Δ𝑡2𝑢_𝑥𝑥𝑖^{𝑛+1,𝑠+1}+𝑢_𝑥𝑥^𝑛_𝑖+𝑓_𝑖^𝑛+1/2, (11) where 𝑓_𝑖^𝑛+1/2=(𝑓_𝑖^𝑛+1+𝑓^𝑛_𝑖)/2 and the second superscript “𝑠” denotes the number of iterations 𝑠=0,…,𝑆 and 𝑖=2,…,𝑁−1.

The deferred correction technique [11] is utilized to approximate the second-order derivatives at higher time levels (𝑢_𝑥𝑥)_𝑖^{𝑛+1,𝑠+1},𝑖=2,…,𝑁−1 by the iterative method 𝑢_𝑥𝑥𝑖^{𝑛+1,𝑠+1}=𝑢^𝑙_𝑥𝑥𝑖^{𝑛+1,𝑠+1}+𝑢^ℎ_𝑥𝑥𝑖^𝑛+1,𝑠−𝑢^𝑙_𝑥𝑥𝑖^𝑛+1,𝑠, (12) where (𝑢^ℎ_𝑥𝑥)_𝑖^𝑛+1,𝑠, 𝑖=2,…,𝑁−1, 𝑠=0,…,𝑆, is high-order approximation on wide stencil and (𝑢^𝑙_𝑥𝑥)_𝑖^𝑛+1,𝑘, 𝑘=𝑠,𝑠+1, 𝑖=2,…,𝑁−1, is the lower-order approximation on compact stencil (usually three-point stencil). The expression in the square brackets of (12) is evaluated explicitly using the values known from the previous iteration. When 𝑠=0 we use the solution from the time level 𝑛 (so 𝑢^𝑛+1,0=𝑢^𝑛 and (𝑢_𝑥𝑥)_𝑖^𝑛+1,0=(𝑢_𝑥𝑥)^𝑛_𝑖). Once the iterations converge, the lower-order approximation terms drop out and the approximation of (𝑢_𝑥𝑥)_𝑖^{𝑛+1,𝑠+1} obtained has the same order of approximation as (𝑢^ℎ_𝑥𝑥)^𝑛+1,𝑆_𝑖. There are no difficulties to construct high-order approximation for interior points.

To preserve a compact three using wide stencil in the finite difference scheme at higher time level (𝑛+1,𝑠+1), we use the central second-order finite difference approximation to approximate the lower-order term in (12): 𝑢^𝑙_𝑥𝑥𝑖^𝑛+1,𝑘=1Δ𝑥²Λ_𝑙𝑢_𝑖^𝑛+1,𝑘Λ,𝑘=𝑠,𝑠+1,𝑖=3,…,𝑁−2,_𝑙𝑢_𝑖^𝑛+1,𝑘=𝑢^𝑛+1,𝑘_𝑖−1−2𝑢_𝑖^𝑛+1,𝑘+𝑢^𝑛+1,𝑘_𝑖+1. (13) For the high-order approximation term in (12), we use a symmetric five-point wide stencil for the inner points to reach the fourth order of approximation: 𝑢^ℎ_𝑥𝑥𝑖^𝑛+1,𝑠=1Δ𝑥²Λ_ℎ𝑢_𝑖^𝑛+1,𝑠Λ,𝑖=3,…,𝑁−2,_ℎ𝑢_𝑖^𝑛+1,𝑠=112−𝑢^𝑛+1,𝑠_𝑖−2+16𝑢^𝑛+1,𝑠_𝑖−1−30𝑢_𝑖^𝑛+1,𝑠+16𝑢^𝑛+1,𝑠_𝑖+1−𝑢^𝑛+1,𝑠_𝑖+2. (14) Case 𝑠=0 in (13) gives the fourth order of approximation to approximate the second-order derivatives at the time level 𝑛.

2.2.1. Stability Analysis

To study the stability of scheme (11)–(14), we use the Von-Neumann stability analysis. For simplicity, we assume that 𝑓_𝑖^𝑛+1/2≡0 in (11) and 𝑢 is periodic in 𝑥.

Let us recast scheme (11) in the following form: 𝐸+𝛼Λ_𝑙𝑢_𝑖^{𝑛+1,𝑠+1}=𝛼Λ_𝑙−Λ_ℎ𝑢_𝑖^𝑛+1,𝑠+𝐸−𝛼Λ_𝑙𝑢^𝑛_𝑖, (15) where 𝛼=𝛽Δ𝑡/(2Δ𝑥²). If we define the following operators 𝐴=𝐸+𝛼Λ_𝑙, 𝐵=𝐸−𝛼Λ_ℎ, and 𝐶=𝐸+𝛼Λ_ℎ, where 𝐸 is the identity operator, then (15) can be rewritten as follows: 𝐴𝑢_𝑖^{𝑛+1,𝑠+1}=(𝐴−𝐶)𝑢_𝑖^𝑛+1,𝑠+𝐵𝑢^𝑛_𝑖. (16) Assuming that the operators commute, (𝐴−𝐶)𝐴=𝐴(𝐴−𝐶) (e.g., in the case of uniform grid), it is easy to demonstrate that if 𝑢^𝑛+1,𝑆_𝑖+1=𝑢_𝑖^𝑛+1 and 𝑢_𝑖^𝑛+1,0=𝑢^𝑛_𝑖 we get 𝐴𝑆₊₁𝑢_𝑖^𝑛+1=_𝑘𝑆₌₀𝐴𝑆_−𝑘(𝐴−𝐶)^𝑘𝐵𝑢^𝑛_𝑖+(𝐴−𝐶)𝑆₊₁𝑢^𝑛_𝑖. (17) Let 𝑢^𝑛_𝑖=𝜉^𝑛𝑒^𝐼Θ𝑖, 𝐼=√−1, be the solution of (11)–(14), where Θ=2𝜋Δ𝑥/𝑙 is the phase angle with wavelength 𝑙. From (17), we can derive an equation for the amplification factor in the form 𝜉=𝜑Θ,𝑆,𝛼, (18) where 𝑆 is the number of iterations, and 𝜑=∑𝑆Θ,𝑆,𝛼^𝑘=0𝐴𝑆_−𝑘(𝐴−𝐶)^𝑘𝐵+(𝐴−𝐶)𝑆₊₁𝑒^𝐼Θ𝑖𝐴𝑆₊₁𝑒^𝐼Θ𝑖. (19) For stability of the method it is necessary that the absolute values of the amplification factor are less than one; that is, 𝜉<1. (20) Calculations are tedious and almost impossible to do by hand without mistake. We have therefore automate all calculations in a computer algebra environment based on REDUCE to obtain an explicit form of |𝜑(Θ,𝑆,𝛼)|. Figure 1 shows the values of |𝜉|² in the polar coordinate system (|𝜉|²,Θ) for 𝑆=1,3,and5. If only one iteration is executed in (11), 𝑆=1, inequality (20) holds if 𝛼<1.5, as can be seen from Figure 1(a)). If 3 iterations are done in (11) (Figure 1(b)), 𝑆=3, the amplification factor remains bounded by one at least for 𝛼≤10. In case of 𝑆=5, the stability criteria hold up to 𝛼=30 as can be seen from Figure 1(c)). It can be seen that increasing the number of internal iterations results in increasing the range of 𝛼 needed for stability. This tendency allows to assume that as 𝑆→∞, our method becomes the unconditionally stable Crank-Nikolson method for the heat equation.

(a)

(b)

(c)

Figure 1 

Variation of amplification factor with Θ. (a) 𝑆=1, dashed line 𝛼=2.0, solid line 𝛼=1.5, dash-doted line 𝛼=1.0, (b) 𝑆=3, dashed line 𝛼=10.0, solid line 𝛼=5.0, dash-dotted line 𝛼=1.0, and (c) 𝑆=5, dashed line 𝛼=30.0, solid line 𝛼=15, dash-doted line 𝛼=5.0, doted line 𝛼=1.0.

2.2.2. Fourth-Order Deferred Correction Scheme

Let us first consider the one-dimensional heat conduction problem with initial data and Dirichlet boundary conditions (1)–(3): 𝑢₁^𝑛+1,𝑘=𝛾₁𝑡^𝑛+1,𝑢_𝑁^𝑛+1,𝑘=𝛾₂𝑡^𝑛+1. (21) The finite difference approximations at 𝑥₂ and 𝑥_𝑁−1, which are the points next to the left and right boundaries, are straightforward: 𝑢^𝑙_𝑥𝑥2^𝑛+1,𝑘=1Δ𝑥²𝛾₁𝑡^𝑛+1−2𝑢₂^𝑛+1,𝑘+𝑢₃^𝑛+1,𝑘,𝑢𝑘=𝑠,𝑠+1,_{ℎ𝑥𝑥2𝑛+1,𝑠}=112Δ𝑥²10𝛾₁𝑡^𝑛+1−15𝑢₂^𝑛+1,𝑠−4𝑢₃^𝑛+1,𝑠+14𝑢₄^𝑛+1,𝑠−6𝑢₅^𝑛+1,𝑠+𝑢₆^𝑛+1,𝑠,𝑢^ℎ_𝑥𝑥^𝑛+1,𝑠_𝑁−1=112Δ𝑥²10𝛾₂𝑡^𝑛+1−15𝑢^𝑛+1,𝑠_𝑁−1−4𝑢^𝑛+1,𝑠_𝑁−2+14𝑢^𝑛+1,𝑠_𝑁−3−6𝑢^𝑛+1,𝑠_𝑁−4+𝑢^𝑛+1,𝑠_𝑁−5,𝑢^𝑙_𝑥𝑥^𝑛+1,𝑘_𝑁−1=1Δ𝑥²𝑢^𝑛+1,𝑘_𝑁−2−2𝑢^𝑛+1,𝑘_𝑁−1+𝛾₂𝑡^𝑛+1,𝑘=𝑠,𝑠+1. (22) Cases 𝑠=0 or 𝑘=0 give formulae to approximate (𝑢^𝑙_𝑥𝑥)^𝑛_𝑖 and (𝑢^ℎ_𝑥𝑥)^𝑛_𝑖. Substituting (13), (14), and (22) into (12) the following fourth-order deferred correction approximations of (𝑢_𝑥𝑥)_𝑖^{𝑛+1,𝑠+1},𝑖=2,…,𝑁−1, are 𝑢_𝑥𝑥2^{𝑛+1,𝑠+1}=56Δ𝑥²𝛾₁𝑡^𝑛+1+1Δ𝑥²−2𝑢₂^{𝑛+1,𝑠+1}+𝑢₃^{𝑛+1,𝑠+1}+112Δ𝑥²9𝑢₂^𝑛+1,𝑠−16𝑢₃^𝑛+1,𝑠+14𝑢₄^𝑛+1,𝑠−6𝑢₅^𝑛+1,𝑠+𝑢₆^𝑛+1,𝑠,𝑢_𝑥𝑥𝑖^{𝑛+1,𝑠+1}=1Δ𝑥²𝑢^{𝑛+1,𝑠+1}_𝑖−1−2𝑢_𝑖^{𝑛+1,𝑠+1}+𝑢^{𝑛+1,𝑠+1}_𝑖+1+112Δ𝑥²−𝑢^𝑛+1,𝑠_𝑖−2+4𝑢^𝑛+1,𝑠_𝑖−1−6𝑢_𝑖^𝑛+1,𝑠+4𝑢^𝑛+1,𝑠_𝑖+1−𝑢^𝑛+1,𝑠_𝑖+2,𝑢𝑖=3,…,𝑁−2,_{𝑥𝑥𝑛+1,𝑠+1𝑁−1}=56Δ𝑥²𝛾₂𝑡^𝑛+1+1Δ𝑥²−2𝑢^{𝑛+1,𝑠+1}_𝑁−1+𝑢_𝑁^{𝑛+1,𝑠+1}+112Δ𝑥²9𝑢^𝑛+1,𝑠_𝑁−1−16𝑢^𝑛+1,𝑠_𝑁−2+14𝑢^𝑛+1,𝑠_𝑁−3−6𝑢^𝑛+1,𝑠_𝑁−4+𝑢^𝑛+1,𝑠_𝑁−5. (23) Schemes (23) can be combined and expressed in the following matrix form: 𝐮^{𝑛+1,𝑠+1}_𝑥𝑥=1Δ𝑥²Λ_𝑙𝐮^{𝑛+1,𝑠+1}+1Δ𝑥²Λ_ℎ−Λ_𝑙𝐮^𝑛+1,𝑠+𝛾𝑡^𝑛+1, (24) where Λ_𝑙 is a tridiagonal (𝑁−2)×(𝑁−2) matrix and Λ_ℎ is the corresponding triangular and sparse (𝑁−2)×(𝑁−2) matrix, 𝐮^𝑛+1,𝑘=𝑢₂^𝑛+1,𝑘,𝑢₃^𝑛+1,𝑘,…,𝑢^𝑛+1,𝑘_𝑁−1^𝑇𝐮,𝑘=𝑠,𝑠+1,_{𝑛+1,𝑠+1𝑥𝑥}=𝑢_𝑥𝑥2^{𝑛+1,𝑠+1},𝑢_𝑥𝑥3^{𝑛+1,𝑠+1},…,𝑢_𝑥𝑥^{𝑛+1,𝑠+1}_𝑁−1^𝑇. (25)(26) Substituting (6), (23) into (11), the formulae can be written into matrix form 𝐸−𝛼Λ_𝑙𝐮^{𝑛+1,𝑠+1}=𝛼Λ_ℎ−Λ_𝑙𝐮^𝑛+1,𝑠+𝐸+𝛼Λ_ℎ𝐮^𝑛+Δ𝑡𝛾𝑡^𝑛+1+𝛾𝑡^𝑛+Δ𝑡𝐟^𝑛+1/2. (27) The above matrix form is called FODs for Dirichlet boundary value problem (1)–(3). Thomas algorithms can be used to compute the solutions of FODs. At each step of time 𝑡^𝑛 and the initial stage, the convergence of FODs requires more iterations to converge to the solution of the FOSs. The order of approximation of FODs is 𝑂(Δ𝑡²,Δ𝑥⁴) which is the same as FOSs in the uniform norm.

2.3. Fourth-Order Compact Scheme

Let us briefly represent the main idea and final formulae of compact schemes. Spatial derivatives in the Crank-Nikolson scheme (4) are evaluated by the fourth-order compact finite differences implicit scheme [5, 7, 8, 13, 14, 17].

In [8, 14], the Dirichlet boundary conditions 𝑢(0,𝑚Δ𝑡)=𝛾₁𝑡^𝑚=𝑢^𝑚₁,𝑢(𝑙,𝑚Δ𝑡)=𝛾₂𝑡^𝑚=𝑢^𝑚_𝑁 (28) are used to derive the following fourth-order schemes 𝑢_𝑥𝑥^𝑚₂+𝜎𝑢_𝑥𝑥^𝑚₃=124Δ𝑥²𝑎₁𝑢^𝑚₁+𝑎₂𝑢^𝑚₂+𝑎₃𝑢^𝑚₃+𝑎₄𝑢^𝑚₄+𝑎₅𝑢^𝑚₅+𝑎₆𝑢^𝑚₆=𝑎₁24Δ𝑥²𝛾₁𝑡^𝑚+124Δ𝑥²𝑎₂𝑢^𝑚₂+𝑎₃𝑢^𝑚₃+𝑎₄𝑢^𝑚₄+𝑎₅𝑢^𝑚₅+𝑎₆𝑢^𝑚₆,𝑢_𝑥𝑥^𝑚_𝑖−1+10𝑢_𝑥𝑥^𝑚_𝑖+𝑢_𝑥𝑥^𝑚_𝑖+1=2Δ𝑥²6𝑢^𝑚_𝑖−1−12𝑢^𝑚_𝑖+6𝑢^𝑚_𝑖+1,𝑢𝑖=2,…,𝑁−1,_{𝑥𝑥𝑚𝑁−1}+𝜎𝑢_𝑥𝑥^𝑚_𝑁−2=124Δ𝑥²𝑎₁𝑢^𝑚_𝑁+𝑎₂𝑢^𝑚_𝑁−1+𝑎₃𝑢^𝑚_𝑁−2+𝑎₄𝑢^𝑚_𝑁−3+𝑎₅𝑢^𝑚_𝑁−4+𝑎₆𝑢^𝑚_𝑁−5=𝑎₁24Δ𝑥²𝛾₂𝑡^𝑚+124Δ𝑥²𝑎₂𝑢^𝑚_𝑁−1+𝑎₃𝑢^𝑚_𝑁−2+𝑎₄𝑢^𝑚_𝑁−3+𝑎₅𝑢^𝑚_𝑁−4+𝑎₆𝑢^𝑚_𝑁−5, (29) where the coefficients can be found by matching the Taylor series expansion of left-hand-side terms up to order 𝑂(Δ𝑥⁴)𝑢⁽⁶⁾ which gives the following values of coefficients [8]: 1𝜎=2,𝑎₁=19,𝑎₂=−14,𝑎₃𝑎=−38,₄=44,𝑎₅=−13,𝑎₆=2. (30) Then all derivatives in (4) are approximated by the fourth-order compact formula; we can write 𝐴𝐮^𝑚_𝑥𝑥=1Δ𝑥²𝐵𝐮^𝑚+𝛾^𝑚,𝑚=𝑛,𝑛+1, (31) where 𝐴 and 𝐵 are the corresponding triangular and sparse (𝑁−2)×(𝑁−2) matrices, 𝐮^𝑚_𝑥𝑥=((𝑢_𝑥𝑥)^𝑚₂,(𝑢_𝑥𝑥)^𝑚₃,…,(𝑢_𝑥𝑥)^𝑚_𝑁−1)^𝑇, 𝐮^𝑚=(𝑢^𝑚₂,𝑢^𝑚₃,…,𝑢^𝑚_𝑁−1)^𝑇 and 𝛾^𝑚 = (𝛾₁(𝑡^𝑚),0,…,0,𝛾₂(𝑡^𝑚))^𝑇, 𝑚 = 𝑛,𝑛+1. Schemes (4) and (29) can be combined and expressed in the following matrix form: (𝐴−𝛼𝐵)𝐮^𝑛+1=(𝐴+𝛼𝐵)𝐮^𝑛+Δ𝑡𝛾𝑡^𝑛+1+𝛾𝑡^𝑛+Δ𝑡𝐟^𝑛+1/2. (32) This scheme is called FOCs for Dirichlet boundary value problem (1)–(3). We like to mention that the above scheme has truncation error 𝑂(Δ𝑡²,Δ𝑥⁴). Note that the triangular and sparse (𝑁−2)×(𝑁−2) coefficient matrices in FOCs are time independent; hence, we have to store the inverse of the coefficient matrix 𝐴−𝛼𝐵 before the time marching in the implementation of computational efficiency.

3. Numerical Examples

In this section, three numerical examples are carried out. The first two are linear heat-conducting problem, with Dirichlet boundary conditions, which are used to confirm our theoretical analysis. Then we apply the FODS to the Burgers equation. For simplicity, we fix our problem domain Ω={𝑥∣0≤𝑥≤1}. In all computations, we used Δ𝑡=Δ𝑥²/4 and 𝜖=10⁻¹⁰. The following stopping criterion is used: max_{1≤𝑖≤𝑁}𝑢^𝑛+1,𝑆_𝑖+1−𝑢^𝑛+1,𝑆_𝑖<𝜖,𝑠=0,…,𝑆, (33) where “𝑆” denotes the number of the last iteration.

The computations are performed using uniform grids of 11, 21, 41, 81, and 161 nodes. The initial and boundary conditions are obtained based on the exact solutions. For the testing purpose only, all computations are performed for 0≤t≤1.

Example 1 (the homogeneous heat equation with the homogeneous Dirichlet boundary conditions). One has 𝑢_𝑡=𝑢_𝑥𝑥𝑢()()()(),0≤𝑥≤1,𝑡>0,𝑥,0=sin𝜋𝑥,𝑢0,𝑡=0,𝑢1,𝑡=0. (34) The exact solution is 𝑢(𝑥,𝑡)=𝑒^−𝜋2𝑡sin(𝜋𝑥). The results of performance over the time interval 𝑡∈[0,1] for the FOCs, FODs, and FOSs are represented in Table 1, where the maximum error and the rate of convergence at time instant 𝑡=1 are shown.

Table 1 

Maximum absolute error, order of convergence, and CPU time in seconds of the FOCs, FODs, and FOSs for test problem (34) at time instant 𝑡=1.

Example 2 (the nonhomogeneous heat equation with non-homogeneous Dirichlet boundary conditions). One has 𝑢_𝑡=𝑢_𝑥𝑥+𝜋²−1𝑒^−𝑡()𝑢()()cos𝜋𝑥+4𝑥−2,0≤𝑥≤1,𝑡>0,𝑥,0=cos𝜋𝑥+𝑥²,𝑢(0,𝑡)=𝑒^−𝑡,𝑢(1,𝑡)=−𝑒^−𝑡+4𝑡+1. (35) The exact solution is 𝑢(𝑥,𝑡)=𝑒^−𝑡cos(𝜋𝑥)+𝑥²+4𝑥𝑡. The results of performance over the time domain 𝑡∈[0,1] for the FOC, FOD, and FOS schemes are represented in Table 2, where the maximum error and the rate of convergence at time instant 𝑡=1 are shown.

Table 2 

Absolute error, the rate of convergence, and CPU time in seconds of the FOCs, FODs, and FOSs for the test problem (35) at time instant 𝑡=1.

The last two columns of Tables 1 and 2 demonstrate the average number of iterations in FODs at one time step and the CPU time required to obtain the solution at time instant 𝑡=1. The average number of iterations means the total number of iterations divided by the number of time steps. As a rule, at the initial stage the convergence of deferred correction requires more iterations. For larger instants of time, the convergence occurs after 2~7 iterations as can be seen from Tables 1 and 2. All of schemes are seen to be the fourth order of accuracy, as the error is reduced approximately by factor four when the mesh is refined by half. The maximum error of the FODs and FOCs is almost the same, since the iterative scheme FODs is constructed by applying the deferred correction technique on the FOSs. It can be stated that when the iterations converge, the solution of FODs, therefore, converges to the solution of FOSs in each step of time. As shown in Tables 1 and 2, there is hardly a difference in the computational efficiency between FODs and FOSs. Both schemes are more efficient than FODs. An explanation is due to the iteration needed for the convergence of solutions on each step of time.

Although the FODs use more computational time as compared with FOCs and FOSs, it is recommended that the construction of FODs can be easily implemented. Moreover, the scheme does not need to store the inverse of coefficient matrices as required in FOCs and FOSs. Therefore, the method is easily extended to multidimensional cases.

It is suggested that the differed correction technique can solve problems which need high accuracy of computational methods. Also this technique can be easily implemented and extended for solving problem with Neumann boundary conditions. In addition, such technique can be easily used to create standard code and applied in case of nonuniform grids.

Considering Burgers equation 𝑢_𝑡=𝛽𝑢_𝑥𝑥−𝑢𝑢_𝑥,0≤𝑥≤1,𝑡>0, (36) with the exact solution [21] is given by 𝑢()=𝑒𝑥,𝑡𝜉+𝜂+𝜂−𝜉^𝜌1+𝑒^𝜌, (37) where 𝜌=𝜉(𝑥−𝜂𝑡−𝜐)/𝛽. The initial and Dirichlet boundary conditions are considered to be in agreement with the exact solution proposed here. For Burgers equation (36), we solve it by the following fourth-order deferred correction scheme: 𝑢_𝑖^{𝑛+1,𝑠+1}−𝑢^𝑛_𝑖=𝛽Δ𝑡2𝑢_𝑥𝑥𝑖^{𝑛+1,𝑠+1}+𝑢_𝑥𝑥^𝑛_𝑖+𝑓^𝑛_𝑖, (38) where 𝑓^𝑛_𝑖=−[(𝑢²/2)_𝑥]^𝑛_𝑖. The nonlinear term 𝑓^𝑛_𝑖 is approximated with the fourth-order approximation and all the second-derivative terms in (38) are approximated by the fourth-order formula (6) and the fourth-order deferred correction schemes (23). The scheme (38) can be combined and expressed in the following matrix form: 𝐸−𝛼Λ_𝑙𝐮^{𝑛+1,𝑠+1}=𝛼Λ_ℎ−Λ_𝑙𝐮^𝑛+1,𝑠+𝐸+𝛼Λ_ℎ𝐮^𝑛+Δ𝑡𝛾𝑡^𝑛+1+𝛾𝑡^𝑛+Δ𝑡𝐟^𝑛, (39) where 𝐸 is identity matrix, Λ_𝑙 is tridiagonal (𝑁−2)×(𝑁−2) matrix, and Λ_ℎ is the corresponding triangular and sparse (𝑁−2)×(𝑁−2) matrix and can be solved by using Thomas algorithm.

Example 3 (the Burgers equation (36) and the constant values 𝜈=0.125, 𝜉=0.6, 𝜂=0.4, and 𝛽=1 with appropriate initial and Dirichlet boundary condition in agreement with exact solution (37)). This problem was solved using different time step and mesh sizes over the time interval 0<𝑡≤1. The results of performance over the time interval 𝑡∈[0,1] for the FODs are represented in Tables 3 and 4, where the maximum error and the rate of convergence at time instant 𝑡=1 are shown.

Table 3 

Maximum absolute error, order of convergence, and CPU time in seconds for Example 3 at time instant 𝑡=1 with fixed mesh size Δ𝑥=0.05.

Table 4 

Maximum absolute error, order of convergence, and CPU time in seconds for Example 3 at time instant 𝑡=1 with time step size Δ𝑡=Δ𝑥⁴.

In order to analyze the results found in application to the Burgers equation (36), Table 3 demonstrates rate of convergence, average number of iteration at each time step, and CPU time required to obtain the solution of Example 3 by using FODs at time instant 𝑡=1 when Δ𝑥=0.05 with various time step sizes. Table 4 shows the rate of convergence, average number of iteration at each time step, and CPU time required to obtain the solution of Example 3 at time instant 𝑡=1 and using uniform grids of 11, 21, and 41 with time step sizes Δ𝑡=Δ𝑥⁴ and 𝜖=10⁻¹⁰.

It can be seen from Tables 3 and 4 that numerical results are in good agreement with the exact solution. We only observe 𝑂(Δ𝑡) convergence rate and the error is dominated by time error. An explanation for this phenomenon is due to the nonlinear term, which is approximated at time level 𝑛, instead of at time level 𝑛+1/2 for the FODs (38).

4. Conclusion

In this paper, a new set of fourth-order schemes for the one-dimensional heat conduction problem with Dirichlet boundary conditions is constructed using a deferred correction technique. The construction of high-order deferred correction schemes requires only a regular three-point stencil at higher time level which is similar to the standard second-order Crank-Nikolson method. The greatest significance of FODs, compared with FOCs and FOSs, is the easier development and that it can be solved by using Thomas algorithms. Numerical examples confirm the order of accuracy. We also implement our algorithms to nonlinear problems. However, theoretical analysis for nonlinear problems needs further investigation. Posterior idea for this project is to use another way to make 𝑢𝑢_𝑥 term as follows [21, 22]: 𝑢^𝑛+1𝑢_𝑥𝑖^𝑛+1≈𝑢^𝑛+1𝑢_𝑥^𝑛_𝑖+𝑢^𝑛𝑢_𝑥𝑖^𝑛+1−𝑢^𝑛𝑢_𝑥^𝑛_𝑖, (40) where better results are expected to be found. The first two terms on the right-hand side of above equation make the coefficient matrices of FOCs, FODs, and FOSs vary with time. That is, the inverse coefficient matrices of FOCs and FOSs have to be stored on each step of time while FODs have no need. For this reason, the FODs is simple to implement although FODs need more iterations for the convergence of solution on each step of time.

Acknowledgments

This work is financially supported by the Commission on Higher Education (CHE), the Thailand Research Fund (TRF), the University of Phayao (UP), Project MRG5580014. The author would like to express their deep appreciation to Professor Sergey Meleshko for the kind assistance and valuable advice on the REDUCE calculations.