A Moving Mesh Method for Singularly Perturbed Problems

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: 𝑃_𝜀∶𝜕−𝜀²𝜕𝑥²𝑢_𝜀()+𝜕𝑥,𝑡𝑏()𝑢𝜕𝑥𝑥,𝑡_𝜀()+𝜕𝑥,𝑡𝑢𝜕𝑡^𝜀()(),()𝑢𝑥,𝑡=𝑓𝑥,𝑡𝑥,𝑡∈Ω,_𝜀(0,𝑡)=𝑢_𝜀()𝑢1,𝑡=0,_𝜀(𝑥,0)=𝑔(𝑥), (1) with 0<𝜀<1, 𝜆≥𝑏_𝑥≥𝜆>0 and 𝛽≥𝑏≥𝛽>0 for all (𝑥,𝑡)∈Ω=Ω∪Γ_𝐵∪Γ_𝐿∪Γ_𝑅 where 𝜆, 𝜆, 𝛽, 𝛽 are constants, ()×(],ΓΩ=0,10,𝒯_𝐿=)},Γ{(0,𝑡∣0≤𝑡≤𝒯_𝐵=)},Γ{(𝑥,0∣0≤𝑥≤1_𝑅={(1,𝑡)∣0≤𝑡≤𝒯},Γ=Γ_𝐿∪Γ_𝑅∪Γ_𝐵. (2) 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 𝜓(𝑥,𝑡)≥0, for all (𝑥,𝑡)∈Γ. Then 𝐿_𝜀𝜓(𝑥,𝑡)≥0 for all (𝑥,𝑡)∈Ω, this implies that 𝜓(𝑥,𝑡)≥0 for all (𝑥,𝑡)∈Ω.

Proof. Assume that there exists 𝐫=(𝑥_𝑟,𝑡_𝑟)∈Ω such that 𝜓(𝐫)=min_Ω𝜓<0; then 𝐫∉Γ since 𝜓(𝑥,𝑡)≥0 for all (𝑥,𝑡)∈Γ; hence 𝐫∈Ω. Let 𝜓^∗(𝑥,𝑡)=𝜓(𝑥,𝑡)𝑒^{(𝛽/𝜀)(1−𝑥)} (3) for all (𝑥,𝑡)∈Ω. Then 𝜓^∗(𝑥,𝑡)≥0 for all (𝑥,𝑡)∈Γ, and 𝜓(𝐫)<0; thus the minimum of 𝜓^∗ must be also negative. Let 𝐪=(𝑥_𝑞,𝑡_𝑞)∈Ω such that 𝜓^∗(𝐪)=min_Ω𝜓^∗<0. (4) Applying the differential operator 𝐿_𝜀 to 𝜓 gives 𝐿_𝜀𝜓(𝑥,𝑡)=𝐿_𝜀𝜓^∗𝑒^{−(𝛽/𝜀)(1−𝑥)}𝜕=−𝜀²𝜕𝑥²𝜓^∗𝑒^{−(𝛽/𝜀)(1−𝑥)}+𝜕𝜕𝑥𝑏𝜓^∗𝑒^{−(𝛽/𝜀)(1−𝑥)}+𝜕𝜓𝜕𝑡^∗𝑒^{−(𝛽/𝜀)(1−𝑥)}=−𝜀𝜓^∗_𝑥𝑥+𝑏_𝑥𝜓^∗+𝑏𝜓^∗_𝑥+𝜓^∗_𝑡𝑒^{−(𝛽/𝜀)(1−𝑥)}−2𝛽𝜓^∗_𝑥𝑒^{−(𝛽/𝜀)(1−𝑥)}+−𝛽²𝜀+𝑏𝛽𝜀𝜓^∗𝑒^{−(𝛽/𝜀)(1−𝑥)} (5) which can be written as 𝐿_𝜀𝜓(𝑥,𝑡)=𝜂₁+𝜂₂+𝜂₃𝑒^{−(𝛽/𝜀)(1−𝑥)}, (6) where 𝜂₁=−𝜀𝜓^∗_𝑥𝑥+𝑏_𝑥𝜓^∗+𝑏𝜓^∗_𝑥+𝜓^∗_𝑡,𝜂₂=−𝛽²𝜀+𝑏𝛽𝜀𝜓^∗,𝜂₃=−2𝛽𝜓^∗_𝑥. (7) The argument now divides into two cases depending on the position of 𝐪, 𝐪∈Γ_𝑡={(𝑥,𝒯)∣0<𝑥<1} or 𝐪∉Γ_𝑡. If 𝐪∉Γ_𝑡, we have that 𝜓^∗_𝑥𝑥(𝐪)>0,𝜓^∗_𝑥(𝐪)=𝜓^∗_𝑡(𝐪)=0. (8) It can be seen that 𝜂₁<0, 𝜂₂<0 and 𝜂₃=0 for all 𝑥∈(0,1). This leads to the following inequality: 𝐿_𝜀𝜓(𝐪)<0. If 𝐪∈Γ_𝑡, we have that 𝜓^∗_𝑥𝑥(𝐪)>0,𝜓^∗_𝑥(𝐪)=0,𝜓^∗_𝑡(𝐪)≤0. (9) It also follows that 𝜂₁<0, 𝜂₂<0, and 𝜂₃=0 for all (𝑥,𝑡)∈Ω and also leads to the following inequality: 𝐿_𝜀𝜓(𝐪)<0. 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 𝑢_𝜀≤𝑔+𝑓𝛽, (10) for all (𝑥,𝑡)∈Ω.

Proof. Consider the barrier functions 𝜓^±()=𝑔𝑓𝑥,𝑡+𝑥𝛽±𝑢_𝜀(),𝜓𝑥,𝑡_±(0,𝑡)=𝑔±𝑢_𝜀()=𝑔𝜓0,𝑡≥0,_±()=𝑔+𝑓1,𝑡𝛽±𝑢_𝜀()=𝑔+𝑓1,𝑡𝛽𝜓≥0,_±()=𝑔𝑓𝑥,0+𝑥𝛽±𝑢_𝜀(𝑥,0)≥0. (11) For (𝑥,𝑡)∈Ω, 𝐿_𝜀𝜓^±()𝜕𝑥,𝑡=−𝜀²𝜕𝑥²𝑔𝑓+𝑥𝛽±𝑢_𝜀()+𝜕𝑥,𝑡𝑏𝑔𝑓𝜕𝑥+𝑥𝛽±𝑢_𝜀()+𝜕𝑥,𝑡𝑔𝑓𝜕𝑡+𝑥𝛽±𝑢_𝜀(𝑥,𝑡)=𝑏_𝑥𝑔𝑓+𝑥𝛽𝑓+𝑏𝛽±𝑓(𝑥,𝑡)≥0. (12) The maximum principle now applies, and we have 𝜓^±(𝑥,𝑡)≥0 for all (𝑥,𝑡)∈Ω from which we have the required result.

Lemma 3. Let 𝑢_𝜀(𝑥,𝑡) be the solution of 𝑃_𝜀; then the spatial derivatives 𝜕^𝑘𝑢_𝜀(𝑥,𝑡)/𝜕𝑥^𝑘 satisfy the bounds 𝜕^𝑘𝑢_𝜀(𝑥,𝑡)𝜕𝑥^𝑘≤𝐶𝜀^−𝑘max𝑓,𝑢_𝜀𝜕,𝑓𝑜𝑟𝑘=1,2,₃𝑢_𝜀(𝑥,𝑡)𝜕𝑥³≤𝐶𝜀⁻³max𝑓,𝑓^′,𝑢_𝜀. (13)

Proof. Note that ¹_𝑥𝜕𝑏()𝑢𝜕𝑠𝑠,𝑡_𝜀(𝑠,𝑡)𝑑𝑠=𝑏(𝑠,𝑡)𝑢_𝜀(𝑠,𝑡)¹_𝑠=𝑥 (14) which gives ¹_𝑥𝑓()−𝜕𝑠,𝑡𝑏()𝑢𝜕𝑠𝑠,𝑡_𝜀(𝑠,𝑡)𝑑𝑠≤𝑓+𝐶𝑢_𝜀 (15) for all (𝑥,𝑡)∈Ω where 𝐶 depends on ‖𝑏‖. From the mean-value theorem, there exists a point 𝑧∈(1−𝜀,1) such that 𝜕𝑢_𝜀(𝑧,𝑡)=𝑢𝜕𝑥^𝜀(1,𝑡)−𝑢_𝜀(1−𝜀,𝑡)𝜀; (16) hence 𝜀𝜕𝑢_𝜀(𝑧,𝑡)≤𝑢𝜕𝑥^𝜀. (17) Integrating 𝑃_𝜀 with respect to 𝑥, we obtain −𝜀𝜕𝑢_𝜀(1,𝑡)−𝜕𝑥𝜕𝑢_𝜀(𝑥,𝑡)=𝜕𝑥_1𝑥𝑓()−𝜕𝑠,𝑡𝑏()𝑢𝜕𝑠𝑠,𝑡_𝜀(𝑠,𝑡)−𝜕𝑢_𝜀(𝑠,𝑡)𝜕𝑡𝑑𝑠. (18) Using (18) with 𝑥=𝑧 and combining with (15), it follows that 𝜀𝜕𝑢_𝜀(1,𝑡)≤𝑓𝑢𝜕𝑥+𝐶_𝜀. (19) Equation (15) can also be used to give the following bound: 𝜀𝜕𝑢_𝜀(𝑥,𝑡)≤𝑓𝑢𝜕𝑥+𝐶_𝜀, (20) for all (𝑥,𝑡)∈Ω. This proves the result for 𝑘=1. To obtain bounds for the higher derivatives, rewriting (10) 𝜕−𝜀²𝑢_𝜀(𝑥,𝑡)𝜕𝑥²()−𝜕=𝑓𝑥,𝑡𝑏()𝑢𝜕𝑥𝑥,𝑡_𝜀(𝑥,𝑡)−𝜕𝑢_𝜀(𝑥,𝑡),𝜕𝑡 (21) this gives the second derivative bound where 𝐶 depends on ‖𝑔‖, ‖𝑏‖, and ‖𝑏_𝑥‖. Differentiating 𝑃_𝜀 with respect to 𝑥𝜕−𝜀³𝑢_𝜀(𝑥,𝑡)𝜕𝑥³=𝜕𝑓()−𝜕𝜕𝑥𝑥,𝑡𝑏()𝑢𝜕𝑥𝑥,𝑡_𝜀()−𝑥,𝑡𝜕𝑢_𝜀(𝑥,𝑡)𝜕𝑡 (22) 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) where 𝑣_𝜀(𝑥,𝑡) is the solution to problem𝐿_𝜀𝑣_𝜀()(),()𝑣𝑥,𝑡=𝑓𝑥,𝑡𝑥,𝑡∈Ω,_𝜀(𝑥,𝑡)=𝑢_𝜀(𝑥,𝑡),(𝑥,𝑡)∈Γ_𝐵∪Γ_𝐿, (24a)(24b)and the singular component 𝑤_𝜀(𝑥,𝑡) is the solution of the homogenous problem 𝐿_𝜀𝑤_𝜀()()𝑤𝑥,𝑡=0,𝑥,𝑡∈Ω,_𝜀(𝑥,𝑡)=0,(𝑥,𝑡)∈Γ_𝐵∪Γ_𝐿,𝑤_𝜀(𝑥,𝑡)=𝑢_𝜀(𝑥,𝑡)−𝑣_𝜀(𝑥,𝑡),(𝑥,𝑡)∈Γ_𝑅. (25a)(25b)(25c)

Theorem 4. The solution 𝑢_𝜀(𝑥,𝑡) of the continuous problem 𝑃_𝜀 can be decomposed as a sum of the smooth and layer functions 𝑢_𝜀(𝑥,𝑡)=𝑣_𝜀(𝑥,𝑡)+𝑤_𝜀(𝑥,𝑡),(𝑥,𝑡)∈Ω, (26) where for all 𝑘, 0≤𝑘≤3, the smooth component 𝑣_𝜀(𝑥,𝑡) satisfies 𝜕^𝑘𝑣_𝜀(𝑥,𝑡)𝜕𝑥^𝑘≤𝐶1+𝜀^{−(𝑘−2)}; (27) the singular component 𝑤_𝜀(𝑥,𝑡) satisfies 𝑤_𝜀(𝑥,𝑡)≤𝐶𝑒^{−(𝛽/𝜀)(1−𝑥)},∀𝑥∈𝜕Ω,_𝑘𝑤_𝜀(𝑥,𝑡)𝜕𝑥^𝑘≤𝐶𝜀^−𝑘𝑒^{−(𝛽/𝜀)(1−𝑥)},𝑘=1,2,3 (28)(29) for some constant 𝐶 independent of 𝜀.

Proof. To find these bounds, we rewrite the smooth component as 𝑣_𝜀(𝑥,𝑡)=𝑣₀+𝜀𝑣₁+𝜀²𝑣₂(𝑥,𝑡),(𝑥,𝑡)∈Ω, (30) where 𝑣₀(𝑥,𝑡) is the solution to the reduced problem 𝑃₀, 𝑣₁(𝑥,𝑡) satisfies 𝜕𝑏()𝑣𝜕𝑥𝑥,𝑡₁(𝑥,𝑡)+𝜕𝑣₁(𝑥,𝑡)=𝜕𝜕𝑡²𝑣₀(𝑥,𝑡)𝜕𝑥²,()𝑣𝑥,𝑡∈Ω,₁(𝑥,𝑡)=0,(𝑥,𝑡)∈Γ_𝐵∪Γ_𝐿, (31) and 𝑣₂(𝑥,𝑡) satisfies 𝐿_𝜀𝑣₂()=𝜕𝑥,𝑡²𝑣₁(𝑥,𝑡)𝜕𝑥²,()𝑣𝑥,𝑡∈Ω,₂(𝑥,𝑡)=0,(𝑥,𝑡)∈Γ. (32) We clearly have 𝐿_𝜀𝑣_𝜀(𝑥,𝑡)=𝑓(𝑥) in Ω with 𝑣_𝜀(𝑥,𝑡)=𝑣₀(𝑥,𝑡)+𝜀𝑣₁(𝑥,𝑡) on Γ_𝑅. It can be easily seen that 𝑣₀(𝑥,𝑡) and 𝑣₁(𝑥,𝑡) are all bounded by a constant independent of 𝜀 and 𝜕^𝑖+𝑗𝑣₁(𝑥,𝑡)𝜕𝑥^𝑖𝜕𝑡^𝑗≤𝐶,(𝑥,𝑡)∈Ω (33) for 0≤𝑖+2𝑗≤4. Therefore 𝑣₂(𝑥,𝑡) is the solution of a problem similar to 𝑃_𝜀; hence from Lemma 2, we obtain that 𝑣₂≤𝐶⟹𝑣_𝜀≤𝐶1+𝜀+𝜀². (34) 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 0≤𝑘≤3𝜕^𝑘𝑣_𝜀(𝑥,𝑡)𝜕𝑥^𝑘≤𝐶1+𝜀^{−(𝑘−2)}. (35)
To find bounds for the singular component, we consider the following mesh functions: 𝜓^±(𝑥,𝑡)=𝐶𝑒^{−(𝛽/𝜀)(1−𝑥)}±𝑤_𝜀(𝑥,𝑡), (36) where 𝐶=‖𝑣₀‖+‖𝑣₁‖. 𝜓^±(0,𝑡)=𝐶𝑒^{−(𝛽/𝜀)}±𝑤_𝜀(0,𝑡)=𝐶𝑒^{−(𝛽/𝜀)}𝜓≥0,_±(1,𝑡)=𝐶±𝑤_𝜀()𝜓1,𝑡≥0,_±(𝑥,0)=𝐶𝑒^{−(𝛽/𝜀)(1−𝑥)}±𝑤_𝜀(𝑥,0)=𝐶𝑒^{−(𝛽/𝜀)(1−𝑥)}≥0. (37) Applying the differential operator 𝐿_𝜀, we obtain that 𝐿_𝜀𝜓^±(𝑥,𝑡)=𝐶𝑏𝛽𝜀−𝛽²𝜀+𝑏_𝑥𝑒^{−(𝛽/𝜀)(1−𝑥)}≥0. (38) Thus from the maximum principle, we can say that 𝜓^±(𝑥,𝑡)≥0; hence 𝑤_𝜀(𝑥,𝑡)≤𝐶𝑒^{−(𝛽/𝜀)(1−𝑥)}. (39) To establish the bounds for the derivatives of 𝑤_𝜀(𝑥,𝑡), we start by noting that ¹_𝑥𝜕𝑏()𝑤𝜕𝑠𝑠,𝑡_𝜀(𝑠,𝑡)𝑑𝑠≤𝐶𝑤_𝜀(𝑥,𝑡)≤𝐶𝑒^{−(𝛽/𝜀)(1−𝑥)}, (40) where 𝐶 depends on ‖𝑏‖. From the mean value theorem, there exists a point 𝑧∈(1−𝜀,1) such that 𝜕𝑏()𝑤𝜕𝑥𝑧,𝑡_𝜀()=𝑏()𝑤𝑧,𝑡1,𝑡_𝜀(1,𝑡)−𝑏(1−𝜀,𝑡)𝑤_𝜀(1−𝜀,𝑡)𝜀. (41) Using the triangle inequality, it can easily be seen that 𝜕𝑤_𝜀(𝑧,𝑡)≤1𝜕𝑥|𝑏(×()𝑧,𝑡)|𝜕𝑏𝑧,𝑡𝑤𝜕𝑥^𝜀(+𝑏()𝑤𝑧,𝑡)1,𝑡_𝜀(1,𝑡)−𝑏(1−𝜀,𝑡)𝑤_𝜀(1−𝜀,𝑡)𝜀. (42) Using (39) at the point 𝑥=1−𝜀,1 and 𝑥=𝑧, we obtain that 𝜕𝑤_𝜀(𝑧,𝑡)𝜕𝑥≤𝐶𝜀⁻¹𝑒^{−(𝛽/𝜀)(1−𝑧)}. (43) Integrating (25a), (25b), and (25c) with respect to 𝑥, we obtain that 𝜀𝜕𝑤_𝜀(1,𝑡)𝜕𝑥=𝜀𝜕𝑤_𝜀(𝑥,𝑡)+𝜕𝑥_1𝑥𝜕𝑏()𝑤𝜕𝑠𝑠,𝑡_𝜀(𝑠,𝑡)+𝜕𝑤_𝜀(𝑠,𝑡)𝜕𝑡𝑑𝑠 (44) at the point 𝑥=𝑧𝜀𝜕𝑤_𝜀(1,𝑡)=𝜀𝜕𝑥𝜕𝑤_𝜀(𝑧,𝑡)+𝜕𝑥_1𝑧𝜕𝑏()𝑤𝜕𝑠𝑠,𝑡_𝜀(𝑠,𝑡)+𝜕𝑤_𝜀(𝑠,𝑡).𝜕𝑡𝑑𝑠 (45) Hence 𝜀𝜕𝑤_𝜀(1,𝑡)𝜕𝑥≤𝐶𝑒^{−(𝛽/𝜀)(1−𝑧)}. (46) From (44) and (46), we obtain that 𝜀𝜕𝑤_𝜀(𝑥,𝑡)𝜕𝑥≤𝐶𝑒^{−(𝛽/𝜀)(1−𝑥)},∀(𝑥,𝑡)∈Ω (47) which is the required bound for the first derivative. From (25a), (25b), and (25c), it can be seen that 𝜀𝜕²𝑤_𝜀(𝑥,𝑡)𝜕𝑥²=𝜕𝑏()𝑤𝜕𝑥𝑥,𝑡_𝜀(𝑥,𝑡)+𝜕𝑤_𝜀(𝑥,𝑡),𝜀𝜕𝜕𝑡₃𝑤_𝜀(𝑥,𝑡)𝜕𝑥³=𝜕𝜕𝜕𝑥𝑏()𝑤𝜕𝑥𝑥,𝑡_𝜀(𝑥,𝑡)+𝜕𝑤_𝜀(𝑧,𝑡),𝜕𝑡 (48) 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 𝑑𝑢=𝑑𝑡𝑑𝑥𝑑𝑡𝜕𝑢+𝜕𝜕𝑥()𝜕𝑥𝒜𝑢+𝑓. (49) Discretizing (49) using an implicit scheme, 𝑢_𝑖,𝑗+𝛿𝑡𝑥̇_{𝑖,𝑗−1}𝐷⁻_𝑥+𝐷⁻_𝑥𝒜𝑢_𝑖,𝑗=𝑢_{𝑖,𝑗−1}+𝛿𝑡𝑓𝑥_𝑖,𝑗 (50) which can be written as a system 𝑇𝑢=Ψ. (51)𝑢_𝑖,𝑗 represents the solution at the point (𝑥_𝑖,𝑡_𝑗), 𝐷⁻_𝑥𝑢_𝑖,𝑗=𝑢_𝑖,𝑗−𝑢_{𝑖−1,𝑗}ℎ_𝑖,𝑗,𝐷⁺_𝑥𝑢_𝑖,𝑗=𝑢_𝑖+1,𝑗−𝑢_𝑖,𝑗ℎ_𝑖+1,𝑗. (52)ℎ_𝑖,𝑗=𝑥_𝑖,𝑗−𝑥_{𝑖−1,𝑗}, 𝑢_{𝑖,𝑗−1} 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] 𝑀_𝑖+1,𝑗+𝑀_𝑖,𝑗𝑥̇_𝑖+1,𝑗−𝑥̇_𝑖,𝑗−𝑀_𝑖,𝑗+𝑀_{𝑖−1,𝑗}𝑥̇_𝑖,𝑗−𝑥̇_{𝑖−1,𝑗}𝐸=−^{𝑖,𝑗−1}𝜏, (53) where 𝐸_{𝑖,𝑗−1}=𝑀_{𝑖+1,𝑗−1}+𝑀_{𝑖,𝑗−1}𝑥_{𝑖+1,𝑗−1}−𝑥_{𝑖,𝑗−1}−𝑀_{𝑖,𝑗−1}+𝑀_{𝑖−1,𝑗−1}𝑥_{𝑖,𝑗−1}−𝑥_{𝑖−1,𝑗−1}. (54)𝜏 is a user defined parameter, which determines how fast the grid moves towards the equidistributed grid, and 𝑀 is the arc length monitor function 𝑀=1+𝜕𝑢𝜕𝑥². (55)

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, 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_𝑖,𝑗|𝑥̇_𝑖,𝑗|<𝐶₁, since Ψ_𝑖,0=𝑢_𝑖,0+𝛿𝑡𝑓𝑥_𝑖,1⟹Ψ_𝑖,0≤𝐶+𝛿𝑡𝑓=𝐶₂<∞. (56) Assume that this result is true for some time 𝑡=𝑡_𝑛−1; hence |Ψ_{𝑖,𝑛−1}|≤𝐶₂<∞. Since the matrix 𝑇 is invertible at 𝑡=𝑡_𝑛=𝒯, it follows that 𝑢_𝑖,𝑛 exists for all 𝑖 and is bounded by some constant 𝐶, max_𝑖,𝑛|𝑢_𝑖,𝑛|≤𝐶.

