# Fourth-Order Deferred Correction Scheme for Solving Heat Conduction Problem

Authors: D. Yambangwai; N. P. Moshkin
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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 orderp, the stencil is p+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:(1)ut=βuxx+f(x,t),0<x<l,t>0,(2)u(x,0)=u0(x),0<x<l,(3)DirichletBC:u(0,t)=γ1(t),u(l,t)=γ2(t),t>0,
where the diffusion coefficient β is positive, u(x,t) represents the temperature at point (x,t), and f(x,t), γ1(t), γ2(t) are sufficiently smooth functions.The rest of this paper is organized as follows. Section2 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Δt denote the temporal mesh size. For simplicity, we consider a uniform mesh consisting of N points: x1,x2,…,xN where xi=(i-1)Δx and the mesh size is Δx=l/(N-1). Below we use the notations uim and (uxx)im to represent the numerical approximations of u(xi,tm) and uxx(xi,tm) where tm=mΔt and u(p) is the value of the pth derivative of the given function u.


### 2.1. Fourth-Order Semi-Implicit Scheme

The application to the well-known Crank-Nikolson scheme to (1) results in the following expression:
(4)uin+1-uinΔt=β2[(uxx)in+1+(uxx)in]+fin+1/2,
where fin+1/2=(fin+1+fin)/2,i=2,…,N-1. The Dirichlet boundary conditions
(5)u(0,mΔt)=γ1(tm)=u1m,u(l,mΔt)=γ2(tm)=uNm
are used to derive the following fourth-order approximation of second derivative terms:
(6)(uxx)2m=112Δx2(a1u1m+a2u2m+a3u3m+a4u4m+a5u5m+a6u6m)=a112Δx2γ1(tm)+112Δx2×(a2u2m+a3u3m+a4u4m+a5u5m+a6u6m),(uxx)im=112Δx2×(-ui-2m+16ui-1m-30uim+16ui+1m-ui+2m),i=3,…,N-2,(uxx)N-1m=112Δx2(a1uNm+a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m)=a112Δx2γ2(tm)+112Δx2×(a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m),
where the coefficients can be found by matching the Taylor series expansion of left-hand-side terms up to order O(Δx4)u(6) which gives the following values of coefficients:
(7)a1=10,a2=-15,a3=-4,a4=14,a5=-6,a6=1.
Schemes (6) can be combined and expressed in the following matrix form:
(8)uxxm=1Δx2Λhum+γ(tm),
where Λh is the corresponding triangular and sparse (N-2)×(N-2) matrix,
(9)uxxm=((uxx)2m,(uxx)3m,…,(uxx)N-1m)T,um=(u2m,u3m,…,uN-1m)T,γ(tm)=(γ1(tm),0,…,0,γ2(tm))T.
Substituting (6) into (4) gives us the following matrix form:
(10)(E-αΛh)un+1=(E+αΛh)un+Δt[γ(tn+1)+γ(tn)]+Δtfn+1/2,
where α=βΔt/(2Δx2), fn+1/2=(f2n+1/2,f3n+1/2,…,fN-1n+1/2)T, and E denote the (N-2)×(N-2) identity matrix. The scheme (10) is FOSs for the heat-conducting problem with Dirichlet boundary condition. The order of approximation is O(Δt2,Δx4) in the uniform norm. The triangular and sparse (N-2)×(N-2) coefficient matrix in FOSs are time independent; hence, we have to store the inverse of the coefficient matrix E-αΛh 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:(11)uin+1,s+1-uinΔt=β2[(uxx)in+1,s+1+(uxx)in]+fin+1/2,
where fin+1/2=(fin+1+fin)/2 and the second superscript “s” denotes the number of iterations s=0,…,S^ and i=2,…,N-1.The deferred correction technique [11] is utilized to approximate the second-order derivatives at higher time levels (uxx)in+1,s+1,i=2,…,N-1 by the iterative method
(12)(uxx)in+1,s+1=(uxxl)in+1,s+1+[(uxxh)in+1,s-(uxxl)in+1,s],
where (uxxh)in+1,s, i=2,…,N-1, s=0,…,S^, is high-order approximation on wide stencil  and (uxxl)in+1,k, k=s,s+1, i=2,…,N-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 s=0 we use the solution from the time level n (so un+1,0=un and (uxx)in+1,0=(uxx)in). Once the iterations converge, the lower-order approximation terms drop out and the approximation of (uxx)in+1,s+1 obtained has the same order of approximation as (uxxh)in+1,S^. 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(n+1,s+1), we use the central second-order finite difference approximation to approximate the lower-order term in (12):
(13)(uxxl)in+1,k=1Δx2Λluin+1,k,k=s,s+1,i=3,…,N-2,Λluin+1,k=ui-1n+1,k-2uin+1,k+ui+1n+1,k.
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:
(14)(uxxh)in+1,s=1Δx2Λhuin+1,s,i=3,…,N-2,Λhuin+1,s=112(-ui-2n+1,s+16ui-1n+1,s-30uin+1,s+16ui+1n+1,s-ui+2n+1,s).
Case s=0 in (13) gives the fourth order of approximation to approximate the second-order derivatives at the time level n.


