HIGH-ORDER EXPONENTIALLY FITTED DIFFERENCE SCHEMES FOR SINGULARLY PERTURBED TWO-POINT BOUNDARY VALUE PROBLEMS.

1. Introduction. In this paper we study the numerical solution of the singularly perturbed two-point boundary value problem for the ODE:

(1.1) [mathematical expression not reproducible]

(1.2) u(0) = [[alpha].sub.0], u(1)=[[alpha].sub.1].

We assume that the perturbation parameter [epsilon] and the coefficient c satisfy

(1.3) [epsilon] > 0 and c(x) [less than or equal to] 0 forx [member of] [0,1].

It is known that classical methods fail for this problem when [epsilon] is small relative to the mesh width. Various numerical approaches have been proposed, for example, difference schemes [5, 13, 15] and collocation methods [2, 6, 8, 9, 12]. An extensive review of methods is given by Kadalbayo and Patidar in .

To derive a difference scheme for the above problem, we consider functions that satisfy a differential equation of the form

(D - [[rho].sub.1])(D - [[rho].sub.2]) * * * (D - [rho]k)s = 0

for some given real numbers [[rho].sub.1], [[rho].sub.2], . . ., [rho]k. Here, D stands for the differentiation operator (Df = f = df/dx for some function f). For mutually distinct [[rho].sub.i], the solution of the above equation is an exponential sum

(1.4) [mathematical expression not reproducible]

Generally, if some [[rho].sub.i] are equal, we do not only have exponentials exp([[rho].sub.i]x) in the sum but also functions of the form x exp([[rho].sub.i]x), x2 exp([[rho].sub.i]x),... The resulting expressions are often called extended exponential sums  to be distinguished from the proper exponential sums (1.4). In this paper we consider extended exponential sums, and we will simply refer to them as exponential sums. The so-defined method belongs to a class of exponentially fitted methods. Exponential fitting is a well-known approach widely used for singularly perturbed problems .

The novelty of this paper is a new approach for the derivation of difference schemes. Using interpolation formulae, we approximate the derivatives in the differential equation (1.1) and obtain a difference scheme. The consistency of the schemes follows directly from the estimate for the interpolation error of the exponential sum. The presented results give an algorithm for the derivation of consistent schemes of arbitrary order. Unfortunately, we do not have a general result of stability (and convergence) for the proposed schemes. So, in this paper, we consider stability of three-point schemes and prove that they are first-order convergent in the singular perturbation case.

In the forthcoming sections we present basic facts about exponential sums, explain the derivation of difference schemes, prove consistency of the schemes, and analyze the convergence for three-point schemes. Finally, we exemplify our methods at several singularly perturbed problems.

2. Exponentially fitted differentiation formulae. As an approximation for u(x), we consider interpolation by particular classes of exponential sums. Since a solution of (1.1)-(1.2) asymptotically behaves as an exponential function on the boundary layer, we just use one exponential term. Hence, we use an exponential sum from the null-space of the differential operator

(2.1) [D.sup.k-1] - (D+[rho])s = 0,

and our exponential sum is of the form

(2.2) [mathematical expression not reproducible]

Let be given a real number [rho], an increasing sequence of points [mathematical expression not reproducible] ([t.sub.0] < t1 < ... < [t.sub.k-1]), and a sequence of real numbers (yi)0 . Then, there is a unique exponential sum (2.2) satisfying the interpolation conditions

(2.3) s([t.sub.i]) = [y.sub.i], i = 0,..., k - 1.

Interpolation by exponential sums is widely investigated in  and the references therein. is more convenient to write an exponential sum in Lagrange form. Let [mathematical expression not reproducible] be an exponential sum from the null-space of the differential operator (2.1) satisfying

(2.4) [mathematical expression not reproducible]

for j = 0,..., k - 1, i = -m,... ,k - 1 - m, and for some m (0 [less than or equal to] m [less than or equal to] k - 1). For equidistant points [mathematical expression not reproducible] ([t.sub.i+1] =[t.sub.i] + h,h> 0), the exponential sum

(2.5) [mathematical expression not reproducible]

satisfies the interpolation conditions (2.3). Note that the so-defined s is not from the null-space of (2.1) but from the null-space of

[D.sup.k-1] - (D+[rho]/h)s = 0.

The interpolation error for exponential sums is studied in  and . Here we reveal results for the error bound applied to the sum (2.5).

LEMM[A.sup.2].1 (). Let u [member of] [C.sup.k] (a,b), and let s be an exponential sum (2.5) of order k that interpolates u at an equidistant sequence of points [mathematical expression not reproducible], (a [less than or equal to] [t.sub.0] < . . . < [t.sub.k-1] [less than or equal to] b, h = [t.sub.i+1] - ti). Then for all x [member of] [a, b] and for r = 0,1,..., k - 1, the following estimates hold:

