High Order Schemes for Reaction-Diffusion Singularly Perturbed Systems

In this paper we are interested in solving efficiently a singularly perturbed linear system of differential equations of reaction-diffusion type. Firstly, a non-monotone finite difference scheme of HODIE type is constructed on a Shishkinmesh. The previous method is modified at the transition points such that an inverse monotone scheme is obtained.We prove that if the diffusion parameters are equal it is a third order uniformly convergent method. If the diffusion parameters are different some numerical evidence is presented to suggest that an uniformly convergent scheme of order greater than two is obtained. Nevertheless, the uniform errors are bigger and the orders of uniform convergence are less than in the case corresponding to equal diffusion parameters.


Introduction
In this work we consider the singularly perturbed boundary value problem given by the linear reaction-diffusion system L ε u = f , x ∈ Ω = (0, 1), u(0) = u 0 , u(1) = u 1 , where the differential operator L ε is defined by a 11 (x) a 12 (x) a 21 (x) a 22 (x) .
We will assume that the diffusion parameters 0 < ε 1 ≤ ε 2 ≤ 1, can take arbitrary small values having, in general, different order of magnitude, that the data of problem (1) are sufficiently smooth functions and also that the coefficients of the coupling reaction term satisfy i.e., the reaction matrix is an M -matrix. First order uniform convergence of the central finite difference scheme constructed on a Shishkin mesh was proved in [7] and in [4] this was improved to almost second order. Linß and Madden [6] extended this result to the case of an arbitrary number of equations, when the reaction coefficient matrix A satisfy another type of conditions, which include these ones given in (2) for the case that the coupled system has only two equations as problem (1) here considered. Also, in [3] precise information of the asymptotic nature of the solution and its derivatives, for a problem having n equations with n diffusion parameters, has been recently established by means of an appropriate decomposition of the solution, revealing that the solution exhibits overlapping boundary layers with a width O(ε −1/2 i ), i = 1, 2, . . . , n at both endpoints x = 0 and x = 1. It was also proved that the central finite difference scheme constructed on a piecewise uniform mesh of Shishkin type, is first order uniformly convergent in the maximum norm.
High order convergence schemes are very interesting in practice because they provide accurate numerical approximations with a low computational cost. Nevertheless, at the moment, in the literature there are not numerical methods for problem (1) with this desirable property. The aim of this work is to see how the HODIE technique permits to obtain a uniformly convergent method having order bigger than two. In some cases the proof of the uniform convergence of the method is fulfilled, but in general we only dispose of numerical evidences showing the efficiency of the HODIE method.
Henceforth, C denotes a generic positive constant independent of the diffusion parameters, and also of the discretization parameter.

The numerical method
To construct the numerical method we first define a piecewise uniform Shishkin mesh. Following [7] the mesh points are and σ 0 ≥ 4. If τ ε1 = 1/8 and τ ε2 = 1/4, we modify slightly the mesh points; now they are given by . Below we denote the local step sizes by h j = x j − x j−1 , j = 1, . . . , N . On this mesh we impose that the local error be zero on the set of vector polynomials of the form Following the construction made in [1] for the scalar case, we write the finite difference scheme L N = (L N 1 , L N 2 ) T in the form Then, it is not difficult to prove that for j = 1, · · · , N − 1, i = 1, 2 the coefficients r s of the scheme are given, in function of the coefficients q s, by and also that it holds The value of the free parameter q 3 i,j is taken equal to the one obtained for the scalar case in [1] and we will see that this choice is also appropriate for the case of systems. This value depends on the location of the mesh points and also on the ratio between the step sizes of the Shishkin mesh and the diffusion parameters. Concretely, for j = 1, · · · , N/8 − 1, 7N/8 + 1, · · · , N − 1, i.e., x j ∈ (0, τ ε1 ) (1 − τ ε1 , 1) and i = 1, 2, the coefficients q s are For j = N/4, · · · , 3N/4, i.e., x j ∈ [τ ε2 , 1 − τ ε2 ], and i = 1, 2 we distinguish two cases: first, if 2H 2 a ii ∞ /3 ≤ ε i , then the coefficients are defined again by (6); in the other case, when 2H 2 a ii ∞ /3 > ε i , they are given by corresponding to the classical discretization of central differences. Note that in this case (5) does not hold. Last case is when j = N/8, ]. Now, for the second equation, i = 2, the coefficients are again given by (6). Nevertheless for the first equation, i = 1, again we must distinguish two cases; first, when 2h 2 ε2 a 11 ∞ /3 ≤ ε 1 , the coefficients are given by (6); in the other case, 2h 2 ε2 a 11 ∞ /3 > ε 1 , they are given by (7). Note that, in general, the coefficients defined in (6) are not positive and then the associated matrix to the numerical scheme is not an M-matrix. Nevertheless, we will see the efficiency of this scheme. As an example, we solve the particular problem (see [7]) setting by with u 0 = u 1 = 0. For this problem the exact solution is unknown and therefore to approximate the pointwise errors |(U − u)(x j )|, j = 0, · · · , N , we use a variant of the double mesh principle. So, we calculate a numerical approximation U to u given by the scheme (3) on the mesh {x j } that contains the mesh points of the original piecewise Shishkin mesh and their midpoints, i.e., the mesh points are defined byx 2j = x j , j = 0, . . . , N,x 2j+1 = (x j + x j+1 )/2, j = 0, . . . , N − 1. Then, at the original mesh points x j , j = 0, 1, · · · , N , the maximum errors and the uniform errors are approximated by where, in order to permit the stabilization of the errors, we take S as the set From these estimates of the pointwise errors we obtain the corresponding orders of convergence and the uniform orders of convergence in a standard way, by using In all cases we take the constant σ 0 = 4; in practice if this constant is smaller, the desired order of uniform convergence is not achieved. On the other hand, if it is greater than 4 the numerical errors are bigger but the order of uniform convergence is preserved. Table 1 displays the results obtained with the HODIE scheme; from these results we observe that the order of uniform convergence is four except by a logarithmic factor, as it is usual on Shishkin meshes. Nevertheless, the discrete operator of this scheme is not of positive type and we do not have the proof of the uniform (l ∞ , l ∞ )-stability. In [6] this uniform stability was proved without using the inverse monotonicity of the discrete operator, but unfortunately so far we have not been able to apply this technique to the HODIE operator. In [5] a non-monotone FEM scheme was used to solve a scalar reaction-diffusion problem, proving also its uniform stability in the maximum norm. Therefore we propose a slight modification of scheme (3) to have a new scheme satisfying the discrete maximum principle. We clearly see that only the discretization associated with the transition points does not give the correct coefficients sign pattern to have an M-matrix. Then, we change the discretization corresponding to the indexes j = N/8, N/4, 3N/4, 7N/8, such that r − i,j < 0, r + i,j < 0, r − i,j +r + i,j +r c i,j > 0 and q 1 i,j , q 2 i,j , q 3 i,j be positive. It is straightforward to obtain that the coefficients q s are given by It is easy to proof that the discrete operator is of positive type and therefore it satisfies the discrete maximum principle.

