# A Moving Mesh Method for Singularly Perturbed Problems

Authors: Stephen T. Sikwila; Stanford Shateyi
Journal: Abstract and Applied Analysis
 (2013)

---

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


---

## Body

## 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:(1)Pε:{-ε∂2∂x2uε(x,t)+∂∂x(b(x,t)uε(x,t))+∂∂tuε(x,t)=f(x,t),(x,t)∈Ω,uε(0,t)=uε(1,t)=0,uε(x,0)=g(x),
with 0<ε<1, λ-≥bx≥λ>0 and β-≥b≥β>0 for all (x,t)∈Ω-=Ω∪ΓB∪ΓL∪ΓR where λ, λ-, β, β- are constants,
(2)Ω=(0,1)×(0,𝒯],ΓL={(0,t)∣0≤t≤𝒯},ΓB={(x,0)∣0≤x≤1},ΓR={(1,t)∣0≤t≤𝒯},Γ=ΓL∪ΓR∪ΓB.
The functions f and b are assumed to be sufficiently smooth. In general the solution uε will be smooth on Ω for all values of t. 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 Pε. 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 (x,t)∈Ω, ∥dx(t)/dt∥≤C, for some constant C and at t=𝒯, ∥∂uε(x,t)/∂t∥~O(ε) and ∥dx(t)/dt∥~O(ε).


## 2. The Continuous Problem

