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

A novel rational spectral collocation method is presented combined with the singularity-separated technique for a system of singularly perturbed boundary value problems. *e solution is expressed as u � w + v, where w is the solution of the corresponding auxiliary boundary value problem and v is a singular correction with explicit expressions. *e rational spectral collocation method in barycentric form with the sinh transformation is applied to solve the auxiliary third boundary problem.*e 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.


Introduction
e 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. ere is a vast literature dealing with SPPs with smooth coefficients and source term for single equation [1][2][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]. e 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 of m ≥ 2 singularly perturbed linear equations are considered in the unknown vector function u � (u 1 , . . . , u m ) T : and it yields u(0) � b 0 , where A � (a ij ) and B � (b ij ) are m × m matrices and 0 < ε ≪ 1 is a small parameter whose presence makes the problem singularly perturbed. We assumed that f � (f 1 , . . . , f m ) T ∈ C 2 [0, 1] m and both b 0 � (b 01 , . . . , b 0m ) T and b 1 � (b 11 , . . . , b 1m ) T are constant vectors. ese 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 � 0(or x � 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. e 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. e 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. is paper is organized as follows. e asymptotic analysis of coupled system is outlined in Section 2. e 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.
e 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.

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).

Weakly Coupled System of Reaction-Diffusion Problems.
If A ≠ 0, B ≡ 0 in (1), it is convenient to introduce the notation to denote the reaction-diffusion operator: Assume that A has positive diagonal entries and nonpositive off-diagonal entries, i.e., and A is also strictly diagonally dominant i.e., m k�1 k≠i |a ik | < a ii and m j�1 a ij > α 2 > 0, i � 1, . . . , m. ese hypotheses ensure that problem (1) has a unique solution. Proof. Let A have an eigenvalue of λ � 0, according to the Gerschgorin disc theorem, ∃i 0 ,  The direct application of the maximum principle is the following stability result.

Lemma 3. If u(x)
is the solution of (2), then we have the stability bound inequality: Shishkin decomposition is used for asymptotic analysis of problem (3) [14], which splits the solution into regular and layer components: where the regular component w(x) and the layer component v(x) are the solutions of the two problems: then w(x) and v(x) have the following estimates.

Mathematical Problems in Engineering
Lemma 4 (see [15]). Let u � w + v be the solution of problem (3) and w and v be the solutions of (7) and (8), respectively. en, we have where j � 0, 1, 2.

Theorem 1.
e solution of (3) and (2) has the following asymptotic expansion: where The above theorem suggests that the boundary layer regions of the solution u(x) 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

Strongly Coupled System of Convection-Diffusion
Problem. If B ≠ 0, the notation L C is introduced to denote the convection-diffusion operator: 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 e above assumption enables us to predict that the layer in u(x) will be at x � 0 [16].

Theorem 2.
Under above assumptions, the solution of (12) and (2) has the following asymptotic expansion: For the convection-diffusion problems, there is a boundary layer region [0, x * ] in the solution u of (12) and (2).
at is, the boundary layers are located at the left endpoint, and its width is

Weakly Coupled System of Reaction-Diffusion Problem
e general solution of the homogenous equation L R u � 0 can be expressed in the form where λ j is an eigenvalue of the A, the vector V j is the associated eigenvector of A, and C j and C j ′ 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: System (16) can be written in the following matrix form: us, e eigenvalue λ of D can be found from the auxiliary equation: is not the eigenvalue of D, and using the general transform, we can obtain therefore, we can get Mathematical Problems in Engineering 3 where λ j is an eigenvalue of A. e eigenvalue of matrix D is the square root of λ j , so we can obtain all the eigenvalues of matrix D: erefore, the general solution of the homogenous equation L R u � 0 can be expressed in the form (15).

Corollary 1.
e general solution of (3) can be expressed as follows: where w is a special solution, λ j is an eigenvalue of A, the vector V j is the associated eigenvector of A, and C j and C j ′ are the constant vectors, j � 1, 2, . . . , m. So, the solution u(x) is composed of two parts, u ≈ u ss � w + v, in which w(x) is the regular term and v(x) is the singular term, and it has the following explicit expression: The parameters C j and C j ′ , j � 1, 2, . . . , m, in (25) can be defined by the boundary conditions. For any f(x), how to construct the regular solution w is key.
We construct an auxiliary third boundary value problem to solve the special solution of (3): 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 C j and C j ′ in explicit singular function (15):

