# The Novel Rational Spectral Collocation Method for a Singularly Perturbed Reaction-Diffusion Problem with Nonsmooth Data

Authors: Lufeng Yang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

A singularly perturbed reaction-diffusion problem with a discontinuous source term is considered. A novel rational spectral collocation combined with a singularity-separated method for this problem is presented. The solution is expressed asu=w+v, where w is the solution of corresponding auxiliary boundary value problem and v is a singular correction with direct expressions. The rational spectral collocation method combined with a sinh transformation is applied to solve the weakened singularly boundary value problem. According to the asymptotic analysis, the sinh transformation parameters can be determined by the width and position of the boundary layers. The parameters in the singular correction can be determined by the boundary conditions of the original problem. Numerical experiment supports theoretical results and shows that compared with previous research results, the novel method has the advantages of a high computational accuracy in singularly perturbed reaction-diffusion problems with nonsmooth data.


---

## Body

## 1. Introduction

Singular perturbation problems arise in mathematical modeling of physical and engineering problems, such as the boundary layer of fluid mechanics, the turning point of quantum mechanics, and the flow of large Reynolds numbers. In recent decades, the singular perturbation problem has received extensive attention. Detailed theories and analysis of the singular perturbation problem can be found in the literature [1, 2]. These problems have steep gradients in narrow layers, which is a serious obstacle in calculating classical numerical methods.Many studies of singularly perturbed problems with nonsmooth data have been researched. The first-order Schwarz method with uniform parameters on Shishkin’s mesh was proposed for a singularly perturbed reaction-diffusion problem with a nonsmooth source term [3].Chandru and Shanthi applied a boundary value technique and hybrid difference scheme for singularly perturbed boundary value problem of reaction-diffusion type with discontinuous source term in [4, 5], respectively. The second-order finite element method was presented on a Bakhvalov–Shishkin’s mesh for singularly perturbed problems with the interior layer [6]. A uniformly convergent difference scheme was employed for singularly perturbed semilinear problems with a discontinuous source term and proved that this scheme was the first-order convergence [7]. A uniformly numerical convergence scheme based on piecewise-uniform Shishkin’s meshes was developed for a reaction-diffusion equation with a discontinuous diffusion coefficient and proved to have second-order convergence [8]. A hybrid difference scheme on a Shishkin’s mesh was introduced for a singularly perturbed reaction-diffusion problem with one and two parameters [9] with interior layers.The second-order singularly perturbed reaction-diffusion boundary value problem with discontinuous source term is considered:(1)Lux≔−εu″+cu=fx,x∈Ω−∪Ω+,u0=α,u1=β,where 0<ε≪1 is a small positive parameter, c>0. Assuming the function fx is sufficiently smooth in Ω, a jump discontinuity is expected at d∈Ω denoted by wd=wd+−wd−. In general, this discontinuity of f gives rise to interior layers in the first derivative of the exact solution. Here, Ω−=0,d,Ω+=d,1.The rational spectral collocation method was proposed in the literature [10]. A conformal map maps the collocation points clustered near the poles of −1,1, which means more points are near the poles. The parameters of the mapping depend on the position and width of the boundary layer. A rational spectral collocation in a barycentric form with sinh transform (RSC-sinh method) is applied to solve a coupled system of singularly perturbed problems and third-order singularly perturbed problems [11, 12].To weaken the singularity and improve the accuracy of numerical simulation, the singularity-separated method (SSM) for the singular perturbation problem with constant coefficients was proposed by Chen and Yang [13]. Finite element methods with SSM were used to solve a singular perturbation problem with a single boundary layer. A novel rational spectral collocation method is presented combined with the singularity-separated technique for a system of singularly perturbed boundary value problems [14].A novel numerical method based on the rational spectral collocation in a barycentric form with a singularly-separated method (RSC-SSM) is proposed for solving the second-order singularly perturbed reaction-diffusion problem with nonsmooth source term.This paper is organized as follows. The asymptotic analysis and sinh transform is outlined in Section2. The algorithmic details of RSC-SSM for second-order singularly perturbed boundary value problem are considered, and the error estimates for the method are discussed in Section 3. The singularly perturbed reaction-diffusion problem with discontinuous source term is solved in Section 4, which supports theoretical results and provides a favorable comparison with existing methods. Finally, the conclusions are drawn in Section 5.


## 2. Preliminaries

