# A Fitted Mesh Cubic Spline in Tension Method for Singularly Perturbed Problems with Two Parameters

Authors: Tariku Birabasa Mekonnen; Gemechis File Duressa
Journal: International Journal of Mathematics and Mathematical Sciences
 (2022)

---

## Abstract

A numerical treatment via a difference scheme constructed by the Crank–Nicolson scheme for the time derivative and cubic spline in tension for the spatial derivatives on a layer resolving nonuniform Bakhvalov-type mesh for a singularly perturbed unsteady-state initial-boundary-value problem with two small parameters is presented. Error analysis of the constructed scheme is discussed and shown to be parameter-uniformly convergent with second-order convergence. Numerical experimentation is taken to confirm the theoretical findings.


---

## Body

## 1. Introduction

The following class of unsteady-state singularly perturbed one-dimensional parabolic convection-diffusion problems(1)Lε,μu≔∂u∂t−ε∂2u∂x2−μax,t∂u∂x+bx,tu=fx,t,on the domain D=Ω×0,T,Ω=0,1, together with prescribed initial and boundary conditions of Dirichlet-type,(2)ux,0=sx,x∈Ω¯,u0,t=0=u1,t,t∈0,T,where 0<ε≪1 and 0≤μ≪1 are two small positive parameters considered. The functions ax,t,bx,t,fx,t, and sx are assumed to be sufficiently smooth satisfying ax,t≥α>0,bx,t≥β>0,∀x,t∈D¯. Furthermore, we assume the compatibility condition at the two corner points, i.e., u0,0=s0,u1,0=s1. Within this frame, a unique solution ux,t exists for the problem and exhibits a boundary layer on both the left and right lateral boundaries of the spatial domain Ω¯ with significantly different width depending on the relation between μ2 and ε [1–6].Singularly perturbed problems (SPPs) widely appear in modeling physical problems such as fluid dynamics, heat and mass transfer in chemical engineering, quantum mechanics [7], elasticity, the theory of plates and shells, oil and gas reservoir simulation, and magnetic-hydrodynamic flow in general [2] and two-parameter singularly perturbed problems (TPSPPs) occur in chemical flow reactor theory [5] as well as in the case of boundary layers controlled by suction (or blowing) of some fluid [8] in particular.Several articles dealt with the solution of an SPP and showed that it is revealed to have a sudden change in the boundary layer region and behaves normally in the outer region. Due to this multiscale character, the classical numerical methods cease to give satisfactory numerical results and are not competitive computationally. Therefore, to overcome the shortcomings, researchers are forced to construct numerical methods that are flexible enough to provide a solution whose accuracy is independent of the perturbation parameter(s), and its singular nature can easily be captured. The sounding numerical methods on hand are the specially designed (fitted) ones. Just to list some of these, the authors of [9–15] used the fitted operator method, the authors of [16–18] implemented the fitted mesh method, Gowrisankar and Natesan [19] constructed a layer-adapted scheme through the equidistribution of a positive integrable monitor function, and Khan et al. [20] applied variable mesh in which more mesh points are in the layer region. On the other part, the authors of [21, 22] used unfitted methods and the results yielded are not ε−uniform.In this paper, we construct a numerical scheme for the parabolic convection-diffusion problem (1-2) through the fitted mesh method in which Crank–Nicolson discretization is used for the time variable and cubic spline in tension on layer-adapted meshes of Bakhvalov-type for the spatial variable discretization. The scheme constructed from this composition is in its first use for the problem under consideration.


## 2. Properties of the Continuous Problem