Strongly Coupled System of Convection-Diffusion Problem
Theorem 4. If A ≡ 0, the general solution of (12) can be expressed as follows: where w is a special solution, λ j and V j are the distinct eigenvalues and eigenvectors of B, and C j is the constant which is determined by boundary condition (2). (12) can be transformed into an equivalent problem including a system of two ODEs as follows: System (29) can be written in the following matrix form: where E m and O m are identity matrix and all-zeros matrix, respectively. en, we have Let λ is an eigenvalue of D, and e eigenvalue λ � s 0 � 0 has multiplicity m, λ � s j � − (λ j /ε), (j � 1, 2, . . . , m).
Next, we need to show that φ j (x) � e − (λ j /ε) V j is the solution of equation L C u � 0: So, the general solution of (12) can be expressed in the form (28).
Let the singular correction function be In the same way, the regular term w is a solution of an auxiliary third boundary value problem: According to the boundary conditions, we can obtain the undetermined coefficients in an explicit singular function: can be expressed as [17] u where ω k N k�0 are barycentric weights. For Chebyshev-Gauss-Lobatto points x k � − cos(kπ/N), barycentric weights are chosen as [17] ω 0 � 1 2 , 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: where p N (x) 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.
e nth order derivative of polynomial p N (x) evaluated at the point x j , and can be expressed in the form where D (n) jk is the entry of the nth order differentiation matrix.
e entries of the first-and second-order differentiation matrices, are given by [17] 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. erefore, the conformal map g could be chosen to enlarge the ellipse of analyticity of u ∘ g. us, 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.

e sinh Transform.
In order to approximate the rapid changes in the boundary layer region. Tee and Trefethen have constructed the conformal map [19]: where λ and μ are the location and width of the boundary layers, respectively. e transformed Chebyshev points g λ,μ (x k ) N k�0 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:

Mathematical Problems in Engineering
In order to better distinguish the singular perturbation problem with two boundary layers, Tee proposed the combined sinh transform as All derivatives of the piecewise map g 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 μ � − 2(ln

Theorem 5. Let u(x i ) and u N (x i ) be the exact solution and numerical solution of the original problem
where 9 A is the spectral radius of A.
Proof. Let z � u − u N , then the error z satisfies According to the proof of eorem 5, the singular correct function v(x) is the solution of the homogeneous equations Lv � 0. So, set R � u − u SS : Simultaneously, obtained from (27) According to the Lemma 3, we can get And, using the convergence of rational spectral collocation (39), we obtain 6 Mathematical Problems in Engineering Joint (55) and (56), we have the error estimation Similarly, for the convection-diffusion type, we have where 9 B is the spectral radius of B. So, I can get the convergence of the RSC-SSM that can almost achieve spectral accuracy.

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. e relative maximum errors of the solution are given by where u N and u E 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: e exact solution of this problem is where is problem has two boundary layer regions, i.e., the boundary layer are located at the two endpoints of the underlying interval [0, 1]. In our method, we use the transform (44), and the parameter is chosen as α � 0.99. e eigenvalues of A/ε are λ 1 � 1/ε and λ 2 � 3/ε, and V 1 � (1, 1) T and V 2 � (1, − 1) T 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). e undetermined coefficients of singular term can be obtained by (27). Figure 1 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 − 2 and ε � 10 − 6 .
ere 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.
Example 2. Consider the convection-diffusion problem with constant coefficients:   Figure 4 and Table 2, and both verify the high accuracy and efficiency of our method. Figure 4 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 − 2 and ε � 10 − 6 . ere 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.

Concluding Remarks
In this paper, a new numerical method RSC-SSM has been proposed to solve a system of singularly perturbed boundary value problems. e 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) e RSC-SSM is easy to implement, and it enjoys computational efficiency, accuracy, and stability over some popular numerical approaches. (2) e accuracy of numerical approximation depends not only on the number of the grid nodes, but also on parameter ε. e 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.
e 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.

Data Availability
e data used to support the findings of this study are included within the article Conflicts of Interest e author declares that there are conflicts of interest. Mathematical Problems in Engineering 11