In this section, some useful lemmas are given, which provide information about the boundary and interior layers occurring in the solution of singularly perturbed problem (1).Theorem 1.
The problem (1) has a solution u∈C1Ω∩C2Ω−∪Ω+.Proof.
This theorem is proved constructively. Letu1 and u2 be particular solutions of the following differential equations, respectively:(2)−εu1″+bxu1=f−x,x∈Ω−,−εu2″+bxu2=f+x,x∈Ω+.
Consider the function(3)ux=u1x+u0−u10ϕ1x+Aϕ2x,x∈Ω−,u2x+Bϕ1x+u1−u21ϕ2x,x∈Ω+,where ϕ1x and ϕ2x are the solutions of the boundary value problems:(4)−εϕ1″+bxϕ1=0,x∈Ω,ϕ10=1,ϕ11=0,−εϕ2″+bxϕ2=0,x∈Ω,ϕ20=0,ϕ21=1,and A and B are constants to be chosen so that u∈C1Ω. Using the fact that the maximum principle states that the maximum and minimum values taken on by ϕi∈Ω¯,i=1,2 must occur on ∂Ω, then we have 0<ϕi<1,i=1,2, and (5)ϕ1′<0,ϕ2′>0,x∈Ω.
We need to choose the constantsA and B so that u∈C1Ω; then, we impose(6)ud−=ud+,u′d−=u′d+.
Additionally, the constants A and B satisfy(7)Aϕ2d−−Bϕ1d+=g0,Aϕ2′d−−Bϕ1′d+=g1,where(8)gi=u2id+−u1id−+u1−u21ϕ2id+−u0−u10ϕ1id−,i=0,1.
Since(9)ϕ2d−−ϕ1d+ϕ2′d−−ϕ1′d+=−ϕ1′d+ϕ2d−+ϕ1d+ϕ2′d−>0,the system of linear equation (7) has a solution.Lemma 1.
Let us suppose that a functionux satisfies u0≥0,u1≥0,Lux≥0,∀x∈Ω−∪Ω+ and ud=0,u′d≤0, then ux≥0,∀x∈Ω¯.Proof.
Letx0 be such that ux0=minΩ¯ux. If ux0>0, the conclusion is clearly correct. Suppose that ux0<0; it is clear that either x0∈Ω−∪Ω+ or x0=d. In the first case, u′x0=0,u″x0≥0 and(10)Lux0=−εu″x0+cux0<0,which is a contradiction. In the second case, the proof depends on whether or not u is differentiable at d. If u is not differentiable at d, then u′d≠0. Since u′d−≥0,u′d+≤0, it is clear that u′d>0, which is a contradiction. If u is differentiable at d, according to the mean value theorem, u∈C1Ω¯. Since ud<0, it follows that there exists a neighbourhood Ωh=d−h,d such that ux<0,∀x∈Ωh. Choose a point x1≠d,x1∈Ωh such that ux1>ud. Then, the mean value theorem indicates that ∃x2∈x1,d:(11)u′x2=ud−ux1d−x1<0,and similarly, there also exists x3∈Ωh:(12)u″x3=u′d−u′x2d−x2=−u′x2d−x2>0.
Thus, we have(13)Lux3=−εu″x3+cux3<0,which is a contradiction.
The direct application of the maximum principle is the following stability result.Theorem 2.
Letux be the solution of (1), which can be estimated as(14)uΩ¯≤maxα,β,fxb.Proof.
SettingM=maxα,β,|fx|/b, we can construct v=M±u satisfying(15)Lv=−εv″±bv=bM±Lu=bM±f≥0,v0=M±α≥0,v1=M±β≥0.
Furthermore, the solutionu of (1) has first-order continuous differentiability, that is,(16)vd=ud=0,v′d=u′d=0.
It follows from the maximum principle thatv≥0 for all x∈Ω¯; thus, we obtain(17)uΩ¯≤maxα,β,fxb.Theorem 3.
The problem (1) has a unique solution u∈C1Ω∩C2Ω−∪Ω+.Proof.
The existence of problem (1) has been proved in Theorem 1. Let u1x,u2x be the solution of problem (1), and u1x≠u2x,∀x∈0,1; set hx=u1x−u2x. Then, we have(18)Lhx=Lu1x−Lu2x=fx−fx=0,h0=h1=0,hd=h′d=0.
It follows at once from Theorem2, and we can obtain hx=0.
Thus, problem (1) has a unique solution.


### 2.1. Asymptotic Expansion Approximation