#### 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 fin+1/2≡0 in (11) and u is periodic in x.Let us recast scheme (11) in the following form:
(15)(E+αΛl)uin+1,s+1=α(Λl-Λh)uin+1,s+(E-αΛl)uin,
where α=βΔt/(2Δx2). If we define the following operators A=E+αΛl, B=E-αΛh, and C=E+αΛh, where E is the identity operator, then (15) can be rewritten as follows:
(16)Auin+1,s+1=(A-C)uin+1,s+Buin.
Assuming that the operators commute, (A-C)A=A(A-C) (e.g.,  in the case of uniform grid), it is easy to demonstrate that if uin+1,S^+1=uin+1 and uin+1,0=uin we get
(17)AS^+1uin+1=(∑k=0S^‍AS^-k(A-C)k)Buin+(A-C)S^+1uin.
Let uin=ξneIΘi, I=-1, be the solution of (11)–(14), where Θ=2πΔx/l is the phase angle with wavelength l. From (17), we can derive an equation for the amplification factor in the form
(18)|ξ|=|φ(Θ,S^,α)|,
where S^ is the number of iterations, and
(19)|φ(Θ,S^,α)|=|[(∑k=0S^AS^-k(A-C)k)B+(A-C)S^+1]eIΘi||AS^+1eIΘi|.
For stability of the method it is necessary that the absolute values of the amplification factor are less than one; that is,
(20)|ξ|<1.
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 |φ(Θ,S^,α)|. Figure 1 shows the values of |ξ|2 in the polar coordinate system (|ξ|2,Θ) for S^=1,3,and5. If only one iteration is executed in (11), S^=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)), S^=3, the amplification factor remains bounded by one at least for α≤10. In case of S^=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 S^→∞, our method becomes the unconditionally stable Crank-Nikolson method for the heat equation.Variation of amplification factor withΘ. (a) S^=1, dashed line α=2.0, solid line α=1.5, dash-doted line α=1.0, (b) S^=3, dashed line α=10.0, solid line α=5.0, dash-dotted line α=1.0, and (c) S^=5, dashed line α=30.0, solid line α=15, dash-doted line α=5.0, doted line α=1.0.
(a)
(b)
(c)