In this section, we consider some properties of the continuous problem through the maximum principle and bounds on the solution and on its derivatives that enable us to see a priori estimates for the solutionux,t and its derivatives.Lemma 1.
LetΠx,t∈C2,1D¯. If Πx,t≥0, for all x,t on the boundary ∂D of the domain D, and Lε,μΠ≥0, for all x,t∈D, then Πx,t≥0, for all x,t∈D¯.Proof.
Assume to the contrary that there exists a pointx∗,t∗∈D¯ such that(3)Πx∗,t∗=minx,t∈D¯Πx,t<0.
It is obvious thatx∗,t∗∉∂D and implies x∗,t∗∈D. Applying the differential operator Lε,μ in equation (1) on Πx,t at the critical point x∗,t∗ yields(4)Lε,μΠx∗,t∗≔∂Π∂tx∗,t∗−ε∂2Π∂x2x∗,t∗−μax∗,t∗∂Π∂xx∗,t∗+bx∗,t∗Πx∗,t∗.
By the partial derivative test, we have∂2Π/∂x2x∗,t∗≥0,∂Π/∂tx∗,t∗=0, and ∂Π/∂xx∗,t∗=0. Incorporating our assumptions, bx∗,t∗≥β>0 with Πx∗,t∗<0, we can arrive at Lε,μΠx∗,t∗≤0. This contradicts the given statement. Therefore, we can conclude that the minimum of Πx,t is nonnegative.Lemma 2.
A solutionux,t of equation (1) satisfies the estimate(5)u≤β−1f+maxsx,where . is the maximum norm.Proof.
For the proof, one can refer [23].Lemma 3.
The derivatives of the solutionux,t satisfy the following bound for all nonnegative integers i,j such that i+3j≤4.(6)∂i+ju∂xi∂tjD¯≤ε−i/2ifμ2≤Cε,μεiifμ2≥Cε,where the constant C is independent of ε and μ and depends only on the bounded derivatives of the coefficients and the source term.Proof.
The proof is given in [2].


## 3. Discretization and Mesh Generation

Consider a uniform mesh on the time domain0,T with M mesh intervals spaced by k=T/M and given by DtM=tm=m−1k,m=1,2,…,M+1. Applying the Crank–Nicolson method for the time variable in equations (1) and (2), we obtain the following two-level semidiscrete scheme:(7)ux,0=sx,Lε,μMu˜≡−εu˜xxx−μaxu˜xx+qxu˜x=gx,u˜0=0=u˜1,where qx=bx,tm+2/k,u˜x=ux,tm,ax=ax,tm,px=bx,tm−2/k,Hx=fx,tm+fx,tm−1 and gx=εuxxx,tm−1+μax,tm−1uxx,tm−1−pxux,tm−1+Hx.A characteristic equation corresponding to equation (7) that is helpful to describe the boundary layers occurring at the two end points of Ω¯ is(8)−εr2x−μaxrx+qx=0.This defines two continuous functions(9)r1x=−μax−μax2+4εqx2ε,r2x=−μax+μax2+4εqx2ε.Let(10)η1=−maxx∈Ω¯r1xandη2=minx∈Ω¯r2x.Lemma 4.
For a fixed number0<ρ<1, a certain order δ, and the solution u˜x of equation (3), the following holds:(11)∂iu˜∂xi≤C1+η1ie−ρη1ix+η2ie−ρη2i1−x,for0≤i≤δ,where C is a constant independent of ε,μ.Proof.
For the proof, readers can refer [24].Lemma 5.
The error estimate in the temporal direction satisfies(12)TEt≤Ck2.Proof.
The proof is given in [11].


### 3.1. Graded Mesh

