# Rational Spectral Collocation Combined with the Singularity Separated Method for a System of Singularly Perturbed Boundary Value Problems

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

---

## Abstract

A novel rational spectral collocation method is presented combined with the singularity-separated technique for a system of singularly perturbed boundary value problems. The solution is expressed asu=w+v, where w is the solution of the corresponding auxiliary boundary value problem and v is a singular correction with explicit expressions. The rational spectral collocation method in barycentric form with the sinh transformation is applied to solve the auxiliary third boundary problem. The parameters of the singular correction can be determined by the boundary conditions of the original problem. Numerical experiments are carried out to support theoretical results and provide a favorable comparison with research results of other work.


---

## Body

## 1. Introduction

The singular perturbation problems (SPPs) such as fluid mechanics boundary layers, quantum mechanics turning points, and flow of large Reynolds numbers arise in the mathematical modeling physical and engineering problems. In the past decades, singular perturbation has received extensive attention. Singular perturbation boundary value problems have steep gradients in the narrow layers, which was a serious obstacle in the calculation of classical numerical methods. There is a vast literature dealing with SPPs with smooth coefficients and source term for single equation [1–3]. Compared with single equation singularly perturbed problems, coupled systems can simulate more complicated physical phenomena. Only a few authors have developed numerical methods to deal with the problem of coupled singular perturbation system on smooth data. A finite difference method on a piecewise uniform mesh for the reaction-diffusion type was proposed [4]. The result of Shishkin’s method in terms of convergence, stability, and error estimation can be found in Madden and Stynes [5], Linss and Madden [6, 7], and Matthews et al. [8]. An upwinding finite difference scheme on piecewise-uniform Shishkin meshes for the convection-diffusion type was presented, and it was solved by Jacobi iteration to compute the solution [9].In this paper, a coupled system ofm≥2 singularly perturbed linear equations are considered in the unknown vector function u=u1,…,umT:(1)Lu:=−εu″−Bu′+Au=f,x∈Ω=0,1,and it yields(2)u0=b0,u1=b1,where A=aij and B=bij are m×m matrices and 0<ε≪1 is a small parameter whose presence makes the problem singularly perturbed. We assumed that f=f1,…,fmT∈C20,1m and both b0=b01,…,b0mT and b1=b11,…,b1mT are constant vectors. These layer behaviors can be examined in two different cases:Case 1. B≡0, the system is called reaction-diffusion type. Generally speaking, the solution of this type has two boundary layers with width Oε at x=0 and x=1 under some assumptions.Case 2. B≠0, the system is called convection-diffusion type, the solution of this type has a single boundary layer with width Oε at x=0orx=1 under proper assumptions.If (1) is coupled through its convective terms, we say it is strongly coupled; otherwise, if B≡0 or B is just a nonzero diagonal matrix, it is said to be weakly coupled.The rational spectral collocation method was proposed in the literature [10]. A conformal map which maps the collocation points clustered near the poles of −1,1 into a new set of collocation points. The parameters of the mapping are determined by the position and width of the boundary layer. Chen and Wang [11, 12] applied a rational spectral collocation in barycentric form with the sinh transform (RSCAT) method to solve a coupled system of singularly perturbed problems and third-order singularly perturbed problems.To weaken the singularity and improve the accuracy of numerical simulation, the singularity-separated technique (SST) for singular perturbation problem with constant coefficients was proposed by Chen and Yang [13], and finite element methods with SST were used to solve a singular perturbation problem with a single boundary layer.Here, we present a novel numerical method based on rational spectral collocation in barycentric form with the singularly-separated method (RSC-SSM) to solve singularly perturbed boundary value problem in various types, both weakly coupled and strongly coupled.This paper is organized as follows. The asymptotic analysis of coupled system is outlined in Section2. The algorithmic details of the RSC-SSM for a coupled system of singularly perturbed problems are provided and the error estimates for the method are discussed in Section 3. The coupled system of singularly perturbed problems is solved in Section 4, which supports theoretical results and provides a favorable comparison with existing methods. Finally, we present some concluding remarks in Section 5.


## 2. Preliminaries

For the construction of the RSC-SSM, it is necessary to understand the properties of exact solution, especially the position and width of the boundary layer. Some useful lemmas are listed, including the maximum principle, stability result, the error estimation of the solution, and its derivatives are established for the boundary value problem (1).


### 2.1. Weakly Coupled System of Reaction-Diffusion Problems

