A Moving Mesh Method for Singularly Perturbed Problems

and Applied Analysis 3 Proof. Consider the barrier functions ψ ± (x, t) = 󵄩 󵄩 󵄩 󵄩 g 󵄩 󵄩 󵄩 󵄩 + x 󵄩 󵄩 󵄩 󵄩 f 󵄩 󵄩 󵄩 󵄩


Introduction
The use of adapted meshes [1][2][3] in the numerical solution of differential equations has become a popular technique for improving existing approximation schemes.When considering an adaptive mesh algorithm for the solution of time dependent differential equations [4][5][6], the techniques which underpin the grid movement are often found in the literature [4,6,7] for the generation of adapted grids for the numerical solution of steady problems.One such technique is equidistribution, first introduced by de Boor [8], involving locating mesh points such that some measure of the solution geometry or error is equalized over each subinterval; a typical example is redistributing the arc length of the solution.To approximate the solution accurately in these regions, it is necessary to generate a mesh that is dense where the solution is changing rapidly and to remove unneeded points from regions where the solution is becoming smoother.Thus the mesh must have a dynamic behaviour in much the same way as the solution.This problem has been addressed by Huang et al. [9], who proposed a general adaptive mesh method known as the moving mesh method.
For the moving mesh methods, the number of grid points is fixed.The mesh points move continuously in the space time domain and concentrate in regions where the solution is steep.The movement of the mesh is governed by a mesh equation which moves the mesh around in an orderly fashion.Huang et al. [9] developed several forms of the mesh equations known as the moving mesh partial differential equations (MMPDEs).Here a simple equidistribution relation in one spatial dimension is differentiated with respect to time in order to derive equations prescribing the correct velocities of nodes in order to preserve the equidistribution principle as the solution and grid evolve.The mesh equation and the original differential equation can be solved simultaneously or decoupled to get the physical solution and mesh.How the MMPDEs are formulated and solved [10,11] is crucial to the efficiency and robustness of the method.Zhou et al. [12] applied a difference scheme to a singularly perturbed problem.The study used two algorithms on moving mesh methods by using Richardson extrapolation which can improve the accuracy of numerical solution.Yang [13] considered a kind of nonconservative singularly perturbed two-point value problems in fluid dynamics.Cen [14] examined a class of delay differential equations with a perturbation parameter .More recently, Gowrisankar and Natesan [5] numerically studied singularly perturbed parabolic convection-diffusion problems exhibiting regular boundary layers.
In order to obtain a robust moving mesh method which can solve a wide range of problems, we are going to adopt the equidistribution principle moving mesh strategy with the arc length as the monitor function.
In general the solution   will be smooth on Ω for all values of .Boundary and interior layers [15][16][17] are normally present in the solutions of problems involving such equations.These layers are thin regions in the domain where the gradient of the solution steepens as the singular perturbation parameter  tends to zero.In problems, in which large solution variations are common, the choice of a nonuniform mesh cannot only retain the accuracy but also improve the efficiency of an existing method by concentrating mesh points in regions of interest.If the regions of high spatial activity are moving in time, then techniques that also adapt the grid in time are needed.The moving mesh method [6] will be used to solve   .The drawback of this strategy is that, with the introduction of the mesh equations which govern mesh movement, the system becomes nonlinear for any linear problem; hence very little theoretical analysis [1,7,18,19] has been carried out to explain the convergence behaviour of the method.The following assumptions will be made: for all (, ) ∈ Ω, ‖()/‖ ≤ , for some constant  and at  = T, ‖  (, )/‖ ∼ () and ‖()/‖ ∼ ().

The Continuous Problem
The differential operator   for   satisfies the following maximum principle.
An immediate consequence of this is the following bound on the solution of any problem from   .

