Exponentially Fitted Finite Difference Schemes for Reaction-Diffusion Equations

The purpose of this work is to introduce a new kind of finite difference formulation inspired from Fourier analysis, for reaction-diffusion equations. Compared to classical schemes, the proposed scheme is much more accurate and has interesting stability properties. Convergence properties and stability of the scheme are discussed. Numerical examples are provided to show better performance of the method, compared with other existing methods in the literature.


Introduction
The movement of many individuals such as basic particles in physics or molecules can change under the effect of two important processes: chemical reactions in which the substances are transformed into each other, and diffusion mechanism which causes the substances to spread over a media.Such a process is governed by the following reaction-diffusion system where D is a matrix of diffusion coefficients, R is a vector valued function whose elements are the densities of the substances and represents reaction process.We are concerned with simulating the dynamics of single or many substances numerically in one dimension.The analytical and numerical treatments of (1) attract the attention of many researchers with diverse backgrounds.Numerical analysis of several finite differences (FD) schemes for reaction-diffusion systems is given in the research monograph [1].Apart from the classical approach, Ramirez et al. [2] derived nonstandard FD schemes for steady state reaction-diffusion equations based on Green's integral formulation of the exact solution.In this work, we will derive so-called exponentially fitted finite difference schemes that can be obtained by only changing denominators in discrete derivative formulas.

Exponentially Fitted Finite Difference Formulation
Finite difference discretizations of ODEs and PDEs are widely used by the practitioners due to their easy analysis and implementation issues.However, many drawbacks, such as violation of positivity requirements and strict conditions on step sizes are encountered in the applications.To overcome these difficulties, some variants of FD schemes are developed.The very well known modification to FD schemes is nonstandard finite differences (NSFD) proposed by Mickens [3].The rules of NSFD are based on complicated denominators and nonlocal representation of nonlinear terms in the equation.The other version of FD is exponentially fitted (EF) methods.In the last few decades, EF methods have become so popular that hundreds of papers have been published in theoretical and application-oriented journals.EF methods for differential equations date back to Gautschi [4], who proposed a method to integrate exactly appropriate trigonometric functions, just as classical methods integrate exactly algebraic polynomials.Paternoster [5] presents the overview of recent literature of EF methods.As for mathematical chemistry and physics, Anastassi and Simos [6] proposed EF methods for the numerical integration of the Schrödinger equation and related initial value problems with oscillating solution.The similarities and differences between NSFD and EF methods are discussed by Erdogan [7] and Gurski [8].Other versions of EF methods are employed to solve singular perturbation problems which are numerically challenging.A historical review of these types of methods is given in [9].Recent studies such as [10,11] present mathematical analysis and computational aspects of EF methods for singular perturbation problems.
EF methods include a fitting parameter, the so-called frequency that affects the convergence and stability properties of the method.For ordinary differential equations, the problem of how to choose fitting parameters are discussed in [12,13].However, the usage EF method in the full discretization of PDE's has not been well studied.In the recent work of the author [14], an exponentially fitted method is designed for solving advection type of initial boundary value problems by employing time and space frequencies in finite difference formulas.This approach does not only yield more accurate numerical results but also preserves the qualitative behavior of the equation under consideration.
In the construction of EF schemes, we will need the following EF based finite difference derivative formulas [15] that are exact for the reference sets {1, e ωx } and {1, e Ikx , e −Ikx } respectively for the first and second derivatives.In the next chapters, these derivative approximations are inserted into the differential equation, and strategies for parameter selections are discussed.

Pure Diffusion
For the sake of clarity, we start with considering the diffusion equation with initial condition u(x, 0) = u 0 (x) and periodic boundary conditions.The separation of variables method for (3) gives assuming the solution is of the form u(x, t) = T(t)X(x).It is clear that the time and spatial components of the solution are e ωt and Asin(kx) + Bcos(kx).The relation between time and spatial components is called dispersion relation [16].However, classical finite difference methods never take this issue into account explicitly.Our contribution is to consider the dispersion relation through the construction of finite difference schemes.Fortunately, exponentially fitted methods and nonstandard finite differences provide us with a very useful tool: denominators in terms of non-polynomial functions [17,18].We propose forward time derivative and the second space derivative in (2) to obtain an explicit EF scheme.The resulting scheme becomes The assignments of ω and k and the relation between these parameters are of great importance.Ehrhardt and Mickens [18] obtained a similar scheme by sub-equation method.However, their motivation and procedure to determine ω and k were quite different.The main contribution of the present work is to describe a parameter selection strategy based on dispersion relation and local error formula.The local error of the scheme ( 6) is given by where the derivatives (u t ) j i , (u x ) j i etc. are computed at the point (x i , t j ).The selections , that vanish the first two terms in the local error formula and increase the time and space order.Furthermore, one can realize that these selections reveal the dispersion relation ( 5) by considering In the implementation of ( 6), the classical finite difference approximations of u xxxx and u xx might be employed as follows for computing k numerically, These approximations do not vanish the local error terms exactly.Nevertheless, this approach reduces the local error significantly.It is clear that k is updated at each (x i , t j ).
k becomes constant within the time evolution and exact solution is obtained as long as initial profile is composed of only one Fourier component.However, in case of multi-frequency initial profile, k will not be constant any longer.It is also possible to take average values of parameters over time and space indices.