IfA≠0,B≡0 in (1), it is convenient to introduce the notation to denote the reaction-diffusion operator:(3)LRu:=−εu″+Au=f.Assume thatA has positive diagonal entries and nonpositive off-diagonal entries, i.e.,(4)aij=>0,i=j,≤0,i≠j,and A is also strictly diagonally dominant i.e., ∑k=1k≠imaik<aii and ∑j=1maij>α2>0,i=1,…,m. These hypotheses ensure that problem (1) has a unique solution.Lemma 1.
Under above assumptions, all eigenvalues of matrixA are positive.Proof.
LetA have an eigenvalue of λ=0, according to the Gerschgorin disc theorem, ∃i0,ai0i0≤Ri0A, which is a contradiction. So, there is no zero eigenvalue in A, that is, A≠0.
Setfx=D+xB,x∈0,1, where D=diagA,B=A−D, then fx is a continuous function on 0,1 and f0=D=∏i=1maii>0; assume f1=A<0, then we can obtained ∃ξ∈0,1, such that fx=D+ξB=0. Moreover, A is strictly diagonally dominant, then aii>RiA≥∑j=1,j≠inξaij, that is, D+ξB is also strictly diagonally dominant. Contradictions between fξ=D+ξB≠0 and fξ=0 can be obtained. Therefore, f1=A>0, each of the main determinants of the strictly diagonally dominant matrix is the strictly diagonally dominant, then we have the A is positive definite matrix. So, all eigenvalues of matrix A are positive.Lemma 2.
Let us suppose that a functionux satisfies u0≥0,u1≥0,Lux≤0,∀x∈Ω, then ux≥0;∀x∈Ω¯.
The direct application of the maximum principle is the following stability result.Lemma 3.
Ifux is the solution of (2), then we have the stability bound inequality:(5)u≤M=maxb0,b1+Cmaxfx,u′x≤CMε.
Shishkin decomposition is used for asymptotic analysis of problem (3) [14], which splits the solution into regular and layer components:(6)ux=wx+vx,where the regular component wx and the layer component vx are the solutions of the two problems:(7)LRw=−εw″+Aw=f,w0=A−1f0,w1=A−1f1,(8)LRv=−εw″+Av=f,v0=b0−w0,v1=b1−w1,then wx and vx have the following estimates.Lemma 4.
(see [15]). Let u=w+v be the solution of problem (3) and wandv be the solutions of (7) and (8), respectively. Then, we have(9)wj≤C1+ε2−j/2,vj≤Cε−j/2e−αx/ε+e−α1−x/ε,where j=0,1,2.Theorem 1.
The solution of (3) and (2) has the following asymptotic expansion:(10)ux=wx+vx¯+Oε,where(11)vx¯=vx,x∈0,x∗∪1−x∗,1,0,x∈x∗,1−x∗,0<x∗=−lnεαε≪1.
The above theorem suggests that the boundary layer regions of the solutionux to the reaction-diffusion problems are 0,x∗ and 1−x∗,1, that is, the boundary layers are located at the two endpoints of the underlying interval 0,1 and each of their width is x∗=−lnε/αε.


### 2.2. Strongly Coupled System of Convection-Diffusion Problem

IfB≠0, the notation LC is introduced to denote the convection-diffusion operator:(12)LCu:=−εu″−Bu′+Au=f.To ensure that there is a unique solution to this problem, the following assumptions need to be satisfied:B≠0 is strictly diagonally dominant and hence invertible and define β=minbii,i=1,…,m.The above assumption enables us to predict that the layer inux will be at x=0 [16].Theorem 2.
Under above assumptions, the solution of (12) and (2) has the following asymptotic expansion:(13)ux=wx+vx¯+Oε,where wx¯ is the layer correction function, i.e.,(14)vx¯=vx,x∈0,x∗,0,x∈x∗,1,0<x∗=−lnεβ≪1.
For the convection-diffusion problems, there is a boundary layer region0,x∗ in the solution u of (12) and (2). That is, the boundary layers are located at the left endpoint, and its width is x∗=−lnε/βε.


## 2.1. Weakly Coupled System of Reaction-Diffusion Problems

IfA≠0,B≡0 in (1), it is convenient to introduce the notation to denote the reaction-diffusion operator:(3)LRu:=−εu″+Au=f.Assume thatA has positive diagonal entries and nonpositive off-diagonal entries, i.e.,(4)aij=>0,i=j,≤0,i≠j,and A is also strictly diagonally dominant i.e., ∑k=1k≠imaik<aii and ∑j=1maij>α2>0,i=1,…,m. These hypotheses ensure that problem (1) has a unique solution.Lemma 1.
Under above assumptions, all eigenvalues of matrixA are positive.Proof.
LetA have an eigenvalue of λ=0, according to the Gerschgorin disc theorem, ∃i0,ai0i0≤Ri0A, which is a contradiction. So, there is no zero eigenvalue in A, that is, A≠0.
Setfx=D+xB,x∈0,1, where D=diagA,B=A−D, then fx is a continuous function on 0,1 and f0=D=∏i=1maii>0; assume f1=A<0, then we can obtained ∃ξ∈0,1, such that fx=D+ξB=0. Moreover, A is strictly diagonally dominant, then aii>RiA≥∑j=1,j≠inξaij, that is, D+ξB is also strictly diagonally dominant. Contradictions between fξ=D+ξB≠0 and fξ=0 can be obtained. Therefore, f1=A>0, each of the main determinants of the strictly diagonally dominant matrix is the strictly diagonally dominant, then we have the A is positive definite matrix. So, all eigenvalues of matrix A are positive.Lemma 2.
Let us suppose that a functionux satisfies u0≥0,u1≥0,Lux≤0,∀x∈Ω, then ux≥0;∀x∈Ω¯.
The direct application of the maximum principle is the following stability result.Lemma 3.
Ifux is the solution of (2), then we have the stability bound inequality:(5)u≤M=maxb0,b1+Cmaxfx,u′x≤CMε.
Shishkin decomposition is used for asymptotic analysis of problem (3) [14], which splits the solution into regular and layer components:(6)ux=wx+vx,where the regular component wx and the layer component vx are the solutions of the two problems:(7)LRw=−εw″+Aw=f,w0=A−1f0,w1=A−1f1,(8)LRv=−εw″+Av=f,v0=b0−w0,v1=b1−w1,then wx and vx have the following estimates.Lemma 4.
(see [15]). Let u=w+v be the solution of problem (3) and wandv be the solutions of (7) and (8), respectively. Then, we have(9)wj≤C1+ε2−j/2,vj≤Cε−j/2e−αx/ε+e−α1−x/ε,where j=0,1,2.Theorem 1.
The solution of (3) and (2) has the following asymptotic expansion:(10)ux=wx+vx¯+Oε,where(11)vx¯=vx,x∈0,x∗∪1−x∗,1,0,x∈x∗,1−x∗,0<x∗=−lnεαε≪1.
The above theorem suggests that the boundary layer regions of the solutionux to the reaction-diffusion problems are 0,x∗ and 1−x∗,1, that is, the boundary layers are located at the two endpoints of the underlying interval 0,1 and each of their width is x∗=−lnε/αε.