The case of equal diffusion parameters
To find a theoretical proof of the uniform convergence of the method, we begin with the case where both diffusion parameters take the same value, ε 1 = ε 1 = ε. Note that in this case really there are only two transition points in the Shishkin mesh and the transition parameter is defined by τ = min 1/4, σ 0 √ ε 2 ln N . Following the idea of extending the domain introduced by Shishkin in [8], which was also used in [2] to find a decomposition of the exact solution of a two dimensional scalar equation of reaction-diffusion type, it can be proved the following result showing the asymptotic behavior of the exact solution.
Lemma 1. Let assume a ij , f ∈ C 4 (Ω), i, j = 1, 2. Then, for ε 1 = ε 2 = ε, the exact solution of (1) can be decomposed as and w Note that we have appropriate bounds of the regular and singular components and their derivatives up to sixth order, which we will need in the analysis of the truncation error.
Theorem 1. Let u be the solution of continuous problems (1) and U the solution of the discrete operator (3)- (7) and (10) defined on the previous Shishkin mesh, when ε 1 = ε 2 = ε. Then, the error satisfies Proof. In the case σ = 1/4, using that ε −1/2 ≤ C ln N and the crude bounds u (i) ∞ ≤ Cε −i/2 , i = 0, · · · , 6, a classical analysis proves that U − u ∞ ≤ CN −4 ln 4 N . When σ < 1/4, first we study the error for the regular component. Then, if ε ≤ 2H 2 min i a ii ∞ /3, taking Taylor expansions the local error can be bounded as in [1] for a single equation, and therefore we have Then, the discrete maximum principle proves that On the other hand, if ε > 2H 2 min i a ii ∞ /3, we can obtain At the transition points, using that for any z ∈ C 4 (Ω) it holds we deduce Defining the barrier function Z( where θ is the piecewise linear function using the maximum principle, it can be proven that and taking into account that ε 1/2 ln N ≤ 1, it follows For the singular component we distinguish two cases depending on the location of the mesh point. For x j ∈ [τ, 1 − τ ], using the exponential character of this component, it is not difficult to deduce In the second case, x j ∈ (0, τ ) ∪ (1 − τ, 1), the local error is bounded by Applying again the maximum principle, now on From (13)-(15) the result follows.
For the same example as before, with ε = 2 0 , 2 −2 , . . . , 2 −50 , Table 2 displays the results obtained; from it we clearly observe that the order of uniform convergence is similar to that for the unmodified HODIE scheme. Note that the numerical results indicate an order of uniform convergence higher than this one proven in Theorem 1. In the general case, when the diffusion parameters can be different, the theoretical question is more complicated. An important question is related with the decomposition of the exact solution. In this case it is possible to find a decomposition into a regular and singular part (see [4,6,7] for instance), but it is not clear how it is possible to obtain the bounds (11) for the regular component of the solution; note that we need the bounds of the derivatives up to sixth order, to find appropriate bounds for the local error associate to the scheme. Nevertheless, for us it is interesting to confirm in practice that this scheme gives an order of uniform convergence bigger than two. Table 3 displays the results obtained with the new scheme when the diffusion parameters are not equal. From this table we observe that the method gives an almost third order uniformly convergent method, which is less than the order obtained in the case of equal diffusion parameters. Table 3. Uniform errors and orders of convergence for the modified scheme given by (3)-(7)  Nevertheless, this method improves both the maximum errors and the numerical order of uniform convergence with respect to central finite difference scheme. Because the modified finite difference scheme satisfies the maximum principle, having appropriate bounds for the derivatives of the regular and singular part of the solution, would allow us carry out the analysis of the local error, and therefore prove the desired uniform convergence.