For the construction of the RSC-SSM method, it is necessary to use the exact solution of the problem to analyze the position and width of the boundary layer. The problem (1) is considered as the coupling of the following two subproblems: (19)PL:L−ux=−εu″+cu=f−x,x∈Ω−,u0=α,ud−=γ,(20)PR:L+ux=−εu″+cu=f+x,x∈Ω+,ud+=γ,u1=β,where γ is the value of u at d, which is determined by u′d−=u′d+.The Shishkin decomposition method splits the solution of the singular perturbation problem (1) into regular ω and layer v components. The regular component ω is defined as the solution of(21)Lωx=fx,x∈Ω−∪Ω+,ω0=f−0c,ωd−=f−d−c,ωd+=f+d+c,ω1=f+1c.The singular componentv is given by (22)Lvx=fx,x∈Ω−∪Ω+,vd=−ωd,v′d=−ω′d,v0=u0−ω0,v1=u1−ω1.Lemma 2.
For each integerk satisfying 0≤k≤4, the regular term ω and the singular term v satisfy the bounds(23)ωx≤C1+ε1−k/2E1x,x∈Ω−,C1+ε1−k/2E2x,x∈Ω+,vx≤Cε−k/2E1x,x∈Ω−,Cε−k/2E2x,x∈Ω+,where C is a constant independent of ε and(24)E1x=e−xc/ε+e−d−xc/ε,E2x=e−x−dc/ε+e−1−xc/ε.Proof.
This is found using the stability result in Theorem2 and the techniques in [3].Lemma 3.
(see [12]). Considering the singularly perturbed reaction-diffusion(25)−εu″+cu=fx,x∈0,1,u0=α,u1=β,where fx is sufficiently smooth in 0,1, the solution u of this problem has the following asymptotic expansion:(26)ux=wx+v¯x+Oε,where(27)v¯x=vx,x∈0,τ0∪1−τ0,1,0,x∈τ0,1−τ0,0<τ0=−εclnε≪1,and the following inequality is satisfied:(28)u−w+v¯≤Cε,where C is a generic constant.
Lemma3 indicates that the boundary layer region of the singularly perturbed reaction-diffusion problem with smooth source term is 0,τ0 and 1−τ0,1. In other words, the position of the boundary layer is at the two endpoints in the interval 0,1, and its width is δ=τ0.


### 2.2. Rational Spectral Collocation Method in the Barycentric Form

A rational interpolation functionpNx in the barycentric form of a function ux can be expressed as follows [15]:(29)ux≈pNx=∑k=0Nωk/x−xkuxk∑k=0Nωk/x−xk,where ωkk=0N are barycentric weights and xkk=0N are distinct interpolation points. Most important in practice are the so-called Chebyshev–Gauss–Lobatto points xk=−coskπ/N,k=0,1,…,N, and in this case, the barycentric weights are as follows [15]:(30)ω0=12,ωk=−1k,k=1,2,…,N−1,ωN=−1N2.Theorem 4.
(see [16]). Let D1 and D2 be domains in C containing J=−1,1 and a real interval I, respectively. Let g:D1⟶D2 be a conformal map such that gJ=I. If f:D2⟶C is such that f∘g:D1⟶C is analytic inside and on an ellipse with foci ±1, semimajor axis length L, and semiminor axis length l, then the rational function pNx, which interpolates f at the transformed Chebyshev points x^k=gxk satisfies(31)pNx−fx=OL+l−N,uniformly for all x∈−1,1.
As suggested in (31), the convergence rate of the rational spectral collocation method mainly depends on the analytic region of u in the complex plane. Thus, the conformal map g could be chosen to enlarge the ellipse of analyticity of u∘g. Then, compared with the Chebyshev spectral method, a better approximation of u could be obtained.
An advantage of the rational function in barycentric form is that its derivatives can be calculated directly using differentiation formulae instead of repeatedly using the differential quotient rule. Thenth derivative of the rational interpolating function pNx evaluated at the point xj can be expressed in the following form:(32)pNnxj=∑k=0Ndndxnωk/x−xkuxk∑l=0Nωl/x−xlx=xj=∑k=0NDjknuxk,j=0,1,…,N,where Djkn is the entry of the nth order differentiation matrix. The entries of the first- and second-order differentiation matrices are given as follows [15]:(33)Djk1=ωkωjxj−xk,j≠k,−∑i≠kDji1,j=k,Djk2=2Djk1Djj1−1xj−xk,j≠k,−∑i≠kDji2,j=k.
Note that differentiation matrices (33) only rely on weights ωk and the points xk, which is the reason why the underlying equation does not need to be converted to new coordinates after maps.


### 2.3. The Sinh Transform

To approximate the rapid changes in the boundary layer region, Tee and Trefethen have constructed the conformal map [17]:(34)gλ,δx=λ+μsinhsinh−11−λμ+sinh−11+λμx−12+sinh−11−λμ,where λ,μ are the location and width of the boundary layers, respectively. The transformed Chebyshev points gλ,δxkk=0N are clustered near the location of boundary layer x=λ, and the density is determined by the width of the boundary layer.To better distinguish the singular perturbation problem with two boundary layers, Tee proposed the combined sinh transform as(35)g˜μx=12g−1,μ2x+1−1,x∈−1,0,12g1,μ2x−1+1,x∈0,1.All derivatives of the piecewise mapg˜ at x=0 are continuous to preserve the spectral accuracy. For the reaction-diffusion type, the parameter in (35) should be chosen as μ=2τ0.