[mathematical expression not reproducible]

and

[mathematical expression not reproducible]

for r = 0,1,..., k - 1. The function Ck,r ([rho], y) satisfies

[mathematical expression not reproducible]

Here, by || ||[infinity] we denoted the maximum norm on the segment [a, b].

For [rho] = 0, the exponential sum is a polynomial of order k. In this case, the interpolation error is bounded by [h.sup.k]C,||[u.sup.(k)]||[infinity]. When [rho] [right arrow] [infinity], the exponential sum tends to a polynomial of order k - 1, and the bound reads [h.sup.k-1]C||[u.sup.(k-1)]||[infinity]. This is in concordance with results for the interpolation by polynomial functions.

However, when u is a solution of the singularly perturbed problem (1.1)-(1.2), due to the exponential behavior at the boundary layer, ||[u.sup.(k)] ||[infinity] is not bounded for small [epsilon]. For b(x) < 0, the asymptotic behavior is described by (cf. )

[mathematical expression not reproducible]

where

(2.6) [mathematical expression not reproducible]

and C is a constant independent of [epsilon] and x. When b(x) > 0, an analogous result holds (just substitute x [RIGHT ARROW] 1 - x). Further, a solution u may be bounded in the similar manner when b = 0.

The interpolation error when u is a solution of a singularly perturbed problem is studied in ; the main result is reported in the following lemma.

LEMMA 2.2 (). Let u be a solution ofa singularly perturbed boundary value problem (1.1)-(1.2), and let [mathematical expression not reproducible] be an equidistant sequence of points (0 = [x.sub.0] < [x.sub.1] < ... < [x.sub.n] = 1, h = [x.sub.i] - [x.sub.i]-1)for an arbitrarily chosen integer n satisfying

(2.7) h [greater than or equal to] 4(k - 1)[epsilon]ln(1/[epsilon])/[b.sub.min],

where the constant [b.sub.min] is defined by (2.6). Further, assume that

1. The functions b, c, and f are sufficiently smooth such that u [member of] [C.sup.k](0,1);

2. b(x) = 0 and c(x) [less than or equal to] 0 for all x [member of] [0,1];

3. The parameter [rho]from the exponential part of the exponential sum is of the same sign as the function b and satisfies K [less than or equal to] h/(|[rho]|[epsilon]) [mathematical expression not reproducible] for some positive constant K.

Then the exponential sum s of order k (k [greater than or equal to] 2) that interpolates the solution u at k consecutive mesh points (a subsequence of [mathematical expression not reproducible] denoted by [mathematical expression not reproducible] satisfies

[mathematical expression not reproducible]

for all x [member of] [t1, [t.sub.k-1]] when b(x) < 0 or all x [member of] [t.sub.2], [t.sub.k]] when b(x) > 0. The constant R is independent ofh and [epsilon]. Further,

[mathematical expression not reproducible]

for all x [member of] [t1, [t.sub.k]) when b(x) < 0 or all x [member of] (t1, [t.sub.k]] when b(x) > 0.

Note that Lemma 2.2 assumes b(x) = 0. A similar result may be obtained for the self-adjoint problem (b = 0).

3. A derivation of k-points difference schemes. To discretize the differential equation (1.1)-(1.2) we divide the interval [0,1] into n equal subintervals

0 = [x.sub.0] < [x.sub.1] < * * * < [x.sub.n-1] < [x.sub.n] = 1, h = 1/n, [x.sub.i] = ih.

A finite difference method comprises a discretization of the differential equation using the grid points [x.sub.i], where the unknown function u is approximated by an exponential sum of order k in the neighborhood of the knot [x.sub.i].

For given k (k [greater than or equal to] 3) and arbitrary i [member of] {0,..., n}, we choose [m.sub.i] [member of] {0,1,..., k - 1} satisfying 0 [less than or equal to] i - [m.sub.i] [less than or equal to] n - k + 1. Let [mathematical expression not reproducible] denote an exponential sum defined over k consecutive grid points [x.sub.i-mi],..., [x.sub.i+k-mi-1]. By taking advantage of the equidistant mesh, we may apply the representation of exponential sums given by (2.4) and (2.5):

[mathematical expression not reproducible]

The unknowns [u.sub.i] are determined by the approximation of u([x.sub.i]), u'([x.sub.i]), and u"([x.sub.i]) in (1.1) with [mathematical expression not reproducible], and [mathematical expression not reproducible], respectively:

[mathematical expression not reproducible]

