Dynamically Adapted Mesh Construction for the Efficient Numerical Solution of a Singular Perturbed Reaction-diffusion-advection Equation

This work develops a theory of the asymptotic-numerical investigation of the moving fronts in reaction-diffusion-advection models. By considering the numerical solution of the singularly perturbed Burgers’s equation we discuss a method of dynamically adapted mesh construction that is able to significantly improve the numerical solution of this type of equations. For the construction we use a priori information that is based on the asymptotic analysis of the problem. In particular, we take into account the information about the speed of the transition layer, its width and structure. Our algorithms are able to reduce significantly complexity and enhance stability of the numerical calculations in comparison with classical approaches for solving this class of problems. The numerical experiment is presented to demonstrate the effectiveness of the proposed method. The article is published in the authors’ wording.


Introduction
This work is concerned with analytic-numerical methods for singularly perturbed parabolic equations which may be interpreted as models of reaction-diffusion-advection processes in various applications.Such problems often feature narrow boundary and interior layers (stationary or moving fronts), and are extremely difficult for a numerical treatment.Recently, asymptotic analysis is successfully used to investigate such problems (see [8,15,14,13,18,16,17,10,11,12], where special mesh methods used).Note that the special numerical methods for singularly perturbed parabolic equations were practically studied mostly for problems with stationary boundary and interior layers.There are some specific features in constructing special difference schemes in the case of moving interior layers.For singularly perturbed parabolic problems with such type of solutions, special schemes were constructed in [10,6,7,5].Such schemes are fairly complicated, that calls the necessity to develop simpler schemes and alternative numerical methods (based on a posteriori adapted grids).In the presented paper we propose an effective combined analytic-numerical approach, which use some a priori information from asymptotic analysis of the moving front type solution, in order to simplify numerical calculations.Numerical methods for nonlinear singularly perturbed interior layer problems using an approximate layer location are developed, particularly, in the works [10,9] (see, also, the references therein), where different type of refined meshes were proposed.
The motivation to combine asymptotic and numerical methods is that if we use numerical method for a singularly perturbed problem, two opposite phenomena connected with small parameter arise: on the one hand, the smaller this parameter, the more rough and unstable numerical solution we obtain; on the other hand, the smaller this parameter, the more precise a priori information about exact solution we are able to get by using asymptotic methods.This fact gives the possibility for a productive combination of asymptotic and numerical approaches.Another reason is that using asymptotic analysis of singularly perturbed problem we can prove the existence of the exact solution with boundary or internal layers (moving fronts) and can decrease the spatial dimension for numerical calculations: the dimension of spatial variables in the equation for the moving front location is lower per unit then the original problem or this equation is not partial, but ODE (or it is not differential equation at all).
The main purpose of this paper is to suggest an algorithm of numerical solution of the singularly perturbed problems with the moving internal layers (fronts) which based on the asymptotic analysis of the problem.Using the rigorous asymptotic analysis, which also state the existence of the moving front type solution and gives the asymptotic expansion of the front location, we construct so called dynamically adapted mesh (DAM).
Our DAM is constructed by the following way.We use the basic (coarse) mesh with the mesh interval choosing in accordance with the width of the interior layer.The information about the location and width of the interior layer (front) we get by the asymptotic analysis.Inside the moving front region two course mesh intervals are divided into some additional subintervals (so inside the moving front region the mesh is a version of a Shishkin mesh).In order to justify the number of these additional intervals some method of the a posteriori estimates (e.g.Richardson extrapolation) can be applied.The coarse mesh (and, therefore, the front region intervals) is chosen to guarantee that the front is located inside these intervals.
Note, that our variant of the coarse mesh is ε-dependent.For the problems with moving internal layers this choice of the coarse mesh simplifies the interpolation procedure for the moving refined mesh and allows to save computing resources.Of course, all presented ideas can be realized on the mesh when the total number of nodes (coarse and refined meshes together) not depends of ε (Shishkin mesh), but in this case nodes of the coarse mesh are not fixed because of refined mesh is moving with the front.In this paper we did not investigate a priori the convergence of the numerical method, but our variant of the mesh choice allows to simplify the procedure of a posteriori error estimates because not requires to do interpolation to the new coarse mesh nodes on each time step.
We also have to note that the problem to define the layer location in our example is explicitly solvable.In more general case we have to create numerical method to solve this problem which is done at the end of Section 1..The paper is structured as follows.In Section 1. we describe methods which allow to get a priori information that will be used for the DAM constructing.In Section 2. we briefly describe the main ideas that are used for constructing the DAM.In Section 3. we perform a numerical experiment for the example and explain some nuances of the mesh constructing.Finally, in Section 4. we discuss some perspectives of our approach.

