# A Higher-Order Finite Difference Scheme for Singularly Perturbed Parabolic Problem

Authors: Shifang Tian; Xiaowei Liu; Ran An
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, we deal with a singularly perturbed parabolic convection-diffusion problem. Shishkin mesh and a hybrid third-order finite difference scheme are adopted for the spatial discretization. Uniform mesh and the backward Euler scheme are used for the temporal discretization. Furthermore, a preconditioning approach is also used to ensure uniform convergence. Numerical experiments show that the method is first-order accuracy in time and almost third-order accuracy in space.


---

## Body

## 1. Introduction

We consider the singularly perturbed parabolic problem posed on the domainG=Ωx×Ωt=0,1×0,T as follows:(1)Lu+∂u∂t≔−ε∂2u∂x2−ax,t∂u∂x+cx,tu+∂u∂t=fx,t,x,t∈G,u0,t=u1,t=0,x∈Ω¯t,ux,0=u0x,x∈Ω¯x,where Ω¯t=0,T, Ω¯x=0,1, and ε is a small positive perturbation parameter satisfying(2)0<ε≤ε∗≪1,with a positive constant ε∗, and u0x is the initial value when t=0. We assume that the functions ax,t, cx,t, and fx,t are sufficiently smooth, and that ax,t and cx,t satisfy(3)ax,t≥β>0,cx,t>0,∀x,t∈G¯,where β is a positive constant.With these conditions, there exists a uniqueC4,2G¯-solution u of problem (1) (see, for instance, [1]). The C4,2G¯-solutions u of problem (1) satisfied from its original function u to its fourth partial derivative for spatial variable uxxxx are continuous, and solutions from its original function u to its second partial derivative for temporal variable utt are also continuous.Singular perturbation problems play an important role in many areas, such as astronomy, mechanics, and fluid dynamics. It also has a broad background and important applications in control systems with different time scales [2–4]. It is especially important to find a uniform and effective approximate solution when the exact solution cannot be obtained. There are many methods for solving singular perturbation problems. Recent convergence analysis of the finite element method is referred to [5–14]. Except the finite element method, the finite difference method is the most widely used one at present. Nowadays, more and more people begin to study higher-order finite difference schemes for solving singular perturbation problems. In 1988, Vulanović in [15] proposed a third-order hybrid finite difference scheme and showed numerical results on the Shishkin mesh. Afterwards, Vulanovic and Nhan in [16] improved on what had already been done and proposed a new uniformly convergent numerical scheme. Both the methods proposed in [15, 17] have been analyzed on a piecewise-uniform Shishkin mesh and were proved to be almost third-order accuracy in space. Comparing the previous scheme, the difference is that when ε is large enough, the accuracy of the new scheme is better than the that of the previous scheme. However, when ε is small enough, there is no difference between these two methods.In this paper, our primary aim is to propose and analyze a higher-order hybrid finite difference scheme for problem (1). This is accomplished by discretizing the domain Ωx using the Shishkin mesh and by considering the uniform mesh in the temporal direction. In order to obtain the fully discrete scheme, we adopt the two-stage discretization process. The first stage consists of discretizing the time derivative with the backward difference scheme on the uniform mesh. In the second stage, discretize in the spatial direction by utilizing a hybrid finite difference scheme on the Shishkin mesh.The ultimate goal of numerical methods for problem (1) is to obtain a series of discrete solutions so as to achieve a numerical approximation of the continuous solution. Such that its error converges to 0 uniformly as N⟶+∞, where N is the number of discretization on the spatial mesh. Apart from this, in numerical experiments, we also need to illustrate that the proposed scheme is almost third-order accurate in space.The rest of the paper is as follows. In Section2, we define the meshes for temporal and spatial discretization and introduce some special difference operators. In Section 3, we define some difference operators and the final finite difference scheme. In Section 4, we give the linear equations needed to solve the problem and get the coefficient matrix and the right end term. In Section 5, we preprocess the coefficient matrix. In Section 6, we give the pseudo code to solve the problem. In Section 7, we give the results of numerical experiments. In Section 8, some final conclusions are given.


## 2. The Mesh