LetN+1 be the number of graded mesh points, 0=x1<x2<…<xN+1=1, of the interval Ω¯ that are arbitrarily spaced by hn=xn+1−xn,n=1,…,N yielding a space mesh DxN=x1,x2,…,xN+1. Now, divide the domain Ω¯ into three subintervals, namely, 0,τ1,τ1,1−τ2, and 1−τ2,1, such that the second subinterval is with equidistant N/2 mesh subintervals and the rest two are with gradually graded N/4 mesh subintervals. The two transition points τ1 and τ2 are chosen such that exp−η1x≤η1−2 in τ1,1 and exp−η21−x≤η2−2 in 0,1−τ2 and given by(13)τ1=min14,2η1logη1andτ2=min14,2η2logη2.Since we are using nonuniform meshes in the spatial direction, we do not assumeτ1 and τ2 to be 1/4. Furthermore, the relation ηj−1logηj≤CN−1,j=1,2 is given in [24].For the two boundary layers, following the technique in [3], we have two generating functions(14)ψ1ξ=−log1−41−1/η1ξandξψ2=−log1−41−1/η21−ξ,which satisfy ψ10=0,ψ11/4=logη1,ψ23/4=logη2,andψ21=0. Analogous to [24], the mesh points in the spatial direction can be defined by(15)xn=−2η1log1−4N1−1η1n−1,n=1,2,…,N4+1,τ1+21−τ1−τ2n−1/N−1/4,n=N4+2,…,3N4+1,1+2η2log1−41−1η21−n−1N,n=3N4+2,…,N+1.This gives condensed and refined meshes in the layer regions and points where the mesh switches from fine to coarse and vice versa are chosen to bexN/4+1=τ1 and x3N/4+1=1−τ2. The step sizes are related by(16)h1<h2<…hN/4,hN/4+1=hN/4+2=…=h3N/4,h3N/4+1>h3N/4+2>…>hN.Lemma 6.
The mesh sizeshn=xn+1−xn,n=1,2,…,N satisfy(17)hn≤Cη1−1,n=1,2,…,N/4,CN−1,n=N/4+1,…,3N/4,Cη2−1,n=3N/4+1,…,N.Proof.
Forn=1,2,…,N/4,(18)hn=xn+1−xn=−2η1log1−4N1−1η1n−−2η1log1−4N1−1η1n−1,=−2η1ψ1ξn+2η1ψ1ξn−1=−2η1ψ1ξn−ψ1ξn−1.
However, from the mean value theorem, there existsςn∈ξn−1,ξn such that(19)ψ1′ςn=ψ1ξn−ψ1ξn−1/ξn−ξn−1.
This giveshn≤Cη1−1,n=1,2,…,N/4.
Forn=N/4+1,…,3N/4, from equation (15), we have hn=xn+1−xn=2/N1−τ1−τ2≤CN−1. The proof in the third subinterval can be shown in a similar way as in the first subinterval.


### 3.2. Total Difference Scheme

LetQx,tm interpolate the solution Uxn,tm=Unm and pass through the points xn,tm,Unm and xn+1,tm,Un+1m such that(20)Qx,tm=−hn2θ2sinθsinθx−xnhnYn+1m+sinθxn+1−xhnYnm+hn2θ2x−xnhnYn+1m+hn2θ2Un+1m+xn+1−xhnYnm+θ2hn2Unm,where Ynm=U'′nm and θ is a free parameter used to maintain consistency. The first derivatives of Qx,tm from the right and left of x=xn are given, respectively, by(21)Q′xn+,tm=Un+1m−Unmhn+hnθ21−θcscθYn+1m−1−θcotθYnm,Q′xn−,tm=Unm−Un−1mhn−1+hn−1θ21−θcotθYnm−1−θcscθYn−1m.Then, from the continuity condition, we haveQ′xn+,tm=Q′xn−,tm and obtain(22)λ1hn−1Yn−1m+λ2hn−1+hnYnm+λ1hnYn+1m=Un+1m−Unmhn−Unm−Un−1mhn−1,where λ1=θcscθ−1/θ2,λ2=1−θcotθ/θ2. Application of the L’Hopital’s rule gives(23)λ1+λ2=1/2.Atx=xn−1,xn,xn+1 using the notation Wnm=U'′nm+U'′nm−1, from equation (7), we can have(24)Wn−1m=1ε−μan−1mU′n−1m+an−1m−1U′n−1m−1+qn−1mUn−1m+pn−1m−1Un−1m−1−Hn−1m,Wnm=1ε−μanmU′nm+anm−1U′nm−1+qnmUnm+pnm−1Unm−1−Hnm,Wn+1m=1ε−μan+1mU′n+1m+an+1m−1U′n+1m−1+qn+1mUn+1m+pn+1m−1Un+1m−1−Hn+1m,where, from [25],(25)U′n−1m=1hn−1hnhn−1+hn−hn2+2hn−1hnUn−1m+hn−1+hn2Unm−hn−12Un+1m,U′nm=1hn−1hnhn−1+hn−hn2Un−1m+hn2−hn−12Unm+hn−12Un+1m,U′n+1m=1hn−1hnhn−1+hnhn2Un−1m−hn−1+hn2Unm+hn−12+2hn−1hnUn+1m.Now, adopting the cubic spline in tension used in [25] yields(26)λ1hn−1Wn−1m+λ2hn−1+hnWnm+λ1hnWn+1m=Un+1m−Unmhn−Unm−Un−1mhn−1+Un+1m−1−Unm−1hn−Unm−1−Un−1m−1hn−1.Substituting equations (24) and (25) into equation (26) and writing it in compact form, we arrive at(27)Lε,μM,NUnm≔An−Un−1m+An0Unm+An+Un+1m+Bn−Un−1m−1+Bn0Unm−1+Bn+Un+1m−1=Fnm,where n=2,3,…,N,m=2,3,…,M,(28)An−=μλ1hn+2hn−1hn−1+hnan−1m+λ2hnhn−1anm−λ1hn2hn−1hn−1+hnan+1m+λ1hn−1qn−1m−εhn−1,An0=μ−λ1hn+hn−1hnan−1m−λ2hn2−hn−12hn−1hnanm+λ1hn−1+hnhn−1an+1m+λ2hn−1+hnqnm+εhn+hn−1hn−1hn,An+=μλ1hn−12hnhn−1+hnan−1m−λ2hn−1hnanm−λ1hn−1+2hnhn−1+hnan+1m+λ1hnqn+1m−εhn,Bn−=μλ1hn+2hn−1hn−1+hnan−1m−1+λ2hnhn−1anm−1−λ1hn2hn−1hn−1+hnan+1m−1+λ1hn−1pn−1m−1−εhn−1,Bn0=μ−λ1hn+hn−1hnan−1m−1−λ2hn2−hn−12hn−1hnanm−1+λ1hn−1+hnhn−1an+1m−1+λ2hn−1+hnpnm−1+εhn+hn−1hn−1hn,Bn+=μλ1hn−12hnhn−1+hnan−1m−1−λ2hn−1hnanm−1−λ1hn−1+2hnhn−1+hnan+1m−1+λ1hnpn+1m−1−εhn,Fnm=λ1hn−1Hn−1m+λ2hn+hn−1Hnm+λ1hnHn+1m.The system (27) together with the given conditions is uniquely solvable. Furthermore, the matrix associated with this system of equations satisfies the property of an M-matrix that is given in [26].