#### 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):
(21)u1n+1,k=γ1(tn+1),uNn+1,k=γ2(tn+1).
The finite difference approximations at x2 and xN-1, which are the points next to the left and right boundaries, are straightforward:
(22)(uxxl)2n+1,k=1Δx2(γ1(tn+1)-2u2n+1,k+u3n+1,k),k=s,s+1,(uxxh)2n+1,s=112Δx2(10γ1(tn+1)-15u2n+1,s-4u3n+1,s+14u4n+1,s-6u5n+1,s+u6n+1,s),(uxxh)N-1n+1,s=112Δx2(10γ2(tn+1)-15uN-1n+1,s-4uN-2n+1,s+14uN-3n+1,s-6uN-4n+1,s+uN-5n+1,s),(uxxl)N-1n+1,k=1Δx2(uN-2n+1,k-2uN-1n+1,k+γ2(tn+1)),k=s,s+1.
Cases s=0 or k=0 give formulae to approximate (uxxl)in and (uxxh)in. Substituting (13), (14), and (22) into (12) the following fourth-order deferred correction approximations of (uxx)in+1,s+1,i=2,…,N-1, are
(23)(uxx)2n+1,s+1=56Δx2γ1(tn+1)+1Δx2(-2u2n+1,s+1+u3n+1,s+1)+112Δx2(9u2n+1,s-16u3n+1,s+14u4n+1,s-6u5n+1,s+u6n+1,s),(uxx)in+1,s+1=1Δx2(ui-1n+1,s+1-2uin+1,s+1+ui+1n+1,s+1)+112Δx2(-ui-2n+1,s+4ui-1n+1,s-6uin+1,s+4ui+1n+1,s-ui+2n+1,s),i=3,…,N-2,(uxx)N-1n+1,s+1=56Δx2γ2(tn+1)+1Δx2(-2uN-1n+1,s+1+uNn+1,s+1)+112Δx2(9uN-1n+1,s-16uN-2n+1,s+14uN-3n+1,s-6uN-4n+1,s+uN-5n+1,s).
Schemes (23) can be combined and expressed in the following matrix form:
(24)uxxn+1,s+1=1Δx2Λlun+1,s+1+1Δx2(Λh-Λl)un+1,s+γ(tn+1),
where Λl is a tridiagonal (N-2)×(N-2) matrix and Λh is the corresponding triangular and sparse (N-2)×(N-2) matrix,
(25)un+1,k=(u2n+1,k,u3n+1,k,…,uN-1n+1,k)T,k=s,s+1,(26)uxxn+1,s+1=((uxx)2n+1,s+1,(uxx)3n+1,s+1,…,(uxx)N-1n+1,s+1)T.
Substituting (6), (23) into (11), the formulae can be written into matrix form
(27)(E-αΛl)un+1,s+1=α(Λh-Λl)un+1,s+(E+αΛh)un+Δt[γ(tn+1)+γ(tn)]+Δtfn+1/2.
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 tn 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 O(Δt2,Δx4) 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
(28)u(0,mΔt)=γ1(tm)=u1m,u(l,mΔt)=γ2(tm)=uNm
are used to derive the following fourth-order schemes
(29)(uxx)2m+σ(uxx)3m=124Δx2(a1u1m+a2u2m+a3u3m+a4u4m+a5u5m+a6u6m)=a124Δx2γ1(tm)+124Δx2(a2u2m+a3u3m+a4u4m+a5u5m+a6u6m),(uxx)i-1m+10(uxx)im+(uxx)i+1m=2Δx2(6ui-1m-12uim+6ui+1m),i=2,…,N-1,(uxx)N-1m+σ(uxx)N-2m=124Δx2(a1uNm+a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m)=a124Δx2γ2(tm)+124Δx2(a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m),
where the coefficients can be found by matching the Taylor series expansion of left-hand-side terms up to order O(Δx4)u(6) which gives the following values of coefficients [8]:
(30)σ=12,a1=19,a2=-14,a3=-38,a4=44,a5=-13,a6=2.
Then all derivatives in (4) are approximated by the fourth-order compact formula; we can write
(31)Auxxm=1Δx2Bum+γm,m=n,n+1,
where A and B are the corresponding triangular and sparse (N-2)×(N-2) matrices, uxxm=((uxx)2m,(uxx)3m,…,(uxx)N-1m)T, um=(u2m,u3m,…,uN-1m)T and γm = (γ1(tm),0,…,0,γ2(tm))T, m = n,n+1. Schemes (4) and (29) can be combined and expressed in the following matrix form:
(32)(A-αB)un+1=(A+αB)un+Δt[γ(tn+1)+γ(tn)]+Δtfn+1/2.
This scheme is called FOCs for Dirichlet boundary value problem (1)–(3). We like to mention that the above scheme has truncation error O(Δt2,Δx4). Note that the triangular and sparse (N-2)×(N-2) coefficient matrices in FOCs are time independent; hence, we have to store the inverse of the coefficient matrix A-αB before the time marching in the implementation of computational efficiency.


## 2.1. Fourth-Order Semi-Implicit Scheme

