Research Article  Open Access
Tamás Mona, István Lagzi, Ágnes Havasi, "Solving ReactionDiffusion and Advection Problems with Richardson Extrapolation", Journal of Chemistry, vol. 2015, Article ID 350362, 9 pages, 2015. https://doi.org/10.1155/2015/350362
Solving ReactionDiffusion and Advection Problems with Richardson Extrapolation
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 onedimensional numerical (reactiondiffusion) problems.
1. Introduction
Phenomena occurring in nature and in laboratories can be understood and described by mathematical models. Most of those models are timedependent 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 reactiondiffusion systems. Generally, reactiondiffusion systems are mathematical models that describe the spatial and temporal variations of concentrations of chemical substances. From the mathematical point of view, the reactiondiffusion 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 realtime 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 reactiondiffusion and advection systems. Four different problems are investigated and examined, including (i) simple diffusion of a chemically inert compound, (ii) diffusion with first and secondorder 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 reactiondiffusion and advection problems.
3. Model Applications
In the following sections we investigate the efficiency of the passive and active Richardson extrapolations on several onedimensional model problems. First the problems were discretized in space using secondorder finite differences and a small spatial step size. Then the obtained systems of timedependent 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 secondorder 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 L7501n3 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 onedimensional 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 spacetime domain of the problem. Then, after a standard spatial discretization of the secondorder 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 firstorder method we should obtain values of around 1, while for a secondorder formula we should obtain values of around 2. One can see that the forward Euler method behaves like a firstorder 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 16fold 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
Onedimensional 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 Neumanntype boundary condition was used at the inflow boundary. The firstorder 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.