## 3.1. Graded Mesh

LetN+1 be the number of graded mesh points, 0=x1<x2<…<xN+1=1, of the interval Ω¯ that are arbitrarily spaced by hn=xn+1−xn,n=1,…,N yielding a space mesh DxN=x1,x2,…,xN+1. Now, divide the domain Ω¯ into three subintervals, namely, 0,τ1,τ1,1−τ2, and 1−τ2,1, such that the second subinterval is with equidistant N/2 mesh subintervals and the rest two are with gradually graded N/4 mesh subintervals. The two transition points τ1 and τ2 are chosen such that exp−η1x≤η1−2 in τ1,1 and exp−η21−x≤η2−2 in 0,1−τ2 and given by(13)τ1=min14,2η1logη1andτ2=min14,2η2logη2.Since we are using nonuniform meshes in the spatial direction, we do not assumeτ1 and τ2 to be 1/4. Furthermore, the relation ηj−1logηj≤CN−1,j=1,2 is given in [24].For the two boundary layers, following the technique in [3], we have two generating functions(14)ψ1ξ=−log1−41−1/η1ξandξψ2=−log1−41−1/η21−ξ,which satisfy ψ10=0,ψ11/4=logη1,ψ23/4=logη2,andψ21=0. Analogous to [24], the mesh points in the spatial direction can be defined by(15)xn=−2η1log1−4N1−1η1n−1,n=1,2,…,N4+1,τ1+21−τ1−τ2n−1/N−1/4,n=N4+2,…,3N4+1,1+2η2log1−41−1η21−n−1N,n=3N4+2,…,N+1.This gives condensed and refined meshes in the layer regions and points where the mesh switches from fine to coarse and vice versa are chosen to bexN/4+1=τ1 and x3N/4+1=1−τ2. The step sizes are related by(16)h1<h2<…hN/4,hN/4+1=hN/4+2=…=h3N/4,h3N/4+1>h3N/4+2>…>hN.Lemma 6.
The mesh sizeshn=xn+1−xn,n=1,2,…,N satisfy(17)hn≤Cη1−1,n=1,2,…,N/4,CN−1,n=N/4+1,…,3N/4,Cη2−1,n=3N/4+1,…,N.Proof.
Forn=1,2,…,N/4,(18)hn=xn+1−xn=−2η1log1−4N1−1η1n−−2η1log1−4N1−1η1n−1,=−2η1ψ1ξn+2η1ψ1ξn−1=−2η1ψ1ξn−ψ1ξn−1.
However, from the mean value theorem, there existsςn∈ξn−1,ξn such that(19)ψ1′ςn=ψ1ξn−ψ1ξn−1/ξn−ξn−1.
This giveshn≤Cη1−1,n=1,2,…,N/4.
Forn=N/4+1,…,3N/4, from equation (15), we have hn=xn+1−xn=2/N1−τ1−τ2≤CN−1. The proof in the third subinterval can be shown in a similar way as in the first subinterval.


