Numerical and Soft Computing Methods for Characteristic Value Problems of ODE and ODEs SystemsView this Special Issue
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-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.
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 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 . 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: 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 and the mesh size is . 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: where . The Dirichlet boundary conditions 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 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 matrix, Substituting (6) into (4) gives us the following matrix form: where , , and denote the 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 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: where and the second superscript “” denotes the number of iterations and .
The deferred correction technique  is utilized to approximate the second-order derivatives at higher time levels by the iterative method where , , , is high-order approximation on wide stencil and , , , 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 we use the solution from the time level (so and ). Once the iterations converge, the lower-order approximation terms drop out and the approximation of obtained has the same order of approximation as . 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 , 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 in (13) gives the fourth order of approximation to approximate the second-order derivatives at the time level .
2.2.1. Stability Analysis
Let us recast scheme (11) in the following form: where . 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 and we get Let , , be the solution of (11)–(14), where 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 . Figure 1 shows the values of in the polar coordinate system for . If only one iteration is executed in (11), , inequality (20) holds if , as can be seen from Figure 1(a)). If 3 iterations are done in (11) (Figure 1(b)), , the amplification factor remains bounded by one at least for . In case of , the stability criteria hold up to 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.
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): The finite difference approximations at and , which are the points next to the left and right boundaries, are straightforward: Cases or give formulae to approximate and . Substituting (13), (14), and (22) into (12) the following fourth-order deferred correction approximations of , are Schemes (23) can be combined and expressed in the following matrix form: where is a tridiagonal matrix and is the corresponding triangular and sparse matrix, Substituting (6), (23) into (11), the formulae can be written into matrix form 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 are used to derive the following fourth-order schemes 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 : Then all derivatives in (4) are approximated by the fourth-order compact formula; we can write where and are the corresponding triangular and sparse matrices, , and , . Schemes (4) and (29) can be combined and expressed in the following matrix form: 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 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 . In all computations, we used and . 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 .
Example 1 (the homogeneous heat equation with the homogeneous Dirichlet boundary conditions). One has The exact solution is . The results of performance over the time interval for the FOCs, FODs, and FOSs are represented in Table 1, where the maximum error and the rate of convergence at time instant are shown.
Example 2 (the nonhomogeneous heat equation with non-homogeneous Dirichlet boundary conditions). One has The exact solution is . The results of performance over the time domain for the FOC, FOD, and FOS schemes are represented in Table 2, where the maximum error and the rate of convergence at time instant 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 . 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  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 matrix, and is the corresponding triangular and sparse matrix and can be solved by using Thomas algorithm.
Example 3 (the Burgers equation (36) and the constant values , , , and 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 . The results of performance over the time interval for the FODs are represented in Tables 3 and 4, where the maximum error and the rate of convergence at time instant 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 when 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 and using uniform grids of 11, 21, and 41 with time step sizes and .
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 for the FODs (38).
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]: 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.
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.
I. Christie, “Upwind compact finite difference schemes,” Journal of Computational Physics, vol. 59, no. 3, pp. 353–368, 1985.View at: Google Scholar
J. H. Ferziger and M. Peric, Computational Methods in Fluid Dynamics, Springer, Berlin, Germany, 2002.
J. C. Kalita, D. C. Dalal, and A. K. Dass, “A class of higher order compact schemes for the unsteady two-dimensional convection-diffusion equation with variable convection coefficients,” International Journal for Numerical Methods in Fluids, vol. 38, no. 12, pp. 1111–1131, 2002.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
J. Zhao, W. Dai, and S. Zhang, “Fourth-order compact schemes for solving multidimensional heat problems with Neumann boundary conditions,” Numerical Methods for Partial Differential Equations, vol. 24, no. 1, pp. 165–178, 2008.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
M. C. Miriane, M. G. Paulo, and C. R. Estaner, “Error analysis in the solution of unsteady nonlinear convection-diffusion problems,” Research Journal of Physical and Applied Science, vol. 1, no. 1, pp. 20–22, 2012.View at: Google Scholar
B. Jiang, The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics, Scientific Computation, Springer, Berlin, Germany, 1998.View at: MathSciNet