Lemma 2. Let
The maximum principle now applies, and we have  ± (, ) ≥ 0 for all (, ) ∈ Ω from which we have the required result.
Equation ( 15) can also be used to give the following bound: for all (, ) ∈ Ω.This proves the result for  = 1.To obtain bounds for the higher derivatives, rewriting ( 10) this gives the second derivative bound where  depends on ‖‖, ‖‖, and and using the idea from (16), this gives the bound for the third derivative.
Consider the following decomposition of the solution into the smooth and singular components: where V  (, ) is the solution to problem and the singular component   (, ) is the solution of the homogenous problem Theorem 4. The solution   (, ) of the continuous problem   can be decomposed as a sum of the smooth and layer functions where for all , 0 ≤  ≤ 3, the smooth component V  (, ) satisfies for some constant  independent of .
Proof.To find these bounds, we rewrite the smooth component as where V 0 (, ) is the solution to the reduced problem and V 2 (, ) satisfies We clearly have   V  (, ) = () in Ω with V  (, ) = V 0 (, ) + V 1 (, ) on Γ  .It can be easily seen that V 0 (, ) and V 1 (, ) are all bounded by a constant independent of  and for 0 ≤  + 2 ≤ 4. Therefore V 2 (, ) is the solution of a problem similar to   ; hence from Lemma 2, we obtain that To get the bounds for the spatial derivatives, we only consider V 2 (, ) since V 0 (, ) and V 1 (, ) are independent of .As previously stated that the problem for V 2 is similar to the problem for   (, ), we can use Lemma 3 from which we obtain for 0 ≤  ≤ 3 To find bounds for the singular component, we consider the following mesh functions: where Applying the differential operator   , we obtain that Thus from the maximum principle, we can say that  ± (, ) ≥ 0; hence To establish the bounds for the derivatives of   (, ), we start by noting that          ∫ which yields the required estimates for the second and third derivatives for the singular component   .

The Discretized Problem
In this section the discretization process for   is considered.By changing the time derivative into the Lagrangian form, this enables the introduction of node velocities into the system.Setting A = (/) − ,   can be written as Discretizing (49) using an implicit scheme, which can be written as a system where is a user defined parameter, which determines how fast the grid moves towards the equidistributed grid, and  is the arc length monitor function First it will be shown that the solution from (50) is bounded.An inductive proof is used to show that such constants exist and are indeed finite.At  =  0 = 0, we have a uniform mesh.The initial condition  ,0 is assumed not to be identically zero over the whole domain and to be bounded max  | ,0 | ≤  by some constant 0 <  < ∞.The system for the MMPDE can be written in such a way that we get an -matrix on the left hand side.This means that a solution for the system exists, and we can assume that max where for all 1 ≤  ≤  − 1.The discrete differential operator    in (57) satisfies the following discrete maximum principle.
Theorem 5 (discrete maximum principle).Let  , be any mesh function defined on Ω  .If  , ≥ 0 for all (  ,   ) ∈ Γ  and     , ≥ 0 for all (  ,   ) ∈ Ω  , then  , ≥ 0 for all Proof.Assume that there exists r = (  ,   ) ∈ Ω  such that then r ∉ Γ  , which implies that r ∈ Ω  , and we also know that We have to show that  −  ( ,  , ) ≤ 0, and we proceed by contradiction, suppose that Since r is an arbitrary point, we have that  0,  0, < 0, but  0,  0, ≥ 0 which is a contradiction.Our supposition that  −  ( ,  , ) > 0 is false, so we have that  −  ( ,  , ) ≤ 0. All the arguments given above imply that     , ≤ 0. Therefore we must have     , = 0, but we have that − 2   , ≤ 0 and  −  ( ,  , ) ≤ 0; hence Using the same argument as before, this implies that  0, < 0, but  0, ≥ 0 which is a contradiction.So     , < 0, which is a contradiction.Thus our original assumption must be false, and we conclude that the minimum of  , is nonnegative.

