Absolutely stable difference scheme for a general class of singular perturbation problems

This paper presents an absolutely stable noniterative difference scheme for solving a general class of singular perturbation problems having left, right, internal, or twin boundary layers. The original two-point second-order singular perturbation problem is approximated by a first-order delay differential equation with a variable deviating argument. This delay differential equation is transformed into a three-term difference equation that can be solved using the Thomas algorithm. The uniqueness and stability analysis are discussed, showing that the method is absolutely stable. An optimal estimate for the deviating argument is obtained to take advantage of the second-order accuracy of the central finite difference method in addition to the absolute stability property. Several problems having left, right, interior, or twin boundary layers are considered to validate and illustrate the method. The numerical results confirm that the deviating argument can stabilize the unstable discretized differential equation and that the new approach is effective in solving the considered class of singular perturbation problems.


Introduction
Singular perturbation problems (SPPs), also known as stiff, arise very frequently in various areas of applied science, such as fluid mechanics, heat transfer, materials science, chemical and electrical engineering, superconductor theory, and chemical reactor systems [1][2][3][4][5][6][7][8][9]. It is a well-known fact that SPPs exhibit stiff internal or boundary layers as the singular perturbation parameter approaches zero . As this happens, the stiffness ratio increases and the classical numerical methods become inadequate for solving such problems over standard locally uniform meshes, where nonphysical oscillations emerge in the numerical solution due to the stability restriction on the step-size. An interesting collection of examples of exotic numerical solutions of singular perturbation problems obtained using inadequate numerical discretization schemes can be found in [10]. The construction of special-purpose numerical methods for ensuring convergence regardless of the value of the perturbation parameter was first raised in [11,12] using fitting factor (stabilization term) and fitting mesh (layer adapted mesh) techniques. It is well known that depending on the behavior of the differential operator coefficients in SPPs, the solutions exhibit a boundary layer, twin boundary-layers, or interior layer [3,4,[13][14][15][16][17][18]. A variety of numerical and semianalytical numerical methods for solving different classes of SPPs exhibiting a boundary layer, twin boundary-layers, or interior layer can be found in . In all of these SPP classes, a priori information about the location of the layers is used to construct appropriate numerical methods for SPPs. And consequently, the numerical methods developed for SPPs have different numerical integration algorithms depending on the location of the layers. The purpose of this paper is to present a noniterative absolutely stable difference scheme for solving a general class of SPPs having left, right, internal, or twin boundary layers. The method does not depend on asymptotic expansions and needs no prior information about the location of the layers, and so it is designed to be suitable for the practicing engineers and applied mathematicians who need a practical tool for solving SPPs (see, for example, [15,38] for asymptotic techniques, designed for a similar purpose). The method depends on replacing the original two-point second-order SPP by a first-order approximate delay differential equation with a variable deviating argument. This delay differential equation is transformed into a three-term difference equation that can be solved using the Thomas algorithm [50]. The stability analysis of the method is discussed showing that the method is absolutely stable under a certain condition on the deviating argument whereas there is no stability restriction on the step-size. An optimal estimate for the deviating argument is obtained to take advantage of the second-order accuracy of the central finite difference method [48,50] in addition to the absolute stability property [19,28,[48][49][50]. Several problems having left, right, interior, or twin boundary layers are considered to validate and illustrate the method. To analyze the effect of the deviating argument on the solution accuracy, the maximum absolute solution error is presented in tables and figures for constant and variable deviating argument values. The numerical results confirm that the deviating argument can stabilize the unstable discretized differential equation and that the new approach is effective in solving the considered class of singular perturbation problems.

Description of the method
Consider the following linear SPP: with BCs where 0 < ε 1, a, b, A and B are given constants, and p(x), q(x) and r( Moreover, we assume that q(x) ≥ 0 for all x ∈ [a, b]. Under these assumptions, problem (1) has a unique solution with boundary or interior layers . The interval [a, b] is divided into N subintervals [x i , x i+1 ], i = 0, 1, 2, . . . , N -1, each of length h, i.e., h = (ba)/N and x i = a + ih. For the sake of simplicity, we use (1) by -ε and letting ω = 1/ε results in Let δ i = δ(x i ) be a positive variable deviating argument (0 < δ i 1). Then by using the Taylor expansion about the point x i , we have By substituting Eq. (3) into Eq. (2), we have Now, applying standard forward and centered finite difference formulas for y i in Eq. (4) results in Again, using the Taylor expansion and the standard backward finite difference, we obtain: By substituting Eq. (6) into Eq. (5), we get the following difference equation: where The difference Eq. (7) and the two BCs in (1) results in a tridiagonal system that can be easily solved using the Thomas algorithm to obtain the unknowns y 1 to y N-1 .

Thomas algorithm
A brief description of the Thomas algorithm [50], is presented as follows.
In this algorithm, the solution of the difference Eq. (7) can be written as where W i and T i are to be determined. From Eq. (12) we have and substituting (12) and (13) into (7), we get Comparing (14) and (12) results in Starting with initial conditions W 0 = 0 and T 0 = y 0 = A, the values of W i and T i for i = 1, 2, . . . , N -1 are computed in a forward process, from Eqs. (15)- (16), and then y i is computed in a backward process, from Eq. (12), using the BC y N = B.

Uniqueness and stability analysis
In this section, the uniqueness and stability condition of the proposed algorithm for the difference Eq. (7) is analyzed.

Theorem 1
The numerical scheme E i y i-1 + F i y i + G i y i+1 ∼ = H i defined in (7)-(11) results in a unique stable solution provided Proof The existence of a unique solution to the tridiagonal system results from establishing that the tridiagonal coefficient matrix of the system is diagonally dominant. It is clear that E i in (8) and G i in (10) are nonnegative when 2 hδ i + ωp i 2h ≥ 0 and 2 hδ i -ωp i 2h ≥ 0, respectively. These inequalities are verified if the condition δ i ≤ 4 ω|p i | is verified. Following this condition, F i is negative and verifies the inequality |F i | ≥ |G i | + |E i | = 4 hδ i . Thus, under condition (17), the tridiagonal coefficient matrix is diagonally dominant, and the numerical scheme has a unique solution.
For stability analysis, suppose that a small numerical error τ i-1 was introduced in the calculation of W i-1 so that it yields an approximate valueW i-1 such that Then, from (18) and (15), we have Thus, the numerical scheme defined by (7)-(11) results in a unique stable solution under condition (17).
Now, if we chose δ i = 4 ω|p i | , then for the left-end boundary layers, p i < 0, we have E i = 0, and the present scheme is reduced to a stable forward integration scheme. Also, on the other hand, for the right-end boundary layers, p i > 0, we have G i = 0, and the present scheme is reduced to a stable backward integration scheme. Moreover, in the above two cases, the propagation error in (20) has vanished.
Moreover, at p i = 0, for twin-boundary layers or internal layers, or at h ≤ 2 ω p i ∞ , we set δ i = 2h ≤ 4 ω|p i | , and the present scheme is reduced to the standard central finite difference scheme (CFD) in [48,50]. Thus, to take advantage of the second-order accuracy of the standard central finite difference scheme in addition to the absolute stability property, δ i can be selected according to The stability in CFD is controlled by the step-size restriction whereas the present scheme is absolutely stable without any restrictions on the step-size. These details are combined in the following algorithm:

Algorithm steps
Step Step III: Do Steps (i) and (ii) for index2 = 1, 2, 3, . . . , N -1 Step IV: Compute the solution This algorithm is easily adaptable in any mathematical environment. We present the Matlab code in Appendix A.

Numerical results
To demonstrate the applicability of the method, we consider in the following several problems having left, right, interior, or twin boundary layers. These problems have been discussed in the literature and their approximate solutions are available for comparison.

Figure 2
The solution of Example 2 for ε = 10 -4 using the ASCD and CFD schemes Table 3 The solution errors of Example 3 for ε = 10 -3 and ε = 10 -4 at h = 0.001 The results in Tables 1, 2 and 3 show that at δ i = 4/ω|p i | or δ i = δ op = 4/ω|p i |, there is no propagation error with the Thomas algorithm (20). Also, at ε = 10 -3 the ASCD is reduced to the CFD scheme, where δ op = 2h, and the results are identical. Although the ASCD stability is verified at δ i = 4/ω|p i |, selecting δ i = δ op takes advantage of the second-order accuracy of the CFD in addition to the absolute stability property depending on the stepsize value. The solution obtained at ε = 10 -4 using the ASCD is more accurate than that obtained by CFD due to the absolute stability property of the ASCD, while the stability restriction of CFD results in nonphysical oscillations in the numerical solution as illustrated in Figs. 1, 2 and 3. Figure 4 shows the solution error of Example 3 at ε = 10 -4 using different values of the deviating argument δ i = k/ω|p i |, k = 4, 3, 2, 1. The results in Fig. 4 reveal that a more accurate solution is obtained at δ i = 4/ω|p i |. Moreover, Fig. 5 depicts the solution error of Example 3 at ε = 10 -4 using constant and variable deviating arguments δ i = 4/ω|p i | ∞ , δ i = 4/ω|p i |, respectively. Figure 5 reveals that a more accurate solution is obtained using a variable deviating argument δ i .
The results in Tables 4 and 5 and Figs. 6 and 7 show that the ASCD can handle problems with right-end boundary layers effectively as well as those with left-end boundary layers.

Internal or twin-boundary-layer problems
Example 6 Consider the following SPBVP [49] given by Figure 6 The solution of Example 4 for ε = 10 -4 using the ASCD and CFD schemes Table 5 Solution errors of Example 5 for ε = 10 -3 and ε = 10 -4 at h = 0.001 with BCs y(-1) = 0 and y(1) = 2. The exact solution is given by The solution errors are listed in Table 6 for ε = 10 -3 and ε = 10 -4 at h = 0.001. The numerical solutions using the ASCD and CFD schemes for ε = 10 -4 are presented in Fig. 8.

Figure 7
The solution of Example 5 for ε = 10 -4 using the ASCD and CFD schemes Table 6 The solution errors of Example 6 for ε = 10 -3 and ε = 10 -4 at h = 0.001 The solution errors are listed in Table 7 for ε = 10 -3 and ε = 10 -4 at h = 0.001. The numerical solutions using the ASCD and CFD methods for ε = 10 -4 are presented in Fig. 9.
The results in Table 6 and Fig. 8 show that the ASCD and CFD lead to stable and accurate results for the internal layer problem (26). The reason for absenting the stability restriction of CFD in solving internal layer problems is that the internal layer occurs at a turning point x = x in at which p(x in ) = 0. The results in Table 7 show that the solution obtained for the twin boundary layer problem (28) at ε = 10 -4 using the ASCD is more accurate than those of the CFD and thus is due to the absolute stability property of the ASCD, while for the CFD the stability restriction results in nonphysical oscillations in the numerical solution as shown in Fig. 9. Figure 10 depicts the solution error of Example 7 for ε = 10 -4 using different values of the deviating argument δ i = k/ω|p i |, k = 4, 3, 2, 1. The results in Fig. 10 show that a more accurate solution is obtained at δ i = 4/ω|p i |. Figure 11 illustrates the solution error of Exam- Figure 8 The solution of Example 6 for ε = 10 -4 using the ASCD and CFD schemes Table 7 The solution errors of Example 7 for ε = 10 -3 and ε = 10 -4 at h = 0.001 ple 7 for ε = 10 -4 when adopting constant and variable deviating arguments δ i = 4/ω|p i | ∞ , δ i = 4/ω|p i |, respectively. Figure 11 reveals that a more accurate solution is obtained using a variable deviating argument.

Conclusions
This paper presented an absolutely stable noniterative central difference scheme for solving a general class of SPPs having left, right, internal, or twin boundary layers. The original two-point second-order SPPs is approximated by a first-order delay differential equation with a variable deviating argument. This delay differential equation was transformed into a three-term difference equation that was solved using the Thomas algorithm. The stability analysis of the method is discussed showing that the method is absolutely stable under Figure 9 The solution of Example 7 for ε = 10 -4 using the ASCD and CFD schemes Figure 10 The solution errors of Example 7 at ε = 10 -4 and different values of δ i Figure 11 Solution errors of Example 7 at ε = 10 -4 using constant and variable deviating argument δ i a certain condition on the deviating argument without stability restriction on the stepsize. An optimal estimation for the deviating argument is obtained taking advantage of the second-order accuracy of the central finite difference method in addition to the abso-lute stability property. Several examples having left, right, interior, or twin boundary layers were considered to illustrate the method. To analyze the effect of the deviating argument on the solution accuracy, the maximum absolute error is reported in several tables and figures for constant and variable deviating argument values. The results confirm that the variable deviating argument results in more accurate results than the constant one. Moreover, the results reveal that the deviating argument can stabilize the unstable discretized differential equation and that the present method is effective in solving the considered class of SPPs.