## 2.2. Strongly Coupled System of Convection-Diffusion Problem

IfB≠0, the notation LC is introduced to denote the convection-diffusion operator:(12)LCu:=−εu″−Bu′+Au=f.To ensure that there is a unique solution to this problem, the following assumptions need to be satisfied:B≠0 is strictly diagonally dominant and hence invertible and define β=minbii,i=1,…,m.The above assumption enables us to predict that the layer inux will be at x=0 [16].Theorem 2.
Under above assumptions, the solution of (12) and (2) has the following asymptotic expansion:(13)ux=wx+vx¯+Oε,where wx¯ is the layer correction function, i.e.,(14)vx¯=vx,x∈0,x∗,0,x∈x∗,1,0<x∗=−lnεβ≪1.
For the convection-diffusion problems, there is a boundary layer region0,x∗ in the solution u of (12) and (2). That is, the boundary layers are located at the left endpoint, and its width is x∗=−lnε/βε.


## 3. The Singularity-Separated Technology

Note that ifε=10−k, then k≥3,e−1/ε=10−0.43×10k≤10−430≈0 can be neglected.


### 3.1. Weakly Coupled System of Reaction-Diffusion Problem

Theorem 3.
The general solution of the homogenous equationLRu=0 can be expressed in the form(15)u=∑j=1mCjeλj/εx−1−Cj′e−λj/εxVj,where λj is an eigenvalue of the A, the vector Vj is the associated eigenvector of A, and CjandCj′ are the constant vectors, j=1,2,…,m.Proof.
We can write higher order differential equations as a system with a very simple change of variable:(16)u′=v,−εv′+Au=0.
System (16) can be written in the following matrix form:(17)EmOmOm−εEmuv′=OmEm−AOmuv.
Thus,(18)uv′=OmEmAεOmuv,LetD=OmEmAεOm.
The eigenvalueλ of D can be found from the auxiliary equation:(19)λE2m−D=λEm−Em−AελEm=0.
Ifλ=0, we have λE2m−D=−1m2/εmA≠0. So, λ=0 is not the eigenvalue of D, and using the general transform, we can obtain(20)λEm−Em−AελEm=λEmOm−AελEm−Aλε=λEm−Aλε=0,therefore, we can get(21)λ2=λjε,j=1,2,…,m.where λj is an eigenvalue of A. The eigenvalue of matrix D is the square root of λj, so we can obtain all the eigenvalues of matrix D:(22)sj=λjε,sj′=−λjε,j=1,2,…,m.
Next, we need to show thatφjx=esjx−1Vjandψjx=e−sjxVj are the solutions of equation LRu=0:(23)LRφj=−εφj″x+Aφjx=−εsj2esjx−1Vj+Aesjx−1Vj=−λjVj+AVjesjx−1=0,that is, φjx is the solution of equation LRu=0. Similarly, it can be proved that ψjx is the solution of equation LRu=0.
Therefore, the general solution of the homogenous equationLRu=0 can be expressed in the form (15).Corollary 1.
The general solution of (3) can be expressed as follows:(24)u=w+∑j=1mCjeλj/εx−1−Cj′e−λj/εxVj,where w is a special solution, λj is an eigenvalue of A, the vector Vj is the associated eigenvector of A, and CjandCj′ are the constant vectors, j=1,2,…,m.
So, the solutionux is composed of two parts, u≈uss=w+v, in which wx is the regular term and vx is the singular term, and it has the following explicit expression:(25)v=∑j=1mCjϕjx+Cj′ψjx.
The parametersCjandCj′,j=1,2,…,m, in (25) can be defined by the boundary conditions. For any fx, how to construct the regular solution w is key.
We construct an auxiliary third boundary value problem to solve the special solution of (3):(26)Lw=−εw″+Aw=fx,Aw0=f0,Aw1=f1,then εw″0=0,εw″1=0, and the singularities of the solution are weakened. Based on the boundary condition (2), we can obtain the undetermined coefficients CjandCj′ in explicit singular function (15):(27)∑j=0mCj′Vj=b0−w0,∑j=0mCjVj=b1−w1.
Due to the linear independence ofVjj=1m, the matrix V=V1,V2,…,Vm is invertible. Equation (27) has a unique solution.