Here, in this section, we describe the uniform mesh for the temporal discretization of the domainΩt and the Shishkin mesh for the spatial discretization of the domain Ωx.We will often use the assumption that(4)ε≤C∗N−1,where C∗ is a sufficiently small positive constant which is independent of both ε and N. All constants, independent of ε and N, are denoted generically by C.


### 2.1. The Uniform Mesh

For the time domain0,T, we use a uniform mesh with time step Δt, such that(5)ΩtM=tn=nΔtn=0,…,MΔt=TM,where M is the number of mesh points in the t-direction on the interval 0,T.


### 2.2. Shishkin Mesh

Since problem (1) has a boundary layer along the side x=1, the mesh should be condensing in the neighborhood of x=1. Define(6)σ=min12,αεβlnN,with α≥4 (the proof of range of α was given in [16]). To define the piecewise-uniform mesh, we divide the domain 0,1 into two subdomains, such that 0,1=0,1−σ∪1−σ,1 and then divide each of the subdomains into N/2 equal intervals. Set(7)h≔hi=1−σN,fori=0,1,2,…,N/2,H≔hi=σN,fori=N/2+1,…,N.Now, we denote the spatial grids by(8)ΩxN=0=x0,x1,…,xN/2=1−σ,…,xN=1,where(9)xi=ih,fori=0,1,2,…,N/2,1−σ+i−N2H,fori=N/2+1,…,N,and N≥4 be a positive even integer. Here, the transition point 1−σ separates the coarse and fine portions of the mesh.Moreover, define(10)xi+z=xi+zhi+1,ifz∈0,1,xi+zhi,ifz∈−1,0,where z is some fixed constant.


## 2.1. The Uniform Mesh

For the time domain0,T, we use a uniform mesh with time step Δt, such that(5)ΩtM=tn=nΔtn=0,…,MΔt=TM,where M is the number of mesh points in the t-direction on the interval 0,T.


## 2.2. Shishkin Mesh

Since problem (1) has a boundary layer along the side x=1, the mesh should be condensing in the neighborhood of x=1. Define(6)σ=min12,αεβlnN,with α≥4 (the proof of range of α was given in [16]). To define the piecewise-uniform mesh, we divide the domain 0,1 into two subdomains, such that 0,1=0,1−σ∪1−σ,1 and then divide each of the subdomains into N/2 equal intervals. Set(7)h≔hi=1−σN,fori=0,1,2,…,N/2,H≔hi=σN,fori=N/2+1,…,N.Now, we denote the spatial grids by(8)ΩxN=0=x0,x1,…,xN/2=1−σ,…,xN=1,where(9)xi=ih,fori=0,1,2,…,N/2,1−σ+i−N2H,fori=N/2+1,…,N,and N≥4 be a positive even integer. Here, the transition point 1−σ separates the coarse and fine portions of the mesh.Moreover, define(10)xi+z=xi+zhi+1,ifz∈0,1,xi+zhi,ifz∈−1,0,where z is some fixed constant.


## 3. Discretization