## 2.1. Asymptotic Expansion Approximation

For the construction of the RSC-SSM method, it is necessary to use the exact solution of the problem to analyze the position and width of the boundary layer. The problem (1) is considered as the coupling of the following two subproblems: (19)PL:L−ux=−εu″+cu=f−x,x∈Ω−,u0=α,ud−=γ,(20)PR:L+ux=−εu″+cu=f+x,x∈Ω+,ud+=γ,u1=β,where γ is the value of u at d, which is determined by u′d−=u′d+.The Shishkin decomposition method splits the solution of the singular perturbation problem (1) into regular ω and layer v components. The regular component ω is defined as the solution of(21)Lωx=fx,x∈Ω−∪Ω+,ω0=f−0c,ωd−=f−d−c,ωd+=f+d+c,ω1=f+1c.The singular componentv is given by (22)Lvx=fx,x∈Ω−∪Ω+,vd=−ωd,v′d=−ω′d,v0=u0−ω0,v1=u1−ω1.Lemma 2.
For each integerk satisfying 0≤k≤4, the regular term ω and the singular term v satisfy the bounds(23)ωx≤C1+ε1−k/2E1x,x∈Ω−,C1+ε1−k/2E2x,x∈Ω+,vx≤Cε−k/2E1x,x∈Ω−,Cε−k/2E2x,x∈Ω+,where C is a constant independent of ε and(24)E1x=e−xc/ε+e−d−xc/ε,E2x=e−x−dc/ε+e−1−xc/ε.Proof.
This is found using the stability result in Theorem2 and the techniques in [3].Lemma 3.
(see [12]). Considering the singularly perturbed reaction-diffusion(25)−εu″+cu=fx,x∈0,1,u0=α,u1=β,where fx is sufficiently smooth in 0,1, the solution u of this problem has the following asymptotic expansion:(26)ux=wx+v¯x+Oε,where(27)v¯x=vx,x∈0,τ0∪1−τ0,1,0,x∈τ0,1−τ0,0<τ0=−εclnε≪1,and the following inequality is satisfied:(28)u−w+v¯≤Cε,where C is a generic constant.
Lemma3 indicates that the boundary layer region of the singularly perturbed reaction-diffusion problem with smooth source term is 0,τ0 and 1−τ0,1. In other words, the position of the boundary layer is at the two endpoints in the interval 0,1, and its width is δ=τ0.


## 2.2. Rational Spectral Collocation Method in the Barycentric Form

A rational interpolation functionpNx in the barycentric form of a function ux can be expressed as follows [15]:(29)ux≈pNx=∑k=0Nωk/x−xkuxk∑k=0Nωk/x−xk,where ωkk=0N are barycentric weights and xkk=0N are distinct interpolation points. Most important in practice are the so-called Chebyshev–Gauss–Lobatto points xk=−coskπ/N,k=0,1,…,N, and in this case, the barycentric weights are as follows [15]:(30)ω0=12,ωk=−1k,k=1,2,…,N−1,ωN=−1N2.Theorem 4.
(see [16]). Let D1 and D2 be domains in C containing J=−1,1 and a real interval I, respectively. Let g:D1⟶D2 be a conformal map such that gJ=I. If f:D2⟶C is such that f∘g:D1⟶C is analytic inside and on an ellipse with foci ±1, semimajor axis length L, and semiminor axis length l, then the rational function pNx, which interpolates f at the transformed Chebyshev points x^k=gxk satisfies(31)pNx−fx=OL+l−N,uniformly for all x∈−1,1.
As suggested in (31), the convergence rate of the rational spectral collocation method mainly depends on the analytic region of u in the complex plane. Thus, the conformal map g could be chosen to enlarge the ellipse of analyticity of u∘g. Then, compared with the Chebyshev spectral method, a better approximation of u could be obtained.
An advantage of the rational function in barycentric form is that its derivatives can be calculated directly using differentiation formulae instead of repeatedly using the differential quotient rule. Thenth derivative of the rational interpolating function pNx evaluated at the point xj can be expressed in the following form:(32)pNnxj=∑k=0Ndndxnωk/x−xkuxk∑l=0Nωl/x−xlx=xj=∑k=0NDjknuxk,j=0,1,…,N,where Djkn is the entry of the nth order differentiation matrix. The entries of the first- and second-order differentiation matrices are given as follows [15]:(33)Djk1=ωkωjxj−xk,j≠k,−∑i≠kDji1,j=k,Djk2=2Djk1Djj1−1xj−xk,j≠k,−∑i≠kDji2,j=k.
Note that differentiation matrices (33) only rely on weights ωk and the points xk, which is the reason why the underlying equation does not need to be converted to new coordinates after maps.


## 2.3. The Sinh Transform