### 3.2. Strongly Coupled System of Convection-Diffusion Problem

Theorem 4.
IfA≡0, the general solution of (12) can be expressed as follows:(28)u=w+∑j=1mCje−λj/εxVj,where w is a special solution, λjandVj are the distinct eigenvalues and eigenvectors of B, and Cj is the constant which is determined by boundary condition (2).Proof.
Letu′=v, then (12) can be transformed into an equivalent problem including a system of two ODEs as follows:(29)u′=v,−εv′−Bv=0.
System (29) can be written in the following matrix form:(30)EmOmOm−εEmuv′=OmEmOmBuv,where EmandOm are identity matrix and all-zeros matrix, respectively. Then, we have(31)uv′=OmEmOm−Bεuv,setD=OmEmOm−Bε.
Letλ is an eigenvalue of D, and(32)λE2m−D=λEm−EmOmλEm+Bε=0.
The eigenvalueλ=s0=0 has multiplicity m, λ=sj=−λj/ε,j=1,2,…,m.
Next, we need to show thatφjx=e−λj/εVj is the solution of equation LCu=0:(33)LCφj=−εφj″x−Bφj′x=−λj2εe−λj/εxVj+Bλjεe−λj/εxVj=−λjVj+BVjλjεe−λj/εx=0.
So, the general solution of (12) can be expressed in the form (28).
Let the singular correction function be(34)vx=∑j=1mCje−λj/εxVj.
In the same way, the regular termw is a solution of an auxiliary third boundary value problem:(35)Lw=−εw″+Bw′=fx,Bw′0=f0,w1=b1.
According to the boundary conditions, we can obtain the undetermined coefficients in an explicit singular function:(36)∑j=1mCj′Vj=b0−w0.


## 3.1. Weakly Coupled System of Reaction-Diffusion Problem

Theorem 3.
The general solution of the homogenous equationLRu=0 can be expressed in the form(15)u=∑j=1mCjeλj/εx−1−Cj′e−λj/εxVj,where λj is an eigenvalue of the A, the vector Vj is the associated eigenvector of A, and CjandCj′ are the constant vectors, j=1,2,…,m.Proof.
We can write higher order differential equations as a system with a very simple change of variable:(16)u′=v,−εv′+Au=0.
System (16) can be written in the following matrix form:(17)EmOmOm−εEmuv′=OmEm−AOmuv.
Thus,(18)uv′=OmEmAεOmuv,LetD=OmEmAεOm.
The eigenvalueλ of D can be found from the auxiliary equation:(19)λE2m−D=λEm−Em−AελEm=0.
Ifλ=0, we have λE2m−D=−1m2/εmA≠0. So, λ=0 is not the eigenvalue of D, and using the general transform, we can obtain(20)λEm−Em−AελEm=λEmOm−AελEm−Aλε=λEm−Aλε=0,therefore, we can get(21)λ2=λjε,j=1,2,…,m.where λj is an eigenvalue of A. The eigenvalue of matrix D is the square root of λj, so we can obtain all the eigenvalues of matrix D:(22)sj=λjε,sj′=−λjε,j=1,2,…,m.
Next, we need to show thatφjx=esjx−1Vjandψjx=e−sjxVj are the solutions of equation LRu=0:(23)LRφj=−εφj″x+Aφjx=−εsj2esjx−1Vj+Aesjx−1Vj=−λjVj+AVjesjx−1=0,that is, φjx is the solution of equation LRu=0. Similarly, it can be proved that ψjx is the solution of equation LRu=0.
Therefore, the general solution of the homogenous equationLRu=0 can be expressed in the form (15).Corollary 1.
The general solution of (3) can be expressed as follows:(24)u=w+∑j=1mCjeλj/εx−1−Cj′e−λj/εxVj,where w is a special solution, λj is an eigenvalue of A, the vector Vj is the associated eigenvector of A, and CjandCj′ are the constant vectors, j=1,2,…,m.
So, the solutionux is composed of two parts, u≈uss=w+v, in which wx is the regular term and vx is the singular term, and it has the following explicit expression:(25)v=∑j=1mCjϕjx+Cj′ψjx.
The parametersCjandCj′,j=1,2,…,m, in (25) can be defined by the boundary conditions. For any fx, how to construct the regular solution w is key.
We construct an auxiliary third boundary value problem to solve the special solution of (3):(26)Lw=−εw″+Aw=fx,Aw0=f0,Aw1=f1,then εw″0=0,εw″1=0, and the singularities of the solution are weakened. Based on the boundary condition (2), we can obtain the undetermined coefficients CjandCj′ in explicit singular function (15):(27)∑j=0mCj′Vj=b0−w0,∑j=0mCjVj=b1−w1.
Due to the linear independence ofVjj=1m, the matrix V=V1,V2,…,Vm is invertible. Equation (27) has a unique solution.


## 3.2. Strongly Coupled System of Convection-Diffusion Problem