In this section, we will give different difference operators corresponding to different points on the Shishkin mesh and combine these difference operators to form the final numerical scheme.Firstly, we denote byUij the approximation of u at point xi,tj and set fi+zj=fxi+z,tj. Then, we use Dχ,z0Uij, Dχ,z1Uij, and Dχ,z2Uij as the approximations of uxi+z,tj, uxxi+z,tj, and uxxxi+z,tj, respectively. They are defined by the following equation [16]:(11)Dχ,z0Uij=12zz−1Ui−1j+21−z2Uij+zz+1Ui+1j,Dχ,z1Uij=16χ−3z2+6z−2Ui−1j+33z2−4z−1Uij+3−3z2+2z+2Ui+1j+3z2−1Ui+2j,Dχ,z2Uij=1χ21−zUi−1j+3z−2Uij+1−3zUi+1j+zUi+2j,where χ is the step size of a uniform mesh and z is a constant satisfying z∈−1,1. We set ei+zn as the truncation error between the numerical solution and the exact solution.(12)ei+zn=Dχ,znUij−uxi+z,tj,n=0,1,2.Firstly,ei+z0=Dχ,z0Uij−uxi+z,tj.Lemma 1.
Suppose that(13)∂m∂xmux,t≤C,x,t∈Ω¯x×Ω¯t,m≥0.The truncation error associated toei+z0 satisfies(14)ei+z0∞≤Cχ3,where χ is the step size.Proof.
We substituteuxi−1,tj,uxi,tj,anduxi+1,tj for Ui−1j,Uij,andUi+1j in operator Dχ,z0Uij and apply the Taylor expansion to obtain(15)ei+z0=12zz−1uxi,tj−χu′xi,tj+12χ2u″xi,tj−16χ3u3xi,tj+⋯+21−z2uxi,tj+zz+1uxi,tj−χu′xi,tj+12χ2u″xi,tj−16χ3u3xi,tj+⋯−uxi,tj−χu′xi,tj+12χ2u″xi,tj−16χ3u3xi,tj+⋯=16zχ3u3xi,tj−16z3χ3u3xi,tj+Oχ4.
Thus,(16)ei+z0∞≤Cχ3.
Secondly,ei+z1=Dχ,z1Uij−uxi+z,tj.Lemma 2.
Suppose that(17)∂m∂xmux,t≤C,x,t∈Ω¯x×Ω¯t,m≥0.The truncation error associated toei+z1 satisfies(18)ei+z1∞≤Cχ3,where χ is the step size.Proof.
Similar to above, we substituteuxi−1,tj,uxi,tj,uxi+1,tj,anduxi+2,tj for Ui−1j,Uij,Ui+1j,andUi+2j in operator Dχ,z1Uij and again apply the Taylor expansion to obtain(19)ei+z1=16χ−3z2+6z−2uxi,tj−χu′xi,tj+12χ2u′′xi,tj−16χ3u3xi,tj+⋯+33z2−4z−1uxi,tj+3−3z2+2z+2uxx,tj−χu′xi,tj+12χ2u″xi,tj−16χ3u3xi,tj+⋯+3z2−1uxx,tj−χu′xi,tj+122χ2u″xi,tj−162χ3u3xi,tj+⋯−uxi,tj−χzu′xi,tj+12χz2u′′xi,tj−16χz3u3xi,tj+⋯=16zχ3u3xi,tj−16z3χ3u4xi,tj+Oχ4.
Thus,(20)ei+z1∞≤Cχ3.
In conclusion, bothDχ,z0Uij and Dχ,z1Uij are third-order accurate with respect to the spatial variable x for any value of z; if z=1/3, Dχ,z1Uij is transformed into the classical three-point scheme. And in the same way, operator D″Uij (27), D¯0UN/2j (29), D^″UN/2j (30), and time difference operator Dt−Uij can all be proven.
Moreover,ei+z2=Dχ,z2Uij−uxi+z,tj.Lemma 3.
Assume that(21)∂m∂xmux,t≤C,x,t∈Ω¯x×Ω¯t,m≥0.The truncation error associated toei+z2 satisfies that if z=3−15/6,(22)ei+z2∞≤Cχ3,else(23)ei+z2∞≤Cχ2,where χ is the step size.Proof.
Once more, substitutinguxi−1,tj,uxi,tj,uxi+1,tj,anduxi+2,tj for Ui−1j,Uij,Ui+1j,andUi+2j in operator Dχ,z2Uij and applying the Taylor expansion results in(24)ei+z2=1χ21−zuxi,tj−χu′xi,tj+12χ2u″xi,tj−16χ3u3xi,tj+124χ4u4xi,tj−1120χ5u5xi,tj+⋯+3z−2uxi,tj+1−3zuxi,tj−χu′xi,tj+12χ2u″xi,tj−16χ3u3xi,tj+124χ4u4xi,tj−1120χ5u5xi,tj+⋯+zuxi,tj−χu′xi,tj+122χ2u″xi,tj−162χ3u3xi,tj+1242χ4u4xi,tj−11202χ5u5xi,tj+⋯−u″xi,tj+zχu3xi,tj+12zχ2u4xi,tj+16zχ3u5xi,tj+⋯=34−112z−z22χ2u4xi,tj+14+130z+z36χ3u5xi,tj+Oχ4,we can found that the operator Dχ,z2Uij in general is second-order accurate, and if z=3−15/6, it is third-order accurate.
These schemes can be used to create the following difference operatorΛχ,z:(25)Λχ,zUij=−εDχ,z2Uij−axi+z,tjDχ,z1Uij+cxi+z,tjDχ,z0Uij,i=1,…,N2−2,where z=3−15/6 and χ=h. The operator Λχ,z is only used as part of the discretization on the Shishkin grid because the Shishkin grid is not uniform in the entire computational domain. More specifically, the difference operator Λχ,z cannot be applied at xN/2−1 and xN−1.Remark 1.
Because schemeΛχ,z has point xi+2,tj in it, the Shishkin grid used is divided into two intervals, 0,xN/2 and xN/2+1,xN. Since the step size is different and the point xi+2 spans two intervals, the difference operator Λχ,z cannot be applied at either xN/2−1 or xN−1.
Now, we introduce the difference operator as follows:(26)Λ˜χ,zUij=−εD″Uij+axi+z,tjDχ,z1Uij+cxi+z,tjDχ,z0Uij,i=N2+1,…,N−1,with(27)D″Uij=2hij+hi+1jUi+1j−Uijhi+1j−Uij−Ui−1jhij,where z=1/3 and χ=H, and it is also third-order accurate.
Then, we give the scheme at pointxN/2 by means of one-side difference schemes as follows:(28)Λ^H,1−zUN/2j=−εD^″UN/2j+axN/2+1−z,tjDH,1−z1UN/2j+cxN/2+1−z,tjD¯0UN/2j,with(29)D¯0UN/2j=16+z2UN/2j+23UN/2+1j+16−z2UN/2+2j,and(30)D^″UN/2j=1H2UN/2j−2UN/2+1j+UN/2+2j,where z=1−1/3 and χ=H, and both D¯0UN/2j and D^″UN/2j have third-order accuracy in space.
In addition, about thet-direction, the discretization of utxi,tj by the backward Euler scheme is defined by(31)Dt−Uij=Uij−Uij−1Δt≈utxi,tj,with the time step Δt, and it is first-order accurate with respect to the temporal variable t.
Finally, we combineDt−Uij with three difference operators Λχ,zUij, Λ˜χ,zUij, and Λ^χ,zUij at different points, respectively, and finally propose the following numerical scheme:(32)Λh,zUij+Dt−Uij=fi+zjz=3−156,for1≤i≤N2−2,Λ˜h,zUij+Dt−Uij=fi+zjz=13,fori=N2−1,Λ^H,zUij+Dt−Uij=fi+zjz=1−13,fori=N2,Λ˜H,zUij+Dt−Uij=fi+zjz=13,forN2+1≤i≤N−1,with j=1,…,M.