Here we used the abbreviations [b.sub.i] = b([x.sub.i]), [c.sub.i] = c([x.sub.i]), and fi = f([x.sub.i]). Now,

[mathematical expression not reproducible]

(3.1) [mathematical expression not reproducible]

Since the boundary conditions (1.2) give two equations

(3.2) [u.sub.0] = [[alpha].sub.0], [u.sub.n] = [[alpha].sub.1],

we need n - 1 equations of the form (3.1) to determine the n + 1 unknowns [u.sub.j].

Not all of the choices for such equations will lead to a regular system of equations or to a stable method. However, all methods are consistent as we prove in the next section. A reasonable strategy for the choice of equations is to select mutually different [x.sub.i]'s. Later in this paper we construct and examine three convergent three-point methods to illustrate our approach.

Up to this point we did not discuss the choice of the parameter [rho] from the exponential sum si. We determine [rho] from the condition that the scheme is exact when the solution of the problem (1.1)-(1.2) with constant coefficients is an exponential function. In this case, [rho] is the largest (by the absolute value) root of the quadratic equation

(3.3) [mathematical expression not reproducible]

We denote by [[rho].sub.i] the parameter [rho] associated with the exponential sum [mathematical expression not reproducible]

In the case b(x) [not equal to] 0, [[rho].sub.i] is given by

(3.4) [mathematical expression not reproducible]

This definition is used for the choice of the tension parameters in the collocation by tension splines . Another possibility is to choose [rho] according to an asymptotic expansion of the solution in the boundary layer. This approach is widely used for exponentially fitted difference schemes. Such a choice leads to

(3.5) [mathematical expression not reproducible]

Both choices (3.4) and (3.5) are hardly distinguishable for small perturbation parameters [epsilon]. For [b.sub.i] > 0 we have

[mathematical expression not reproducible]

If c [equivalent to] 0, the choices (3.4) and (3.5) are equal.

When b [equivalent to] 0, the solution of the problem (1.1)-(1.2) has two boundary layers, and [rho] is calculated according to

[mathematical expression not reproducible]

4. Consistency of the difference schemes. Here we study the application of the exact solution u of the problem (1.1)-(1.2) to the scheme given by (3.1). Let us define

(4.1) [mathematical expression not reproducible]

where i and [m.sub.i] are arbitrarily chosen.

There is an unique exponential sum [S.sub.i], of order k that interpolates u(x) at the grid points xj, j = i - mi,..., i + k - [m.sub.i] - 1 (u([x.sub.j]) = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE]. For exponential sums of order k, the difference scheme (3.1) is exact:

[mathematical expression not reproducible]

Further, since u satisfies the differential equation (1.1), i.e., Lu(x) = f(x), [tau]i reads

[mathematical expression not reproducible]

Now, the approximation error for u'(x) is bounded by (cf. Lemma 2.1)

(4.2) [mathematical expression not reproducible]

and the approximation error for u"(x) is bounded by

(4.3) [mathematical expression not reproducible]

where

(4.4) [mathematical expression not reproducible]

Note that [mathematical expression not reproducible] interpolates u at the nodes [x.sub.i-mi],..., [x.sub.i-mi+k-1] and that the second argument of the functions [mathematical expression not reproducible]. Using (4.2) and (4.3) we obtain

(4.5) [mathematical expression not reproducible]

where [mathematical expression not reproducible] is defined by

(4.6) [mathematical expression not reproducible]

Lemma 2.2 and (4.5) imply that

[mathematical expression not reproducible]

for some constant R. The previous findings can be summarized in the next theorem.

THEOREM 4.1. Let n [greater than or equal to] k and k [greater than or equal to] 3 be arbitrary, and let u [member of] [C.sup.k](0,1) be a solution of (1.1)-(1.2). Then for any i [member of] {0,..., n} and [m.sub.i] (0 [less than or equal to] [m.sub.i] [less than or equal to] k - 1, 0 [less than or equal to] i - [m.sub.i] [less than or equal to] n - k + 1), the difference scheme defined by (3.1) is consistent, and [tau]i defined by (4.1) satisfies

[mathematical expression not reproducible]

and

[mathematical expression not reproducible]

where h = 1/n. The function [mathematical expression not reproducible]([[rho].sub.i]) is defined by (4.6). Further, if b(x) = 0 on [0,1], then there exists a constant C, independent on h, such that

[mathematical expression not reproducible]

for b(x) < 0 and [m.sub.i] < k - 1, or b(x) > 0 and [m.sub.i] > 0.

To study stability (and convergence) of difference schemes, it is more convenient to use a matrix notation. Our goal is to bound the error u([x.sub.i]) -ui, where u is the exact solution of the problem (1.1)-(1.2). Let us define the vectors