Theorem 4.
IfA≡0, the general solution of (12) can be expressed as follows:(28)u=w+∑j=1mCje−λj/εxVj,where w is a special solution, λjandVj are the distinct eigenvalues and eigenvectors of B, and Cj is the constant which is determined by boundary condition (2).Proof.
Letu′=v, then (12) can be transformed into an equivalent problem including a system of two ODEs as follows:(29)u′=v,−εv′−Bv=0.
System (29) can be written in the following matrix form:(30)EmOmOm−εEmuv′=OmEmOmBuv,where EmandOm are identity matrix and all-zeros matrix, respectively. Then, we have(31)uv′=OmEmOm−Bεuv,setD=OmEmOm−Bε.
Letλ is an eigenvalue of D, and(32)λE2m−D=λEm−EmOmλEm+Bε=0.
The eigenvalueλ=s0=0 has multiplicity m, λ=sj=−λj/ε,j=1,2,…,m.
Next, we need to show thatφjx=e−λj/εVj is the solution of equation LCu=0:(33)LCφj=−εφj″x−Bφj′x=−λj2εe−λj/εxVj+Bλjεe−λj/εxVj=−λjVj+BVjλjεe−λj/εx=0.
So, the general solution of (12) can be expressed in the form (28).
Let the singular correction function be(34)vx=∑j=1mCje−λj/εxVj.
In the same way, the regular termw is a solution of an auxiliary third boundary value problem:(35)Lw=−εw″+Bw′=fx,Bw′0=f0,w1=b1.
According to the boundary conditions, we can obtain the undetermined coefficients in an explicit singular function:(36)∑j=1mCj′Vj=b0−w0.


## 4. The Rational Spectral Collocation Method Combined with Singularity-Separated Method

### 4.1. The Barycentric Form of Rational Interpolation

Rational functionpNx in barycentric form, which interpolates function ux at N+1 distinct points xkk=0N can be expressed as [17](37)ux≈pNx=∑k=0Nωk/x−xkuxk∑k=0Nωk/x−xk,where ωkk=0N are barycentric weights. For Chebyshev–Gauss–Lobatto points xk=−coskπ/N, barycentric weights are chosen as [17](38)ω0=12,ωk=−1k,k=1,2,…,N−1,ωN=−1N2.According to the theorem in Baltensperger et al. [18], we have the convergence analysis of the rational interpolation polynomial in barycentric form with transformed Chebyshev points as follows:(39)pNx−ux=Oρ−N.∀x∈−1,1,where pNx is the rational interpolating polynomial in barycentric form of u, g is a conformal map, and ρ>1 is the sum of its major and minor axes of an ellipse.An advantage of the rational interpolation in barycentric form is that its derivatives can be calculated directly using differentiation formulae, instead of using the differential quotient rule repeatedly. Thenth order derivative of polynomial pNx evaluated at the point xj, and can be expressed in the form(40)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 by [17](41)Djk1=ωkωjxj−xk,j≠k,−∑i≠kDji1,j=k,Djk2=2Djk1Djj1−1xj−xk,j≠k.−∑i≠kDji2,j=k.As suggested in (39), the convergence rate of the rational spectral collocation method mainly depends on the analytic region of u in the complex plane. Therefore, the conformal map g could be chosen to enlarge the ellipse of analyticity of u∘g. Thus, compared with the Chebyshev spectral method, a better approximation of u could be obtained.Note that differentiation matrices (41) only rely on weights ωk and new points x^k, which is why the underlying equation does not need to be converted to new coordinates after maps.


### 4.2. The sinh Transform

In order to approximate the rapid changes in the boundary layer region. Tee and Trefethen have constructed the conformal map [19]:(42)gx=λ+μsinhsinh−11−λμ+sinh−11+λμx−12+sinh−11−λμ,where λandμ 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 boundary layer. For the convection-diffusion type, parameters in (42) should be chosen as follows:(43)λ=−1,μ=−2lnεβε.In order to better distinguish the singular perturbation problem with two boundary layers, Tee proposed the combined sinh transform as(44)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 so that the spectral accuracy could be preserved. For the reaction-difffusion type, the parameter in (44) should be chosen as μ=−2lnε/αε.


### 4.3. RSC-SSM

In this part we will show the RSC-SSM in a detailed algorithm.For reaction-diffusion, we solve auxiliary third boundary value problem (26) using a rational spectral collocation method with the sinh transformation. First, by introducing the transformation x=y+1/2 and defining u^y=ux=uy+1/2, then u′x=2u^′yandu″x=4u^″y. Evaluating equations in (26) at points yk,k=0,1,…,N, yields(45)−4εEm⊗D2+A⊗EmW=F,where(46)W=w^1Y⋮w^mY,w^iY=w^iy0⋮w^iyN,F=f^1Y⋮f^mY,f^iY=f^iy0⋮f^iyN,for i=1,…,m, ⊗ is Kronecker product.Boundary conditions in (26) suggest that(47)A⊗Em1W=f^iy0,A⊗EmN+1W=f^iyN,i=1,2,…,m.Solving the linear algebra system composed of (45) and (47), we can obtain the numerical solution of (26), which has a boundary layer which occurs at y=1 with the width μ=2x∗. Let the transformed Chebyshev collocation points be(48)Y=ykk=0N=g˜μ−coskπNk=0N,where g˜μ is expressed by (44) in which μ=2x∗.As for the convection-diffusion type, evaluating equations in (35) at points yk,k=0,1,…,N, yields(49)−4εEm⊗D2−2B⊗D1W=F.Boundary conditions in (35) suggest that(50)−2B⊗D11W=f^iy0,w^iyN=b1i,i=1,2,…,m.Solving the linear algebra system composed of (49) and (50), we can obtain the numerical solution of (35).The numerical solution of original problem (1) and (2) can be obtained as uN=W+V.