To approximate the rapid changes in the boundary layer region, Tee and Trefethen have constructed the conformal map [17]:(34)gλ,δx=λ+μsinhsinh−11−λμ+sinh−11+λμx−12+sinh−11−λμ,where λ,μ are the location and width of the boundary layers, respectively. The transformed Chebyshev points gλ,δxkk=0N are clustered near the location of boundary layer x=λ, and the density is determined by the width of the boundary layer.To better distinguish the singular perturbation problem with two boundary layers, Tee proposed the combined sinh transform as(35)g˜μx=12g−1,μ2x+1−1,x∈−1,0,12g1,μ2x−1+1,x∈0,1.All derivatives of the piecewise mapg˜ at x=0 are continuous to preserve the spectral accuracy. For the reaction-diffusion type, the parameter in (35) should be chosen as μ=2τ0.


## 3. The Rational Spectral Collocation with a Singularly-Separated Method

The implementation and error estimation of RSC-SSM for a singularly perturbed reaction-diffusion problem with a discontinuous source term is elaborated in this section.


### 3.1. The Singularity-Separated Technique

Consider the problem (1). Taking the left subproblem (19) as an example, the efficient numerical method can be constructed. The homogenous equation L−u=0 has two eigenvalues:(36)λ1=cε,λ2=−cε.Letu0 be a special solution of L−u=f; its general solution is(37)ux=u0x+C1ϕ1x+C2ϕ2x,ϕ1x=eλ1x−d,ϕ2x=eλ2x.Note that ifε=10−k and k≥4,ϕ11=ϕ20=eλ2≈e−1/ϵ≥e−102≈10−44=0 can be neglected.Thus, the solutionulx of (19) can be decomposed into two parts: ulx=ωlx+vlx, in which ωlx=u0x is the regular term and vlx=C1ϕ1x+C2ϕ2x is the singular term. The regular term ωlx is the solution of an auxiliary boundary value problem:(38)L−ωl=−εωl″+cωl=f−x,ωl0=f−0c,ωld=f−dc.It is known from the boundary conditions of equation (38) that εωl′′0=εωl′′d=0. The singularities of the solution are weakened.Similarly, the solutionurx of (20) can also be decomposed into two parts: urx=ωrx+vrx. The regular term ωrx is the solution of an auxiliary boundary value problem:(39)L+ωr=−εωr″+cωr=f+x,ωrd=f+dc,ωr1=f+1c.In the same way,εωr′′d=εωr′′1=0.Theorem 5.
The singular componentsvlx,vrx of (19) and (20) are(40)vlx=C1le−c/εx+C2lec/εx−d,vrx=C1re−c/εx−d+C2rec/εx−1,where(41)C1l=α−f−0c,C2l=12cε/cf′d+fd,C1r=12cε/cf′d−fd,C2r=β−f+1c,where wlx and wrx are the solutions of the auxiliary in (38) and (39).Proof.
According to the boundary condition of (1), it is known that(42)α=wl0+C1l,β=wr1+C2r.
We can obtainC1l=α−f−0/c,C2r=β−f+1/c. On the contrary, the solution of equation (1) ux∈C1Ω¯, which is ud−=d+,u′d−=u′d+. Thus, we have(43)wld−+C2l=wrd++C1r,wl′d−+C2lc/ε=wr′d+−C1rc/ε.
Furthermore,(44)wrd+−wld−=f+d+c−f−d−c=fdc,wr′d+−wl′d−=f′+d+c−f′−d−c=f′dc.
Thus, the parameters can be obtained:(45)C2l=12cε/cf′d+fd,C1r=12cε/cf′d−fd.
Thus, the solution of the singularly-separated technique can be obtained as(46)ussx=ωlx+vlx,x∈Ω−,ωrx+vrx,x∈Ω+,where the regular terms wlx and wrx are the solutions of the auxiliary in (38) and (39) and the singular terms vlx and vrx are defined as (40).


### 3.2. RSC-SSM Method