[mathematical expression not reproducible]

The equations (3.2) together with the additional n - 1 equations of the form (3.1) define a system matrix that we denote with Lk,n, and [u.sub.n] is a solution of the equation

(4.7) [L.sub.k,n][u.sub.n] = f .

We define a vector [tau] as

[tau] = [L.sub.k,n][u.sub.n]- [L.sub.k,n][u.sub.n] = [L.sub.k,n](u - [u.sub.n]).

If the matrix [L.sub.k,n] is regular, then

[mathematical expression not reproducible]

A bound for [tau] is given in Theorem 4.1, and thus the next corollary follows.

COROLLARY 4.2. Let the integers n [greater than or equal to] k and k [greater than or equal to] 3 be arbitrary, and let u [member of] C (0,1) be a solution of (1.1)-(1.2). If n - 1 equations of the form (3.1) together with (3.2) define a regular matrix [L.sub.k,n], then [u.sub.n] = [[u.sub.0], ..., [[u.sub.n]].sup. T], a solution of (4.7), satisfies

[mathematical expression not reproducible]

and

[mathematical expression not reproducible]

where h = 1/n and j = 0,1,... ,n. If b(x) = 0 on [0,1], then there exists a constant C, independent on h, such that

[mathematical expression not reproducible]

if the limit [mathematical expression not reproducible] exists and when b(x) < 0 and [m.sub.i] < k - 1 or b(x) > 0 and [m.sub.i] > 0.

There is no general approach to bound [mathematical expression not reproducible]. Instead, we have to consider each method separately.

5. Three-point difference schemes. Here we exemplify our approach in the simplest case, the three-point difference schemes. In the analysis of stability (and convergence), in addition to assumptions (1.3), we assume that the coefficient c satisfies

[mathematical expression not reproducible]

Further, the coefficient b is restricted to b(x) = 0 for all x, i.e., we exclude the so-called turning point case. We propose three three-point difference schemes (k = 3) with a uniform choice of mi, [M.sub.I] = m, m = 0,1,2.

The system matrix (4.7), here denoted by [L.sub.3,m,n], is a tridiagonal matrix

[mathematical expression not reproducible]

Since the approximation satisfies the boundary conditions, it holds that [tau]0 = 0 and [tau]n = 0. Therefore, we may omit the first and the last equation from the system, and we can consider the simplified matrix

[mathematical expression not reproducible]

Theorem 4.1 states that the methods are consistent. But, to prove convergence of the methods we have to consider stability, i.e., boundedness of the matrix [mathematical expression not reproducible]. This will be done in the forthcoming section.

5.1. A difference scheme for m = 0. The general scheme (3.1) for k = 3 and m = 0 reads:

[mathematical expression not reproducible]

with the end conditions [u.sub.0] = [[alpha].sub.0] and [u.sub.n] = [[alpha].sub.1]. A straightforward calculation results in the coefficients

[mathematical expression not reproducible]

for i = 0,..., n - 2.

LEMMA 5.1. Let -c(x) [greater than or equal to] [c.sub.min] > 0 and b(x) > 0 for x [MEMBER OF] [0,1]. If the [[rho].sub.i]'s are chosen according to (3.4), then the matrix L3,0,n is diagonally dominant, and it holds that

[mathematical expression not reproducible]

Proof First, note that from (3.3) it easily follows that

[mathematical expression not reproducible]

Therefore, the entries of the matrix L3,0,n satisfy

[mathematical expression not reproducible]

The coefficient [mathematical expression not reproducible] is given by

[mathematical expression not reproducible]

This inequality holds because of the positive signs of b and [[rho].sub.i]. Since

[mathematical expression not reproducible]

it follows that

[mathematical expression not reproducible]

Now we easily obtain that

[mathematical expression not reproducible]

and the bound from Lemma 5.1 is proved.

One may note that for the choice of the [[rho].sub.i]'s according to (3.5), [mathematical expression not reproducible] is positive for small [[rho].sub.i] and negative for large [[rho].sub.i]. Thus, a simple argument of diagonal dominance does not apply in this case.

5.2. A difference scheme for m = 1. The central-point scheme is given by

(5.1) [mathematical expression not reproducible]

i = 1,... ,n - 1, with two additional equations obtained from the boundary conditions [u.sub.0] = [[alpha].sub.0] and [u.sub.n] = [[alpha].sub.1]. The explicit expressions for the coefficients of the scheme (5.1) are

[mathematical expression not reproducible]

for i = 1,..., n - 1. The parameter [[rho].sub.i] is defined either by (3.4) or (3.5). A bound for the matrix [mathematical expression not reproducible] is given in the following lemma.