### 4.4. Errors Analysis

Theorem 5.
LetuxianduNxi be the exact solution and numerical solution of the original problem (3) and (2), respectively. Then, we have(51)u−uN≤c1e−ϱA/ε+c2ρ−N,where ϱA is the spectral radius of A.Proof.
Letz=u−uN, then the error z satisfies(52)z=u−uN=u−us+us−uN≤u−us+w+v−uN.
According to the proof of Theorem 5, the singular correct functionvx is the solution of the homogeneous equations Lv=0. So, set R=u−uSS:(53)Lu−uS=Lu−w−v=−εu″+Au−−εw″+Au−Lv=0.
Simultaneously, obtained from (27)(54)R0=u0−w0−∑j=1mCj′Vj+∑j=1mCjVje−λj=∑j=1mCjVje−λj,R1=u1−w1−∑j=1mCjVj+∑j=1mCj′Vje−λj=∑j=1mCj′Vje−λj.
According to the Lemma 3, we can get(55)R=u−us≤∑j=1mCj′Vje−λj+∑j=1mCjVje−λj≤c1e−ϱA/ε.
And, using the convergence of rational spectral collocation (39), we obtain(56)w+v−uN≤c2ρ−N.
Joint (55) and (56), we have the error estimation(57)u−uN≤u−us+w+v−uN≤c1e−ϱ/ε+c2ρ−N.
Similarly, for the convection-diffusion type, we have(58)u−uN≤u−us+w+v−uN≤c1e−ϱB/ε+c2ρ−N,where ϱB is the spectral radius of B.
So, I can get the convergence of the RSC-SSM that can almost achieve spectral accuracy.


## 4.1. The Barycentric Form of Rational Interpolation

Rational functionpNx in barycentric form, which interpolates function ux at N+1 distinct points xkk=0N can be expressed as [17](37)ux≈pNx=∑k=0Nωk/x−xkuxk∑k=0Nωk/x−xk,where ωkk=0N are barycentric weights. For Chebyshev–Gauss–Lobatto points xk=−coskπ/N, barycentric weights are chosen as [17](38)ω0=12,ωk=−1k,k=1,2,…,N−1,ωN=−1N2.According to the theorem in Baltensperger et al. [18], we have the convergence analysis of the rational interpolation polynomial in barycentric form with transformed Chebyshev points as follows:(39)pNx−ux=Oρ−N.∀x∈−1,1,where pNx is the rational interpolating polynomial in barycentric form of u, g is a conformal map, and ρ>1 is the sum of its major and minor axes of an ellipse.An advantage of the rational interpolation in barycentric form is that its derivatives can be calculated directly using differentiation formulae, instead of using the differential quotient rule repeatedly. Thenth order derivative of polynomial pNx evaluated at the point xj, and can be expressed in the form(40)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 by [17](41)Djk1=ωkωjxj−xk,j≠k,−∑i≠kDji1,j=k,Djk2=2Djk1Djj1−1xj−xk,j≠k.−∑i≠kDji2,j=k.As suggested in (39), the convergence rate of the rational spectral collocation method mainly depends on the analytic region of u in the complex plane. Therefore, the conformal map g could be chosen to enlarge the ellipse of analyticity of u∘g. Thus, compared with the Chebyshev spectral method, a better approximation of u could be obtained.Note that differentiation matrices (41) only rely on weights ωk and new points x^k, which is why the underlying equation does not need to be converted to new coordinates after maps.


## 4.2. The sinh Transform

In order to approximate the rapid changes in the boundary layer region. Tee and Trefethen have constructed the conformal map [19]:(42)gx=λ+μsinhsinh−11−λμ+sinh−11+λμx−12+sinh−11−λμ,where λandμ 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 boundary layer. For the convection-diffusion type, parameters in (42) should be chosen as follows:(43)λ=−1,μ=−2lnεβε.In order to better distinguish the singular perturbation problem with two boundary layers, Tee proposed the combined sinh transform as(44)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 so that the spectral accuracy could be preserved. For the reaction-difffusion type, the parameter in (44) should be chosen as μ=−2lnε/αε.


## 4.3. RSC-SSM