Define the discrete operator 𝐿^𝑁_𝜀 as 𝐿^𝑁_𝜀𝑢_𝑖,𝑗≡−𝜀𝛿²_𝑥+𝐷⁻_𝑥𝑏_𝑖,𝑗+𝐷⁻_𝑡−𝑥̇_𝑖,𝑗𝐷⁻_𝑥𝑢_𝑖,𝑗, (57) where 𝛿²_𝑥𝑢_𝑖,𝑗=2ℎ_𝑖+1,𝑗+ℎ_𝑖,𝑗𝐷⁺_𝑥𝑢_𝑖,𝑗−𝐷⁻_𝑥𝑢_𝑖,𝑗 (58) 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 𝐫=(𝑥_𝑖,𝑡_𝑗)∈Ω^𝑁 such that 𝜓_𝑖,𝑗=𝜓(𝐫)=min_Ω𝑁𝜓<0; (59) then 𝐫∉Γ^𝑁, which implies that 𝐫∈Ω^𝑁, and we also know that 𝛿²_𝑥𝜓_𝑖,𝑗≥0,𝐷⁻_𝑥𝜓_𝑖,𝑗≤0,𝐷⁻_𝑡𝜓_𝑖,𝑗≤0. (60) We have to show that 𝐷⁻_𝑥(𝑏_𝑖,𝑗𝜓_𝑖,𝑗)≤0, and we proceed by contradiction, suppose that 𝐷⁻_𝑥𝑏_𝑖,𝑗𝜓_𝑖,𝑗>0⟹𝑏_𝑖,𝑗𝜓_𝑖,𝑗−𝑏_{𝑖−1,𝑗}𝜓_{𝑖−1,𝑗}>0⟹𝑏_𝑖,𝑗𝜓_𝑖,𝑗>𝑏_{𝑖−1,𝑗}𝜓_{𝑖−1,𝑗}. (61) Since 𝐫 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 −𝜀𝛿²_𝑥𝜓_𝑖,𝑗≤0 and 𝐷⁻_𝑥(𝑏_𝑖,𝑗𝜓_𝑖,𝑗)≤0; hence 𝜓_𝑖+1,𝑗=𝜓_𝑖,𝑗=𝜓_{𝑖−1,𝑗}<0. (62) 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.

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 𝑈_𝜀𝑥_𝑖,𝑡_𝑗=𝑉_𝜀𝑥_𝑖,𝑡_𝑗+𝑊_𝜀𝑥_𝑖,𝑡_𝑗,𝑥_𝑖,𝑡_𝑗∈Ω^𝑁, (63) where 𝑉_𝜀(𝑥_𝑖,𝑡_𝑗) and 𝑊_𝜀(𝑥_𝑖,𝑡_𝑗) are solutions of the respective equations𝐿^𝑁_𝜀𝑉_𝜀𝑥_𝑖,𝑡_𝑗=𝑓𝑥_𝑖,𝑡_𝑗,𝑥_𝑖,𝑡_𝑗∈Ω^𝑁𝑉_𝜀,𝑥_𝑖,𝑡_𝑗=𝑣_𝜀𝑥_𝑖,𝑡_𝑗,𝑥_𝑖,𝑡_𝑗∈Γ^𝑁, (64a)(64b)𝐿^𝑁_𝜀𝑊_𝜀𝑥_𝑖,𝑡_𝑗=0,𝑥_𝑖,𝑡_𝑗∈Ω^𝑁𝑊_𝜀,𝑥_𝑖,𝑡_𝑗=𝑤_𝜀𝑥_𝑖,𝑡_𝑗,𝑥_𝑖,𝑡_𝑗∈Γ^𝑁, (65a)(65b)where 𝑊_𝜀 is the singular component and 𝑉_𝜀 is the smooth component. The error for the numerical solution will be decomposed as 𝑈_𝜀−𝑢_𝜀𝑥_𝑖,𝑡_𝑗=𝑉_𝜀−𝑣_𝜀+𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗,𝑥_𝑖,𝑡_𝑗∈Ω^𝑁. (66) 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 𝜓(𝑥,𝑡)∈𝐶²(Ω), 𝐷⁻_𝑥−𝜕𝑏𝑥𝜕𝑥^𝑖,𝑡_𝑗𝜓𝑥_𝑖,𝑡_𝑗≤12ℎ_𝑖,𝑗𝜕²𝜕𝑥²𝑏(𝑥,𝑡)𝜓(𝑥,𝑡). (67)