The application to the well-known Crank-Nikolson scheme to (1) results in the following expression:
(4)uin+1-uinΔt=β2[(uxx)in+1+(uxx)in]+fin+1/2,
where fin+1/2=(fin+1+fin)/2,i=2,…,N-1. The Dirichlet boundary conditions
(5)u(0,mΔt)=γ1(tm)=u1m,u(l,mΔt)=γ2(tm)=uNm
are used to derive the following fourth-order approximation of second derivative terms:
(6)(uxx)2m=112Δx2(a1u1m+a2u2m+a3u3m+a4u4m+a5u5m+a6u6m)=a112Δx2γ1(tm)+112Δx2×(a2u2m+a3u3m+a4u4m+a5u5m+a6u6m),(uxx)im=112Δx2×(-ui-2m+16ui-1m-30uim+16ui+1m-ui+2m),i=3,…,N-2,(uxx)N-1m=112Δx2(a1uNm+a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m)=a112Δx2γ2(tm)+112Δx2×(a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m),
where the coefficients can be found by matching the Taylor series expansion of left-hand-side terms up to order O(Δx4)u(6) which gives the following values of coefficients:
(7)a1=10,a2=-15,a3=-4,a4=14,a5=-6,a6=1.
Schemes (6) can be combined and expressed in the following matrix form:
(8)uxxm=1Δx2Λhum+γ(tm),
where Λh is the corresponding triangular and sparse (N-2)×(N-2) matrix,
(9)uxxm=((uxx)2m,(uxx)3m,…,(uxx)N-1m)T,um=(u2m,u3m,…,uN-1m)T,γ(tm)=(γ1(tm),0,…,0,γ2(tm))T.
Substituting (6) into (4) gives us the following matrix form:
(10)(E-αΛh)un+1=(E+αΛh)un+Δt[γ(tn+1)+γ(tn)]+Δtfn+1/2,
where α=βΔt/(2Δx2), fn+1/2=(f2n+1/2,f3n+1/2,…,fN-1n+1/2)T, and E denote the (N-2)×(N-2) identity matrix. The scheme (10) is FOSs for the heat-conducting problem with Dirichlet boundary condition. The order of approximation is O(Δt2,Δx4) in the uniform norm. The triangular and sparse (N-2)×(N-2) coefficient matrix in FOSs are time independent; hence, we have to store the inverse of the coefficient matrix E-αΛh 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:(11)uin+1,s+1-uinΔt=β2[(uxx)in+1,s+1+(uxx)in]+fin+1/2,
where fin+1/2=(fin+1+fin)/2 and the second superscript “s” denotes the number of iterations s=0,…,S^ and i=2,…,N-1.The deferred correction technique [11] is utilized to approximate the second-order derivatives at higher time levels (uxx)in+1,s+1,i=2,…,N-1 by the iterative method
(12)(uxx)in+1,s+1=(uxxl)in+1,s+1+[(uxxh)in+1,s-(uxxl)in+1,s],
where (uxxh)in+1,s, i=2,…,N-1, s=0,…,S^, is high-order approximation on wide stencil  and (uxxl)in+1,k, k=s,s+1, i=2,…,N-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 s=0 we use the solution from the time level n (so un+1,0=un and (uxx)in+1,0=(uxx)in). Once the iterations converge, the lower-order approximation terms drop out and the approximation of (uxx)in+1,s+1 obtained has the same order of approximation as (uxxh)in+1,S^. 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(n+1,s+1), we use the central second-order finite difference approximation to approximate the lower-order term in (12):
(13)(uxxl)in+1,k=1Δx2Λluin+1,k,k=s,s+1,i=3,…,N-2,Λluin+1,k=ui-1n+1,k-2uin+1,k+ui+1n+1,k.
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:
(14)(uxxh)in+1,s=1Δx2Λhuin+1,s,i=3,…,N-2,Λhuin+1,s=112(-ui-2n+1,s+16ui-1n+1,s-30uin+1,s+16ui+1n+1,s-ui+2n+1,s).
Case s=0 in (13) gives the fourth order of approximation to approximate the second-order derivatives at the time level n.