## 3.2. Total Difference Scheme

LetQx,tm interpolate the solution Uxn,tm=Unm and pass through the points xn,tm,Unm and xn+1,tm,Un+1m such that(20)Qx,tm=−hn2θ2sinθsinθx−xnhnYn+1m+sinθxn+1−xhnYnm+hn2θ2x−xnhnYn+1m+hn2θ2Un+1m+xn+1−xhnYnm+θ2hn2Unm,where Ynm=U'′nm and θ is a free parameter used to maintain consistency. The first derivatives of Qx,tm from the right and left of x=xn are given, respectively, by(21)Q′xn+,tm=Un+1m−Unmhn+hnθ21−θcscθYn+1m−1−θcotθYnm,Q′xn−,tm=Unm−Un−1mhn−1+hn−1θ21−θcotθYnm−1−θcscθYn−1m.Then, from the continuity condition, we haveQ′xn+,tm=Q′xn−,tm and obtain(22)λ1hn−1Yn−1m+λ2hn−1+hnYnm+λ1hnYn+1m=Un+1m−Unmhn−Unm−Un−1mhn−1,where λ1=θcscθ−1/θ2,λ2=1−θcotθ/θ2. Application of the L’Hopital’s rule gives(23)λ1+λ2=1/2.Atx=xn−1,xn,xn+1 using the notation Wnm=U'′nm+U'′nm−1, from equation (7), we can have(24)Wn−1m=1ε−μan−1mU′n−1m+an−1m−1U′n−1m−1+qn−1mUn−1m+pn−1m−1Un−1m−1−Hn−1m,Wnm=1ε−μanmU′nm+anm−1U′nm−1+qnmUnm+pnm−1Unm−1−Hnm,Wn+1m=1ε−μan+1mU′n+1m+an+1m−1U′n+1m−1+qn+1mUn+1m+pn+1m−1Un+1m−1−Hn+1m,where, from [25],(25)U′n−1m=1hn−1hnhn−1+hn−hn2+2hn−1hnUn−1m+hn−1+hn2Unm−hn−12Un+1m,U′nm=1hn−1hnhn−1+hn−hn2Un−1m+hn2−hn−12Unm+hn−12Un+1m,U′n+1m=1hn−1hnhn−1+hnhn2Un−1m−hn−1+hn2Unm+hn−12+2hn−1hnUn+1m.Now, adopting the cubic spline in tension used in [25] yields(26)λ1hn−1Wn−1m+λ2hn−1+hnWnm+λ1hnWn+1m=Un+1m−Unmhn−Unm−Un−1mhn−1+Un+1m−1−Unm−1hn−Unm−1−Un−1m−1hn−1.Substituting equations (24) and (25) into equation (26) and writing it in compact form, we arrive at(27)Lε,μM,NUnm≔An−Un−1m+An0Unm+An+Un+1m+Bn−Un−1m−1+Bn0Unm−1+Bn+Un+1m−1=Fnm,where n=2,3,…,N,m=2,3,…,M,(28)An−=μλ1hn+2hn−1hn−1+hnan−1m+λ2hnhn−1anm−λ1hn2hn−1hn−1+hnan+1m+λ1hn−1qn−1m−εhn−1,An0=μ−λ1hn+hn−1hnan−1m−λ2hn2−hn−12hn−1hnanm+λ1hn−1+hnhn−1an+1m+λ2hn−1+hnqnm+εhn+hn−1hn−1hn,An+=μλ1hn−12hnhn−1+hnan−1m−λ2hn−1hnanm−λ1hn−1+2hnhn−1+hnan+1m+λ1hnqn+1m−εhn,Bn−=μλ1hn+2hn−1hn−1+hnan−1m−1+λ2hnhn−1anm−1−λ1hn2hn−1hn−1+hnan+1m−1+λ1hn−1pn−1m−1−εhn−1,Bn0=μ−λ1hn+hn−1hnan−1m−1−λ2hn2−hn−12hn−1hnanm−1+λ1hn−1+hnhn−1an+1m−1+λ2hn−1+hnpnm−1+εhn+hn−1hn−1hn,Bn+=μλ1hn−12hnhn−1+hnan−1m−1−λ2hn−1hnanm−1−λ1hn−1+2hnhn−1+hnan+1m−1+λ1hnpn+1m−1−εhn,Fnm=λ1hn−1Hn−1m+λ2hn+hn−1Hnm+λ1hnHn+1m.The system (27) together with the given conditions is uniquely solvable. Furthermore, the matrix associated with this system of equations satisfies the property of an M-matrix that is given in [26].