LEMMA 5.2. Let - c(x) [greater than or equal to] [c.sub.min] >0 andb(x) =0 for x [member of] [0,1]. If the [[rho].sub.i]'s are chosen according to (3.4) or (3.5), then the matrix L3,1,n is diagonally dominant, and

[mathematical expression not reproducible]

Proof First, we prove that the coefficients [mathematical expression not reproducible], and [mathematical expression not reproducible], i = 1,..., n - 1, satisfy [mathematical expression not reproducible]

From the definition of [[rho].sub.i] (3.4), we have

[mathematical expression not reproducible]

Hence,

(5.2) [mathematical expression not reproducible]

If [[rho].sub.i] is defined according to (3.5), then

(5.3) [mathematical expression not reproducible]

For r1 we have, using (5.2) or (5.3),

[mathematical expression not reproducible]

because [b.sub.i] and [[rho].sub.i] (and therefore [b.sub.i] and [e[rho].sub.i] -1) have the same sign. For ti we have, using (5.2) or (5.3), that

[mathematical expression not reproducible]

where we again used that [[rho].sub.i] and [b.sub.i] have the same sign. From the definition of [r.sub.i], [s.sub.i], and ti it is clear that [s.sub.i] + [r.sub.i] + [t.sub.i] = [c.sub.i], and

[mathematical expression not reproducible]

Since

[mathematical expression not reproducible]

we conclude that L3 1 n is diagonally dominant, and

[mathematical expression not reproducible]

5.3. A difference scheme for m = 2. When m = 2, the three-point scheme reads

(5.4) [mathematical expression not reproducible]

with the end conditions [u.sub.0] = [[alpha].sub.0] and [u.sub.n] = [[alpha].sub.1]. The coefficients from the scheme (5.4) are

[mathematical expression not reproducible]

[mathematical expression not reproducible]

[mathematical expression not reproducible]

for i = 2,..., n.

LEMMA 5.3. Let -c(x) [greater than or equal to] [c.sub.min] > 0 and b(x) > 0 for x G [0,1]. If the [[rho].sub.i] 's are chosen according to (3.4) or (3.5), then the matrix L3,2,n is diagonally dominant, and

[mathematical expression not reproducible]

for h = 1/n satisfying

(5.5) [mathematical expression not reproducible]

Proof. Because of

[mathematical expression not reproducible]

for [[rho].sub.i] [greater than or equal to] 0, the coefficients [mathematical expression not reproducible] and [mathematical expression not reproducible] satisfy

[mathematical expression not reproducible]

Also, for positive [[rho].sub.i],

[mathematical expression not reproducible]

Now,

[mathematical expression not reproducible]

The last inequality follows from the condition (5.5). As in the previous cases,

[mathematical expression not reproducible]

which proves the diagonal dominance and the bound of the lemma.

We summarize all previous results in the following theorem.

THEOREM5.4. Letu [member of] C (0,1) be the solution of (1.1)-(1.2), where - c(x) [greater than or equal to] [c.sub.min] > 0 and b(x) = 0 on [0,1]. Let [u.sub.n] = [[[u.sub.0], ..., [u.sub.n]].sup.T] denote a solution of the difference scheme (4.7), for k = 3, m [member of] {0,1, 2}, and h = 1/n. In addition, let the [[rho].sub.i]'s be chosen according to (3.4), and let h satisfy

[mathematical expression not reproducible]

in the cases m = 0 for b < 0, and m = 2 for b > 0. Then, the error at the nodes [x.sub.j], j = 0,..., n, is bounded by

[mathematical expression not reproducible]

and

[mathematical expression not reproducible]

The function Cm is defined by (4.6), while the norm [mathematical expression not reproducible] is defined by (4.4). Further, there exists a constant C, independent of h, such that

[mathematical expression not reproducible]

when b < 0 and m<2orb>0 and m > 0.

Proof. In the case m = 1, the assertion of the theorem follows directly from Corollary 4.2 by a direct substitution of [mathematical expression not reproducible] with the bound from Lemma 5.2.

Bounds for [mathematical expression not reproducible] and [mathematical expression not reproducible] are given by Lemma 5.3 and Lemma 5.1 in the case when b(x) > 0. When b(x) < 0, the simple substitution t = 1 - x transforms the problem (1.1)-(1.2) into the case -b(x) > 0. With such a substitution, the scheme for m = 0 becomes the scheme for m = 2 and vice versa. We may also verify that substituting [[rho].sub.i] into -[[rho].sub.i] and [b.sub.i] into - [b.sub.i] transforms the scheme m = 2 into the scheme for m = 0, as well as the scheme m = 0 into the scheme for m = 2. (Note that for [b.sub.i] > 0 the parameter [[rho].sub.i] is positive, while for [b.sub.i] < 0 the parameter [[rho].sub.i] is negative.) Thus, [mathematical expression not reproducible] and [mathematical expression not reproducible] are bounded for b(x) = 0, and the error estimate follows directly from Corollary 4.2.