Proof. Using integration by parts to reduce the order of differentiation in the integral, (𝐷⁻_𝑥−𝜕/𝜕𝑥)(𝑏(𝑥_𝑖,𝑡_𝑗)𝜓(𝑥_𝑖,𝑡_𝑗)) can be expressed as 𝐷⁻_𝑥−𝜕𝑏𝑥𝜕𝑥^𝑖,𝑡_𝑗𝜓𝑥_𝑖,𝑡_𝑗=1ℎ_𝑖,𝑗^{𝑥𝑖,𝑗}_{𝑥𝑖−1,𝑗}𝑥_{𝑖−1,𝑗}𝜕−𝑠²𝜕𝑠²𝑏(𝑠,𝑡)𝜓(𝑠,𝑡)𝑑𝑠. (68) It follows that 𝐷⁻_𝑥−𝜕𝑥𝜕𝑥𝑏𝜓_𝑖,𝑡_𝑗≤1ℎ_𝑖,𝑗𝜕²𝜕𝑥²𝑏(𝑥,𝑡)𝜓(𝑥,𝑡)^{𝑥𝑖,𝑗}_{𝑥𝑖−1,𝑗}𝑠−𝑥_{𝑖−1,𝑗}≤1𝑑𝑠2ℎ_𝑖,𝑗𝜕²𝜕𝑥²𝑏(𝑥,𝑡)𝜓(𝑥,𝑡). (69)