## 4. Convergence and Stability Analysis

Obtaining stability conditions for the Crank–Nicolson by using either the Fourier analysis or the matrix method is subjected to severe restrictions. This is worse for problems with variable coefficients like in our case. Thus, the maximum principle analysis can be taken as an alternative means [27]. So, to establish the parameter uniform convergence of the proposed scheme, let us see the following lemmas.Lemma 7.
Assume that(29a)μhn+hn2k≤Cε,1≤n≤N4,(29b)hn2k≤Cε+μhn,3N4+1≤n≤N.Then, under these assumptions, the matrix associated with equation (27) at each time level is an M-matrix. Hence, the operator Lε,μM,N satisfies the following discrete maximum principle.Lemma 8.
(discrete maximum principle). Let the discrete functionΨnm satisfy Ψ1m≥0 and ΨN+1m≥0. Then, Lε,μM,NΨnm≥0,n=2,…,N implies Ψnm≥0,∀n∈1,2,…,N+1.Proof.
To follow the proof by contradiction, let there exist a mesh pointl,m for l∈1,2,…,N+1 such that(30)Ψlm=min1≤n≤N+1Ψnm,and suppose Ψlm≤0. Then, this indicates l≠1,N+1, and for l∈2,3,…,N, we obtain(31)Lε,μM,NΨlm≤0,which contradicts the assumption Lε,μM,NΨnm≥0,n=2,…,N. This follows that Ψnm≥0,∀n∈1,2,…,N+1.Lemma 9.
For any discrete functionΦ defined on the tensor product ΩxN×ΩtM such that Φxn,tm=0 for all xn,tm on the boundary of ΩxN×ΩtM, then(32)Φxn,tm≤1q∗max1≤n≤N+1,1≤m≤M+1Lε,μN,MΦxn,tm,for qnm≥q∗>0.Proof.
The proof follows from the application of the discrete maximum principle to a mesh function given by(33)Z±xn,tm=1q∗max1≤n≤N+1,1≤m≤M+1Lε,μN,MΦxn,tm±Φxn,tm.Lemma 10.
The truncation error in the spatial direction, denoted byTEx, satisfies(34)TEx≤CN−2.Proof.
Let the truncation error in the spatial variable with fixed time be given by(35)TEx=An−Un−1m+An0Unm+An+Un+1m+Bn−Un−1m−1+Bn0Unm−1+Bn+Un+1m−1−Fnm.
The Taylor series expansion of each term in space up to the third order derivative and collectingUnm,U′nm,U'′nm,U'′'nm,Unm−1,U′nm−1,U'′nm−1, and U'′'nm−1 give(36)TEx=ζ0,nmUnm+ζ1,nmU′nm+ζ2,nmU″nm+ζ3,nmU′″nm+ζ0,nm−1Unm−1+ζ1,nm−1U′nm−1+ζ2,nm−1U‴nm−1+ζ3,nm−1U‴nm−1,where upon simplification and restricting the expansion of the coefficients to their first term yield(37)ζ0,nm=ζ1,nm=ζ0,nm−1=ζ1,nm−1=0,ζ2,nm=ζ2,nm−1=ε−12+λ1+λ2hn−1+hn.
Using the relationηj−1logηj≤CN−1,j=1,2, Lemma 4 and equations (17) and (23), it can be further simplified to(38)TEx≤CN−2.
As the operatorLε,μM,N satisfies the discrete maximum principle on the tensor product DtM×DxN and the error is estimated, the current method is uniformly stable in the maximum norm. Applying the triangular inequality from equations (12) and (38), we can have the following theorem.Theorem 1.
Letux,t,x,t∈D¯ be the solution of equations (1) and (2) and Uxn,tm, xn,tm∈DxN×DtM be the solution of equation (13). Then, the error estimate for the fully discrete scheme is given by(39)maxn,muxn,tm−Uxn,tm≤CN−2+k2.Therefore, the presented method is convergent independent of the perturbation parameters and its rate of convergence is two.