Theorem 5.4 may be applied to the choice of [[rho].sub.i] 's given by (3.5) with the exception of the cases m = 0 for b > 0 and m = 2 for b < 0.

6. Numerical examples. In this section we illustrate our methods with examples having a known analytical expression for the exact solution.

EXAMPLE 6.1. Consider:

(6.1) [mathematical expression not reproducible]

(6.2) [mathematical expression not reproducible]

The exact solution of problem (6.1)-(6.2) is

(6.3) [mathematical expression not reproducible]

We apply difference schemes to this problem using 2n subintervals, for n = 2,..., 24. The perturbation parameter [epsilon] is set to [10.sup.-5]. The maximum error at the grid points in dependency on the number of subintervals is shown in Figure 6.1.

As could be expected, the central point scheme (m = 1) gives the best approximation. Further, we compare our methods to those proposed by Reddy and Chakravarthy  (denoted by RP) and Awoke and Reddy  (denoted by AR). Both methods are exponentially fitted, and the AR method is an improvement of the RP method. When [epsilon] is smaller than the mesh width h, the AR method is slightly better for this example than our method for m = 1, while the RP method gives a solution that is indistinguishable to our method. But both methods (AR and RP) fail to converge to the solution when h [right arrow] 0. This lack of convergence is not observed in some other examples (not shown here). If the inner solution is an exponential function, then the proposed schemes give the exact solution since the formulae used for the first and second derivative are exact for exponential functions. Therefore, the approximation error decreases in the case h < 4(k - 1)[epsilon] ln(1/[epsilon])/[b.sub.min], too. The difference schemes denoted by AR  and RP  are not exact for exponential functions and this may be the reason for their behavior.

In Figure 6.1 we also present the numerical order of convergence ([r.sub.n]). It is calculated according to

[mathematical expression not reproducible]

where [e.sub.n] is error for the partition with [2.sup.n] subintervals. For [epsilon] < h, the convergence of the methods is linear as was shown in Theorem 5.4. More precisely, for small [epsilon], the error is proportional to [epsilon] + h. In this case, [epsilon] may be neglected and the convergence is linear.

The same theorem also provides a bound for relatively large [epsilon]. Now, the error is proportional to [epsilon]h + h?, i.e., the convergence is at least linear. From Figure 6.1, it is evident that when h is relatively small with respect to [epsilon], then the method with m = 1 converges quadratically. An analysis of the method in the case h [much less than] 0 may explain this behavior.

In this example we used the choice of [[rho].sub.i] given by (3.4). We also tried the choice in (3.5), but there was no difference in the results.

Further, we illustrate some properties of the proposed methods that are not discussed in this paper. First, we applied high-order methods to the problem (6.1)-(6.2). All the used methods appeared to be convergent (Figure 6.2). The order of convergence is k - 2 for the k-point scheme.

EXAMPLE 6.2. An important property of methods for singularly perturbed problems is the independence of the convergence of the perturbation parameter [epsilon]. Such convergence is called [epsilon]-uniform convergence. A method is [epsilon]-uniformly convergent if there exists constants C and m, independent of [epsilon], such that the solution u of the problem (1.1)-(1.2) and its approximation at the mesh points [u.sub.i] satisfy

[mathematical expression not reproducible]

for all i.

Although Example 6.1 suggests [epsilon]-uniform convergence of the proposed methods, it is known that exponentially fitted methods are not [epsilon]-uniformly convergent . In  we analyzed the interpolation error by exponential sums for the solution of singularly perturbed boundary value problems, and we demonstrated that the interpolation error is not independent of [epsilon].

Here we consider the equation 

(6.4) [epsilon]u" + 2(1 - x)u' - 2u = ([epsilon] - 2x)[e.sup.x],

with boundary conditions

(6.5) u(0) = 2 and u(1) = [e.sup.-1/[epsilon]] + e.

The exact solution is given by

u(x) = [e.sup.-(2x-[x.sup.2])/[epsilon]] + [member of] [e.sup.x].

This solution has a boundary layer near the point x = 0.

