Journal of Chemistry
Volume 2015 (2015), Article ID 350362, 9 pages
http://dx.doi.org/10.1155/2015/350362

Research Article

Solving Reaction-Diffusion and Advection Problems with Richardson Extrapolation

Tamás Mona,¹ István Lagzi,² and Ágnes Havasi³

¹Department of Meteorology, Eötvös Loránd University, Pázmány Péter sétány 1/A, Budapest 1117, Hungary
²Department of Physics, Budapest University of Technology and Economics, Budafoki út 8, Budapest 1111, Hungary
³Department of Applied Analysis and Computational Mathematics, Eötvös Loránd University and MTA-ELTE Numerical Analysis and Large Networks Research Group, Pázmány Péter sétány 1/C, Budapest 1117, Hungary

Received 9 December 2014; Accepted 11 March 2015

Academic Editor: Henryk Kozlowski

Copyright © 2015 Tamás Mona et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 solution

has 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 equation

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

Figure 1: Discretization scheme of the diffusion problem.

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 formula

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

Table 1: Absolute error of the diffusion problem without and with passive/active Richardson extrapolation (with the estimated convergence order in parentheses).

Time step
Absolute error

Euler
Passive Richardson
Active Richardson

—

—

—

(1.000)

(2.019)

(1.991)

(1.000)

(2.027)

(1.987)

(1.000)

(2.035)

(1.978)

(1.000)

(2.043)

(1.971)

(1.000)

(2.054)

(1.961)

(1.000)

(2.117)

(1.852)

(1.000)

(2.109)

(1.413)

(1.000)

(1.183)

(0.934)

(1.000)

(0.520)

(0.650)

(1.000)

(0.252)

(0.308)

(1.000)

(0.016)

(0.024)

(1.000)

(−0.004)

(−0.001)

(0.999)

(0.000)

(0.000)

(0.997)

(−0.001)

(0.001)

(0.991)

(0.000)

(0.000)

(0.975)

(0.002)

(0.000)

Figure 2: Absolute error as a function of time step for the Euler method without and with passive and active Richardson extrapolation, by grid step and coefficient . The reference solution was calculated using .

Figure 3: Relative error as a function of time step for the Euler method without and with passive and active Richardson extrapolation, by grid step and coefficient . The reference solution was calculated using .

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.

Table 2: Computational times of the diffusion problem without and with passive/active Richardson extrapolation.

Time step
Computational time (s)

Euler
Passive Richardson
Active Richardson

0.20
0.05
0.12
0.14

0.17
0.07
0.13
0.17

0.14
0.08
0.15
0.20

0.13
0.08
0.18
0.23

0.11
0.09
0.20
0.25

0.10
0.10
0.22
0.29

0.05
0.19
0.43
0.56

0.03
0.29
0.66
0.83

0.03
0.38
0.87
1.11

0.02
0.48
1.08
1.39

0.01
0.96
2.16
2.76

0.002
4.76
10.79
14.36

0.001
9.57
21.64
27.62

0.0002
47.95
106.92
138.50

0.0001
94.58
213.67
275.30

0.00002
474.86
1077.50
1381.50

0.00001
953.78
2153.56
2885.96

Figure 4: Relative error as a function of computational time for the Euler method without and with passive and active Richardson extrapolation, by grid step and coefficient . The reference solution was calculated using .

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.

Table 3: Absolute error of the diffusion problem with linear reaction term without and with passive/active Richardson extrapolation (with the estimated convergence order in parentheses).

Time step
Absolute error

Euler
Passive Richardson
Active Richardson

—

—

—

(1.000)

(2.012)

(1.990)

(1.000)

(2.016)

(1.986)

(1.000)

(2.021)

(1.982)

(1.000)

(2.027)

(1.976)

(1.000)

(2.034)

(1.970)

(1.000)

(2.083)

(1.926)

(1.000)

(2.247)

(1.786)

(1.000)

(2.402)

(1.561)

(1.000)

(2.364)

(1.263)

(1.000)

(0.751)

(0.668)

(1.000)

(−0.179)

(0.040)

(1.000)

(−0.015)

(0.001)

(1.000)

(−0.002)

(0.000)

(0.999)

(0.000)

(0.000)

(0.997)

(0.000)

(0.000)

(0.990)

(0.001)

(0.001)

Table 4: Absolute error of the diffusion problem with quadratic reaction term without and with passive/active Richardson extrapolation (with the estimated convergence order in parentheses).

Time step
Absolute error

Euler
Passive Richardson
Active Richardson

—

—

—

(1.027)

(2.134)

(2.596)

(1.019)

(2.138)

(2.469)

(1.015)

(2.135)

(2.388)

(1.012)

(2.128)

(2.330)

(1.011)

(2.121)

(2.288)

(1.007)

(2.091)

(2.185)

(1.003)

(2.056)

(2.098)

(1.002)

(2.040)

(2.065)

(1.002)

(2.030)

(2.047)

(1.001)

(2.014)

(2.019)

(1.000)

(1.873)

(1.719)

(1.000)

(1.035)

(0.618)

(1.000)

(0.171)

(0.086)

(1.001)

(0.011)

(0.006)

(1.002)

(0.002)

(0.000)

(1.007)

(−0.007)

(−0.005)

Figure 5: Absolute error as a function of time step in case of a quadratic reaction term for the Euler method without and with passive and active Richardson extrapolation, by grid step and coefficient . The reference solution was calculated using .

6. Advection

One-dimensional advection is described by the equation

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

Figure 6: Discretization scheme for the advection problem.

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.

Table 5: Absolute error of the advection problem without and with passive/active Richardson extrapolation (with the estimated convergence order in parentheses).