The key to the success of the RSC-SSM is obtaining a high-precision numerical solution for weakened singular perturbation problems (38). The rational spectral collocation with a sinh transform is selected.First, by introducing the linear transformx=d/2t+1,x∈0,d,t∈−1,1 and defining a new function u^t=ux=ud/2t+1, then u′x=2/du^′t,u″x=4/d2u^′′t, and we can rewrite (38) as(47)Lω^=−4d2εω^l″t+cω^lt=f^−t,t∈−1,1,ω^l−1=f^−−1c,ω^l1=f^−1c.The boundary layers are located at the two endpoints of the underlying interval−1,1, and its width is δ=2/dτ0.Let the transformed Chebyshev collocation points be(48)x^kk=0N=gδ−coskπNk=0N,where g is shown in (35).Evaluating (47) at points x^kk=0N yields(49)−4d2εD2Wl+cEN+1Wl=F,where(50)Wl=w0,w1,…,wNT,withwkevaluatingwtk,F=f^0,f^1,…,f^NT,withf^kevaluatingf^tk,and EN+1 is a unit matrix of N+1-order.LetA=−4/d2εD2+cEN+1, then equation (49) can be rewritten as(51)AWl=F.The boundary conditions in (47) suggest that(52)Wl1=f^−−1c,WlN+1=f^−1c.The Dirichlet boundary conditions are imposed on the first and last rows ofA. We substitute elements in these rows with zeros. The first zero in the first row and the last zero in the last row are replaced by 1. Simultaneously, the first and last elements of vector F are replaced by f−−1/c and f−1/c.By solving the linear algebra system of (51), we can obtain the numerical solution Wl, which is the regular term in (35). Similarly, the numerical solution Wr, which is the solution of (47), can be obtained.According to the boundary condition of the original problem (1) and the value of the regular components ωl and ωr, the parameters of singular terms (40) are obtained.The numerical solution of equation (1) can be obtained as(53)uNx=Wlx+vlx,x∈Ω−,Wrx+vrx,x∈Ω+.


### 3.3. Error Analysis

In this section, we establish the error estimates of the RSC-SSM (53).Theorem 6.
Letuxi and uNxi be the exact solution and numerical solution of the original problem (1), respectively. Then, we have(54)u−uN≤c1e−c/ε+c2ρ−N.Proof.
Letz=u−uN; then, the error z satisfies(55)z=u−uN=u−uSS+uSS−uN≤u−uSS+uSS−uN=u−uSS+maxωl−Wl,ωr−Wr.
Since the singular correct functionvx is the solution of the homogenous equations Lv=0, set R=u−uSS, then(56)LR=Lu−uSS=Lu−ω−v=−εu″+cu−−εω″+cω−Lv=0,and, according to the proof of Theorem 5, we have(57)Rd≈ud+−wrd+−C1r−ud−−wld−−C2l=0,R′d≈u′d+−wr′d++cεC1r−u′d−−wl′d−−cεC2l=0.
Simultaneously,(58)R0=u0−ωl0−vl0=α−f−0c−C1l−C2le−dc/ε=−C2le−dc/ε,R1=u1−ωr1−vr1=β−f+1c−C1re−1−dc/ε−C2r=−C1red−1c/ε.
According to Lemma2, we can obtain(59)R=u−uSS≤maxR0,R1=maxC2le−dc/ε,C1red−1c/ε≤12cfd+f′de−c/ε.
A combination of the convergence of rational spectral collocation (31) leads to(60)uSS−uN=maxωl−Wl,ωr−Wr≤c2ρ−N.
By considering (59) and (60), it leads to(61)u−uN≤u−uSS+uSS−uN≤c1e−c/ε+c2ρ−N.
Thus, we can establish a rational spectral method with a singularity-separated technique that can almost achieve spectral accuracy.


## 3.1. The Singularity-Separated Technique

Consider the problem (1). Taking the left subproblem (19) as an example, the efficient numerical method can be constructed. The homogenous equation L−u=0 has two eigenvalues:(36)λ1=cε,λ2=−cε.Letu0 be a special solution of L−u=f; its general solution is(37)ux=u0x+C1ϕ1x+C2ϕ2x,ϕ1x=eλ1x−d,ϕ2x=eλ2x.Note that ifε=10−k and k≥4,ϕ11=ϕ20=eλ2≈e−1/ϵ≥e−102≈10−44=0 can be neglected.Thus, the solutionulx of (19) can be decomposed into two parts: ulx=ωlx+vlx, in which ωlx=u0x is the regular term and vlx=C1ϕ1x+C2ϕ2x is the singular term. The regular term ωlx is the solution of an auxiliary boundary value problem:(38)L−ωl=−εωl″+cωl=f−x,ωl0=f−0c,ωld=f−dc.It is known from the boundary conditions of equation (38) that εωl′′0=εωl′′d=0. The singularities of the solution are weakened.Similarly, the solutionurx of (20) can also be decomposed into two parts: urx=ωrx+vrx. The regular term ωrx is the solution of an auxiliary boundary value problem:(39)L+ωr=−εωr″+cωr=f+x,ωrd=f+dc,ωr1=f+1c.In the same way,εωr′′d=εωr′′1=0.Theorem 5.
The singular componentsvlx,vrx of (19) and (20) are(40)vlx=C1le−c/εx+C2lec/εx−d,vrx=C1re−c/εx−d+C2rec/εx−1,where(41)C1l=α−f−0c,C2l=12cε/cf′d+fd,C1r=12cε/cf′d−fd,C2r=β−f+1c,where wlx and wrx are the solutions of the auxiliary in (38) and (39).Proof.
According to the boundary condition of (1), it is known that(42)α=wl0+C1l,β=wr1+C2r.
We can obtainC1l=α−f−0/c,C2r=β−f+1/c. On the contrary, the solution of equation (1) ux∈C1Ω¯, which is ud−=d+,u′d−=u′d+. Thus, we have(43)wld−+C2l=wrd++C1r,wl′d−+C2lc/ε=wr′d+−C1rc/ε.
Furthermore,(44)wrd+−wld−=f+d+c−f−d−c=fdc,wr′d+−wl′d−=f′+d+c−f′−d−c=f′dc.
Thus, the parameters can be obtained:(45)C2l=12cε/cf′d+fd,C1r=12cε/cf′d−fd.
Thus, the solution of the singularly-separated technique can be obtained as(46)ussx=ωlx+vlx,x∈Ω−,ωrx+vrx,x∈Ω+,where the regular terms wlx and wrx are the solutions of the auxiliary in (38) and (39) and the singular terms vlx and vrx are defined as (40).