The dependence of the error on the number of intervals for three-, four-, and five-point methods is presented in Figure 6.3. Figure 6.3 clearly demonstrates that the proposed methods are not [epsilon]-uniformly convergent. For h [greater than or equal to] 4(k - 1)[epsilon] ln(1/[epsilon])/[b.sub.min] the behavior of the error is described by Lemma 2.2 (i.e., Theorem 4.1): [mathematical expression not reproducible] (blue points in the figure). Experimental results in  and the figure suggest that the condition (2.7) can be extended to h [greater than or equal to] (k - 1)[epsilon] ln(1/[epsilon])/[b.sub.min] (red points in the figure). After that, the error starts to increase (except for the three-point method). For h [less than or equal to] [10.sup.-5] = [epsilon] the error again decreases, which is in accordance with the bound from Theorem 4.1 (bound for 'large' [rho]). Although the error for the three-point method seems not to be affected by [epsilon], in  we have shown that the interpolation error is not [epsilon]-uniform convergent for three-point methods. Numerical experiments from  suggest that an implementation of a dense mesh in the boundary layer results in [epsilon]-uniform convergence.

EXAMPLE 6.3. Next, we applied the methods to the a[d.sub.v]ection-diffusion problem 

(6.6) [epsilon]y" + y = - [e.sub.1-x],

(6.7) y(0) = 0, y(1) = 0,

with the exact solution

[mathematical expression not reproducible]

The dependence of the error on the number of intervals is displayed in Figure 6.4 (for [epsilon] = [10.sup.-5]). The behavior of the methods is the same as in the case c(x) = 0 (problem (6.1)-(6.2), Figure 6.1). The convergence of the methods is linear for h > [epsilon], which is in agreement with the theoretical results. The order of the convergence is the same in the case h < [epsilon] for m = 0 and m = 2. The method defined by m = 1 converges quadratically for small values of h.

EXAMPLE 6.4. In the last example we demonstrate the applicability of the proposed schemes to nonlinear problems. The problem

(6.8) [mathematical expression not reproducible]

with homogeneous boundary conditions y(0) = 0, y(1) = 0, and f [equivalent to] 0 has an approximate asymptotic solution

(6.9) [mathematical expression not reproducible]

In the example we impose

[mathematical expression not reproducible]

and boundary conditions

[mathematical expression not reproducible]

Now, the function (6.9) is the exact solution of equation (6.8). We solve this equation with the three-point method (k = 3, m = 1) and the four-point method (k = 4, m = 2). The nonlinear system is solved by five iterations of a Newton-Raphson method. The results are presented in Figure 6.5. Again, the behavior of the error is the same as in Examples 6.1 and 6.3.

It is noteworthy that condition (1.3) is important for the convergence of the method. In a semilinear problem of the form [epsilon] y" + by' + a(x, y) = f(x), the assumption [mathematical expression not reproducible] [less than or equal to] 0 guarantees convergence. For example, the method did not converge for the problem [epsilon] y" + 2y' + [e.sup.y] = f(x).

7. Concluding remarks. In this paper we introduced a family of exponentially fitted difference schemes for singularly perturbed two-point boundary value problems. The schemes are derived from interpolating formulae for exponential sums. The consistency of the considered difference schemes follows from the error bounds given in . Although this paper deals with exponential sums with a basis of the form 1, x,..., [x.sup.k-2], exp([rho]x), interpolating formulae may be applied to a wider class of exponential sums used in the literature, for example, sums with a basis 1, x,..., [x.sup.k-2], exp([rho]x), x exp([rho]x), ..., [x.sup.l] exp([rho]x) or 1, x,..., [x.sup.k-2], exp([rho]x), exp(-[rho]x). Consistency may follow from a generalization of the results in .

We illustrated this approach for three three-point difference schemes. All three methods are stable and convergent. As expected, the examples show that the central difference scheme gives the best approximation. The proposed methods are comparable to other exponentially fitted three-point methods such as, for example, the methods proposed in [1, 14]. Still, our methods have better properties than these two methods.

For exponential fitting, we use the fitting parameter [rho] in (3.4) proposed for collocation methods by tension spline [9, 12]. All three methods are convergent for such a choice of [rho] as well as for the choice (3.5) commonly used in similar schemes. Moreover, the performed numerical tests indicate no difference between the approximations obtained by these two choices of [rho]. Some methods from the literature use a uniform choice of [rho]. In our case, a nonuniform choice guarantees convergence. Numerical tests show that a uniform choice of [rho] does not improve the method. Further, there are methods that use exponential fitting in the boundary layer only (for example ) or that behave well when exponential fitting is applied only in the boundary layer although the methods assume the exponential fitting over a whole interval (for example ). Such approach is in concordance with the fact that the solution of singularly perturbed problem exhibits strong exponential behavior only in the boundary layer. This is not possible in our case since exponential fitting on a whole interval is essential for the stability of the methods. Numerical tests show that the methods did not converge when exponential fitting was applied on the boundary layer alone.