The differential operatorLε for Pε satisfies the following maximum principle.Theorem 1 (maximum principle).
Letψ(x,t) be any function in the domain of Lε and assume that ψ(x,t)≥0, for all (x,t)∈Γ. Then Lεψ(x,t)≥0 for all (x,t)∈Ω, this implies that ψ(x,t)≥0 for all (x,t)∈Ω-.Proof.
Assume that there existsr=(xr,tr)∈Ω- such that ψ(r)=minΩ-ψ<0; then r∉Γ since ψ(x,t)≥0 for all (x,t)∈Γ; hence r∈Ω. Let
(3)ψ*(x,t)=ψ(x,t)e(β/ε)(1-x)
for all (x,t)∈Ω-. Then ψ*(x,t)≥0 for all (x,t)∈Γ, and ψ(r)<0; thus the minimum of ψ* must be also negative. Let q=(xq,tq)∈Ω such that
(4)ψ*(q)=minΩ-ψ*<0.
Applying the differential operator Lε to ψ gives
(5)Lεψ(x,t)=Lε(ψ*e-(β/ε)(1-x))=-ε∂2∂x2(ψ*e-(β/ε)(1-x))+∂∂x(bψ*e-(β/ε)(1-x))+∂∂t(ψ*e-(β/ε)(1-x))=(-εψxx*+bxψ*+bψx*+ψt*)e-(β/ε)(1-x)-2βψx*e-(β/ε)(1-x)+(-β2ε+bβε)ψ*e-(β/ε)(1-x)
which can be written as
(6)Lεψ(x,t)=(η1+η2+η3)e-(β/ε)(1-x),
where
(7)η1=-εψxx*+bxψ*+bψx*+ψt*,η2=(-β2ε+bβε)ψ*,η3=-2βψx*.
The argument now divides into two cases depending on the position of q, q∈Γt={(x,𝒯)∣0<x<1} or q∉Γt. If q∉Γt, we have that
(8)ψxx*(q)>0,ψx*(q)=ψt*(q)=0.
It can be seen that η1<0, η2<0 and η3=0 for all x∈(0,1). This leads to the following inequality: Lεψ(q)<0. If q∈Γt, we have that
(9)ψxx*(q)>0,ψx*(q)=0,ψt*(q)≤0.
It also follows that η1<0, η2<0, and η3=0 for all (x,t)∈Ω and also leads to the following inequality: Lεψ(q)<0. This is a contradiction, and thus our original assumption is false and we can conclude that the minimum of ψ*(x,t) is non-negative.An immediate consequence of this is the following bound on the solution of any problem fromPε.Lemma 2.
Letuε(x,t) be the solution of Pε; then
(10)∥uε∥≤∥g∥+∥f∥β,
for all (x,t)∈Ω-.Proof.
Consider the barrier functions(11)ψ±(x,t)=∥g∥+x∥f∥β±uε(x,t),ψ±(0,t)=∥g∥±uε(0,t)=∥g∥≥0,ψ±(1,t)=∥g∥+∥f∥β±uε(1,t)=∥g∥+∥f∥β≥0,ψ±(x,0)=∥g∥+x∥f∥β±uε(x,0)≥0.
For (x,t)∈Ω,
(12)Lεψ±(x,t)=-ε∂2∂x2(∥g∥+x∥f∥β±uε(x,t))+∂∂x(b(∥g∥+x∥f∥β±uε(x,t)))+∂∂t(∥g∥+x∥f∥β±uε(x,t))=bx(∥g∥+x∥f∥β)+b∥f∥β±f(x,t)≥0.
The maximum principle now applies, and we have ψ±(x,t)≥0 for all (x,t)∈Ω- from which we have the required result.Lemma 3.
Letuε(x,t) be the solution of Pε; then the spatial derivatives ∂kuε(x,t)/∂xk satisfy the bounds
(13)∥∂kuε(x,t)∂xk∥≤Cε-kmax{∥f∥,∥uε∥},fork=1,2,∥∂3uε(x,t)∂x3∥≤Cε-3max{∥f∥,∥f′∥,∥uε∥}.Proof.
Note that(14)∫x1‍∂∂s(b(s,t)uε(s,t))ds=b(s,t)uε(s,t)|s=x1
which gives
(15)|∫x1‍(f(s,t)-∂∂s(b(s,t)uε(s,t)))ds|≤∥f∥+C∥uε∥
for all (x,t)∈Ω where C depends on ∥b∥. From the mean-value theorem, there exists a point z∈(1-ε,1) such that
(16)∂uε(z,t)∂x=uε(1,t)-uε(1-ε,t)ε;
hence
(17)|ε∂uε(z,t)∂x|≤∥uε∥.
Integrating Pε with respect to x, we obtain
(18)-ε(∂uε(1,t)∂x-∂uε(x,t)∂x)00=∫x1‍(f(s,t)-∂∂s(b(s,t)uε(s,t))-∂uε(s,t)∂t)ds.
Using (18) with x=z and combining with (15), it follows that
(19)|ε∂uε(1,t)∂x|≤∥f∥+C∥uε∥.
Equation (15) can also be used to give the following bound:
(20)|ε∂uε(x,t)∂x|≤∥f∥+C∥uε∥,
for all (x,t)∈Ω. This proves the result for k=1. To obtain bounds for the higher derivatives, rewriting (10)
(21)-ε∂2uε(x,t)∂x2=f(x,t)-∂∂x(b(x,t)uε(x,t))-∂uε(x,t)∂t,
this gives the second derivative bound where C depends on ∥g∥, ∥b∥, and ∥bx∥. Differentiating Pε with respect to x(22)-ε∂3uε(x,t)∂x3=∂∂x(f(x,t)-∂∂x(b(x,t)uε(x,t))0000-∂uε(x,t)∂t)
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:(23)uε(x,t)=vε(x,t)+wε(x,t),(x,t)∈Ω-,
where vε(x,t) is the solution to problem(24a)Lεvε(x,t)=f(x,t),(x,t)∈Ω,(24b)vε(x,t)=uε(x,t),(x,t)∈ΓB∪ΓL,and the singular component wε(x,t) is the solution of the homogenous problem (25a)Lεwε(x,t)=0,(x,t)∈Ω,(25b)wε(x,t)=0,(x,t)∈ΓB∪ΓL,(25c)wε(x,t)=uε(x,t)-vε(x,t),(x,t)∈ΓR.Theorem 4.
The solutionuε(x,t) of the continuous problem Pε can be decomposed as a sum of the smooth and layer functions
(26)uε(x,t)=vε(x,t)+wε(x,t),(x,t)∈Ω-,
where for all k, 0≤k≤3, the smooth component vε(x,t) satisfies
(27)∥∂kvε(x,t)∂xk∥≤C(1+ε-(k-2));
the singular component wε(x,t) satisfies
(28)|wε(x,t)|≤Ce-(β/ε)(1-x),∀x∈Ω-,(29)∥∂kwε(x,t)∂xk∥≤Cε-ke-(β/ε)(1-x),k=1,2,3
for some constant C independent of ε.Proof.
To find these bounds, we rewrite the smooth component as(30)vε(x,t)=(v0+εv1+ε2v2)(x,t),(x,t)∈Ω-,
where v0(x,t) is the solution to the reduced problem P0, v1(x,t) satisfies
(31)∂∂x(b(x,t)v1(x,t))+∂v1(x,t)∂t=∂2v0(x,t)∂x2,(x,t)∈Ω,v1(x,t)=0,(x,t)∈ΓB∪ΓL,
and v2(x,t) satisfies
(32)Lεv2(x,t)=∂2v1(x,t)∂x2,(x,t)∈Ω,v2(x,t)=0,(x,t)∈Γ.
We clearly have Lεvε(x,t)=f(x) in Ω with vε(x,t)=v0(x,t)+εv1(x,t) on ΓR. It can be easily seen that v0(x,t) and v1(x,t) are all bounded by a constant independent of ε and
(33)∥∂i+jv1(x,t)∂xi∂tj∥≤C,(x,t)∈Ω-
for 0≤i+2j≤4. Therefore v2(x,t) is the solution of a problem similar to Pε; hence from Lemma 2, we obtain that
(34)∥v2∥≤C⟹∥vε∥≤C(1+ε+ε2).
To get the bounds for the spatial derivatives, we only consider v2(x,t) since v0(x,t) and v1(x,t) are independent of ε. As previously stated that the problem for v2 is similar to the problem for uε(x,t), we can use Lemma 3 from which we obtain for 0≤k≤3(35)∥∂kvε(x,t)∂xk∥≤C(1+ε-(k-2)).
To find bounds for the singular component, we consider the following mesh functions:(36)ψ±(x,t)=Ce-(β/ε)(1-x)±wε(x,t),
where C=∥v0∥+∥v1∥.
(37)ψ±(0,t)=Ce-(β/ε)±wε(0,t)=Ce-(β/ε)≥0,ψ±(1,t)=C±wε(1,t)≥0,ψ±(x,0)=Ce-(β/ε)(1-x)±wε(x,0)=Ce-(β/ε)(1-x)≥0.
Applying the differential operator Lε, we obtain that
(38)Lεψ±(x,t)=C(bβε-β2ε+bx)e-(β/ε)(1-x)≥0.
Thus from the maximum principle, we can say that ψ±(x,t)≥0; hence
(39)|wε(x,t)|≤Ce-(β/ε)(1-x).
To establish the bounds for the derivatives of wε(x,t), we start by noting that
(40)|∫x1‍∂∂s(b(s,t)wε(s,t))ds|≤C∥wε(x,t)∥≤Ce-(β/ε)(1-x),
where C depends on ∥b∥. From the mean value theorem, there exists a point z∈(1-ε,1) such that
(41)∂∂x(b(z,t)wε(z,t))=b(1,t)wε(1,t)-b(1-ε,t)wε(1-ε,t)ε.
Using the triangle inequality, it can easily be seen that
(42)|∂wε(z,t)∂x|22≤1|b(z,t)|2222×(|∂b(z,t)∂xwε(z,t)|2222222+|b(1,t)wε(1,t)-b(1-ε,t)wε(1-ε,t)ε|).
Using (39) at the point x=1-ε,1 and x=z, we obtain that
(43)|∂wε(z,t)∂x|≤Cε-1e-(β/ε)(1-z).
Integrating (25a), (25b), and (25c) with respect to x, we obtain that
(44)ε∂wε(1,t)∂x=ε∂wε(x,t)∂x+∫x1‍(∂∂s(b(s,t)wε(s,t))+∂wε(s,t)∂t)ds
at the point x=z(45)|ε∂wε(1,t)∂x|=|ε∂wε(z,t)∂x0000+∫z1‍(∂∂s(b(s,t)wε(s,t))+∂wε(s,t)∂t)ds|.
Hence
(46)|ε∂wε(1,t)∂x|≤Ce-(β/ε)(1-z).
From (44) and (46), we obtain that
(47)|ε∂wε(x,t)∂x|≤Ce-(β/ε)(1-x),∀(x,t)∈Ω
which is the required bound for the first derivative. From (25a), (25b), and (25c), it can be seen that
(48)ε∂2wε(x,t)∂x2=∂∂x(b(x,t)wε(x,t))+∂wε(x,t)∂t,ε∂3wε(x,t)∂x3=∂∂x(∂∂x(b(x,t)wε(x,t))+∂wε(z,t)∂t),
which yields the required estimates for the second and third derivatives for the singular component wε.