In this part we will show the RSC-SSM in a detailed algorithm.For reaction-diffusion, we solve auxiliary third boundary value problem (26) using a rational spectral collocation method with the sinh transformation. First, by introducing the transformation x=y+1/2 and defining u^y=ux=uy+1/2, then u′x=2u^′yandu″x=4u^″y. Evaluating equations in (26) at points yk,k=0,1,…,N, yields(45)−4εEm⊗D2+A⊗EmW=F,where(46)W=w^1Y⋮w^mY,w^iY=w^iy0⋮w^iyN,F=f^1Y⋮f^mY,f^iY=f^iy0⋮f^iyN,for i=1,…,m, ⊗ is Kronecker product.Boundary conditions in (26) suggest that(47)A⊗Em1W=f^iy0,A⊗EmN+1W=f^iyN,i=1,2,…,m.Solving the linear algebra system composed of (45) and (47), we can obtain the numerical solution of (26), which has a boundary layer which occurs at y=1 with the width μ=2x∗. Let the transformed Chebyshev collocation points be(48)Y=ykk=0N=g˜μ−coskπNk=0N,where g˜μ is expressed by (44) in which μ=2x∗.As for the convection-diffusion type, evaluating equations in (35) at points yk,k=0,1,…,N, yields(49)−4εEm⊗D2−2B⊗D1W=F.Boundary conditions in (35) suggest that(50)−2B⊗D11W=f^iy0,w^iyN=b1i,i=1,2,…,m.Solving the linear algebra system composed of (49) and (50), we can obtain the numerical solution of (35).The numerical solution of original problem (1) and (2) can be obtained as uN=W+V.


## 4.4. Errors Analysis

Theorem 5.
LetuxianduNxi be the exact solution and numerical solution of the original problem (3) and (2), respectively. Then, we have(51)u−uN≤c1e−ϱA/ε+c2ρ−N,where ϱA is the spectral radius of A.Proof.
Letz=u−uN, then the error z satisfies(52)z=u−uN=u−us+us−uN≤u−us+w+v−uN.
According to the proof of Theorem 5, the singular correct functionvx is the solution of the homogeneous equations Lv=0. So, set R=u−uSS:(53)Lu−uS=Lu−w−v=−εu″+Au−−εw″+Au−Lv=0.
Simultaneously, obtained from (27)(54)R0=u0−w0−∑j=1mCj′Vj+∑j=1mCjVje−λj=∑j=1mCjVje−λj,R1=u1−w1−∑j=1mCjVj+∑j=1mCj′Vje−λj=∑j=1mCj′Vje−λj.
According to the Lemma 3, we can get(55)R=u−us≤∑j=1mCj′Vje−λj+∑j=1mCjVje−λj≤c1e−ϱA/ε.
And, using the convergence of rational spectral collocation (39), we obtain(56)w+v−uN≤c2ρ−N.
Joint (55) and (56), we have the error estimation(57)u−uN≤u−us+w+v−uN≤c1e−ϱ/ε+c2ρ−N.
Similarly, for the convection-diffusion type, we have(58)u−uN≤u−us+w+v−uN≤c1e−ϱB/ε+c2ρ−N,where ϱB is the spectral radius of B.
So, I can get the convergence of the RSC-SSM that can almost achieve spectral accuracy.


## 5. Numerical Experiments