Lemma 7. Assuming the bound (27), the error in the smooth component of the numerical solution satisfies the following error bound: 𝑉_𝜀−𝑣_𝜀𝑥_𝑖,𝑡_𝑗≤𝐶ℎ_max,𝑗+𝛿𝑡, (70) where ℎ_max,𝑗=max_𝑖{ℎ_𝑖,𝑗∣1≤𝑖≤𝑁} for some constant 𝐶 independent of 𝜀 and 𝑁.

Proof. We start by considering the local truncation error for the smooth component 𝐿^𝑁_𝜀𝑉_𝜀−𝑣_𝜀𝑥_𝑖,𝑡_𝑗=𝐿_𝜀−𝐿^𝑁_𝜀𝑣_𝜀𝑥_𝑖,𝑡_𝑗=𝜕−𝜀²𝜕𝑥²−𝛿²_𝑥+𝜕𝜕𝑥−𝐷⁻_𝑥𝑏𝑥_𝑖,𝑡_𝑗+𝜕−𝐷𝜕𝑡^−𝑡−𝐷⁻_𝑡𝑥_𝑖,𝑗𝐷⁻_𝑥𝑣_𝜀𝑥_𝑖,𝑡_𝑗. (71) It can be easily shown that 𝐷⁻_𝑥𝑣_𝜀𝑥_𝑖,𝑡_𝑗=1𝑥_𝑖,𝑗−𝑥_{𝑖−1,𝑗}^{𝑥𝑖,𝑗}_{𝑥𝑖−1,𝑗}𝜕𝑣_𝜀(𝑠,𝑡)≤𝜕𝑠𝑑𝑠𝜕𝑣_𝜀(𝑥,𝑡)(),𝐷𝜕𝑥≤𝐶1+𝜀_−𝑡𝑥𝑡_𝑗=1𝑡_𝑗−𝑡_𝑗−1^𝑡𝑗_𝑡𝑗−1𝑑𝑥(𝑠)≤(𝑡)𝑑𝑠𝑑𝑠𝑑𝑥𝑑𝑡≤𝐶. (72) Using the local truncation error estimates from [16], (69) (72), it follows that 𝐿^𝑁_𝜀𝑉_𝜀−𝑣_𝜀𝑥_𝑖,𝑡_𝑗≤𝜀3ℎ_𝑖+1,𝑗+ℎ_{𝑖−1,𝑗}𝜕³𝑣_𝜀(𝑥,𝑡)𝜕𝑥³+12ℎ_𝑖,𝑗𝜕²𝜕𝑥²𝑏𝑣_𝜀(+1𝑥,𝑡)2𝑡_𝑗−𝑡_𝑗−1𝜕²𝑣_𝜀(𝑥,𝑡)𝜕𝑡²+𝜕𝑣_𝜀(𝑥,𝑡)(𝑡)𝜕𝑥𝑑𝑥ℎ𝑑𝑡≤𝐶_max,𝑗+𝛿𝑡. (73) We introduce the two mesh functions 𝜓^±_𝑖,𝑗=𝐶𝛽⁻¹ℎ_max,𝑗+𝛿𝑡𝛽+𝑥_𝑖,𝑗±𝑉_𝜀−𝑣_𝜀𝑥_𝑖,𝑡_𝑗 (74) it can be seen from the mesh functions that 𝜓^±_0,𝑗=𝐶ℎ_max,𝑗𝜓+𝛿𝑡≥0,_±𝑁,𝑗=𝐶𝛽⁻¹ℎ_max,𝑗𝜓+𝛿𝑡𝛽+1≥0,_±𝑖,0=𝐶𝛽⁻¹𝑁⁻¹+𝛿𝑡𝛽+𝑥_𝑖,0𝐿≥0,_𝑁𝜀𝜓^±_𝑖,𝑗=𝐶𝛽⁻¹ℎ_max,𝑗𝛽𝑏+𝛿𝑡_𝑖,𝑗−𝑏_{𝑖−1,𝑗}ℎ_𝑖,𝑗+𝑏_𝑖,𝑗𝑥_𝑖,𝑗−𝑏_{𝑖−1,𝑗}𝑥_{𝑖−1,𝑗}ℎ_𝑖,𝑗±𝐿^𝑁_𝜀𝑉_𝜀−𝑣_𝜀𝑥_𝑖,𝑡_𝑗≥𝐶ℎ_max,𝑗𝑏+𝛿𝑡^𝑖,𝑗−𝑏_{𝑖−1,𝑗}ℎ_𝑖,𝑗+𝑏_𝑖,𝑗𝛽±𝐿^𝑁𝑉_𝜀−𝑣_𝜀𝑥_𝑖,𝑡_𝑗≥0. (75) By the discrete maximum principle, we conclude that 𝜓^±_𝑖,𝑗≥0 and so for all (𝑥_𝑖,𝑡_𝑗)∈Ω^𝑁; hence 𝑉_𝜀−𝑣_𝜀𝑥_𝑖,𝑡_𝑗≤𝐶ℎ_max,𝑗+𝛿𝑡. (76)