## 3. The Discretized Problem

In this section the discretization process forPε is considered. By changing the time derivative into the Lagrangian form, this enables the introduction of node velocities into the system. Setting 𝒜u=ε(∂u/∂x)-bu, Pε can be written as
(49)dudt=dxdt∂u∂x+∂∂x(𝒜u)+f.
Discretizing (49) using an implicit scheme,
(50)ui,j+δt(x˙i,j-1Dx-+Dx-𝒜)ui,j=u_i,j-1+δtf(xi,j)
which can be written as a system
(51)Tu=Ψ.ui,j represents the solution at the point (xi,tj),
(52)Dx-ui,j=ui,j-ui-1,jhi,j,Dx+ui,j=ui+1,j-ui,jhi+1,j.hi,j=xi,j-xi-1,j, u_i,j-1 is the interpolated value of u on the new grid obtained at time t=tj, and x˙ is the node velocity obtained from the moving mesh partial differential equation (MMPDE) (53) derived by Huang et al. [6]
(53)(Mi+1,j+Mi,j)(x˙i+1,j-x˙i,j)-(Mi,j+Mi-1,j)(x˙i,j-x˙i-1,j)=-Ei,j-1τ,
where
(54)Ei,j-1=(Mi+1,j-1+Mi,j-1)(xi+1,j-1-xi,j-1)-(Mi,j-1+Mi-1,j-1)(xi,j-1-xi-1,j-1).τ is a user defined parameter, which determines how fast the grid moves towards the equidistributed grid, and M is the arc length monitor function
(55)M=1+(∂u∂x)2.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 t=t0=0, we have a uniform mesh. The initial condition ui,0 is assumed not to be identically zero over the whole domain and to be bounded maxi|ui,0|≤C- by some constant 0<C-<∞. The system for the MMPDE can be written in such a way that we get an M-matrix on the left hand side. This means that a solution for the system exists, and we can assume that maxi,j|x˙i,j|<C-1, since
(56)Ψi,0=u_i,0+δtf(xi,1)⟹|Ψi,0|≤C-+δt∥f∥=C-2<∞.
Assume that this result is true for some time t=tn-1; hence |Ψi,n-1|≤C-2<∞. Since the matrix T is invertible at t=tn=𝒯, it follows that ui,n exists for all i and is bounded by some constant C, maxi,n|ui,n|≤C.Define the discrete operatorLεN as
(57)LεNui,j≡(-εδx2+Dx-bi,j+Dt--x˙i,jDx-)ui,j,
where
(58)δx2ui,j=2hi+1,j+hi,j(Dx+ui,j-Dx-ui,j)
for all 1≤i≤N-1. The discrete differential operator LεN in (57) satisfies the following discrete maximum principle.Theorem 5 (discrete maximum principle).
Letψi,j be any mesh function defined on Ω-N. If ψi,j≥0 for all (xi,tj)∈ΓN and LεNψi,j≥0 for all (xi,tj)∈Ω-N, then ψi,j≥0 for all (xi,tj)∈Ω-N.Proof.
Assume that there existsr=(xi,tj)∈Ω-N such that
(59)ψi,j=ψ(r)=minΩ-Nψ<0;
then r∉ΓN, which implies that r∈ΩN, and we also know that
(60)δx2ψi,j≥0,Dx-ψi,j≤0,Dt-ψi,j≤0.
We have to show that Dx-(bi,jψi,j)≤0, and we proceed by contradiction, suppose that
(61)Dx-(bi,jψi,j)>0⟹bi,jψi,j-bi-1,jψi-1,j>0⟹bi,jψi,j>bi-1,jψi-1,j.
Since r is an arbitrary point, we have that b0,jψ0,j<0, but b0,jψ0,j≥0 which is a contradiction. Our supposition that Dx-(bi,jψi,j)>0 is false, so we have that Dx-(bi,jψi,j)≤0. All the arguments given above imply that LεNψi,j≤0. Therefore we must have LεNψi,j=0, but we have that -εδx2ψi,j≤0 and Dx-(bi,jψi,j)≤0; hence
(62)ψi+1,j=ψi,j=ψi-1,j<0.
Using the same argument as before, this implies that ψ0,j<0, but ψ0,j≥0 which is a contradiction. So LεNψi,j<0, which is a contradiction. Thus our original assumption must be false, and we conclude that the minimum of ψi,j is nonnegative.