Decomposition of Numerical Solution and Error Estimates
Let   denote the discrete solution, and assume that the discrete solution can be decomposed into the sum where   (  ,   ) and   (  ,   ) are solutions of the respective equations where   is the singular component and   is the smooth component.The error for the numerical solution will be decomposed as The error (  −   )(  ,   ) can now be estimated using the error estimates for the singular and smooth components of the solution.
Proof.We start by considering the local truncation error for the smooth component (73) We introduce the two mesh functions it can be seen from the mesh functions that for some constant  independent of  and .
Proof.We need to consider two separate cases since the role of the boundary layer is crucial.We start with the case when where  max = max{‖‖, ‖  ‖, ‖  ‖}.Consider the following barrier function: and the mesh functions It can be easily seen that Applying the difference operator    to the barrier function, Hence it follows that     ± , ≥ 0 for 1 ≤  ≤  − 1; by the discrete maximum principle, it follows that  ± , ≥ 0; hence To consider the case when  −1 > 2/ ln 1/ and  <  −1 , we will assume that ℎ , ≤ ℎ −1, and that there exists a point   ∈ [0, 1] such that 1−  ≥ / ln 1/.The point   splits the interval into the coarse mesh Ω   and fine mesh Ω   where We give separate proofs in the coarse and fine mesh subintervals.First suppose that   ∈ Ω   , in this subinterval maximum grid size ℎ min, ≥ 1/ since there is no boundary layer, so both   and  are small.Using the triangle inequality, we have To derive a similar bound on   , we introduce the mesh function It can easily be seen that Also note that Hence the required result follows.

Numerical Experiments
Numerical results from this section will show that the moving mesh method produces numerical results which converge uniformly at  = T.The constant  will be set at  = 0.05, and no time step control mechanism or smoothing will be used.
Consider the following problem: where  is chosen such that is the exact solution when   (, )/ ∼ () at  = T = 10.
The solutions for the mesh equation and the physical pde will be obtained separately.This strategy is known as the decoupled approach.Discretizing (105a) and (105b), we obtain a system similar to (50) with  , = 1+ , (1 −  , ).The discretization for the monitor function is given as follows: The tolerance is set at 10 −6 .Table 1 shows the maximum point-wise errors (   −  exact

𝜀
) and the maximum error   max .   is the numerical solution on Ω  .From Table 1, it can be seen that for a fixed value of , as  increases, the error reduces.Also as  becomes smaller, the errors for any  become larger but stabilize after a while.This reflects the uniformity of the error estimates.The values of the rate of convergence given in Table 1 Table 2 shows that the method is almost of first order for a sufficiently large values of .Figures 3 and 5 show the mesh trajectory for  = 0.05 with different  values.The ability of the method to capture the boundary layer can clearly be seen in these figures, the region where the mesh points are concentrated shows the region of high derivatives, and as  increases, the mesh points become concentrated in the neighbourhood of  = 1 which is the boundary layer region.It should also be noted that the regular region is not starved of mesh points.The term 1 in the arc length monitor function plays this role; if a smaller value is used, more and more points will be packed in the boundary layer region, thereby starving the other region of mesh points which may lead to larger errors.Figure 4 shows the effect of increasing the value of , and an increase of this value will lead to a quasiuniform mesh and less points being placed in the boundary layer region.The number of mesh points placed in the layer region might be insufficient, thereby the method might fail to resolve the boundary layer.It is advisable to use small values of  to fully resolve the boundary layer.Figures 1 and 2 show the effect of reducing the value of  on the thickness and steepness on the boundary layer as can be seen in the neighbourhood of  = 1.

Conclusion
Theoretical analysis of the discrete problem was only done when the problem reached its steady state; this is mainly due to the absence of any theoretical analysis of the MMPDEs leading to the derivation of bounds for the nodal velocities.There should be a limit to how the mesh points are redistributed at the next iteration or time step, but up to date no bounds have been given explicitly.Apart from this drawback, this method fully resolves the boundary layer and is computationally efficient as can clearly be seen from the numerical results.
(, ) be the solution of   ; then