In order to demonstrate the accuracy and efficiency of the RSC-SSM, we use RSCAT and RSC-SSM to solve singularly perturbed problems with exact solutions and compare the results. We verify the theoretical results obtained in the previous section through numerical experiments. The relative maximum errors of the solution are given by(59)e=uN−uE∞uE∞,where uNanduE are the numerical solution and the exact solution, respectively.In our computations, all experiments are carried outperformed using MATLAB (version R2014a) on a personal computer with a 2.5 GHz central processing unit (Intel Core i5-2450M), 4.00 GB of memory, and Windows 7 operating system.Example 1.
Consider the reaction-diffusion problem with constant coefficients:(60)B=0,A=2−1−12,f=12,b0=B1=00.
The exact solution of this problem is(61)u=121,1Tu1+12−1,1Tu2,where(62)u1=3ex−2/ε−e−x/ε−ex−1/ε+e−x+1/ε1−e−2/ε+3,u2=13e3x−2/ε−e−3x/ε−e3x−1/ε+e−3x+1/ε1−e−23/ε+13.
This problem has two boundary layer regions, i.e., the boundary layer are located at the two endpoints of the underlying interval0,1. In our method, we use the transform (44), and the parameter is chosen as α=0.99. The eigenvalues of A/ε are λ1=1/εandλ2=3/ε, and V1=1,1TandV2=1,−1T are the corresponding eigenvectors. So, the general solution of this problem can be expressed (24), and the special solution w is a solution of an auxiliary third boundary value problem (26). The undetermined coefficients of singular term can be obtained by (27).
Figure1 plots the relative maximum errors of both RSCAT and RSC-SSM in semilog scale for various values of ε. Compared with RSCAT, the convergence rates have considerably improved in our RSC-SSM on over trend. Furthermore, Table 1 provides the comparison of numerical results obtained by our RSC-SSM and RSCAT with several choices of ε. We perceive that our method gives much better results with significantly fewer number of unknowns. With the increase of N, the relative maximum error of RSCAT method decreases, while quite the opposite is true for the RSC-SSM. Figure 2 displays plots of numerical and exact solutions with ε=10−2andε=10−6. There are more points located in the boundary layer region. Figure 3 shows the point-wise errors of the function in the whole region, it is observed that compared with RSCAT, the smaller the parameter ε is, the more accurate the RSC-SSM result.Figure 1
Relative maximum errors in Example 1. (a)ε = 10−1. (b) ε = 10−4. (c) ε = 10−6. (d) ε = 10−8.
(a)
(b)
(c)
(d)Table 1
Comparison of the relative maximum errors for Example 1.
ε
N
RSC-SSM
RSCAT
u
v
u
v
1e − 2
20
5.19e − 05
4.14e − 05
1.61e − 03
1.29e − 03
1e − 4
40
2.33e − 12
4.10e − 12
6.74e − 06
5.70e − 06
1e − 6
60
7.49e − 14
2.04e − 13
1.05e − 07
1.06e − 07
1e − 8
60
1.50e − 15
8.66e − 16
5.63e − 06
4.77e − 06Figure 2
Numerical solutions vs. exact solutions(u1: solid lines and u2: dashed lines). (a) ε = 10−2 and N = 40. (b) ε = 10−6 and N = 60.
(a)
(b)Figure 3
Pointwise relative errors in Example 1: RSCAT (dashed lines) vs. RSC-SSM (solid lines). (a)ε = 10−2 and N = 20. (b) ε = 10−6 and N = 60.
(a)
(b)Example 2.
Consider the convection-diffusion problem with constant coefficients:(63)B=3−1−1−14−2−1−24,f=−411−7,A=0,b0=−14−1,b1=e−1/ε−2e−4/ε+1e−1/ε+e−4/ε+2e−6/ε−2e−1/ε+e−4/ε−2e−6/ε.
The exact solution of this problem is(64)u=1,1,1Te−x/ε+−2,1,1Te−4x/ε+0,2,−2Te−6x/ε+x,−2x,x−1T.
For the sake of comparison, the parameters in transform (42) are chosen as β=2.55. This problem has rapid changes near the point x=0. The eigenvalues of B/ε are λ1=1/ε,λ2=4/ε,andλ3=6/ε, and the corresponding eigenvectors are V1=1,1,1T,V2=2,−1,−1T,andV3=0,−1,1T. So, the general solution of this problem can be expressed (28), and the special solution w is a solution of auxiliary third boundary value problem (35). The undetermined coefficients of singularity term can be obtained by (36). We compare the maximum relative errors of the RSC-SSM with the RSCAT method in Figure 4 and Table 2, and both verify the high accuracy and efficiency of our method.
Figure4 plots the relative maximum errors of both RSCAT and RSC-SSM in semilog scale for various values of ε. It shows that the convergence rates have greatly improved in our RSC-SSM. Besides, Table 2 illustrates the comparison of numerical results obtained by our RSC-SSM and the method in [12] with several choices of ε. We observe that our method gives much better results with significantly less number of unknowns. Figure 5 displays the plots of numerical and exact solutions for the case with ε=10−2andε=10−6. There are more points located in the boundary layer region. Figure 6 shows the pointwise errors of the function in the whole region; it is observed that compared with RSCAT, we can conclude that with the increase of N, the error of the RSC-SSM remains stable.Figure 4
Relative maximum errors in Example 2: RSCAT (dashed lines) vs. RSC-SSM (solid lines). (a)ε = 10−2. (b) ε = 10−6.
(a)
(b)Table 2
Comparison of the relative maximum errors for Example 2.
ε
N
RSC-SSM
RSCAT
u
v
w
u
v
w
1e − 2
20
3.39e − 15
1.55e − 15
3.11e − 15
1.21e − 03
4.70e − 03
1.76e − 02
1e − 4
40
3.12e − 14
9.77e − 15
3.52e − 14
1.53e − 05
2.58e − 04
1.01e − 03
1e − 6
60
4.33e − 13
1.31e − 13
2.67e − 13
1.38e − 06
2.05e − 05
8.35e − 05
1e − 8
60
6.87e − 12
1.23e − 12
3.18e − 12
3.34e − 04
9.71e − 03
3.81e − 02Figure 5
Numerical solutions vs. exact solutions. RSCAT (dashes lines) vs. RSC-SSM (solid lines). (a)ε = 10−2 and N = 20. (b) ε = 10−6 and N = 60.
(a)
(b)Figure 6
Pointwise errors. RSCAT (dashes lines) vs. RSC-SSM method (solid lines).
(a)
(b)


## 6. Concluding Remarks

In this paper, a new numerical method RSC-SSM has been proposed to solve a system of singularly perturbed boundary value problems. The numerical solution consists of two parts: the solution of weaker singularly auxiliary equation and the singular correction term. Numerical experiments confirm that compared to RSCAT, the present RSC-SSM has the following advantages:(1)
The RSC-SSM is easy to implement, and it enjoys computational efficiency, accuracy, and stability over some popular numerical approaches.(2)
The accuracy of numerical approximation depends not only on the number of the grid nodes, but also on parameterε. The smaller the ε is, the thinner the boundary layer is, and the better results obtained.(3)
In the RSC-SSM, we make good use of the eigenvalue of coefficient matrix, which plays a important role in constructing general solution of singularly perturbed problems.The numerical results demonstrate the spectral accuracy of proposed algorithms and agree with the theoretical analysis very well. And, the theoretical and numerical frameworks presented in this paper are essential for extension to more complicated problems.



---

*Source: 9030565-2019-11-14.xml*