A priori information from the asymptotic analysis
We illustrate our approach by considering the following class of problems: where functions A(u, x) and B(u, x) are sufficiently smooth for (x, u) ∈ [0; 1]×(−∞; +∞).The rigorous treatment of problem (1) was obtained in [3], where the existence and asymptotics of the front type solution were established.The main information, which we use for our approach, is: 1. speed or position of the transition layer; 2. width of the transition layer; 3. structure of the transition layer (we use the fact that transition layer has the exponential behavior).
We will obtain this information using the asymptotic of the solution of problem (1).The main ideas of the asymptotic procedure are briefly represented below.
Let consider two initial value problems: We suppose, that the solutions of problems (2) exist on x ∈ [0, 1], and denote them as ϕ l (x) and ϕ r (x) respectively.Assume, that the following inequalities are fulfilled for all x ∈ [0, 1]: The problem of formation of the moving fronts is discussed in [2] and [1].Particularly, for the reaction-diffusion problem the fast formation stage is described in [2].We assume that the front is already formed at the moment t = 0 and the initial function u init (x) has a thin transition layer between levels ϕ l (x) and ϕ r (x) located at the point x 00 ∈ [0, 1].
For our problem this assumption can model a discontinuous-initial-condition situation.
We define the moving front position by the function x t.p. (t, ε), that is the point of the intersection of the solution of (1) u(x, t, ε) and the level 1 2 (ϕ l (x) + ϕ r (x)).
We put x t.p. (0, ε) = x 00 and will seek the position and speed of the transition point in the form of power series of ε where has solution x(t) such that and the following inequalities are fulfilled a) ϕ r (x) where Condition (6), a) guarantee that (1) has no stationary solutions, and ( 6), b) provides solvability of some equations in the further asymptotic procedure.
We define by Dl T and Dr T the left and right sides with respect to the point x t.p. (t, ε) respectively and will find the solution of (1) separately in these parts The functions u l and u r can be written in the form Here ūl,r (x, ε) -regular functions, which represent the solution far from transition point x t.p. (t, ε); function Q l,r (ξ, t, ε), where ξ = (x − x t.p. (t, ε))/ε, describe transition layer (moving front) near this point (ξ ≤ 0 for the function with the index l and ξ ≥ 0 for the function with the index r).
The functions ūl,r (x, ε) and Q l,r (ξ, t, ε) we seek as power series of ε: Terms of series ( 9) can be built by the standard procedure (see [3]) separately at two sides of the point x t.p. (t, ε).We assume that the functions u l and u r are joint with the continuous first derivatives (C (1) -matching conditions): Substituting ( 9), a) into (1), at zero order of ε we obtain equations (2).So we have ūl 0 (x) = ϕ l (x), ūr 0 (x) = ϕ r (x) and Functions ūl,r k (x), k ≥ 1 are the solutions of the following problems where Āl,r (x) = A(ϕ l,r (x), x), Bl,r (x) = B(ϕ l,r (x), x) and f l,r k (x) is known function.Solutions of (12) can be written explicitly.
In order to obtain the transition layer functions (9), b) we must do the change of variables (x, t) by (ξ, t) in (1), where ξ = (x Substituting these expressions and ( 9), b) into (1) and equating terms with different powers of ε, we obtain equations for functions Q l,r i (ξ, t) , i = 1, 2, . . .For Q l,r 0 (ξ, t) we have the following problem: Using the continuous function we can rewrite (13) as The exponential estimates for the solutions Q l,r 0 (ξ, t) of ( 15) were established in [3]: where C and κ are positive constants.
Using C (1) -matching conditions and explicit formulas for ∂ Q ∂ξ we obtain the problem (5) for zero order term of the moving front position Note that for general cases the initial value problem (17) has to be solved numerically.In the example, which we use at Section 3., we can solve it explicitly.
The solutions of ( 18) can be found explicitly and for all t ∈ [0, T ] the functions where C and κ are some positive constants.Using C (1) -matching conditions (10) at first order of ε and explicit formulas for Q l,r 1 (ξ, t) and , we obtain the problem for first order term of the moving front position where K(x 0 (t)) and G 1 (x 0 (t), t) are known functions (see [3]).First two orders in ε of moving front speed v 0 (t) + εv 1 (t) are defined in (17) and (20).So the location of the moving front at two higher orders of asymptotics is defined as the solutions of Cauchy problems (17) and (20).
The width of the transition layer is defined by the estimates (16).It means that interior layer exponentially tends to functions ϕ r,l (x) and its thickness is The structure of the interior layer is defined by the problem (15) and by estimates (16).
In Section 3. we demonstrate our approach by the following particular case of problem (1), where A(u, x) = −u, B(u, x) = u: For this example we have: and the conditions (3) and ( 6) are satisfied for all x ∈ [0, 1]: From (17) for the zero order term of the moving front speed we have From (20) for the first order term we get and therefore x 1 (t) ≡ 0.
If we are not able to calculate position of the transition point analytically (that is possible in some cases of particular A(u, x) and B(u, x)),we are able to do this numerically.Here we want to describe corresponding numerical algorithm of x t.p. determination.
At the beginning we introduce uniform mesh X N 0 only on x-dimension that has number of nodes N 0 + 1 (that equals to N 0 intervals): . This mesh will be used for calculation of the solution q(x).Then it is possible to solve equations ( 2) and find numerically ϕ l (x) and ϕ r (x) on the mesh X N 0 .For numerical solving of (2) any suitable scheme can be applied.Then right-hand side of the equation ( 5) can be calculated on the same mesh X N 0 .
After that we have to solve the problem (5).At first, we introduce uniform mesh T M 0 only on t-dimension that has number of nodes M 0 + 1 (that equals to M 0 intervals): We suggest to use Rosenbrock scheme with complex coefficient (CROS1), which is monotone, stable and has the order of accuracy O(τ 2 ) (see [4]): .
Note: It is very important to mention that we have to perform interpolation of the tabulated function V (x) from basic mesh X N 0 on each point x(t m ).It is need not to perform interpolation with order of accuracy greater than 3. So, we have dependence of transition point position x on time t as the solution of (5) (see Figure 1(a)).After that we are able to interpolate corresponding function t(x) on the uniform mesh X, which has intervals equal to |ε log ε|/C, C > 1 (where |ε log ε| is thickness of the interior layer).This interpolation is possible to perform because of strongly monotonic of the function x(t) (see Figure 1(b)).So, we have grid T M which is quasi-uniform grid [20].