We derived high-order schemes, higher than the seventh-order scheme given in . The obtained results may be directly applied to schemes of higher order. We already have results for the stability of four-point schemes. However, a stability analysis of such schemes would be more complex.

Further, in the stability analysis we consider a singular perturbation problem with b(x) = 0 and c(x) < 0. Numerical examples show that these schemes are applicable to equations where b(x) = 0 or c(x) = 0. In these cases, stability may be analyzed in the same way.

The proposed schemes are not [epsilon]-uniform convergent. However, convergence is guaranteed up to the point when the mesh width h becomes smaller than 4(k - 1)[epsilon]ln(1/[epsilon])/[b.sub.min]. If higher accuracy is required, then higher-order schemes may be used. The construction of a grid that would result in [epsilon]-uniform convergence is currently under investigation.

REFERENCES

 A. AWOKE AND Y. N. REDDY, An exponentially fitted special second-order finite difference method for solving singular perturbation problems, Appl. Math. Comput., 190 (2007), pp. 1767-1782.

 T. BOSNER, Basis of splines associated with singularly perturbed a[d.sub.v]ection-diffusion problems, Math. Commun., 15 (2010), pp. 1-12.

 D. BRAESS, Nonlinear Approximation Theory, Springer, Berlin, 1986.

 P. P. CHAKRAVARTHY, K. PHANEENDRA, AND Y. N. REDDY, A seventh order numerical method for singular perturbation problems, Appl. Math. Comput., 186 (2007), pp. 860-871.

 P. A. FARRELL, A. F. HEGARTY, J. J. H. MILLER, E. O'RIORDAN, AND G. I. SHISHKIN, Robust Computational Techniques for Boundary Layers, Chapman & Hall, Boca Raton, 2000.

 J. E. FLAHERTY AND W. MATHON, Collocation with polynomial and tension splines for singularly-perturbed boundary value problems, SIAM J. Sci. Statist. Comput., 1 (1980), pp. 260-289.

 M. K. KADALBAJOO AND K. C. PATIDAR, A survey of numerical techniques for solving singularly perturbed ordinary differential equations, Appl. Math. Comput., 130 (2002), pp. 457-510.

 I. KAVCIC, M. ROGINA, AND T. BOSNER, Singularly perturbed a[d.sub.v]ection-diffusion-reaction problems: comparison of operator-fitted methods, Math. Comput. Simulation, 81 (2011), pp. 2215-2224.

 M. MARUSIC, A fourth/second order accurate collocation method for singularly perturbed two-point boundary value problems using tension splines, Numer. Math., 88 (2001), pp. 135-158.

 __, Limit properties of interpolation by exponential sums, in Proceedings of BAIL 2002, Perth, 2002, S. Wang and N. Fowkes, eds., University of Western Australia, Perth, 2002, pp. 183-188.

 __, On e-uniform convergence of exponentially fitted methods, Math. Commun., 19 (2014), pp. 545-559.  M. MARUSIC AND M. ROGINA, A collocation method for singularly perturbed two-point boundary value problems with splines in tension, A[d.sub.v]. Comput. Math., 6 (1996), pp. 65-76.

 J. J. H. MILLER, E. O'RIORDAN, AND G. I. SHISHKIN, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, River Edge, 1996.

 Y N. REDDY AND P. P. CHAKRAVARTHY, An exponentially fitted finite difference method for singular perturbation problems, Appl. Math. Comput., 154 (2004), pp. 83-101.

 H.-G. ROOS, M. STYNES, AND L. TOBISKA, Numerical Methods for Singularly Perturbed Differential Equations, Springer, Berlin, 1996.

 G. I. SHISHKIN, Approximation of solutions of singularly perturbed boundary value problems with a parabolic boundary layer, U.S.S.R. Comput. Math. and Math. Phys., 29 (1989), pp. 1-10.

 J. VIGO-AGUIAR AND S. NATESAN, An efficient numerical method for singular perturbation problems, J. Comput. Appl. Math., 192 (2006), pp. 132-141.

MILJENKO MARUSIC ([dagger])

(*) Received July 25, 2016. Accepted June 4, 2018. Published online on September 11, 2018. Recommended by M. Gander. Work supported by the grant #037-0982913-2762 of the Croatian Ministry of Science, Education and Sports.

([dagger]) Department of Mathematics, University of Zagreb, Bijenicka c. 30, 10 000 Zagreb, Croatia (miljenko.marusic@math.hr).

DOI: 10.1553/etna_vol48s329
COPYRIGHT 2018 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.