## 4. Linear Problem

The corresponding difference schemes ofuxx, ux, and u at point xi+z,tj and ut at point xi,tj are substituted into equation (1). When combined with scheme (32), the following linear equations (33) are obtained:(33)r1Ui−1j+r2Uij+r3Ui+1j+r4Ui+2j=gj,where Uij is the approximation of uxi,tj and r1,r2,r3,r4, and gj are defined as follows: if 1≤i≤N/2−2,(34)r1=−ε1−zΔth2−Δtbxi+z,tj−3z2+6z−26h+12Δtzz−1+12zz−1,r2=−εΔt3z−2h2−Δtbxi+z,tj3q2−4z−12h+Δt1−z2+zz+1,r3=−εΔt1−3zh2−Δtbxi+z,tj−3z2+2z+22h+12Δtzz+1+12zz+1,r4=−εΔtzh2−Δtbxi+z,tj3z2−16h,gj=Δtfxi+z,tj+Dh,z0Uij−1,where Δt is the time steps, h is the space steps defined by (7), z=3−15/6, and Dh,z0Uij−1 is defined by (11). If i=N/2−1,(35)r1=−εΔth2−Δtbxi+z,tj−3z2+6z−26h+12Δtzz−1+12zz−1,r2=εΔth2+εΔth2−Δtbxi+z,tj3z2−4z−12h+Δt1−z2+1−z2,r3=−εΔth2−Δtbxi+z,tj−3z2+2z+22h+12Δtzz+1+12zz−1,r4=0,gj=Δtfxi+z,tj+Dh,z0Uij−1,where Δt is the time steps, h is the space steps defined by (7), z=1/3, and Dh,z0Uij−1 is defined by (11). If i=N/2,(36)r1=0,r2=−εΔtH2−Δtbxi+z,tj3z2−2z−22H+12Δt13+z+1213+z,r3=−εΔtzH2−Δtbxi+z,tj3z2−16H+23Δt+23,r4=−εΔtH2−Δtbxi+z,tj3z2−6z+26H+12Δt13−z+1213−z,gj=Δtfxi+1−z,tj+D¯0UN/2j−1,where Δt is the time steps, H is the space steps defined by (7), z=1−1/3, and D¯0UN/2j−1 is defined by (11). If N/2+1≤i≤N−1,(37)r1=−εΔtH2−Δtbxi+z,tj−3z2+6z−26H+12Δtzz−1+12zz−1,r2=εΔtH2+εΔtH2−Δtbxi+z,tj3z2−4z−12H+Δt1−z2+1−z2,r3=−εΔtH2−Δtbxi+z,tj−3z2+2z+22H+12Δtzz+1+12zz+1,r4=0,gj=Δtfxi+z,tj+DH,z0Uij−1,where Δt is the time steps, H is the space steps defined by (7), z=1/3, and DH,z0Uij−1 is defined by (11).Finally, the linear system for numerical scheme (32) is obtained, that is,(38)A^x=g.Here, the coefficient matrixA^ is defined by(39)row123⋮N2−2N2−1N2N2+1N2+2⋮N−2N−1r2r3r4r1r2r3r4r1r2r3r4⋱⋱⋱⋱r1r2r3r4r1r2r3r2r3r4r1r2r3r1r2r3⋱⋱⋱r1r2r3r1r2,where each unwritten element is 0. The unknown term x is defined by x0,x1,x2,…,xN−1,xNT with x0=0 and xN=1, and the right end term g is defined by g0,g1,g2,…,gN−1,gNT.