## 3.2. RSC-SSM Method

The key to the success of the RSC-SSM is obtaining a high-precision numerical solution for weakened singular perturbation problems (38). The rational spectral collocation with a sinh transform is selected.First, by introducing the linear transformx=d/2t+1,x∈0,d,t∈−1,1 and defining a new function u^t=ux=ud/2t+1, then u′x=2/du^′t,u″x=4/d2u^′′t, and we can rewrite (38) as(47)Lω^=−4d2εω^l″t+cω^lt=f^−t,t∈−1,1,ω^l−1=f^−−1c,ω^l1=f^−1c.The boundary layers are located at the two endpoints of the underlying interval−1,1, and its width is δ=2/dτ0.Let the transformed Chebyshev collocation points be(48)x^kk=0N=gδ−coskπNk=0N,where g is shown in (35).Evaluating (47) at points x^kk=0N yields(49)−4d2εD2Wl+cEN+1Wl=F,where(50)Wl=w0,w1,…,wNT,withwkevaluatingwtk,F=f^0,f^1,…,f^NT,withf^kevaluatingf^tk,and EN+1 is a unit matrix of N+1-order.LetA=−4/d2εD2+cEN+1, then equation (49) can be rewritten as(51)AWl=F.The boundary conditions in (47) suggest that(52)Wl1=f^−−1c,WlN+1=f^−1c.The Dirichlet boundary conditions are imposed on the first and last rows ofA. We substitute elements in these rows with zeros. The first zero in the first row and the last zero in the last row are replaced by 1. Simultaneously, the first and last elements of vector F are replaced by f−−1/c and f−1/c.By solving the linear algebra system of (51), we can obtain the numerical solution Wl, which is the regular term in (35). Similarly, the numerical solution Wr, which is the solution of (47), can be obtained.According to the boundary condition of the original problem (1) and the value of the regular components ωl and ωr, the parameters of singular terms (40) are obtained.The numerical solution of equation (1) can be obtained as(53)uNx=Wlx+vlx,x∈Ω−,Wrx+vrx,x∈Ω+.


## 3.3. Error Analysis

In this section, we establish the error estimates of the RSC-SSM (53).Theorem 6.
Letuxi and uNxi be the exact solution and numerical solution of the original problem (1), respectively. Then, we have(54)u−uN≤c1e−c/ε+c2ρ−N.Proof.
Letz=u−uN; then, the error z satisfies(55)z=u−uN=u−uSS+uSS−uN≤u−uSS+uSS−uN=u−uSS+maxωl−Wl,ωr−Wr.
Since the singular correct functionvx is the solution of the homogenous equations Lv=0, set R=u−uSS, then(56)LR=Lu−uSS=Lu−ω−v=−εu″+cu−−εω″+cω−Lv=0,and, according to the proof of Theorem 5, we have(57)Rd≈ud+−wrd+−C1r−ud−−wld−−C2l=0,R′d≈u′d+−wr′d++cεC1r−u′d−−wl′d−−cεC2l=0.
Simultaneously,(58)R0=u0−ωl0−vl0=α−f−0c−C1l−C2le−dc/ε=−C2le−dc/ε,R1=u1−ωr1−vr1=β−f+1c−C1re−1−dc/ε−C2r=−C1red−1c/ε.
According to Lemma2, we can obtain(59)R=u−uSS≤maxR0,R1=maxC2le−dc/ε,C1red−1c/ε≤12cfd+f′de−c/ε.
A combination of the convergence of rational spectral collocation (31) leads to(60)uSS−uN=maxωl−Wl,ωr−Wr≤c2ρ−N.
By considering (59) and (60), it leads to(61)u−uN≤u−uSS+uSS−uN≤c1e−c/ε+c2ρ−N.
Thus, we can establish a rational spectral method with a singularity-separated technique that can almost achieve spectral accuracy.


