Journal of Chemistry

Volume 2015, Article ID 350362, 9 pages

http://dx.doi.org/10.1155/2015/350362

## Solving Reaction-Diffusion and Advection Problems with Richardson Extrapolation

^{1}Department of Meteorology, Eötvös Loránd University, Pázmány Péter sétány 1/A, Budapest 1117, Hungary^{2}Department of Physics, Budapest University of Technology and Economics, Budafoki út 8, Budapest 1111, Hungary^{3}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 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.