## 5. Preconditioning

In this section, we analyze theN−1×N−1 matrix A^, which corresponds to the first scheme (32) and is acquired in Section 4. We need to assume (4) and that N is sufficiently large and(40)N≥N∗,where N∗ is a positive integer independent of ε. There exist constants N∗ and C∗ such that (4) and (40) are satisfied, and the matrix A^ has the following structure:(41)row123⋮N2−2N2−1N2N2+1N2+2⋮N−2N−1+−+−+−+−+−+⋱⋱⋱⋱−+−+−+−+−+−+−−+−⋱⋱⋱−+−−+,where each unwritten element is 0.Multiply its equations1,2,…,N/2−1 by h/H. We do this to achieve consistency uniform in ε, but at the same time, the coefficient matrix gets preconditioned as well (the preconditioning be described in [18, 19]).Thus, we take the matrix of the preconditioned system as follows:(42)B^≔diagm1,m2,…,mN−1A^,where(43)mi=hH,fori=0,1,2,…,N2−1,1,fori=N2,…,N−1.


## 6. Pseudo Code

In this section, the pseudo code needed to solve problem (1) using numerical scheme (32) in MATLAB will be presented. In general, if mathematical tools are used to solve problem like this, by scheme (32), there are six steps as follows [17]:(1)
Set the uniform mesh for temporal variable (M is the total number of points in t-direction)(2)
Set the Shishkin mesh for spatial variable (N is the total number of points in space)(3)
Write down the coefficient matrixA^ and the right end term g for the linear system as follows:(44)A^=forn=1:Mtemporal points,fori=N2−2,hi,n=xi,n−xi−1,n,r1Ui−1n+r2Uin+r3Ui+1n+r4Ui+2n=gn,end,fori=N2−1,hi,n=xi,n−xi−1,n,r1Ui−1n+r2Uin+r3Ui+1n+r4Ui+2n=gn,end,fori=N2,hi,n=xi,n−xi−1,n,r1Ui−1n+r2Uin+r3Ui+1n+r4Ui+2n=gn,end,fori=N2+1:N−1,hi,n=xi,n−xi−1,n,r1Ui−1n+r2Uin+r3Ui+1n+r4Ui+2n=gn,end,end,whereUin is the approximation of uxi,tn, hi,n is the space steps, r1,r2,r3,andr4 are the elements of the coefficient matrix, and gn is the right end term (it is defined in Section 4)(4)
A new matrixB^ (from Section 5) is obtained by preprocessing the matrix A^(5)
The new matrixB^ and the right end term g are used to solve problem (1)(6)
The maximum pointwise errors and the orders of convergence are calculated