Lemma 8. For all 𝑁≥4 and at each (𝑥_𝑖,𝑡_𝑗)∈Ω^𝑁, the singular component of the error satisfies 𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗≤𝐶𝑁⁻¹(ln𝑁)²+𝛿𝑡 (77) 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 𝜀⁻¹≤2/𝛽ln1/𝜀 and 𝜀≥𝑁⁻¹; in this case 𝜀⁻¹≤2/𝛽ln𝑁, and our mesh is quasiuniform ℎ_min,𝑗∼ℎ_max,𝑗. Using the standard bound for the local truncation error, we can derive an equivalent expression for the truncation error as for the smooth component 𝐿^𝑁𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗≤𝜀3ℎ_𝑖+1,𝑗+ℎ_{𝑖−1,𝑗}𝜕³𝑤_𝜀(𝑥,𝑡)𝜕𝑥³+12ℎ_𝑖,𝑗𝜕²𝜕𝑥²𝑏𝑤_𝜀(+1𝑥,𝑡)2𝑡_𝑗−𝑡_𝑗−1𝜕²𝑤_𝜀(𝑥,𝑡)𝜕𝑡²+𝜕𝑤_𝜀(𝑥,𝑡)(𝑡)𝜕𝑥𝑑𝑥𝑑𝑡≤𝐶𝜀⁻²ℎ_𝑖+1,𝑗+ℎ_𝑖,𝑗𝑒^{−(𝛽/𝜀)(1−𝑥)}+𝐶𝜀⁻²𝛽_maxℎ_𝑖,𝑗𝑒^{−(𝛽/𝜀)(1−𝑥)}+𝐶𝑡_𝑗−𝑡_𝑗−1+𝐶𝜀⁻¹𝑒^{−(𝛽/𝜀)(1−𝑥)}≤𝐶𝜀⁻²ℎ_max,𝑗1+𝛽_max𝑁+𝐶𝛿𝑡≤𝐶⁻¹(ln𝑁)²+𝛿𝑡, (78) where 𝛽_max=max{‖𝑏‖,‖𝑏_𝑥‖,‖𝑏_𝑥𝑥‖}. Consider the following barrier function: 𝜓_𝑖,𝑗=𝐶𝛽⁻¹𝑁⁻¹(ln𝑁)²+𝛿𝑡𝛽+𝑥_𝑖,𝑗 (79) and the mesh functions 𝜙^±_𝑖,𝑗=𝜓_𝑖,𝑗±𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗. (80) It can be easily seen that 𝜙^±_0,𝑗=𝐶𝑁⁻¹(ln𝑁)²𝜙+𝛿𝑡≥0,_±𝑁,𝑗=𝐶𝛽⁻¹𝑁⁻¹(ln𝑁)²𝜙+𝛿𝑡𝛽+1≥0,_±𝑖,0=𝐶𝛽⁻¹𝑁⁻¹(ln𝑁)²+𝛿𝑡𝛽+𝑥_𝑖,0≥0. (81) Applying the difference operator 𝐿^𝑁_𝜀 to the barrier function, 𝐿^𝑁_𝜀𝜓_𝑖,𝑗=𝐶𝛽⁻¹𝑁⁻¹(ln𝑁)²×𝛽𝑏+𝛿𝑡_𝑖,𝑗−𝑏_{𝑖−1,𝑗}ℎ_𝑖,𝑗+𝑏_𝑖,𝑗𝑥_𝑖,𝑗−𝑏_{𝑖−1,𝑗}𝑥_{𝑖−1,𝑗}ℎ_𝑖,𝑗≥𝐶𝑁⁻¹(ln𝑁)²𝑏+𝛿𝑡^𝑖,𝑗−𝑏_{𝑖−1,𝑗}ℎ_𝑖,𝑗+𝑏_𝑖,𝑗𝛽≥𝐿^𝑁_𝜀𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗. (82) Hence it follows that 𝐿^𝑁_𝜀𝜙^±_𝑖,𝑗≥0 for 1≤𝑖≤𝑁−1; by the discrete maximum principle, it follows that 𝜙^±_𝑖,𝑗≥0; hence 𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗≤𝜓_𝑖,𝑗≤𝐶𝑁⁻¹(ln𝑁)²+𝛿𝑡. (83)
To consider the case when 𝜀⁻¹>2/𝛽ln1/𝜀 and 𝜀<𝑁⁻¹, we will assume that ℎ_𝑖,𝑗≤ℎ_{𝑖−1,𝑗} and that there exists a point 𝑥_𝛼∈[0,1] such that 1−𝑥_𝛼≥𝜀/𝛽ln1/𝜀. The point 𝑥_𝛼 splits the interval into the coarse mesh Ω^𝑁_𝛼𝑙 and fine mesh Ω^𝑁_𝛼𝑟 where Ω^𝑁_𝛼𝑙=𝑥_𝑖∣0≤𝑥_𝑖≤𝑥_𝛼,ℎ_𝑖,𝑗≥1𝑁,Ω^𝑁_𝛼𝑟=𝑥_𝑖∣𝑥_𝛼<𝑥_𝑖≤1,ℎ_𝑖,𝑗<1𝑁. (84) 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 𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗≤𝑊_𝜀𝑥_𝑖,𝑡_𝑗+𝑤_𝜀𝑥_𝑖,𝑡_𝑗. (85) It suffices to bound 𝑊_𝜀 and 𝑤_𝜀 separately. Using (28) for 𝑥_𝑖∈Ω^𝑁_𝛼𝑙, we obtain that 𝑤_𝜀𝑥_𝑖,𝑡_𝑗≤𝐶𝑒^{−(𝛽/𝜀)(1−𝑥𝑖)}≤𝐶𝑒^{−(𝛽/𝜀)(1−𝑥𝛼)}≤𝐶𝑒^−ln1/𝜀=𝐶𝜀≤𝐶𝑁⁻¹. (86) To derive a similar bound on 𝑊_𝜀, we introduce the mesh function 𝑌𝑌𝑥_𝑖,𝑡_𝑗=𝑌_𝑖,𝑗=𝐶𝛽𝑁𝛽+𝑥_𝑖,𝑗. (87) It can easily be seen that 𝑌_0,𝑗=𝐶𝑁𝑌≥0,_𝑁,𝑗=𝐶𝑌𝛽𝑁𝛽+1≥0,_𝑖,0=𝐶𝛽𝑁𝛽+𝑥_𝑖,0𝐷≥0,_−𝑥(𝑏𝑌)𝑥_𝑖,𝑡_𝑗=𝐶𝛽𝑏𝛽𝑁^𝑖,𝑗−𝑏_{𝑖−1,𝑗}ℎ_𝑖,𝑗+𝑏_𝑖,𝑗𝑥_𝑖,𝑗−𝑏_{𝑖−1,𝑗}𝑥_{𝑖−1,𝑗}ℎ_𝑖,𝑗≥0. (88) Also note that 𝛿²_𝑥𝑌(𝑥_𝑖,𝑡_𝑗)=0 and 𝐷⁻_𝑡𝑥_𝑖,𝑗−𝐷⁻_𝑡𝑥_𝑖,𝑗𝐷⁻_𝑥𝑥_𝑖,𝑗=0. Applying the operator 𝐿^𝑁_𝜀 to 𝑌, 𝐿^𝑁_𝜀𝑌𝑥_𝑖,𝑡_𝑗=−𝜀𝛿²_𝑥𝑌_𝑖,𝑗+𝐷⁻(𝑏𝑌)_𝑖,𝑗+𝐷⁻_𝑡𝑌_𝑖,𝑗−𝐷⁻_𝑡𝑌_𝑖,𝑗𝐷⁻_𝑥𝑌_𝑖,𝑗≥0. (89) Now considering the mesh function 𝑌_𝜀−𝑊_𝜀, we have 𝐿^𝑁_𝜀𝑌−𝑊_𝜀𝑥_𝑖,𝑡_𝑗=𝐿^𝑁_𝜀𝑌_𝑖,𝑗−𝐿^𝑁𝑊_𝑖,𝑗=𝐿^𝑁_𝜀𝑌_𝑖,𝑗≥0. (90) From the maximum principle, this means that 𝑌−𝑊_𝜀≥0; hence 𝑊_𝑖,𝑗≤𝑌_𝑖,𝑗=𝐶𝛽𝑥_𝑖𝑁≤𝐶𝑁⁻¹. (91)
It only now remains to give a bound for (𝑊_𝜀−𝑤_𝜀) in the interval (𝑥_𝛼,1]. If 𝑥_𝑖∈Ω^𝑁_𝛼𝑟, the expression for the truncation error for the singular component (78) does not change since in the boundary layer region ℎ_max,𝑗≤1/𝑁. From (65a) and (65b), we have that 𝑊_𝜀1,𝑡_𝑗−𝑤_𝜀1,𝑡_𝑗=0. (92) From (86) and (91), it follows that 𝑊_𝜀𝑥_𝛼,𝑡_𝑗−𝑤_𝜀𝑥_𝛼,𝑡_𝑗≤𝑊_𝜀𝑥_𝛼,𝑡_𝑗+𝑤_𝜀𝑥_𝛼,𝑡_𝑗≤𝐶𝑁⁻¹. (93) Consider the barrier functions 𝜓_𝑖,𝑗=𝐶𝑁⁻¹(ln𝑁)²𝑥+𝛿𝑡^𝑖,𝑗𝛽+1 (94) and the mesh function 𝜙^±_𝑖,𝑗=𝜓_𝑖,𝑗±𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗, (95) for 𝑥_𝑖=𝑥_𝛼 and 𝑥_𝑖=1, 𝜙^±𝑥_𝛼,𝑡_𝑗=𝐶𝑁⁻¹(ln𝑁)²𝑥+𝛿𝑡^𝛼𝛽±𝑊+1_𝜀−𝑤_𝜀𝑥_𝛼,𝑡_𝑗𝜙≥0,_±𝑁,𝑗=𝐶𝑁⁻¹(ln𝑁)²1+𝛿𝑡𝛽𝜙+1≥0,_±𝑖,0=𝐶𝑁⁻¹(ln𝑁)²𝑥+𝛿𝑡^𝑖,0𝛽+1≥0. (96) Applying the discrete difference operator 𝐿^𝑁_𝜀 to 𝜓_𝑖,𝑗 and noting that 𝐷⁻_𝑡𝑥_𝑖,𝑗−𝐷⁻_𝑡𝑥_𝑖,𝑗𝐷⁻_𝑥𝑥_𝑖,𝑗=0, 𝐿^𝑁𝜓_𝑖,𝑗=𝐶𝑁⁻¹(ln𝑁)²×𝑏+𝛿𝑡_𝑖,𝑗𝑥_𝑖,𝑗−𝑏_{𝑖−1,𝑗}𝑥_{𝑖−1,𝑗}𝛽ℎ_𝑖,𝑗+𝑏_𝑖,𝑗−𝑏_{𝑖−1,𝑗}ℎ_𝑖,𝑗≥𝐶𝑁⁻¹(ln𝑁)²𝑏+𝛿𝑡^𝑖𝛽+𝑏_𝑖,𝑗−𝑏_{𝑖−1,𝑗}ℎ_𝑖,𝑗≥𝐿^𝑁𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗; (97) hence it follows that 𝐿^𝑁_𝜀𝜙^±_𝑖,𝑗≥0,𝑥_𝑖∈Ω^𝑁_𝛼𝑟. (98) The discrete maximum principle on Ω^𝑁_𝛼𝑟 then gives 𝜙^±_𝑖,𝑗≥0,𝑥_𝑖∈Ω^𝑁_𝛼𝑟, (99) and it follows that 𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗≤𝜓_𝑖,𝑗≤𝐶𝑁⁻¹(ln𝑁)²+𝛿𝑡. (100) Combining the separate estimates in the two subintervals gives 𝑊_𝜀−𝑤_𝜀𝑥_𝑖,𝑡_𝑗≤𝐶𝑁⁻¹(ln𝑁)²+𝛿𝑡,𝑥_𝑖,𝑡_𝑗∈Ω^𝑁. (101)