Stability of the Method
Both classical explicit finite difference scheme and EF scheme for the diffusion equation can be written as where stand for the CFL numbers for classical FD and EF schemes.
After the dispersion relation ω which tells us that classical FD formulation is employed in the implementation in case very small k values should be noticed.
Applying the well known stability arguments for time dependent PDEs [1,16], the stability requirement is found to be λ EF ≤ 1/2.After cumbersome algebra, it is verified that this requirement is satisfied when λ FD ≤ 1 2 under the assumption k ∈ [− π ∆x , π ∆x ], which prevents aliasing.In other words, stability of the method depends not only on step sizes but also on the parameters of the method.In order to overcome this drawback, we take (10) are very small or large in magnitude.Indeed, this precaution does not reduce the accuracy due to the fact that unwanted large and/or small values of parameters are caused by almost zero gradients that can be accurately resolved with classical FD formulas.

Reaction-Diffusion
In the presence of reaction term, the same derivative approximations are used to obtain the scheme that corresponds to the original PDE However, the dispersion relation (5) does not hold any longer.Considering (8), one computes k and ω as For a coupled reaction-diffusion system, it is suggestive to use different space and time parameters for each of substances.Consider the system with appropriate boundary and initial condition where D is 2 × 2 diffusion matrix.An explicit EF scheme for the coupled system can be written as Apart from employing constant parameters in some special cases, local truncation errors for each line suggest that the selections can decrease the error significantly.It is very usual to have stiffness in coupled reaction-diffusion systems.Explicit EF also fails in the case of stiffness.Therefore, an implicit version of EF finite differences should be used to deal with stiffness.By employing backward time difference in (2), one obtains an implicit version of (17).Local error based parameter selection strategy gives Let us point out that parameter computations are performed explicitly although the proposed method is fully implicit.

Consider the diffusion equation with periodic boundary conditions on
The initial condition u(x, 0) = sin(x) leads to the exact solution u(x, t) = e −t sin(x).In Table 1, the EF method with the settings k = 1 and ω = −1, which will be denoted by EF constant , gives an exact solution up to machine accuracy.However, EF method performed by frequency search algorithm explained in the previous section could not catch the exact solution due to fluctuations in k j i .Nevertheless, it is superior to the classical explicit finite difference scheme

Nonlinear Diffusion
The proposed method is applied to nonlinear diffusion directly.Consider the problem with exact solution u(x, t) = x 2 10 − 6t [19].The discretization of second derivative will be same as the linear case The parameters k and ω are determined in terms of derivatives of u 2 2 not u itself.In Table 2, comparison of errors produced by EF and classical explicit finite difference methods is given.EF needs very moderate step sizes to get higher accuracy.

Fisher Equation
In this example, we will see that the proposed EF approach is not limited to linear equations with single dominant frequency.Consider the Fisher equation with exact solution EF solution is computed by Equations ( 15) and ( 13).Table 3 presents numerical solutions at t = 4 and x ∈ [−30, 30] with ∆x = 1, ∆t = 0.1 and for a = 1/24 , b = c = 1 .The wavelet method [20], which is a kind of spectral method, and B spline technique [21], that yields five points implicit stencil, are used for numerical comparison.Numerical time evolution of the initial profile by EF method is given in Figure 1.

Coupled Nonlinear System
Consider the equation . Its exact travelling wave solution [22] is given by u(x, t) = e k(x+ct) where k = √ 6 3 and c = √ 6 6 .This problem is used to test some semi analytical approximation methods such as Adomian Decomposition Method [23].
A numerical experiment is performed for x ∈ [−30, 30] and t ∈ [0, 10].Explicit EF and classical FD schemes (central difference in space, forward difference in time) are employed.In Figure 2, the order plots are given.In Figure 2a, ∆x = 0.25 is fixed and several time steps within the stability region are considered.In Figure 2b, ∆t = 0.001 is fixed and numerical experiments are performed for different ∆x values.Figures indicate that the proposed parameter selection algorithm could not make leading order terms vanish in local error.This is because the exact solution does not consist of trigonometric space functions.However, the parameter selection strategy reduces the local errors considerably at each (x i , t j ).The classical FD scheme can catch the accuracy that is obtained by EF only for very small ∆t and ∆x values.

Stiff Coupled System
Consider the initial-boundary value problem πx) and homogenous Neumann conditions.Its exact solution is u(x, t) = e −t sin 2 (πx) and v(x, t) = e −t cos 2 (πx).This example is used to test some implicit solvers such as those proposed in [24,25].In this example, Implicit EF method with parameter selections in Equation ( 19) is compared with other implicit solvers.Table 4 presents errors in max norm for a fixed time step and various ∆x values.Although B spline method is based on five point stencil and employs Crank-Nicolson for time discretization, it falls behind Implicit EF.
Table 4. Errors for u(x, 1) and v(x, 1) in max norm at final time t = 1 with ∆t = 0.001 for example 5.

Conclusions
A simple modification of a classical explicit finite difference scheme, which is based on change of denominators, is described and convergence properties of the proposed scheme are discussed.This algorithmic approach has improved the accuracy without considerable computational work.A reader who is familiar with spectral methods [26] can realize that we are trying to include only dominant frequency into the finite difference formulation.The success of the proposed method for diffusion equations relies on the smoothing effect of the diffusion operator that damps high frequencies rapidly.
In a similar manner, all kinds of classical FD schemes can be converted to EF schemes by choosing a suitable reference set.A further study can be performed on the parameter selection strategies, possibly with the instruments of Fourier analysis.

Figure 2 .
Figure 2. (a) log log plot of ∆t values vs. max norm of errors at final time t = 10 for fixed ∆x; (b) log log plot of ∆x values vs. max norm of errors at final time t = 10 for fixed ∆t.

Table 2 .
Errors in max norm at final time t = 1 for example 2.

Table 3 .
Comparison of the methods for example 3.