## 7. Numerical Experiments

In this section, we shall present the numerical results obtained by the proposed numerical schemes (32) for the test problem (45) on the piecewise-uniform rectangular mesh G=ΩxN×ΩtM. In both cases, we perform the numerical experiments by choosing the constants α=4 and β=1 in (6) and the time step Δt=1.0/M.For numerical tests, we consider the following singularly perturbed parabolic problem:(45)−ε∂2u∂x2−x+1∂u∂x+u+∂u∂t=fx,t,x,t∈G,u0,t=u1,t=0,t∈0,1,ux,0=0,x∈0,1,where G≔0,1×0,1. We choose the initial data ux,0=0 and the exact solution for problem (45) as follows [16]:(46)ux,t=te−x/ε−ex+e−e−1/εx.As the exact solution of problem (45) is known, we calculate the maximum pointwise error by(47)EεN,Δt=maxui,j−Ui,jN,Δt,for each ε, where ui,j and Ui,jN,Δt denote the exact solution and numerical solution on xi,tj, respectively. The convergence order is calculated by the following formula:(48)RεN,Δt=log2EεN,ΔtEε2N,Δt/2.The maximum pointwise errorsEεN,Δt and the order of convergence RεN,Δt by using schemes (32) are presented in Table 1. In table, we can observe the ε-uniform convergence of the numerical scheme. The order of convergence in Table 1 is first-order due to the effect of time error. In order to justify the spatial order of convergence precisely, we take M=N3 and the order of convergence is defined by(49)RεN,Δt=log2EεN,ΔtEε2N,Δt/8.Table 1
Temporal errors and converge orders of scheme (32), Eε, Rε.
N,Δtε32,1/3264,1/64128,1/128256,1/256512,1/51210−20.013670.0062940.0028640.0012935.788e − 04Eε10−21.121.141.151.15—Rε10−40.020670.0105050.0052850.002650.00132510−40.97100.99100.99660.9989—10−60.020750.010550.0053150.002670.00133510−60.97540.98980.99540.9977—10−80.020750.010550.0053150.0026650.00133510−80.97540.98980.99540.9983—The numerical results are presented in Table2, where the spatial convergence order is almost third-order.Table 2
Spatial errors and converge orders of scheme (32), Eε, Rε.
N,Δtε32,1/32364,1/643128,1/1283256,1/2563512,1/512310−20.00880.00214.538e − 048.644e − 051.453e − 05Eε10−22.052.222.392.57—Rε10−40.00880.00214.538e − 048.644e − 051.453e − 0510−42.052.222.392.57—10−60.00890.00224.715e − 049.397e − 051.739e − 0510−62.032.222.312.44—10−80.00890.00224.714e − 049.397e − 051.7393e − 0510−82.042.202.332.43—


## 8. Conclusion

A hybrid scheme is proposed for obtaining a numerical solution to the singularly perturbed parabolic problem. The idea is based on the methods presented in the existing research study [15–19]. It can be seen from the results of numerical experiments, whether in space or in time, the scheme is robust insomuch the error of the numerical solution does not increase when ε⟶0. On the contrary, the proposed schemes improve as ε diminishes, becoming almost third-order accurate with the spatial variable and first-order accurate with the temporal variable. The numerical results were compared with those from literature [15–18, 20, 21] which showed that all results reach the expected order of convergence. However, so far it is not possible to construct an arbitrary high-order difference scheme for Shishkin grids, meaning further research is needed. It should be noted that parallelization is not discussed in this article, so the reader is encouraged to refer to additional work [20, 21].



---

*Source: 9941692-2021-08-20.xml*