Abstract
Richardson extrapolation is a simple but powerful computational tool to enhance the accuracy of time integration methods. In the past years a few theoretical and partly practical works have been presented on this method. Detailed numerical applications of this method, however, are rarely found in the literature. Therefore, it is worth investigating whether this promising technique lives up to the expectations also in practice. In this paper we investigate the efficiency of the Richardson method in one-dimensional numerical (reaction-diffusion) problems.
1. Introduction
Phenomena occurring in nature and in laboratories can be understood and described by mathematical models. Most of those models are time-dependent ordinary or partial differential equations that cannot be solved analytically due to the highly nonlinear nature of the functions appearing in the equations. To overcome this problem, several numerical solving strategies have been developed in the past few decades, which can provide an accurate numerical solution to the original problem [1, 2].
Partial differential equations (PDEs), containing time and space as independent variables, can describe the spatiotemporal evolution of a set of physical quantities (physical system). This mathematical framework is used to understand and describe physical systems in various fields of science and technology. There are many processes (e.g., mass transport and pattern formation) which can be described and understood by reaction-diffusion systems. Generally, reaction-diffusion systems are mathematical models that describe the spatial and temporal variations of concentrations of chemical substances. From the mathematical point of view, the reaction-diffusion system is a set of parabolic PDEs [3].
The solution of PDEs should be accurate and the computational time should be small. There are several situations where the set of PDEs is used to make predictions such as in the case of weather prediction models, and thus the computational time should be dramatically reduced to use the solution in real-time applications. There are predominantly two strategies that can be used to increase the accuracy of the solution and reduce the computational time. First is developing sophisticated numerical schemes and methods that can increase the accuracy [4, 5]. This involves more computational effort. Second, computational time of mathematical models consisting of PDEs can be efficiently reduced, applying parallelization strategies on supercomputers, clusters, grids, and graphical processing units (video cards) [3, 6–10].
In this study, we show an efficient numerical method called Richardson extrapolation to increase the accuracy of the solution of various reaction-diffusion and advection systems. Four different problems are investigated and examined, including (i) simple diffusion of a chemically inert compound, (ii) diffusion with first- and second-order reaction of a chemical compound, (iii) advection process, and (iv) phase separation (pattern formation) in the wake of a moving diffusion front.
2. Richardson Extrapolation Methods
Richardson extrapolation is a powerful tool to increase the accuracy of any numerical method. It consists in applying the given numerical scheme with different discretization parameters (usually and ) and combining the obtained numerical solutions by properly chosen weights. Namely, if denotes the order of the selected numerical method, denotes the numerical solution obtained by , and denotes that obtained by , then the combined solutionhas order . This method was first extensively used by Richardson, who called it “the deferred approach to the limit” [11]. The Richardson extrapolation is especially widely used for time integration schemes, when, as a rule, the results obtained by two different time step sizes are combined.
The Richardson extrapolation can be implemented in two different ways. During the passive Richardson extrapolation the combined solution is not used in the further computations, while in the active version it serves as an initial value for the next time step. Results concerning the stability and convergence of different numerical methods combined with the Richardson extrapolation can be found in [12–14]. In this paper, we investigate both the passive and the active Richardson methods in the selected reaction-diffusion and advection problems.
3. Model Applications
In the following sections we investigate the efficiency of the passive and active Richardson extrapolations on several one-dimensional model problems. First the problems were discretized in space using second-order finite differences and a small spatial step size. Then the obtained systems of time-dependent ordinary differential equations were solved by the forward Euler method, both with and without Richardson extrapolation. In this manner the errors resulting from the spatial discretization can be assumed to be much smaller than those arising from the time discretization. Since the forward Euler method is of first order, the Richardson extrapolation should result in a second-order time discretization method.
In all the simulations, each model problem was solved on the same time interval using different time step sizes. In each case we measured the running time, which can be considered as the computational time needed by the given method, since the optimal exploitation of the computer memory was not an issue in these simple models. The codes were written in C programming language. The running times were measured by the clock( ) function contained by the time.h header. All the computations were performed on the same Toshiba Satellite L750-1n3 laptop by Ubuntu 12.04 LTS operational system, and the C compiler g++ was used, which is easy to apply on a Linux terminal.
The model results were compared with a reference solution obtained by the forward Euler method with a very small time step size (). The grid step size in all simulation was fixed ().
4. Diffusion
Our basic model problem was the one-dimensional diffusion equationequipped with Neumann boundary condition and an initial condition describing a peak (constant zero function with the exception of the middle point of the spatial domain where the value was 1). The value of the diffusion coefficient was equal to 1 for simplicity. We generated a mesh with step sizes and on the space-time domain of the problem. Then, after a standard spatial discretization of the second-order spatial derivative, we applied the forward Euler method for time integration, which resulted in the formula
Here stands for the approximation of the exact solution at the mesh point and time layer . This finite difference scheme is illustrated in Figure 1.
The numerical solutions obtained by Figure 1 with and without Richardson extrapolation were compared to the reference solution, and the absolute and relative errors were evaluated. Table 1 shows the absolute errors for decreasing time step sizes. Here and in the further tables the values given in parentheses in the th row of the table show the estimated error order computed from the error obtained for the given time step () and the error obtained for the previous time step () according to the formulaFor a first-order method we should obtain values of around 1, while for a second-order formula we should obtain values of around 2. One can see that the forward Euler method behaves like a first-order method, while both Richardson extrapolations increased the order by one, according to the expectations. The higher order can be observed only for larger time steps, and then the decrease of the errors slows down. This can be explained by the fact that the accuracy of the computation is limited by the accuracy of the reference solution, which is a numerical solution obtained by a very small time step and the same spatial grid size. This maximal accuracy is achieved by the Richardson extrapolation rather early, which shows the robustness of this method. Note that the accuracy that can be achieved by the pure forward Euler method by a time step size of can be achieved by the passive Richardson method already by using a time step size of and by the active Richardson method by a time step size of . Figures 2 and 3 also show that both Richardson methods are significantly more accurate than the forward Euler method in itself when the same time step is used. However, it is not worth reducing the time step size beyond all bounds.
Accuracy is not the only point of view that we should consider. It is also important that the results should be achieved within reasonable computational time. Table 2 shows the running times for the studied methods and the step lengths. One can see that the extrapolation methods roughly take twice as much time as the pure Euler method for the same time step; however, as Figure 4 tells us, they are much more accurate.
For example, if a relative error of 0.7 is required, then the forward Euler method without Richardson extrapolation needs 4.76 s, with passive Richardson extrapolation 0.12 s, and with the active one 0.29 s, to achieve this goal. The latter method results in an almost fortyfold acceleration, while the previous one results in a 16-fold acceleration in the computation. The passive extrapolation is the most efficient method from the three. Interestingly, it has a local minimum in the relative error (here the solution was obtained in 10.79 s). The result of this run is more accurate by more than two magnitudes than the most accurate Eulerian solution that we could obtain.
5. Diffusion with Linear or Quadratic Reaction Terms
For a more extensive analysis of the diffusion problem we supplemented the diffusion term with (i) a linear and (ii) a quadratic reaction term; that is, the following equations were considered:where the reaction rate constant was equal to 1. The diffusion coefficient, the meshes, and the supplementary conditions were the same as in the pure diffusion problem studied in the previous section.
The obtained absolute errors are given in Tables 3 and 4. In the case of the extra linear term, similar results have been obtained as in the pure diffusion problem; that is, the extrapolation methods significantly enhanced the efficiency. However, the addition of the quadratic term did cause some slight differences in the appearance of the error curves. Figure 5 shows no stationary section in the domain of the small time step sizes that were observed in the pure diffusion case. However, a thorough analysis reveals that it is not worth using extremely small time steps in this case, either.
6. Advection
One-dimensional advection is described by the equationwith advection velocity , which was chosen to be 1. The initial condition was the same peak as in the diffusion model, and Neumann-type boundary condition was used at the inflow boundary. The first-order spatial derivative was discretized by a forward difference scheme, and then the forward Euler method was applied for time integration. This procedure leads to the relation
Figure 6 illustrates the applied discretization scheme.
Table 5 and Figures 7 and 8 show the absolute and relative errors of the applied methods. Note that both extrapolation methods produced smaller absolute and relative errors than the Euler method. For large time steps the active Richardson extrapolation is only slightly better than the Eulerian scheme, and then the Richardson extrapolation becomes more and more accurate until the errors cannot decrease any further.
Figure 9 shows the relative errors as a function of the running time. Although the active Richardson extrapolation performs better than the pure Euler method, it is far less efficient than the passive Richardson method, which only needs 19.16 s to produce the most accurate numerical solution. This solution is 500 times more accurate than the best Eulerian solution.
To summarize, both extrapolation methods are suitable for solving the advection problem, but the passive method is significantly more efficient.
7. The Cahn-Hilliard Model
Finally, we consider a more complex reaction-diffusion model problem, the so-called Cahn-Hilliard model [15–17]. It is suitable for describing the phenomenon of phase separation, which is used for the simulation of chemical processes as well as for the formation of certain biological patterns and nucleation in meteorology [18, 19].
We chose the above model because it contains a fourth-order derivative of the unknown function. Our aim was to investigate whether the Richardson extrapolation methods show any noticeable difference due to the presence of the higher-order derivative in a complex reaction-diffusion system. We consider a pattern formation (phase separation) in the wake of a reaction front. This front emerges due to a strongly inhomogeneous initial distribution of the chemical species and . The reaction takes place in a domain occupying the half-space and, initially, the inner electrolyte is homogeneously distributed in it . Chemical species with a much higher initial concentration is brought into contact with the domain at . Assuming a chemical reaction between reactants , the phase separation (pattern formation) behind the advancing front can be described by the set of equationswhere is the reaction rate constant and, for simplicity, the diffusion constants of the reagents are taken to be equal (), and , , , and characterize the phase separation process.
Solving the equations above, we can simulate the pattern formation of chemical species in the wake of the reaction front. The reference solution was computed again using the forward Euler method with a very small time step ().
According to the results, the methods are not sensitive to the appearance of the higher-order derivatives. Both extrapolation methods can handle the problem very well and produce much better results than the application of the forward Eulerian method without extrapolation (Figure 10).
8. Conclusion
In this study, we showed a powerful integration method called Richardson extrapolation, which can be efficiently used to solve various reaction-diffusion and advection problems. This method increases the accuracy of the numerical results at the cost of slightly increased computational time. However, the Richardson extrapolation can provide an increase of around two orders of magnitude in the computational time to achieve the same error tolerance compared to a numerical scheme without using the extrapolation method. Combining this cost-efficient method with parallel computing techniques can provide a powerful method to solve partial differential equations.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the Hungarian Research Found (OTKA K104666), the European Union, and the State of Hungary and cofinanced by the European Social Fund in the framework of TÁMOP 4.2.4. A/1-11-1-2012-0001 National Excellence Program. The European Union and the European Social Fund have provided financial support to the project under the Grant Agreement no. TÁMOP 4.2.1./B-09/1/KMR-2010-0003.