Time step
Absolute error

Euler
Passive Richardson
Active Richardson

—

—

—

(1.009)

(2.000)

(2.019)

(1.008)

(2.000)

(2.016)

(1.007)

(2.000)

(2.013)

(1.006)

(2.000)

(2.012)

(1.005)

(2.000)

(2.010)

(1.004)

(2.000)

(2.007)

(1.002)

(1.999)

(2.002)

(1.001)

(1.998)

(1.998)

(1.001)

(1.996)

(1.995)

(1.001)

(1.965)

(1.981)

(1.000)

(1.023)

(1.744)

(1.000)

(0.037)

(0.814)

(1.000)

(0.002)

(0.139)

(1.001)

(0.000)

(0.011)

(1.002)

(0.000)

(0.001)

(1.007)

(0.000)

(0.000)

Figure 7: Absolute error as a function of time step in case of the advection problem for the Euler method without and with passive and active Richardson extrapolation. In the simulation the grid step was and advection velocity . The reference solution was calculated using .

Figure 8: Relative error as a function of time step for the Euler method without and with passive and active Richardson extrapolation. In the simulation the grid step was and the advection velocity . The reference solution was calculated using .

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.

Figure 9: Relative error as a function of computational time for the Euler method without and with passive and active Richardson extrapolation. In the simulation the grid step was and the advection velocity . The reference solution was calculated using .

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 equations

where 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).

Figure 10: Results of the numerical simulations for the Cahn-Hilliard type reaction-diffusion model obtained by the Euler method without and with passive and active Richardson extrapolations. We used the following parameters in the simulations: , , , , and . The reference solution was calculated using .

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.

References

J. D. Lambert, Numerical Methods for Ordinary Differential Systems: The Initial Value Problem, John Wiley & Sons, New York, NY, USA, 1991.

J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, 2008.

F. Molnár Jr., F. Izsák, R. Mészáros, and I. Lagzi, “Simulation of reaction—diffusion processes in three dimensions using CUDA,” Chemometrics and Intelligent Laboratory Systems, vol. 108, no. 1, pp. 76–85, 2011.

P. A. Zegeling and H. P. Kok, “Adaptive moving mesh computations for reaction-diffusion systems,” Journal of Computational and Applied Mathematics, vol. 168, no. 1-2, pp. 519–528, 2004.

I. Lagzi, D. Kármán, T. Turányi, A. S. Tomlin, and L. Haszpra, “Simulation of the dispersion of nuclear contamination using an adaptive Eulerian grid model,” Journal of Environmental Radioactivity, vol. 75, no. 1, pp. 59–82, 2004.

F. Molnár Jr., T. Szakály, R. Mészáros, and I. Lagzi, “Air pollution modelling using a graphics processing unit with CUDA,” Computer Physics Communications, vol. 181, no. 1, pp. 105–112, 2010.

M. Januszewski and M. Kostur, “Accelerating numerical solution of stochastic differential equations with CUDA,” Computer Physics Communications, vol. 181, no. 1, pp. 183–188, 2010.

A. R. Sanderson, M. D. Meyer, R. M. Kirby, and C. R. Johnson, “A framework for exploring numerical solutions of advection-reaction-diffusion equations using a GPU-based approach,” Computing and Visualization in Science, vol. 12, no. 4, pp. 155–170, 2009.

V. N. Alexandrov, W. Owczarz, P. G. Thomson, and Z. Zlatev, “Parallel runs of a large air pollution model on a grid of Sun computers,” Mathematics and Computers in Simulation, vol. 65, no. 6, pp. 557–577, 2004.

C. Németh, G. Dózsa, R. Lovas, and P. Kacsuk, “The p-grade grid portal,” in Computational Science and Its Applications—ICCSA 2004, pp. 10–19, Springer, 2004.

L. F. Richardson and J. A. Gaunt, “The deferred approach to the limit. Part I. Single lattice. Part II. Interpenetrating lattices,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 226, pp. 299–361, 1927.

I. Farag, G. Havasi, and Z. Zlatev, “Efficient implementation of stable Richardson extrapolation algorithms,” Computers & Mathematics with Applications, vol. 60, no. 8, pp. 2309–2325, 2010.

Z. Zlatev, I. Faragó, and Á. Havasi, “Stability of the Richardson Extrapolation applied together with the θ-method,” Journal of Computational and Applied Mathematics, vol. 235, no. 2, pp. 507–517, 2010.

I. Faragó, Á. Havasi, and Z. Zlatev, “The convergence of explicit Runge-Kutta methods combined with Richardson extrapolation,” in Applications of Mathematics, pp. 99–106, GCSE, 2012.

J. W. Cahn and J. E. Hilliard, “Free energy of a nonuniform system. I. Interfacial free energy,” The Journal of Chemical Physics, vol. 28, no. 2, pp. 258–267, 1958.

J. W. Cahn and J. E. Hilliard, “Spinodal decomposition: a reprise,” Acta Metallurgica, vol. 19, no. 2, pp. 151–161, 1971.

C. M. Elliott and D. A. French, “Numerical studies of the Cahn-Hilliard equation for phase separation,” IMA Journal of Applied Mathematics, vol. 38, no. 2, pp. 97–128, 1987.

Z. Rácz, “Formation of Liesegang patterns,” Physica A: Statistical Mechanics and its Applications, vol. 274, no. 1, pp. 50–59, 1999.

T. Antal, M. Droz, J. Magnin, A. Pekalski, and Z. Rácz, “Formation of Liesegang patterns: simulations using a kinetic Ising model,” The Journal of Chemical Physics, vol. 114, no. 8, pp. 3770–3775, 2001.