### 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 fin+1/2≡0 in (11) and u is periodic in x.Let us recast scheme (11) in the following form:
(15)(E+αΛl)uin+1,s+1=α(Λl-Λh)uin+1,s+(E-αΛl)uin,
where α=βΔt/(2Δx2). If we define the following operators A=E+αΛl, B=E-αΛh, and C=E+αΛh, where E is the identity operator, then (15) can be rewritten as follows:
(16)Auin+1,s+1=(A-C)uin+1,s+Buin.
Assuming that the operators commute, (A-C)A=A(A-C) (e.g.,  in the case of uniform grid), it is easy to demonstrate that if uin+1,S^+1=uin+1 and uin+1,0=uin we get
(17)AS^+1uin+1=(∑k=0S^‍AS^-k(A-C)k)Buin+(A-C)S^+1uin.
Let uin=ξneIΘi, I=-1, be the solution of (11)–(14), where Θ=2πΔx/l is the phase angle with wavelength l. From (17), we can derive an equation for the amplification factor in the form
(18)|ξ|=|φ(Θ,S^,α)|,
where S^ is the number of iterations, and
(19)|φ(Θ,S^,α)|=|[(∑k=0S^AS^-k(A-C)k)B+(A-C)S^+1]eIΘi||AS^+1eIΘi|.
For stability of the method it is necessary that the absolute values of the amplification factor are less than one; that is,
(20)|ξ|<1.
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 |φ(Θ,S^,α)|. Figure 1 shows the values of |ξ|2 in the polar coordinate system (|ξ|2,Θ) for S^=1,3,and5. If only one iteration is executed in (11), S^=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)), S^=3, the amplification factor remains bounded by one at least for α≤10. In case of S^=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 S^→∞, our method becomes the unconditionally stable Crank-Nikolson method for the heat equation.Variation of amplification factor withΘ. (a) S^=1, dashed line α=2.0, solid line α=1.5, dash-doted line α=1.0, (b) S^=3, dashed line α=10.0, solid line α=5.0, dash-dotted line α=1.0, and (c) S^=5, dashed line α=30.0, solid line α=15, dash-doted line α=5.0, doted line α=1.0.
(a)
(b)
(c)


### 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):
(21)u1n+1,k=γ1(tn+1),uNn+1,k=γ2(tn+1).
The finite difference approximations at x2 and xN-1, which are the points next to the left and right boundaries, are straightforward:
(22)(uxxl)2n+1,k=1Δx2(γ1(tn+1)-2u2n+1,k+u3n+1,k),k=s,s+1,(uxxh)2n+1,s=112Δx2(10γ1(tn+1)-15u2n+1,s-4u3n+1,s+14u4n+1,s-6u5n+1,s+u6n+1,s),(uxxh)N-1n+1,s=112Δx2(10γ2(tn+1)-15uN-1n+1,s-4uN-2n+1,s+14uN-3n+1,s-6uN-4n+1,s+uN-5n+1,s),(uxxl)N-1n+1,k=1Δx2(uN-2n+1,k-2uN-1n+1,k+γ2(tn+1)),k=s,s+1.
Cases s=0 or k=0 give formulae to approximate (uxxl)in and (uxxh)in. Substituting (13), (14), and (22) into (12) the following fourth-order deferred correction approximations of (uxx)in+1,s+1,i=2,…,N-1, are
(23)(uxx)2n+1,s+1=56Δx2γ1(tn+1)+1Δx2(-2u2n+1,s+1+u3n+1,s+1)+112Δx2(9u2n+1,s-16u3n+1,s+14u4n+1,s-6u5n+1,s+u6n+1,s),(uxx)in+1,s+1=1Δx2(ui-1n+1,s+1-2uin+1,s+1+ui+1n+1,s+1)+112Δx2(-ui-2n+1,s+4ui-1n+1,s-6uin+1,s+4ui+1n+1,s-ui+2n+1,s),i=3,…,N-2,(uxx)N-1n+1,s+1=56Δx2γ2(tn+1)+1Δx2(-2uN-1n+1,s+1+uNn+1,s+1)+112Δx2(9uN-1n+1,s-16uN-2n+1,s+14uN-3n+1,s-6uN-4n+1,s+uN-5n+1,s).
Schemes (23) can be combined and expressed in the following matrix form:
(24)uxxn+1,s+1=1Δx2Λlun+1,s+1+1Δx2(Λh-Λl)un+1,s+γ(tn+1),
where Λl is a tridiagonal (N-2)×(N-2) matrix and Λh is the corresponding triangular and sparse (N-2)×(N-2) matrix,
(25)un+1,k=(u2n+1,k,u3n+1,k,…,uN-1n+1,k)T,k=s,s+1,(26)uxxn+1,s+1=((uxx)2n+1,s+1,(uxx)3n+1,s+1,…,(uxx)N-1n+1,s+1)T.
Substituting (6), (23) into (11), the formulae can be written into matrix form
(27)(E-αΛl)un+1,s+1=α(Λh-Λl)un+1,s+(E+αΛh)un+Δt[γ(tn+1)+γ(tn)]+Δtfn+1/2.
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 tn 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 O(Δt2,Δx4) which is the same as FOSs in the uniform norm.


