A Moving Mesh Method for Singularly Perturbed Problems
A singularly perturbed time dependent convection diffusion problem is solved on a rectangular domain, using the moving mesh method which uses the equidistribution principle. The problem has a boundary at the steady state. It is shown that the numerical approximations generated by the moving mesh method converge uniformly with respect to the singular perturbation parameter. Theoretical results are obtained which are verified using numerical results.
The use of adapted meshes [1–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–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 , 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. , 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.  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.  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  considered a kind of nonconservative singularly perturbed two-point value problems in fluid dynamics. Cen  examined a class of delay differential equations with a perturbation parameter ε. More recently, Gowrisankar and Natesan  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.
Consider the following time dependent convection-diffusion problem: with , and for all where , , , are constants, The functions and are assumed to be sufficiently smooth. In general the solution will be smooth on for all values of . Boundary and interior layers [15–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  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 , and .
2. The Continuous Problem
The differential operator for satisfies the following maximum principle.
Theorem 1 (maximum principle). Let be any function in the domain of and assume that , for all . Then for all , this implies that for all .
Proof. Assume that there exists such that ; then since for all ; hence . Let for all . Then for all , and ; thus the minimum of must be also negative. Let such that Applying the differential operator to gives which can be written as where The argument now divides into two cases depending on the position of , or . If , we have that It can be seen that , and for all . This leads to the following inequality: . If , we have that It also follows that , , and for all and also leads to the following inequality: . This is a contradiction, and thus our original assumption is false and we can conclude that the minimum of is non-negative.
An immediate consequence of this is the following bound on the solution of any problem from .
Lemma 2. Let be the solution of ; then for all .
Proof. Consider the barrier functions For , The maximum principle now applies, and we have for all from which we have the required result.
Lemma 3. Let be the solution of ; then the spatial derivatives satisfy the bounds
Proof. Note that which gives for all where depends on . From the mean-value theorem, there exists a point such that hence Integrating with respect to , we obtain Using (18) with and combining with (15), it follows that Equation (15) can also be used to give the following bound: for all . This proves the result for . To obtain bounds for the higher derivatives, rewriting (10) this gives the second derivative bound where depends on , , and . Differentiating with respect to 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 is the solution to problemand 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 , , the smooth component satisfies the singular component satisfies for some constant independent of .
Proof. To find these bounds, we rewrite the smooth component as
where is the solution to the reduced problem , satisfies
We clearly have in with on . It can be easily seen that and are all bounded by a constant independent of and
for . Therefore 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 since and are independent of . As previously stated that the problem for is similar to the problem for , we can use Lemma 3 from which we obtain for
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 ; hence To establish the bounds for the derivatives of , we start by noting that where depends on . From the mean value theorem, there exists a point such that Using the triangle inequality, it can easily be seen that Using (39) at the point and , we obtain that Integrating (25a), (25b), and (25c) with respect to , we obtain that at the point Hence From (44) and (46), we obtain that which is the required bound for the first derivative. From (25a), (25b), and (25c), it can be seen that which yields the required estimates for the second and third derivatives for the singular component .
3. 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 , can be written as Discretizing (49) using an implicit scheme, which can be written as a system represents the solution at the point , , is the interpolated value of on the new grid obtained at time , and is the node velocity obtained from the moving mesh partial differential equation (MMPDE) (53) derived by Huang et al.  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 , we have a uniform mesh. The initial condition is assumed not to be identically zero over the whole domain and to be bounded by some constant . 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 , since Assume that this result is true for some time ; hence . Since the matrix is invertible at , it follows that exists for all and is bounded by some constant , .
Define the discrete operator as where for all . 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 for all and for all , then for all .
Proof. Assume that there exists such that then , which implies that , and we also know that We have to show that , and we proceed by contradiction, suppose that Since is an arbitrary point, we have that , but which is a contradiction. Our supposition that is false, so we have that . All the arguments given above imply that . Therefore we must have , but we have that and ; hence Using the same argument as before, this implies that , but which is a contradiction. So , which is a contradiction. Thus our original assumption must be false, and we conclude that the minimum of is nonnegative.
4. 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 equationswhere 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.
Lemma 6. Let ; then for any ,
Proof. Using integration by parts to reduce the order of differentiation in the integral, can be expressed as It follows that
Lemma 7. Assuming the bound (27), the error in the smooth component of the numerical solution satisfies the following error bound: where for some constant independent of and .
Proof. We start by considering the local truncation error for the smooth component It can be easily shown that Using the local truncation error estimates from , (69) (72), it follows that We introduce the two mesh functions it can be seen from the mesh functions that By the discrete maximum principle, we conclude that and so for all ; hence
Lemma 8. For all and at each , the singular component of the error satisfies 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 and ; in this case , and our mesh is quasiuniform . Using the standard bound for the local truncation error, we can derive an equivalent expression for the truncation error as for the smooth component
where . 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 for ; by the discrete maximum principle, it follows that ; hence
To consider the case when and , we will assume that and that there exists a point such that . 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 since there is no boundary layer, so both and are small. Using the triangle inequality, we have It suffices to bound and separately. Using (28) for , we obtain that To derive a similar bound on , we introduce the mesh function It can easily be seen that Also note that and . Applying the operator to , Now considering the mesh function , we have From the maximum principle, this means that ; hence
It only now remains to give a bound for in the interval . If , the expression for the truncation error for the singular component (78) does not change since in the boundary layer region . From (65a) and (65b), we have that From (86) and (91), it follows that Consider the barrier functions and the mesh function for and , Applying the discrete difference operator to and noting that , hence it follows that The discrete maximum principle on then gives and it follows that Combining the separate estimates in the two subintervals gives
Theorem 9. If is the solution of and , the corresponding numerical solution using the method outlined in (50) satisfies for where is a constant independent of and .
Proof. We start by noting that We have already shown in the previous lemmas that Hence the required result follows.
5. Numerical Experiments
Numerical results from this section will show that the moving mesh method produces numerical results which converge -uniformly at . The constant will be set at , 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 . 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 . The discretization for the monitor function is given as follows: The tolerance is set at . Table 1 shows the maximum point-wise errors and the maximum error . 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 are obtained using Table 2 and (108) Table 2 also gives the uniform convergence rates in the last row. These values are calculated using Table 1 and (109) 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 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 which is the boundary layer region. It should also be noted that the regular region is not starved of mesh points. The term 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 .
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.
G. Beckett and J. A. Mackenzie, “Convergence analysis of finite difference approximations on equidistributed grids to a singularly perturbed boundary value problem,” Applied Numerical Mathematics, vol. 35, no. 2, pp. 87–109, 2000.View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
W. Huang, Y. Ren, and R. D. Russell, “Moving mesh methods based on moving mesh partial differential equations,” Journal of Computational Physics,