Dynamically adapted mesh construction
The idea of dynamically adapted mesh construction is the following.If we know the width of the transition layer, we can introduce basic uniform mesh with the steps, which are equal to this width, and then refine two intervals that are the nearest to the transition point (see Figure 2-1 0 ),1)).Then if we know the location of the transition layer on each time step, we can check whether the transition layer is located on these intervals or not.If the transition layer will leave the second interval soon, we refine the next basic interval and perform interpolation of the numerical solution on these additional nodes (see Figure 2-2)).Then we remove from following calculations the nodes of the first refined interval.For appropriate interpolation we should use information about The main problem of this process is the possibility to obtain corresponding a priori information.This problem for one type of reaction-diffusion-advection equation was discussed above in Section 1.. Now, in Section 3., we will describe some numerical experiment and explain some nuances of the mesh constructing for the particular example.

Numerical example
As an example of the application of proposed methods we consider the following Burgers's equation: where u lef t = −5, u right = 2, and the initial condition u init (x) has following form (for example, see Figure 4): Method of lines and Rosenbrock scheme with complex coefficient For numerical solution of the equation (25) we apply the stiff method of lines (SMOL) in order to reduce the PDE to the system of ODEs that can be solved by Rosenbrock scheme with complex coefficient, which is significantly efficient for stiff systems of ODEs [4].
At the beginning we introduce piecewise uniform mesh X N only on x-dimension that has number of nodes N + 1 (that equals to N intervals): . So after finite-difference approximations of derivatives with second order of accuracy in (25) we obtain the following system of This system can be rewritten as where T .The vector-function f has the following structure.For n = 1: for n = 2, N − 2: For numerical solution of this system of ODEs (26) we use Rosenbrock scheme with complex coefficient (CROS1), which is monotone and stable and has the order of accuracy O(τ 2 ) (see [4]).In order to apply this scheme we introduce quasi-uniform mesh T M on t-dimension that has number of nodes M +1 (that equal to M intervals): After that we are able to apply the CROS1 scheme for solving of the system (26): Here E is the identity matrix, f u is the Jacobian matrix, where for n = 1: and for n = N − 1: The other components of f u for the considered equation equal to zero.
Dynamically adapted mesh construction Now we explain in details how to construct the dynamically adapted mesh X N (t) ≡ {X N (t m )}, 0 m M , and how to use it for the calculations by the scheme (27). 1) At first we introduce a basic uniform mesh X N 0 on x with step h 0 = (1 − 0)/N 0 that has number of nodes N 0 + 1 (that corresponds to N 0 intervals): X (0) The number of interval N 0 should be chosen in accordance to the a priori information about thickness of the transition layer that can be calculated by formula (21): 2) Then we introduce the family of piece-wise uniform meshes {X n,n+1 N 0 , N int }, 0 n N 0 − 2, which have the uniform refining on the (n + 1)-th and (n + 2)-th intervals of the basic mesh using N int additional intervals and have totally N = N 0 − 2 + 2N int intervals.possibility for productive combinations of asymptotic and numerical approaches.Using the information, based on the rigorous asymptotic analysis of the problem, we propose an analytic-numerical algorithm for a singularly perturbed reaction-diffusion-advection equations that reduces complexity of the numerical calculations.The class of the problems that was considered in this work is given to illustrate our approach.Our method is not restricted by only this class of the problems.We plan to extend this approach for periodic-parabolic problems with interior layers solutions and some classes of systems.We also plan to improve our procedure in order to use ε-independent meshes and investigate the convergence of the numerical method.

Fig 1
Fig 1. a) Function x(t) as a solution of (5).b) Quasiunidform mesh T M on t constructing

Fig 4 .
Fig 4. The example of u init (x) for ε = 10 −2 (refining of the mesh in a neighborhood of the transition point has been performed)