Theorem 9. If 𝑢_𝜀(𝑥,𝑡) is the solution of 𝑃_𝜀 and 𝑈_𝜀(𝑥_𝑖,𝑡_𝑗), the corresponding numerical solution using the method outlined in (50) satisfies sup_0<𝜀<1𝑈_𝜀−𝑢_𝜀≤𝐶𝑁⁻¹(ln𝑁)²+𝛿𝑡 (102) for 𝑁≥4 where 𝐶 is a constant independent of 𝜀 and 𝑁.

Proof. We start by noting that 𝑈_𝜀−𝑢_𝜀≤𝑉_𝜀−𝑣_𝜀+𝑊_𝜀−𝑤_𝜀. (103) We have already shown in the previous lemmas that 𝑉_𝜀−𝑣_𝜀≤𝐶ℎ_max,𝑊+𝛿𝑡_𝜀−𝑤_𝜀≤𝐶𝑁⁻¹(ln𝑁)²+𝛿𝑡. (104) 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 𝜏=0.05, and no time step control mechanism or smoothing will be used.

Consider the following problem:𝛿−𝜀²𝑢_𝜀(𝑥,𝑡)𝛿𝑥²+𝛿((𝑢𝛿𝑥1+𝑥1−𝑥))_𝜀()+𝑥,𝑡𝛿𝑢_𝜀(𝑥,𝑡)()𝑢𝛿𝑡=𝑓,𝑥,𝑡∈Ω,_𝜀(0,𝑡)=𝑢_𝜀(1,𝑡)=0,𝑢_𝜀()1𝑥,0=1−4𝑥−2², (105a)(105b)where 𝑓 is chosen such that 𝑢_𝜀(𝑥,𝑡)=1−𝑒^{−(1−𝑥)/𝜀}1−𝑒^−1/𝜀−cos𝜋𝑥2 (106) is the exact solution when 𝜕𝑢_𝜀(𝑥,𝑡)/𝜕𝑡∼𝑂(𝜀) at 𝑡=𝒯=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: 𝑀_𝑖,𝑗=𝑢1+^𝑖+1,𝑗−𝑢_{𝑖−1,𝑗}𝑥_𝑖+1,𝑗−𝑥_{𝑖−1,𝑗}². (107) The tolerance is set at 10⁻⁶. 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 are obtained using Table 2 and (108) 𝑅^𝑁=log₂𝑢^𝑁_𝜀−𝑢^exact_𝜀𝑢_𝜀^2𝑁−𝑢^exact_𝜀. (108) Table 2 also gives the uniform convergence rates in the last row. These values are calculated using Table 1 and (109) 𝑅^𝑁_max=log₂max_𝜀𝑢^𝑁_𝜀−𝑢^exact_𝜀max_𝜀𝑢_𝜀^2𝑁−𝑢^exact_𝜀. (109) Table 2 shows that the method is almost of first order for a sufficiently large values of 𝑁.

Table 1 

Maximum pointwise errors.

Table 2 

Convergence rates.

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.

Figure 1 

Numerical solutions for 𝜀=2−4 and 𝑁=64.

Figure 2 

Numerical solutions for 𝜀=2−14 and 𝑁=64.

Figure 3 

Mesh trajectory with 𝜀=2−6, 𝜏=0.05, and 𝑁=32.

Figure 4 

Mesh trajectory with 𝜀=2−6, 𝜏=0.5, and 𝑁=32.

Figure 5 

Mesh trajectory for 𝜀=2−20, 𝜏=0.05, and 𝑁=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.