## 4. Decomposition of Numerical Solution and Error Estimates

LetUε denote the discrete solution, and assume that the discrete solution can be decomposed into the sum
(63)Uε(xi,tj)=Vε(xi,tj)+Wε(xi,tj),(xi,tj)∈Ω-N,
where Vε(xi,tj) and Wε(xi,tj) are solutions of the respective equations(64a)LεNVε(xi,tj)=f(xi,tj),(xi,tj)∈ΩN,(64b)Vε(xi,tj)=vε(xi,tj),(xi,tj)∈ΓN,(65a)LεNWε(xi,tj)=0,(xi,tj)∈ΩN,(65b)Wε(xi,tj)=wε(xi,tj),(xi,tj)∈ΓN,where Wε is the singular component and Vε is the smooth component. The error for the numerical solution will be decomposed as
(66)(Uε-uε)(xi,tj)=((Vε-vε)+(Wε-wε))(xi,tj),(xi,tj)∈Ω-N.
The error (Uε-uε)(xi,tj) can now be estimated using the error estimates for the singular and smooth components of the solution.Lemma 6.
Let(xi,tj)∈ΩN; then for any ψ(x,t)∈C2(Ω-),
(67)|(Dx--∂∂x)(b(xi,tj)ψ(xi,tj))|≤12hi,j∥∂2∂x2(b(x,t)ψ(x,t))∥.Proof.
Using integration by parts to reduce the order of differentiation in the integral,(Dx--∂/∂x)(b(xi,tj)ψ(xi,tj)) can be expressed as
(68)(Dx--∂∂x)(b(xi,tj)ψ(xi,tj))=1hi,j∫xi-1,jxi,j‍(xi-1,j-s)∂2∂s2(b(s,t)ψ(s,t))ds.
It follows that
(69)|(Dx--∂∂x)(bψ)(xi,tj)|≤1hi,j∥∂2∂x2(b(x,t)ψ(x,t))∥∫xi-1,jxi,j‍(s-xi-1,j)ds≤12hi,j∥∂2∂x2(b(x,t)ψ(x,t))∥.Lemma 7.
Assuming the bound (27), the error in the smooth component of the numerical solution satisfies the following error bound:
(70)|(Vε-vε)(xi,tj)|≤C(hmax,j+δt),
where hmax,j=maxi{hi,j∣1≤i≤N} for some constant C independent of ε and N.Proof.
We start by considering the local truncation error for the smooth component(71)LεN(Vε-vε)(xi,tj)=(Lε-LεN)vε(xi,tj)=(-ε(∂2∂x2-δx2)+(∂∂x-Dx-)b(xi,tj)+(∂∂t-(Dt--Dt-xi,jDx-)))vε(xi,tj).
It can be easily shown that
(72)|Dx-vε(xi,tj)|=|1xi,j-xi-1,j∫xi-1,jxi,j‍∂vε(s,t)∂sds|≤∥∂vε(x,t)∂x∥≤C(1+ε),|Dt-x(tj)|=|1tj-tj-1∫tj-1tj‍dx(s)dsds|≤∥dx(t)dt∥≤C.
Using the local truncation error estimates from [16], (69) (72), it follows that
(73)|LεN(Vε-vε)(xi,tj)|≤ε3(hi+1,j+hi-1,j)∥∂3vε(x,t)∂x3∥+12hi,j∥∂2∂x2(bvε)(x,t)∥+12(tj-tj-1)∥∂2vε(x,t)∂t2∥+∥∂vε(x,t)∂x∥∥dx(t)dt∥≤C(hmax,j+δt).
We introduce the two mesh functions
(74)ψi,j±=Cβ-1(hmax,j+δt)(β+xi,j)±(Vε-vε)(xi,tj)
it can be seen from the mesh functions that
(75)ψ0,j±=C(hmax,j+δt)≥0,ψN,j±=Cβ-1(hmax,j+δt)(β+1)≥0,ψi,0±=Cβ-1(N-1+δt)(β+xi,0)≥0,LεNψi,j±=Cβ-1(hmax,j+δt)(βbi,j-bi-1,jhi,j+bi,jxi,j-bi-1,jxi-1,jhi,j)±LεN(Vε-vε)(xi,tj)≥C(hmax,j+δt)(bi,j-bi-1,jhi,j+bi,jβ)±LN(Vε-vε)(xi,tj)≥0.
By the discrete maximum principle, we conclude that ψi,j±≥0 and so for all (xi,tj)∈Ω-N; hence
(76)|(Vε-vε)(xi,tj)|≤C(hmax,j+δt).Lemma 8.
For allN≥4 and at each (xi,tj)∈Ω-N, the singular component of the error satisfies
(77)|(Wε-wε)(xi,tj)|≤C(N-1(lnN)2+δt)
for some constant C independent of ε and N.Proof.
We need to consider two separate cases since the role of the boundary layer is crucial. We start with the case whenε-1≤2/βln1/ε and ε≥N-1; in this case ε-1≤2/βlnN, and our mesh is quasiuniform hmin,j~hmax,j. Using the standard bound for the local truncation error, we can derive an equivalent expression for the truncation error as for the smooth component
(78)|LN(Wε-wε)(xi,tj)|≤ε3(hi+1,j+hi-1,j)∥∂3wε(x,t)∂x3∥+12hi,j∥∂2∂x2(bwε)(x,t)∥+12(tj-tj-1)∥∂2wε(x,t)∂t2∥+∥∂wε(x,t)∂x∥∥dx(t)dt∥≤Cε-2(hi+1,j+hi,j)e-(β/ε)(1-x)+Cε-2βmaxhi,je-(β/ε)(1-x)+C(tj-tj-1)+Cε-1e-(β/ε)(1-x)≤Cε-2hmax,j(1+βmax)+Cδt≤C(N-1(lnN)2+δt),
where βmax=max{∥b∥,∥bx∥,∥bxx∥}. Consider the following barrier function:
(79)ψi,j=Cβ-1(N-1(lnN)2+δt)(β+xi,j)
and the mesh functions
(80)ϕi,j±=ψi,j±(Wε-wε)(xi,tj).
It can be easily seen that
(81)ϕ0,j±=C(N-1(lnN)2+δt)≥0,ϕN,j±=Cβ-1(N-1(lnN)2+δt)(β+1)≥0,ϕi,0±=Cβ-1(N-1(lnN)2+δt)(β+xi,0)≥0.
Applying the difference operator LεN to the barrier function,
(82)LεNψi,j=Cβ-1(N-1(lnN)2+δt)×(βbi,j-bi-1,jhi,j+bi,jxi,j-bi-1,jxi-1,jhi,j)≥C(N-1(lnN)2+δt)(bi,j-bi-1,jhi,j+bi,jβ)≥|LεN(Wε-wε)(xi,tj)|.
Hence it follows that LεNϕi,j±≥0 for 1≤i≤N-1; by the discrete maximum principle, it follows that ϕi,j±≥0; hence
(83)|(Wε-wε)(xi,tj)|≤ψi,j≤C(N-1(lnN)2+δt).
To consider the case whenε-1>2/βln1/ε and ε<N-1, we will assume that hi,j≤hi-1,j and that there exists a point xα∈[0,1] such that 1-xα≥ε/βln1/ε. The point xα splits the interval into the coarse mesh ΩαlN and fine mesh ΩαrN where
(84)ΩαlN={xi∣0≤xi≤xα,hi,j≥1N},ΩαrN={xi∣xα<xi≤1,hi,j<1N}.
We give separate proofs in the coarse and fine mesh subintervals. First suppose that xi∈ΩαlN, in this subinterval maximum grid size hmin,j≥1/N since there is no boundary layer, so both Wε and w are small. Using the triangle inequality, we have
(85)|(Wε-wε)(xi,tj)|≤|Wε(xi,tj)|+|wε(xi,tj)|.
It suffices to bound Wε and wε separately. Using (28) for xi∈ΩαlN, we obtain that
(86)|wε(xi,tj)|≤Ce-(β/ε)(1-xi)≤Ce-(β/ε)(1-xα)≤Ce-ln1/ε=Cε≤CN-1.
To derive a similar bound on Wε, we introduce the mesh function Y(87)Y(xi,tj)=Yi,j=CβN(β+xi,j).
It can easily be seen that
(88)Y0,j=CN≥0,YN,j=CβN(β+1)≥0,Yi,0=CβN(β+xi,0)≥0,Dx-(bY)(xi,tj)=CβN(βbi,j-bi-1,jhi,j+bi,jxi,j-bi-1,jxi-1,jhi,j)≥0.
Also note that δx2Y(xi,tj)=0 and Dt-xi,j-Dt-xi,jDx-xi,j=0. Applying the operator LεN to Y,
(89)LεNY(xi,tj)=-εδx2Yi,j+D-(bY)i,j+Dt-Yi,j-Dt-Yi,jDx-Yi,j≥0.
Now considering the mesh function Yε-Wε, we have
(90)LεN(Y-Wε)(xi,tj)=LεNYi,j-LNWi,j=LεNYi,j≥0.
From the maximum principle, this means that Y-Wε≥0; hence
(91)Wi,j≤Yi,j=CβxiN≤CN-1.
It only now remains to give a bound for(Wε-wε) in the interval (xα,1]. If xi∈ΩαrN, the expression for the truncation error for the singular component (78) does not change since in the boundary layer region hmax,j≤1/N. From (65a) and (65b), we have that
(92)|Wε(1,tj)-wε(1,tj)|=0.
From (86) and (91), it follows that
(93)|Wε(xα,tj)-wε(xα,tj)|≤|Wε(xα,tj)|+|wε(xα,tj)|≤CN-1.
Consider the barrier functions
(94)ψi,j=C(N-1(lnN)2+δt)(xi,jβ+1)
and the mesh function
(95)ϕi,j±=ψi,j±(Wε-wε)(xi,tj),
for xi=xα and xi=1,
(96)ϕ±(xα,tj)=C(N-1(lnN)2+δt)(xαβ+1)±(Wε-wε)(xα,tj)≥0,ϕN,j±=C(N-1(lnN)2+δt)(1β+1)≥0,ϕi,0±=C(N-1(lnN)2+δt)(xi,0β+1)≥0.
Applying the discrete difference operator LεN to ψi,j and noting that Dt-xi,j-Dt-xi,jDx-xi,j=0,
(97)LNψi,j=C(N-1(lnN)2+δt)×(bi,jxi,j-bi-1,jxi-1,jβhi,j+bi,j-bi-1,jhi,j)≥C(N-1(lnN)2+δt)(biβ+bi,j-bi-1,jhi,j)≥|LN(Wε-wε)(xi,tj)|;
hence it follows that
(98)LεNϕi,j±≥0,xi∈ΩαrN.
The discrete maximum principle on ΩαrN then gives
(99)ϕi,j±≥0,xi∈ΩαrN,
and it follows that
(100)|(Wε-wε)(xi,tj)|≤ψi,j≤C(N-1(lnN)2+δt).
Combining the separate estimates in the two subintervals gives
(101)|(Wε-wε)(xi,tj)|≤CN-1((lnN)2+δt),(xi,tj)∈Ω-N.Theorem 9.
Ifuε(x,t) is the solution of Pε and Uε(xi,tj), the corresponding numerical solution using the method outlined in (50) satisfies
(102)sup0<ε<1∥Uε-uε∥≤CN-1((lnN)2+δt)
for N≥4 where C is a constant independent of ε and N.Proof.
We start by noting that(103)∥Uε-uε∥≤∥Vε-vε∥+∥Wε-wε∥.
We have already shown in the previous lemmas that
(104)∥Vε-vε∥≤C(hmax+δt),∥Wε-wε∥≤CN-1((lnN)2+δt).
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 t=𝒯. The constant τ will be set at τ=0.05, and no time step control mechanism or smoothing will be used.Consider the following problem:(105a)-εδ2uε(x,t)δx2+δδx((1+x(1-x))uε(x,t))+δuε(x,t)δt=f,(x,t)∈Ω,(105b)uε(0,t)=uε(1,t)=0,uε(x,0)=1-4(x-12)2,where f is chosen such that
(106)uε(x,t)=1-e-(1-x)/ε1-e-1/ε-cosπx2
is the exact solution when ∂uε(x,t)/∂t~O(ε) 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 bi,j=1+xi,j(1-xi,j). The discretization for the monitor function is given as follows:
(107)Mi,j=1+(ui+1,j-ui-1,jxi+1,j-xi-1,j)2.
The tolerance is set at 10-6. Table 1 shows the maximum point-wise errors (uεN-uεexact) and the maximum error EmaxN. uεN is the numerical solution on ΩN. From Table 1, it can be seen that for a fixed value of ε, as N increases, the error reduces. Also as ε becomes smaller, the errors for any N 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)
(108)RN=log2[∥uεN-uεexact∥∥uε2N-uεexact∥].
Table 2 also gives the uniform convergence rates in the last row. These values are calculated using Table 1 and (109)
(109)RmaxN=log2[maxε∥uεN-uεexact∥maxε∥uε2N-uεexact∥].
Table 2 shows that the method is almost of first order for a sufficiently large values of N.Table 1
Maximum pointwise errors.
ε
Number of intervalsN
16
32
64
128
256
512
2
-
2
0.21
D
-
01
0.12
D
-
01
0.59
D
-
02
0.28
D
-
02
0.10
D
-
02
0.87
D
-
03
2
-
3
0.49
D
-
01
0.28
D
-
01
0.15
D
-
01
0.76
D
-
02
0.35
D
-
02
0.11
D
-
02
2
-
4
0.76
D
-
01
0.45
D
-
01
0.25
D
-
01
0.13
D
-
01
0.64
D
-
02
0.26
D
-
02
2
-
5
0.94
D
-
01
0.58
D
-
01
0.33
D
-
01
0.18
D
-
01
0.91
D
-
02
0.42
D
-
02
2
-
6
0.10
D
-
00
0.67
D
-
01
0.40
D
-
01
0.22
D
-
01
0.12
D
-
01
0.53
D
-
02
2
-
7
0.11
D
-
00
0.71
D
-
01
0.43
D
-
01
0.25
D
-
01
0.14
D
-
01
0.73
D
-
02
2
-
8
0.11
D
-
00
0.73
D
-
01
0.45
D
-
01
0.27
D
-
01
0.16
D
-
01
0.87
D
-
02
2
-
9
0.12
D
-
00
0.73
D
-
01
0.47
D
-
01
0.28
D
-
01
0.17
D
-
01
0.94
D
-
02
2
-
10
0.12
D
-
00
0.74
D
-
01
0.47
D
-
01
0.29
D
-
01
0.17
D
-
01
0.98
D
-
02
2
-
11
0.12
D
-
00
0.74
D
-
01
0.47
D
-
01
0.29
D
-
01
0.17
D
-
01
0.10
D
-
01
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
2
-
19
0.12
D
-
00
0.74
D
-
01
0.47
D
-
01
0.29
D
-
01
0.18
D
-
01
0.10
D
-
01
2
-
20
0.12
D
-
00
0.74
D
-
01
0.47
D
-
01
0.29
D
-
01
0.18
D
-
01
0.10
D
-
01
E
N
max
0.12
D
-
00
0.74
D
-
01
0.47
D
-
01
0.29
D
-
01
0.18
D
-
01
0.10
D
-
01Table 2
Convergence rates.
ε
Number of intervalsN
16
32
64
128
256
2
-
2
0.85
0.96
1.06
1.47
0.25
2
-
3
0.78
0.90
1.00
1.14
1.68
2
-
4
0.74
0.87
0.94
1.03
1.29
2
-
5
0.69
0.81
0.89
0.98
1.13
2
-
6
0.64
0.76
0.84
0.93
1.12
2
-
7
0.63
0.71
0.79
0.87
0.91
2
-
8
0.64
0.68
0.74
0.80
0.85
2
-
9
0.66
0.65
0.72
0.77
0.83
2
-
10
0.67
0.64
0.71
0.75
0.80
2
-
11
0.67
0.64
0.71
0.75
0.77
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
2
-
19
0.68
0.64
0.70
0.73
0.75
2
-
20
0.68
0.64
0.70
0.73
0.75
R
N
max
0.68
0.64
0.70
0.73
0.75Figures3 and 5 show the mesh trajectory for τ=0.05 with different N 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 t increases, the mesh points become concentrated in the neighbourhood of x=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 x=1.Figure 1
Numerical solutions forε=2-4 and N=64.Figure 2
Numerical solutions forε=2-14 and N=64.Figure 3
Mesh trajectory withε=2-6, τ=0.05, and N=32.Figure 4
Mesh trajectory withε=2-6, τ=0.5, and N=32.Figure 5
Mesh trajectory forε=2-20, τ=0.05, and N=64.


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



---

*Source: 214505-2013-11-30.xml*