## 4. Numerical Experiments

To demonstrate the accuracy and efficiency of RSC-SSM, we will use the novel method to solve a singularly perturbed reaction-diffusion problem discontinuous source term and compare the results with other existing methods. We verify the theoretical results obtained in the previous section through numerical experiments. The maximum relative errors of the solution are given by(62)e=uN−u∞uE∞,where uN and u are the numerical and exact solutions, respectively.In our computations, all experiments are carried out using MATLAB (version R2014a) on a PC with a 2.5 GHz central processing unit (Intel Core i5-2450M), 4.00 GB memory, and Windows 7 operating system.Example. Consider the following singularly perturbed problem with discontinuous source term:(63)−ε2u″x+ux=fx,x∈Ω−∪Ω+,u0=y1=f0,where(64)fx=−0.5x,0≤x≤0.5,0.5,0.5<x≤1.The exact solution of this problem can be expressed as(65)ux=0.25A+Be−0.5−x/ε−e−0.5+x/ε−0.5x,0≤x≤0.5,0.25A−Be1−2x/2ε−e2x−3/2ε+0.51−e−1−x/ε,0.5<x≤1,where A and B are constants related to ε:(66)A=ε−e−1/2ε1+e−1/ε,B=1.5−e−1/2ε1−e−1/ε.Obviously, there is an interior layer atx=0.5 in the solution.In the RSC-SSM method, the parameters in combined transform of (35) are chosen as μ=−4εlnε.The RSC-SSM method and the different schemes in [9, 18] are used to solve this problem with several choices of N and ε. The maximum relative errors of the three methods are listed in Table 1.Table 1
Maximum pointwise errors for RSC-SSM and finite different schemes [9, 18].
ε/NRSC-SSMFinite difference [9]Hybrid scheme [18]326412832641281e − 21.07e − 133.49e − 041.77e − 04∗∗∗1e − 32.19e − 147.30e − 039.42e − 04∗∗∗1e − 43.61e − 157.30e − 039.42e − 04∗∗∗1e − 63.11e − 157.30e − 039.42e − 04∗∗∗2−46.75e − 02∗∗1.72e − 038.25e − 044.05e − 042−81.66e − 04∗∗6.45e − 042.49e − 041.08e − 042−121.00e − 13∗∗2.69e − 031.73e − 036.89e − 042−164.95e − 14∗∗4.29e − 032.04e − 039.95e − 042−201.78e − 14∗∗7.23e − 032.52e − 031.04e − 03∗References [9, 18] do not compute.Compared with the existing difference schemes [9, 18], the RSC-SSM accuracy has been significantly improved for the case with ε2≤10−3. However, when the parameters ε2≥10−2, the accuracy of the RSC-SSM is similar to that of the difference method. This is because that the discarded item reached e−1/ε=e−10≈4.5×10−5 in singularly-separated technique. In Table 1, we investigate the effects of the singular perturbation parameter ε on the errors. We can see that a decrease in ε reduces the maximum relative errors.Figure1 depicts the numerical solution and the exact solution of the whole region with different ε. The numerical solution and exact solutions in both the interior and the boundary layers are given in Figure 2.Figure 1
Numerical solution vs. exact solution in Example (N=32): (a) ε = 10−1; (b) ε = 10−4.
(a)(b)Figure 2
Numerical solution vs. exact solution in Example (ε=10−4,N=32): (a) interior layer region; (b) boundary layer region.
(a)(b)Figure3 shows the maximum relative errors in the semilog scale with ε=10−1,10−2,10−3, and 10−4, respectively. It demonstrates the influence of points N on the errors. The convergence rates are almost balanced, and even the errors have increased due to the cumulative effects with an increase in N. Figure 4 shows the pointwise errors of the function in the whole region with ε=10−1,10−2,10−3, and 10−4, respectively. The smaller the parameter ε is, the more accurate the RSC-SSM result. Figures 3 and 4 also show that the error mainly comes from the discarded items in the process of singular separation. This confirms the results discussed in the previous section.Figure 3
The convergence rate.Figure 4
The pointwise error (N=32).


## 5. Conclusions

In this paper, a novel numerical method, named RSC-SSM, has been proposed to solve singularly perturbed boundary value problems with discontinuous source term. The solution is composed of the weaker singularity auxiliary solution and a singular correct function. The numerical experiment illustrates that the proposed RSC-SSM is obviously a better choice to other existing numerical methods available for the singularly perturbed boundary value problems with discontinuous source term. The numerical results demonstrate the almost spectral accuracy of the proposed algorithm and coincide with the theoretical analysis. Moreover, the theoretical and numerical frameworks presented in this paper can be extended to more complex problems.



---

*Source: 8844163-2021-02-26.xml*