## 5. Computational Results and Discussion

To verify the theoretical results experimentally, we solve the following examples.Example 1.
For the problem in [28],(40)∂u∂t−ε∂2u∂x2−μ1+x∂u∂x+u=−16x21−x2,subject toux,0=0,u0,t=0=u1,t.Example 2.
Consider the following singular perturbation initial-boundary-value problem in [28],(41)∂u∂t−ε∂2u∂x2−μ1+expx∂u∂x+1+x5u=−10expt2x21−x2,subject toux,0=0,u0,t=0=u1,t.Example 3.
Consider the following singular perturbation initial-boundary-value problem in [2],(42)∂u∂t−ε∂2u∂x2−μ1+x−x2+t2∂u∂x+1+5xtu=x2−xet−1,subject toux,0=0,u0,t=0=u1,t.
As the considered examples have no exact solution, we calculate the absolute maximum errors using the double mesh principle as follows:(43)Eε,μh,k=max1≤n≤N+1,1≤m≤M+1Uhk−Uh/2k/2,where Uhk is an approximate solution obtained using k and h step sizes in the t and x directions, respectively, and Uh/2k/2 is an approximate solution obtained by bisecting the step sizes. As well, the corresponding rate of convergence is defined by(44)R=logEε,μh,k−logEε,μh/2,k/2log2.
The maximum absolute error and rate of convergence for Examples1–3 are given in Tables 1–3. From these results, one can observe that the current method converges independently of the perturbation parameters and gives more accurate numerical results than the existing results available in the literature. These results are also in excellent agreement with the theoretical findings. The solution profiles illustrated in Figures 1–3 indicate the meshes in the layer regions are more condensed and refined than the outer layer. From the log-log plots in Figures 4–6, we can conclude that our method is parameter uniform.Table 1
Comparison of the maximum pointwise error and rate of convergence for Example1 for μ=10−6 and different values of ε.
ε↓N=32N=64N=128N=256N=512M=10M=20M=40M=80M=160Current method10−26.2022e − 041.8723e − 049.8684e − 055.2922e − 053.0726e − 051.72800.92390.89890.784410−42.2942e − 045.7523e − 051.4390e − 057.5990e − 062.7100e − 061.99581.99910.92121.487510−62.2902e − 045.7428e − 051.4367e − 053.5924e − 068.9814e − 071.99561.99901.99971.999910−82.2902e − 045.7427e − 051.4367e − 053.5924e − 068.9813e − 071.99571.99901.99972.000010−102.2902e − 045.7427e − 051.4367e − 053.5924e − 068.9813e − 071.99571.99901.99972.000010−122.2902e − 045.7427e − 051.4367e − 053.5924e − 068.9813e − 071.99571.99901.99972.0000Method in [28]10−28.6053e − 034.4951e − 032.2857e − 031.1438e − 035.7093e − 040.9370.9760.9991.00210−48.6095e − 034.4529e − 032.2631e − 031.1406e − 035.7256e − 040.9510.9760.9880.99410−68.6014e − 034.4513e − 032.2628e − 031.1406e − 035.7356e − 040.9510.9760.9880.99410−88.6006e − 034.4512e − 032.2628e − 031.1406e − 035.7356e − 040.9510.9760.9880.99410−108.6006e − 034.4512e − 032.2628e − 031.1406e − 035.7356e − 040.9510.9760.9880.99410−128.6006e − 034.4512e − 032.2628e − 031.1406e − 035.7356e − 040.9510.9760.9880.994Table 2
Comparison of the maximum pointwise error and rate of convergence for Example2 for μ=10−6 and different values of ε.
ε↓N=32N=64N=128N=256N=512M=10M=20M=40M=80M=160Current method10−23.7393e − 031.0548e − 032.8089e − 047.2507e − 051.8422e − 051.82581.90891.95381.976710−43.9293e − 031.1111e − 032.9595e − 047.6404e − 051.9412e − 051.82231.90861.95361.976710−63.9337e − 031.1117e − 032.9609e − 047.6436e − 051.9423e − 051.82311.90871.95371.976510−83.9337e − 031.1117e − 032.9608e − 047.6439e − 051.9423e − 051.82311.90871.95361.976510−103.9337e − 031.1117e − 032.9608e − 047.6439e − 051.9423e − 051.82311.90881.95361.976510−123.9337e − 031.1117e − 032.9608e − 047.6439e − 051.9423e − 051.82311.90881.95361.9765Method in [28]10−23.6825e − 021.8188e − 028.5040e − 034.2227e − 032.1179e − 031.0181.0971.0100.99510−43.9442e − 021.9359e − 029.5692e − 034.7539e − 032.3691e − 031.0271.0161.0091.00510−63.9402e − 021.9391e − 029.5773e − 034.7594e − 032.3717e − 031.0231.0181.0091.00510−83.9418e − 021.9392e − 029.5791e − 034.7594e − 032.3718e − 031.0231.0171.0091.00510−103.9418e − 021.9392e − 029.5791e − 034.7594e − 032.3718e − 031.0231.0171.0091.00510−123.9418e − 021.9392e − 029.5791e − 034.7594e − 032.3718e − 031.0231.0171.0091.005Table 3
Comparison of the maximum pointwise error and rate of convergence for Example3 for μ=10−7 and different values of ε.
ε↓N=64N=128N=256N=512M=16M=32M=64M=128Current method10−62.6347e − 056.5876e − 061.6474e − 064.1182e − 071.99981.99962.00011.999910−72.6346e − 056.5877e − 061.6474e − 064.1182e − 071.99971.99962.00011.999910−82.6346e − 056.5877e − 061.6474e − 064.1182e − 071.99971.99962.00011.999910−92.6346e − 056.5877e − 061.6474e − 064.1182e − 071.99971.99962.00011.9999Method in [2]10−63.8754E − 51.0214E − 52.6170E − 66.6241E − 71.92381.96461.98211.991010−73.8753E − 53.8753E − 52.6170E − 66.6241E − 71.92381.96461.98211.991010−83.8753E − 51.0214E − 52.6170E − 66.6241E − 71.92381.96461.98211.991010−93.8753E − 51.0214E − 52.5461E − 66.6241E − 71.92381.96461.98211.9910Figure 1
The numerical solution profile for Example1 when N=64,M=64,ε=10−2,andμ=10−10.Figure 2
The numerical solution profile for Example2 when N=64,M=64,ε=10−2,andμ=10−10.Figure 3
The numerical solution profile for Example3 when N=64,M=64,ε=10−2,andμ=10−10.Figure 4
Log-log plot for Example1 when μ=10−6.Figure 5
Log-log plot for Example2 when μ=10−6.Figure 6
Log-log plot for Example3 when μ=10−7.


## 6. Conclusions

The present study provides the numerical solutions for a singularly perturbed unsteady-state initial-boundary-value problem with two small parameters. The method comprises of developing a scheme through discretization of the time variable by the Crank–Nicolson method on a uniform mesh and the space variable by cubic spline in tension with Bakhvalov-type mesh. The method is flexible and successful in capturing the sudden change of the solution behavior of the problem in the boundary layer region.



---

*Source: 5410754-2022-02-15.xml*