Fourth-Order Deferred Correction Scheme for Solving Heat Conduction Problem

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-ordermethods proposed are the easier development and can be solved by usingThomas algorithms.The stability analysis and numerical experiments have been limited to one-dimensional heat-conducting problems with Dirichlet boundary conditions and initial data.


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 highorder 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 secondorder discretization of heat equation.
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 semiimplicit (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 <  < , Dirichlet BC :  (0, ) =  1 () , (, ) =  2 () ,  > 0, where the diffusion coefficient  is positive, (, ) represents the temperature at point (, ), and (, ),  1 (),  2 () 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.

The Fourth-Order Schemes
Let Δ denote the temporal mesh size.For simplicity, we consider a uniform mesh consisting of  points:  1 ,  2 , . . .,   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 .

Fourth-Order
Semi-Implicit Scheme.The application to the well-known Crank-Nikolson scheme to (1) results in the following expression: where are used to derive the following fourth-order approximation of second derivative terms: where the coefficients can be found by matching the Taylor series expansion of left-hand-side terms up to order (Δ 4 ) (6) which gives the following values of coefficients: Schemes ( 6) can be combined and expressed in the following matrix form: where Λ ℎ is the corresponding triangular and sparse (−2)× ( − 2) matrix, Substituting ( 6) into (4) gives us the following matrix form: where  = Δ/(2Δ 2 ), f +1/2 = ( +1/2 2 ,  +1/2 3 , . . .,  +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 (Δ 2 , Δ 4 ) 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.
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): 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: Case  = 0 in (13) gives the fourth order of approximation to approximate the second-order derivatives at the time level .

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: where  = Δ/(2Δ 2 ).If we define the following operators  =  + Λ  ,  =  − Λ ℎ , and  =  + Λ ℎ , where  is the identity operator, then (15) can be rewritten as follows: 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 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 where Ŝ is the number of iterations, and For stability of the method it is necessary that the absolute values of the amplification factor are less than one; that is, 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 |(Θ, Ŝ, )|.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.

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

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 Δ = Δ 2 /4 and  = 10 −10 .The following stopping criterion is used: 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 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.
Example 2 (the nonhomogeneous heat equation with non-homogeneous Dirichlet boundary conditions).One has The exact solution is (, ) =  − cos() +  2 + 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.
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 with the exact solution [21] is given by 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: where .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: 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.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 Δ = Δ 4 and  = 10 −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).

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 secondorder 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]: 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.

Figure 1
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 Figure1(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.

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.

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.

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 Δ = Δ 4 .