## 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 fin+1/2≡0 in (11) and u is periodic in x.Let us recast scheme (11) in the following form:
(15)(E+αΛl)uin+1,s+1=α(Λl-Λh)uin+1,s+(E-αΛl)uin,
where α=βΔt/(2Δx2). If we define the following operators A=E+αΛl, B=E-αΛh, and C=E+αΛh, where E is the identity operator, then (15) can be rewritten as follows:
(16)Auin+1,s+1=(A-C)uin+1,s+Buin.
Assuming that the operators commute, (A-C)A=A(A-C) (e.g.,  in the case of uniform grid), it is easy to demonstrate that if uin+1,S^+1=uin+1 and uin+1,0=uin we get
(17)AS^+1uin+1=(∑k=0S^‍AS^-k(A-C)k)Buin+(A-C)S^+1uin.
Let uin=ξneIΘi, I=-1, be the solution of (11)–(14), where Θ=2πΔx/l is the phase angle with wavelength l. From (17), we can derive an equation for the amplification factor in the form
(18)|ξ|=|φ(Θ,S^,α)|,
where S^ is the number of iterations, and
(19)|φ(Θ,S^,α)|=|[(∑k=0S^AS^-k(A-C)k)B+(A-C)S^+1]eIΘi||AS^+1eIΘi|.
For stability of the method it is necessary that the absolute values of the amplification factor are less than one; that is,
(20)|ξ|<1.
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 |φ(Θ,S^,α)|. Figure 1 shows the values of |ξ|2 in the polar coordinate system (|ξ|2,Θ) for S^=1,3,and5. If only one iteration is executed in (11), S^=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)), S^=3, the amplification factor remains bounded by one at least for α≤10. In case of S^=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 S^→∞, our method becomes the unconditionally stable Crank-Nikolson method for the heat equation.Variation of amplification factor withΘ. (a) S^=1, dashed line α=2.0, solid line α=1.5, dash-doted line α=1.0, (b) S^=3, dashed line α=10.0, solid line α=5.0, dash-dotted line α=1.0, and (c) S^=5, dashed line α=30.0, solid line α=15, dash-doted line α=5.0, doted line α=1.0.
(a)
(b)
(c)


