Abstract and Applied Analysis

Volume 2013 (2013), Article ID 214505, 11 pages

http://dx.doi.org/10.1155/2013/214505

Research Article

A Moving Mesh Method for Singularly Perturbed Problems

Stephen T. Sikwila¹ and Stanford Shateyi²

¹Department of Mathematics, University of Swaziland, Private Bag 4, Kwaluseni M201, Swaziland

²Department of Mathematics and Applied Mathematics, University of Venda, Private Bag X5050, Thohoyandou, South Africa

Received 29 May 2013; Accepted 19 October 2013

Academic Editor: Chun-Gang Zhu

Copyright © 2013 Stephen T. Sikwila and Stanford Shateyi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

1. Introduction

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 [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.

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 [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 , 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 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 , , 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

and 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. [6]

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 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.

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 [16], (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 .

Table 1: Maximum pointwise errors.

ε
Number of intervals

16
32
64
128
256
512

·
·
·
·
·
·
·

·
·
·
·
·
·
·

·
·
·
·
·
·
·

Table 2: Convergence rates.

ε
Number of intervals

16
32
64
128
256

0.85
0.96
1.06
1.47
0.25

0.78
0.90
1.00
1.14
1.68

0.74
0.87
0.94
1.03
1.29

0.69
0.81
0.89
0.98
1.13

0.64
0.76
0.84
0.93
1.12

0.63
0.71
0.79
0.87
0.91

0.64
0.68
0.74
0.80
0.85

0.66
0.65
0.72
0.77
0.83

0.67
0.64
0.71
0.75
0.80

0.67
0.64
0.71
0.75
0.77

·
·
·
·
·
·

·
·
·
·
·
·

·
·
·
·
·
·

0.68
0.64
0.70
0.73
0.75

0.68
0.64
0.70
0.73
0.75

0.68
0.64
0.70
0.73
0.75

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 .

Figure 1: Numerical solutions for and .

Figure 2: Numerical solutions for and .

Figure 3: Mesh trajectory with , , and .

Figure 4: Mesh trajectory with , , and .

Figure 5: Mesh trajectory for , , and .

6. 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.

References

N. Kopteva and M. Stynes, “A robust adaptive method for a quasi-linear one-dimensional convection-diffusion problem,” SIAM Journal on Numerical Analysis, vol. 39, no. 4, pp. 1446–1467, 2001.

T. Linss and T. U. Dresden, “Uniform point wise convergence of finite difference schemes using grid equi-distribution,” Journal of Computing, vol. 66, pp. 27–39, 2001.

J. Ma, Y. Jiang, and K. Xiang, “Numerical simulation of blowup in nonlocal reaction-diffusion equations using a moving mesh method,” Journal of Computational and Applied Mathematics, vol. 230, no. 1, pp. 8–21, 2009.

E. A. Dorfi and L. O. Drury, “Simple adaptive grids for 1-D initial value problems,” Journal of Computational Physics, vol. 69, no. 1, pp. 175–195, 1987.

S. Gowrisankar and S. Natesan, “Uniformly convergent numerical method for singularly perturbed parabolic initial-boundary-value problems with equidistributed grids,” International Journal of Computer Mathematics, 2013.

W. Huang, Y. Ren, and R. D. Russell, “Moving mesh partial differential equations (MMPDES) based on the equidistribution principle,” SIAM Journal on Numerical Analysis, vol. 31, no. 3, pp. 709–730, 1994.

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.

C. de Boor, “Good Approximation by Splines with Variable Notes II,” in Conference on the Numerical Solution of Differential Equations, Lecture Notes Series 363, Springer, Berlin, Germany, 1973.

W. Huang, Y. Ren, and R. D. Russell, “Moving mesh methods based on moving mesh partial differential equations,” Journal of Computational Physics, vol. 113, no. 2, pp. 279–290, 1994.

S. Li, L. Petzold, and Y. Ren, “Stability of moving mesh systems of partial differential equations,” SIAM Journal on Scientific Computing, vol. 20, no. 2, pp. 719–738, 1998.

Y. Ren and R. D. Russell, “Moving mesh techniques based upon equidistribution, and their stability,” SIAM Journal on Scientific Computing, vol. 13, no. 6, pp. 1265–1286, 1992.

Q. Zhou, Y. Chen, and Y. Yang, “Two improved algorithms and implementation for a singularly perturbed problem on moving meshes,” Journal of Systems Science and Complexity, vol. 24, no. 6, pp. 1232–1240, 2011.

J. Yang, “Alternative convergence analysis for a kind of singularly perturbed boundary value problems,” World Academy of Science, vol. 52, pp. 1–4, 2003.

Z. Cen, “A second-order finite difference scheme for a class of singularly perturbed delay differential equation,” International Journal of Computer Mathematics, pp. 254–302, 1998.

P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O'Riordan, and G. I. Shishkin, Parameter Robust Computational Methods for Singularly Perturbed Differential Equations, CRC Press, 2000.

J. J. H. Miller, E. O'Riordan, and G. I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, River Edge, NJ, USA, 1998.

K. W. Morton, Numerical Solution of Convection Diffusion Problems, Chapman and Hall, London, UK, 1995.

S. Gowrisankar and S. Natesan, “The parameter uniform numerical method for singularly perturbed parabolic reactiondiffusion problems on equidistributed grids,” Applied Mathematics Letters, vol. 26, no. 11, pp. 1053–1060, 2013.

J. Ma and Y. Jiang, “Moving mesh techniques based upon equidistribution, and their stability,” Journal of Computational Physics, vol. 228, no. 18, pp. 6977–6990, 2009.