## 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):
(21)u1n+1,k=γ1(tn+1),uNn+1,k=γ2(tn+1).
The finite difference approximations at x2 and xN-1, which are the points next to the left and right boundaries, are straightforward:
(22)(uxxl)2n+1,k=1Δx2(γ1(tn+1)-2u2n+1,k+u3n+1,k),k=s,s+1,(uxxh)2n+1,s=112Δx2(10γ1(tn+1)-15u2n+1,s-4u3n+1,s+14u4n+1,s-6u5n+1,s+u6n+1,s),(uxxh)N-1n+1,s=112Δx2(10γ2(tn+1)-15uN-1n+1,s-4uN-2n+1,s+14uN-3n+1,s-6uN-4n+1,s+uN-5n+1,s),(uxxl)N-1n+1,k=1Δx2(uN-2n+1,k-2uN-1n+1,k+γ2(tn+1)),k=s,s+1.
Cases s=0 or k=0 give formulae to approximate (uxxl)in and (uxxh)in. Substituting (13), (14), and (22) into (12) the following fourth-order deferred correction approximations of (uxx)in+1,s+1,i=2,…,N-1, are
(23)(uxx)2n+1,s+1=56Δx2γ1(tn+1)+1Δx2(-2u2n+1,s+1+u3n+1,s+1)+112Δx2(9u2n+1,s-16u3n+1,s+14u4n+1,s-6u5n+1,s+u6n+1,s),(uxx)in+1,s+1=1Δx2(ui-1n+1,s+1-2uin+1,s+1+ui+1n+1,s+1)+112Δx2(-ui-2n+1,s+4ui-1n+1,s-6uin+1,s+4ui+1n+1,s-ui+2n+1,s),i=3,…,N-2,(uxx)N-1n+1,s+1=56Δx2γ2(tn+1)+1Δx2(-2uN-1n+1,s+1+uNn+1,s+1)+112Δx2(9uN-1n+1,s-16uN-2n+1,s+14uN-3n+1,s-6uN-4n+1,s+uN-5n+1,s).
Schemes (23) can be combined and expressed in the following matrix form:
(24)uxxn+1,s+1=1Δx2Λlun+1,s+1+1Δx2(Λh-Λl)un+1,s+γ(tn+1),
where Λl is a tridiagonal (N-2)×(N-2) matrix and Λh is the corresponding triangular and sparse (N-2)×(N-2) matrix,
(25)un+1,k=(u2n+1,k,u3n+1,k,…,uN-1n+1,k)T,k=s,s+1,(26)uxxn+1,s+1=((uxx)2n+1,s+1,(uxx)3n+1,s+1,…,(uxx)N-1n+1,s+1)T.
Substituting (6), (23) into (11), the formulae can be written into matrix form
(27)(E-αΛl)un+1,s+1=α(Λh-Λl)un+1,s+(E+αΛh)un+Δt[γ(tn+1)+γ(tn)]+Δtfn+1/2.
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 tn 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 O(Δt2,Δx4) 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
(28)u(0,mΔt)=γ1(tm)=u1m,u(l,mΔt)=γ2(tm)=uNm
are used to derive the following fourth-order schemes
(29)(uxx)2m+σ(uxx)3m=124Δx2(a1u1m+a2u2m+a3u3m+a4u4m+a5u5m+a6u6m)=a124Δx2γ1(tm)+124Δx2(a2u2m+a3u3m+a4u4m+a5u5m+a6u6m),(uxx)i-1m+10(uxx)im+(uxx)i+1m=2Δx2(6ui-1m-12uim+6ui+1m),i=2,…,N-1,(uxx)N-1m+σ(uxx)N-2m=124Δx2(a1uNm+a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m)=a124Δx2γ2(tm)+124Δx2(a2uN-1m+a3uN-2m+a4uN-3m+a5uN-4m+a6uN-5m),
where the coefficients can be found by matching the Taylor series expansion of left-hand-side terms up to order O(Δx4)u(6) which gives the following values of coefficients [8]:
(30)σ=12,a1=19,a2=-14,a3=-38,a4=44,a5=-13,a6=2.
Then all derivatives in (4) are approximated by the fourth-order compact formula; we can write
(31)Auxxm=1Δx2Bum+γm,m=n,n+1,
where A and B are the corresponding triangular and sparse (N-2)×(N-2) matrices, uxxm=((uxx)2m,(uxx)3m,…,(uxx)N-1m)T, um=(u2m,u3m,…,uN-1m)T and γm = (γ1(tm),0,…,0,γ2(tm))T, m = n,n+1. Schemes (4) and (29) can be combined and expressed in the following matrix form:
(32)(A-αB)un+1=(A+αB)un+Δt[γ(tn+1)+γ(tn)]+Δtfn+1/2.
This scheme is called FOCs for Dirichlet boundary value problem (1)–(3). We like to mention that the above scheme has truncation error O(Δt2,Δx4). Note that the triangular and sparse (N-2)×(N-2) coefficient matrices in FOCs are time independent; hence, we have to store the inverse of the coefficient matrix A-αB 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Ω={x∣0≤x≤1}. In all computations, we used Δt=Δx2/4 and ϵ=10-10. The following stopping criterion is used:
(33)max1≤i≤N|uin+1,S^+1-uin+1,S^|<ϵ,s=0,…,S^,
where  “S^” 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 for0≤t≤1.Example 1 (the homogeneous heat equation with the homogeneous Dirichlet boundary conditions).
One has(34)ut=uxx,0≤x≤1,t>0,u(x,0)=sin(πx),u(0,t)=0,u(1,t)=0.
The exact solution is u(x,t)=e-π2tsin(πx). The results of performance over the time interval t∈[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 t=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 t=1.
Types of scheme
Grid points
Maximum error
Rate of convergence
Aver. number of iteration
CPU time in sec.
11
3.8687
×
1
0
-
8
—
1
0.006
21
6.0426
×
1
0
-
10
6.0005
1
0.042
FOCs
41
2.2454
×
1
0
-
11
4.7501
1
0.326
81
1.2821
×
1
0
-
12
4.1304
1
2.564
161
8.0164
×
1
0
-
14
3.9994
1
20.437
11
9.9767
×
1
0
-
9
—
4
0.015
21
1.4996
×
1
0
-
9
2.7361
3
0.085
FODs
41
1.1193
×
1
0
-
10
3.7438
2
0.438
81
7.1438
×
1
0
-
12
3.9698
2
3.450
161
4.4797
×
1
0
-
13
4.1875
2
27.495
11
9.9763
×
1
0
-
9
—
1
0.006
21
1.4996
×
1
0
-
9
2.7361
1
0.043
FOSs
41
1.1193
×
1
0
-
10
3.7438
1
0.334
81
7.1440
×
1
0
-
12
3.9698
1
2.623
161
4.4854
×
1
0
-
13
4.1875
1
20.907Example 2 (the nonhomogeneous heat equation with non-homogeneous Dirichlet boundary conditions).
One has(35)ut=uxx+(π2-1)e-tcos(πx)+4x-2,0≤x≤1,t>0,u(x,0)=cos(πx)+x2,u(0,t)=e-t,u(1,t)=-e-t+4t+1.
The exact solution is u(x,t)=e-tcos(πx)+x2+4xt. The results of performance over the time domain t∈[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 t=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 t=1.
Types of scheme
Grid points
Maximum error
Rate of convergence
Aver. number of iteration
CPU time in sec.
11
1.8470
×
1
0
-
5
—
1
0.006
21
3.6901
×
1
0
-
7
5.6454
1
0.046
FOCs
41
7.5595
×
1
0
-
9
5.6092
1
0.353
81
6.6458
×
1
0
-
10
3.5077
1
2.778
161
4.8841
×
1
0
-
11
3.7663
1
22.141
11
1.7132
×
1
0
-
5
—
7
0.016
21
2.6914
×
1
0
-
7
5.9922
7
0.128
FODs
41
2.7910
×
1
0
-
8
3.2655
6
0.851
81
2.0112
×
1
0
-
9
3.7941
5
5.568
161
1.3116
×
1
0
-
10
3.9415
5
44.375
11
1.2895
×
1
0
-
5
—
1
0.006
21
2.8544
×
1
0
-
7
5.9922
1
0.046
FOSs
41
2.7306
×
1
0
-
8
3.2655
1
0.359
81
1.9590
×
1
0
-
9
3.7941
1
2.821
161
1.3130
×
1
0
-
10
3.9415
1
22.484The last two columns of Tables1 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 t=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(36)ut=βuxx-uux,0≤x≤1,t>0,
with the exact solution [21] is given by
(37)u(x,t)=ξ+η+(η-ξ)eρ1+eρ,
where ρ=ξ(x-ηt-υ)/β. 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:
(38)uin+1,s+1-uinΔt=β2[(uxx)in+1,s+1+(uxx)in]+fin,
where fin=-[(u2/2)x]in. The nonlinear term fin 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:
(39)(E-αΛl)un+1,s+1=α(Λh-Λl)un+1,s+(E+αΛh)un+Δt[γ(tn+1)+γ(tn)]+Δtfn,
where E is identity matrix, Λl is tridiagonal (N-2)×(N-2) matrix, and Λh is the corresponding triangular and sparse (N-2)×(N-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 interval0<t≤1. The results of performance over the time interval t∈[0,1] for the FODs are represented in Tables 3 and 4, where the maximum error and the rate of convergence at time instant t=1 are shown.Table 3
Maximum absolute error, order of convergence, and CPU time in seconds for Example3 at time instant t=1 with fixed mesh size Δx=0.05.
Types of scheme
Time step sizes
Maximum error
Rate of convergence
Aver. number of iteration
CPU time in sec.
FODs
1
0
-
2
8.0945
×
1
0
-
6
—
10
0.015
1
0
-
3
8.0942
×
1
0
-
7
1.0000
6
0.109
1
0
-
4
8.1144
×
1
0
-
8
0.9989
3
0.656
1
0
-
5
8.2989
×
1
0
-
9
0.9902
3
6.281Table 4
Maximum absolute error, order of convergence, and CPU time in seconds for Example3 at time instant t=1 with time step size Δt=Δx4.
Types of scheme
Grid points
Maximum error
Rate of convergence
Aver. number of iteration
CPU time in sec.
11
8.3375
×
1
0
-
8
—
3
0.516
FODs
21
5.2632
×
1
0
-
9
3.9812
2
8.594
41
3.3491
×
1
0
-
10
3.9774
2
200.015In 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 t=1 when Δx=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 t=1 and using uniform grids of 11, 21, and 41 with time step sizes Δt=Δx4 and ϵ=10-10.It can be seen from Tables3 and 4 that numerical results are in good agreement with the exact solution. We only observe O(Δt) 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 n, instead of at time level n+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 makeuux term as follows [21, 22]:
(40)un+1(ux)in+1≈un+1(ux)in+un(ux)in+1-un(ux)in,
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.



---

*Source: 574620-2013-03-31.xml*