# Parameter-uniform fitted operator B-spline collocation method for self-adjoint singularly perturbed two-point boundary value problems.

1. Introduction. We consider the following self-adjoint singularly perturbed boundary value problem,

(1.1) Ly [equivalent to] -[epsilon](a(x)y')' + b(x)y(x) = f(x), x [member of] [bar.[OMEGA]] = [0, 1],

with the Dirichlet boundary conditions,

(1.2) y(0) = [alpha], y(1) = [beta],

where [alpha] and [beta] are given constants and [epsilon] is a small positive parameter. The functions a(x), b(x), and f(x) are sufficiently smooth and satisfy

a(x) [greater than or equal to] [a.sup.*] > 0, b(x) [greater than or equal to] [b.sup.*] > 0, a'(x) [greater than or equal to] 0. (1.3)

Under these conditions the operator L admits the maximum principle .

These types of problems, in which a small parameter multiplies the highest derivative, are known as singular perturbation problems. They arise in the mathematical modeling of physical and chemical processes, for instance, reaction diffusion processes, chemical reactor theory, fluid mechanics, quantum mechanics, fluid dynamics, elasticity, etc. Schatz and Wahlbin  and Boglaev  solved this type of problem by using finite element techniques. Miller  gave a sufficient condition for the uniform first order convergence of a general three-point difference scheme, whereas Niijima  gave a uniformly second order accurate difference scheme. Kadalbajoo and Aggarwal  introduced the B-spline collocation method with a fitted mesh technique for self-adjoint singularly perturbed boundary value problems and proved second order uniform convergence of the method. It is well known that the solution of the SPP converges as [epsilon] [right arrow] 0, 0 [less than or equal to] x [less than or equal to] 1, to the solution of the reduced problem obtained by putting [epsilon] = 0 in the original problem.

Parameter-uniform numerical methods [7, 8] are methods, whose numerical approximations UN satisfy error bounds of the form

[parallel][u.sub.[epsilon]] - [U.sup.N][parallel] [less than or equal to] Cv(N), v(N) [right arrow] 0 as N [right arrow] [infinity],

where [u.sub.[epsilon]] is the solution of the continuous problem, [parallel].[parallel] is the maximum pointwise norm, N is the number of mesh points (independent of [epsilon]) used, and C is a positive constant, which is independent of both [epsilon] and N. In other words, the numerical approximations UN converge to [u.sub.[epsilon]] for all values of [epsilon] in the range 0 < [epsilon] [much less than] 1.

There are two well-known approaches to obtain small truncation error inside the boundary layer/layers when the perturbation parameter [epsilon] is very small. The first is a fitted method based on choosing a fine mesh in the boundary layer/layers region [8, 9, 10, 11] and the second is based on fitted operator methods, i.e., a difference formula reflecting the behavior of the solution in the boundary layer/layers. In this paper, we use the second strategy, and the proposed method, based on a suitably designed fitted operator to the interior layer, is shown to converge with v(N) = 1/[N.sup.2].

2. Continuous problem. The bounds for the solutions and its derivatives are given in this section. Furthermore, the bounds for the smooth and singular components and their derivatives are also given. We first give the maximum principle and stability estimates for the solution of the problem (1.1) and (1.2).

LEMMA 2.1 (Continuous maximum principle). Let [phi] [member of] [C.sup.2]([bar.[OMEGA]]), satisfying [phi](0) [greater than or equal to] 0, [phi](1) [greater than or equal to] 0 and L[phi](x) [greater than or equal to] 0 [for all] x [member of]. Then [phi](x) [greater than or equal to] 0 [for all] x [member of] [bar.[OMEGA]].

Proof. The proof is by contradiction. Suppose that there is a point [x.sup.*] [member of] [bar.[OMEGA]], such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. It is clear from the given conditions that [x.sup.*] [not member of] {0, 1}. Therefore [phi]'([x.sup.*]) = 0 and [phi]"([x.sup.*]) [greater than or equal to] 0. Thus, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which is a contradiction. It follows that [phi]([x.sup.*]) [greater than or equal to] 0 and so [phi](x) [greater than or equal to] 0 [for all] x [member of] [bar.[OMEGA]].

LEMMA 2.2 (Stability). Consider the problem (1.1) and (1.2). If [phi](x) is the solution of this problem, then for some positive constant C, we have

[parallel][phi](x)[parallel] [less than or equal to] C(max(|[alpha]|, |[beta]|) + 1/[b.sup.*] [parallel]f[parallel]), [for all] x [member of] [bar.[OMEGA]].

Proof. Consider the barrier functions

[[PSI].sup.[+ or -]](x) = C(max(|[alpha]|, |[beta]|) + 1/[b.sup.*] [parallel]f[parallel]) [+ or -] [phi](x).

Then it can be easily seen that [[PSI].sup.[+ or -]](0) [greater than or equal to] 0, [[PSI].sup.[+ or -]](1) [greater than or equal to] 0 and for all x [member of] [OMEGA], L[[PSI].sup.[+ or -]](x) [greater than or equal to] 0, by a proper choice of C. Therefore, by applying Lemma 2.1, we obtain [[PSI].sup.[+ or -]](x) [greater than or equal to] 0 [for all] x [member of] [bar.[OMEGA]], which gives the required estimate.

3. Reduction to normal form. Rewrite (1.1) as

-[epsilon]a(x)y" - [epsilon]a'(x)y' + b(x)y = f(x),

or

(3.1) y" + P(x)y' + Q(x)y = F(x),

where

P(x) = a'(x)/a(x), Q(x) = - b(x)/[epsilon]a(x), and F(x) = - f(x)/[epsilon]a(x).

Let

U(x) = exp(- 1/2 [[integral].sup.x.sub.0]P([zeta])d[zeta]).

Consider the transformation

(3.2) y(x) = U(x)V (x).

Then (3.1) can be written in the normal form as

(3.3) V" + [??](x)V = [??](x),

with

V(0) = y(0)/U(0) = [[eta].sub.1], V(1) = y(1)/U(1) = [[eta].sub.2], [[eta].sub.1], [[eta].sub.2] [member of] R,

where

[??](x) = Q(x) - 1/2 P'(x) - 1/4 [(P(x)).sup.2], [??](x) = F(x) exp (1/2 [[integral].sup.x.sub.0]P([zeta]) d[zeta]).

Multiplying (3.3) throughout by -[epsilon], we get

(3.4) [??]V [equivalent to] -[epsilon]V" + W(x)V = Z(x),

with

V(0) = [[eta].sub.1], V(1) = [[eta].sub.2], (3.5)

where

W(x) = -[epsilon][??] (x), Z(x) = -[epsilon][??](x), and W(x) [greater than or equal to] [W.sup.*] > 0.

Note that

Z(x) = -[epsilon][??](x) = -[epsilon]F(x) exp(1/2 [[integral].sup.x.sub.0]P([zeta])d[zeta]) = f(x)/a(x) exp (1/2 [[integral].sup.x.sub.0]P([zeta])d[zeta]).

This shows that Z(x) is independent of [epsilon]. However, W(x) may or may not depend on [epsilon].

Roos  constructed a global uniformly convergent (in [epsilon]) scheme for (3.4) and (3.5), by replacing the coefficients by piecewise polynomials and solved the resulting problem exactly. To solve (3.4) and (3.5), Surla and Jerkovic  derived a difference scheme via a spline in tension and obtained error estimates of the form O(h,min(h, [epsilon])). O'Riordan and Stynes [14, 15, 16, 17] introduced the concept of freezing the coefficients by considering the piecewise constants on each subinterval [[x.sub.i-1], [x.sub.i]] as an approximation for the coefficient terms W(x) and Z(x) of the singularly perturbed boundary value problem(3.4) and (3.5). Stojanovic  gave an optimal difference scheme by considering the quadratic interpolating splines instead of piecewise constants on each subinterval [[x.sub.i-1], [x.sub.i]] as an approximation for the coefficient Z(x). We define the fitting problem associated with (3.4) and (3.5) by

[L.sub.[sigma]] [equivalent to] -((x, [epsilon])V" + W(x)V = Z(x), (3.6)

with

V(0) = [[eta].sub.1], V(1) = [[eta].sub.2], (3.7)

where [sigma](x, [epsilon]) is a fitting factor, which is to be determined subsequently.

4. B-Spline collocation method. In this section, we describe a B-spline collocation method to obtain the approximate solution of boundary value problems (3.6) and (3.7). Let P [equivalent to] {0 = [x.sub.0] < [x.sub.1] < [x.sub.2] ... < [x.sub.N-1] < [x.sub.N] = 1} be the partition of [bar.[OMEGA]] with uniform spacing h = 1/N. We include two more points on each side of the partition P as [x.sub.-2] < [x.sub.-1] < [x.sub.0] and [x.sub.N+2] > [x.sub.N+1] > [x.sub.N]. Then P [equivalent to] {[x.sub.-2] < [x.sub.-1] < [x.sub.0] = 0 < [x.sub.1] < [x.sub.2] ... < [x.sub.N-1] < [x.sub.N] = 1 < [x.sub.N+1] < [x.sub.N+2]}. Let [L.sub.2]([bar.[OMEGA]]) be the space of all square integrable functions on [bar.[OMEGA]] and let X be a linear subspace of [L.sub.2]([bar.[OMEGA]]). We use the cubic B-spline basis functions [B.sub.i](x) for i = 0(1)N; see .

It is easy to see that each [B.sub.i](x) is also a piecewise cubic with knots at [[bar.[OMEGA]].sub.N] and [B.sub.i](x) [member of] X. Thus, [B.sub.i](x) is twice continuously differentiable [for all] x [member of] R. The values of [B.sub.i](x), [B'.sub.i](x) and [B".sub.i](x) at the nodal points [x.sub.i]'s are given in Table 4.1.

Let B = {[B.sub.-1], [B.sub.0], [B.sub.1], ...., [B.sub.N+1]} and let [[PHI].sub.3]([[bar.[OMEGA]].sub.N]) = {[[phi].sub.i] : [[phi].sub.i] = [N+1.summation over (i = -1)][k.sub.i][B.sub.i], [k.sub.i] [member of] R}.

The functions [B.sub.i](x) are linearly independent on [bar.[OMEGA]]. Thus, [[PHI].sub.3]([[bar.[OMEGA]].sub.N]) is (N + 3)-dimensional subspace of X. Let [L.sub.[sigma]] be a linear operator with domain X and with range in X. We seek a function S(x) [member of] [[PHI].sub.3]([[bar.[OMEGA]].sub.N]) that approximates the solution of boundary value problem (3.6) and (3.7), represented as

(4.1) S(x) = [N+1.summation over (i = -1)][c.sub.i][B.sub.i](x),

where [c.sub.i] are unknown real coefficients. Here we have introduced two extra cubic B-splines, [B.sub.-1] and [B.sub.N+1] to satisfy the boundary conditions. Therefore, we have

(4.2) [L.sub.[sigma]]S([x.sub.i]) = Z([x.sub.i]), 0 [less than or equal to] i [less than or equal to] N,

and

(4.3) S([x.sub.0]) = [[eta].sub.1], S([x.sub.N]) = [[eta].sub.2].

Using Table 4.1 and equation (4.1), the system of collocation (4.2), together with boundary conditions (4.3), gives a system of (N + 1) linear equations in (N + 3) unknowns,

(4.4) (-6[[sigma].sub.i]/[h.sup.2] + [W.sub.i])[c.sub.i-1] + (12[[sigma].sub.i]/[h.sup.2] + 4[W.sub.i])[c.sub.i] + (-6[[sigma].sub.i]/[h.sup.2] + [W.sub.i])[c.sub.i+1] = [Z.sub.i],

for 0 [less than or equal to] i [less than or equal to] N, where W([x.sub.i]) = [W.sub.i], Z([x.sub.i]) = [Z.sub.i] and [[sigma].sub.i] is a fitting factor which is to be determined. The given boundary conditions become

(4.5) [c.sub.-1] + 4[c.sub.0] + [c.sub.1] = [[eta].sub.1]

and

(4.6) [c.sub.N-1] + 4[c.sub.N] + [c.sub.N+1] = [[eta].sub.2].

Thus, (4.4), (4.5), and (4.6) lead to a (N + 3) x (N + 3) system with (N + 3) unknowns [c.sup.N] = [([c.sub.-1], [c.sub.0], [c.sub.1], ...., [c.sub.N+1]).sup.t]. Now eliminating [c.sub.-1] from the first equation of (4.4) and from (4.5), we get

(4.7) (36[[sigma].sub.0]/[h.sup.2])[c.sub.0] = [Z.sub.0] - [[eta].sub.1](-6[[sigma].sub.0]/[h.sup.2] + [W.sub.0]).

Similarly, eliminating [c.sub.N+1] from the last equation of (4.4) and Eq. (4.6), we obtain

(4.8) (36[[sigma].sub.N]/[h.sup.2])[c.sup.N] = [Z.sub.N] - [[eta].sub.1]2](-6[[sigma].sub.N]/[h.sup.2] + [W.sub.N]).

Taking (4.7) and (4.8) with the second through (N - 1)st equations of (4.4), we are lead to a system of (N + 1) linear equations,

T[c.sup.N] = [d.sup.N], (4.9)

in (N + 1) unknowns [c.sup.N] = [([c.sub.0], [c.sub.1], ..., [c.sup.N]).sup.t] with right hand side

[d.sup.N] = ([Z.sub.0] - [[eta].sub.1](-6[[sigma].sub.0]/[h.sup.2] + [W.sub.0]), [Z.sub.1], ...,[Z.sub.N-1], [Z.sub.N] - [[eta].sub.2](-6[[sigma].sub.N]/[h.sup.2] + [W.sub.N])).

The coefficient matrix T is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Since W(x) > 0, it can be easily seen that the matrix T is strictly diagonally dominant and hence nonsingular. Since T is nonsingular, we can solve the system T[c.sup.N] = [d.sup.N] for [c.sub.0], [c.sub.1], ..., [c.sup.N] and substitute into the boundary condition (4.5) and (4.6) to obtain [c.sub.-1] and [c.sub.N+1]. Hence the method of collocation using a basis of cubic B-splines applied to problem (3.6) and (3.7) has a unique solution S(x) given by (4.1).

5. Determination of a fitting factor. In order to obtain a fitting factor, we shall use the following lemma .

LEMMA 5.1. Let V(x) [member of] [C.sup.4]([bar.[OMEGA]]), and let W'(0) = W'(1) = 0, then the solution of the problem (3.4) and (3.5) is of the form

V(x) = u(x) + v(x) + w(x),

where

u(x) = p exp (-x[square root of (W(0)/[epsilon])]), v(x) = q exp (-(1 - x)[square root of W(1)/[epsulon])]),

where p and q are bounded functions of [epsilon] independent of x and

[absolute value of w(j)(x)] [less than or equal to] C(1 + [[epsilon].sup.(1-j)/2]), j = 0(1)4,

C is a constant independent of [epsilon].

It is clear that the matrix T is inverse monotone if [h.sup.2][W.sub.j]/6[[sigma].sub.j] [less than or equal to] 1, thus taking the fitting factor of the form,

[[sigma].sub.j] = [h.sup.2][W.sub.j]/6 [mu]([[rho].sub.j]),

where [[rho].sub.j] = [square root of ([W.sub.j]/[epsilon])] and [mu]([[rho].sub.j]) is to be determined. We want to find a fitting factor such that the truncation error for the boundary layer function is equal to zero when W(x) is constant.

From the condition T[u.sub.i] = 0, for W(x) = W = constant, we have

(5.1) (-6[sigma]/[h.sup.2] + W)[u.sub.i-1] + (12[sigma]/[h.sup.2] + 4W)[u.sub.i] + (-6[sigma]/[h.sup.2] + W)[u.sub.i+1] = 0,

for 0 [less than or equal to] i [less than or equal to] N, where

[u.sub.i-1] = p exp (-[x.sub.i-1][square root of (W(0)/[epsilon])]) = [u.sub.i] exp (h[square root of (W(0)/[epsilon])]), [u.sub.i+1] = p exp (-[x.sub.i+1][square root of (W(0)/[epsilon])]) = [u.sub.i] exp (-h[square root of (W(0)/[epsilon])]).

Substituting in (5.1), a simple calculation gives

(5.2) [sigma]([rho]) = [h.sup.2]W/6 [mu]([rho]), when W(x) = constant,

(5.3) [sigma]([[rho].sub.j]) = [h.sup.2][W.sub.j]/6 [mu]([[rho].sub.j]), when W(x) [not equal to] constant,

where

[mu]([[rho].sub.j]) = 1+2[cosh.sup.2]([[rho].sub.j]h/2)/2[sinh.sup.2]([[rho].sub.j]h/2).

6. Convergence of the scheme. In order to prove the uniform convergence of the scheme, first we shall prove the following lemma.

LEMMA 6.1. Let [[sigma].sub.j] be the fitting factor determined in Section 5, then [[sigma].sub.j] approximates [epsilon] with an error O([h.sup.2]), i.e.,

|[[sigma].sub.j] - [epsilon]| [less than or equal to] C[h.sup.2],

for some positive constant C.

Proof. It is clear from (5.3) that

0 [less than or equal to] [[sigma].sub.j] [less than or equal to] C[h.sup.2].

There are two cases:

Case I. [epsilon] [less than or equal to] [h.sup.2]. In this case we have

(6.1) |[[sigma].sub.j] - [epsilon]| [less than or equal to] |[[sigma].sub.j]| + |[epsilon]| [less than or equal to] C[h.sup.2].

Case II. [h.sup.2] [less than or equal to] [epsilon]. In this case we have

This gives

(6.2) |[[sigma].sub.j] - [epsilon]| [less than or equal to] C[h.sup.2].

Therefore for all [epsilon], we have from (6.1) and (6.2)

|[[sigma].sub.j] - [epsilon]| [less than or equal to] C[h.sup.2].

LEMMA 6.2. The B-splines [B.sub.-1], [B.sub.0], ..., [B.sub.N+1] satisfy the inequality

[N+1.summation over (i = -1)]|[B.sub.i](x)| [less than or equal to] 10, 0 [less than or equal to] x [less than or equal to] 1.

Proof. The lemma is shown in . It has previously been applied in .

THEOREM 6.3. Let W(x),Z(x) [member of] [C.sup.2][0, 1] and W(x) [greater than or equal to] [W.sup.*] > 0. The collocation approximation S(x) from the space [[phi].sub.3]([OMEGA]) to the solution V(x) of the boundary value problem (3.6) and (3.7) exists and the error bound is given by the following relation

||V - S|| [less than or equal to] C[h.sup.2],

for sufficiently small h and a positive generic constant C.

Proof. To estimate the error ||V(x) - S(x)||, let [Y.sub.n] be the unique spline interpolant from [[PHI].sub.3]([[bar.[OMEGA]].sub.N]) to the solution V(x) of our boundary value problem (3.6) and (3.7). Since Z(x) [member of] [C.sup.2]([bar.[OMEGA]]), V (x) [member of] C4([bar.[OMEGA]]) and it follows from de Boor-Hall error estimates that

(6.3) ||[D.sup.j](V(x) - [Y.sub.n])|| [less than or equal to] [[lambda].sub.j][h.sup.4-j], j = 0, 1, 2,

where h is the uniform mesh spacing and C is a constant independent of h and N. Let

[Y.sub.n](x) = [N+1.summation over (i = -1)][b.sub.i][B.sub.i](x).

It follows immediately from the estimates (6.3) that

(6.4) |[L.sub.[sigma]]S([x.sub.i]) - [L.sub.[sigma]][Y.sub.n]([x.sub.i])| [less than or equal to] [omega][h.sup.2],

where [omega] = [[??][[lambda].sub.2] + [[lambda].sub.0]||W||[h.sup.2]], with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Let [L.sub.[sigma]][Y.sub.n]([x.sub.i]) = [[??].sub.n]([x.sub.i]) [for all] i and [[??].sub.n] = ([[??].sub.n]([x.sub.0]), [[??].sub.n][([x.sub.1]), ..., [[??].sub.n]([x.sub.N])).sup.t]. Now it is clear from (6.4) that the ith coordinate [[T([c.sup.N] - [b.sup.N])].sub.i] of T ([c.sup.N] - [b.sup.N]), where [b.sup.N] = [([b.sub.0], [b.sub.1], ..., [b.sub.N]).sup.t], satisfies the inequality

(6.5) |[[T([c.sup.N] - [b.sup.N])].sub.i]| [less than or equal to] [omega][h.sup.2].

Since [(T[c.sup.N]).sub.i] = Z([x.sub.i]) and [(T[b.sup.N]).sub.i] = [[??].sub.n]([x.sub.i]) [for all] i = 0(1)N. The ith coordinate of [T([c.sup.N] - [b.sup.N])] is the ith equation,

(6.6) (-6[[sigma].sub.i]/[h.sup.2] + [W.sub.i])[[delta].sub.i-1] + (12[[sigma].sub.i]/[h.sup.2] + 4[W.sub.i])[[delta].sub.i] + (-6[[sigma].sub.i]/[h.sup.2] + [W.sub.i])[[delta].sub.i+1] = [[xi].sub.i],

for 1 [less than or equal to] i [less than or equal to] N - 1, where

[[delta].sub.i] = [b.sub.i] - [c.sub.i], -1 [less than or equal to] i [less than or equal to] N + 1 and [[xi].sub.i] = Z([x.sub.i]) - [[??].sub.n]([x.sub.i]), 1 [less than or equal to] i [less than or equal to] N - 1.

Obviously, |[[xi].sub.i]| [less than or equal to] [omega][h.sup.2]. Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Also consider [delta] = [([[delta].sub.-1], [[delta].sub.0], ..., [[delta].sub.N+1]).sup.t] and define [e.sub.i] = |[[xi].sub.i]| and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Taking absolute value for h sufficiently small, and using the fact that 0 < [W.sup.*] [less than or equal to] W(x) [for all] x [member of] [bar.[OMEGA]], we have from (6.6),

(12[[sigma].sub.i]/[h.sup.2] + 4[W.sup.*])[??] [less than or equal to] [??] + 2[??](6[[sigma].sub.i]/[h.sup.2] + [W.sup.*]).

This gives

(6.7) [??] [less than or equal to] [omega][h.sup.2]/2[W.sup.*] .

Now to estimate [e.sub.-1], [e.sub.0], [e.sub.N] and [e.sub.N+1], we observe that the first and last equation of the the system T([c.sup.N] - [b.sup.N]) = ([Z.sup.n] - [[??].sub.n]) where [Z.sup.n] = ([Z.sub.0], [Z.sub.1], ..., [Z.sub.N]), gives

[e.sub.0] [less than or equal to] [omega][h.sup.4]/36[[sigma].sub.0] and [e.sub.N] [less than or equal to] [omega][h.sup.4]/36[[sigma].sub.N].

Now [e.sub.-1] and [e.sub.N+1] can be evaluated using boundary conditions (4.5) and (4.6) as

[e.sub.-1] [less than or equal to] [omega][h.sup.4]/9[[sigma].sub.0] + [omega][h.sup.2]/2[W.sup.*] and [e.sub.N+1] [less than or equal to] [omega][h.sup.4]/9[[sigma].sub.N] + [omega][h.sup.2]/2[W.sup.*].

Therefore, using the value [omega] = [[??][[lambda].sub.2] + [[lambda].sub.0]||W||[h.sup.2]], we get

(6.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], provided [h.sup.6]/[[sigma].sub.0] ~ 0, [h.sup.6]/[[sigma].sub.N] ~ 0.

Now we have

(6.9) ||S(x) - [Y.sub.n](x) = [N+1.summation over (i = -1)]([c.sub.i] - [b.sub.i])[B.sub.i](x).

Thus using (6.8) and Lemma 6.2, we get

(6.10) ||S - [Y.sub.n]|| [less than or equal to] C[h.sup.2].

Now using (6.3) and (6.10), we obtain

||V - S|| [less than or equal to] C[h.sup.2]

as required.

7. Test examples and numerical results. To demonstrate the efficiency of the method, several examples having boundary layer at one or both end points has been analyzed. For each [epsilon] and N, the maximum absolute errors at nodal points are approximated by the formula

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where y([x.sub.j]) and [y.sub.j] are the exact and computed solution of the given problem and nodal points [x.sub.j]. Also for each N the [epsilon]-uniform error at nodal point is approximated by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

EXAMPLE 7.1. First we consider the problem ,

-[epsilon]y" + y = -[cos.sup.2]([pi]x) - 2[epsilon][[pi].sup.2] cos(2[pi]x), y(0) = 0, y(1) = 0.

The exact solution is given by

exp(-(1 - x)/[square root of [epsilon]])+exp(-x/[square root of [epsilon]]) / 1+exp(-1/[square root of [epsilon]]) - [cos.sup.2]([pi]x).

EXAMPLE 7.2. Now consider the following singular perturbation problem ,

-[epsilon]y" + 4/[(x+1).sup.4] [1 + [square root of [epsilon]](x + 1)]y = f(x), y(0) = 2, y(1) = -1,

where

f(x) = - 4/(x + 1)4 x [{1 + [square root of [epsilon]](x + 1) + 4[[pi].sup.2][epsilon]} cos(4[pi]x/x+1) - 2[pi][epsilon](x + 1) sin(4[pi]x/x+1) + 3{1+[square root of [epsilon]](x+1)} / exp(1/[square root of [epsilon]])-1].

Its exact solution is given by

y(x) = -cos(4[pi]x/x+1) + 3{exp(-2x/[square root of [epsilon]](x+1))-exp(-1/[square root of [epsilon]])} / 1-exp(-1/[square root of [epsilon]]).

EXAMPLE 7.3. Finally consider the problem ,

-[epsilon][(1 + [x.sub.2])y']' + [cosx/[(3-x).sup.3]] y = 4(3[x.sup.2] - 3x + 1)[[([x - 1]/2).sup.2] + 2], y(0) = -1, y(1) = 0.

The exact solution for this problem is not available. Since the exact solution does not exist, the rate of convergence and maximum absolute error based on the double-mesh principle  has been given in Tables 7.7 and 7.8.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] denotes the value of [y.sub.j] for the mesh length h/[2.sup.k].

[FIGURE 7.1 OMITTED]

[FIGURE 7.2 OMITTED]

8. Discussions and conclusions. We have described a B-spline collocation method for the solution of self-adjoint singularly perturbed two-point boundary value problems. It is a practical method and easily can be implemented on a computer. The method has been analyzed for parameter-uniform convergence. For Example 7.3, we have computed the rate of convergence (see Table 7.8) which shows uniform second-order convergence as predicted by the theory.

Graphs have been plotted for Examples 7.1 and 7.2 for values of x [member of] [0, 1] versus the computed solution obtained at different values of x for a fixed [epsilon]. For each plot, we took N = 64 and [epsilon] = [2.sup.-20]. For Example 7.3, the exact solution is not available. We show two graphs, with and without fitting factor. For each graph, we use [epsilon] = [2.sup.-20] and N = 32, N = 64.

[FIGURE 7.3 OMITTED]

[FIGURE 7.4 OMITTED]

[FIGURE 7.5 OMITTED]

[FIGURE 7.6 OMITTED]

For Examples 7.1 and 7.2, it can be seen in Figures 7.1 and 7.3 that the exact and approximate solutions without using fitting factors deviate from each other in the boundary layer regions for smaller [epsilon]. To control these fluctuations, we used the fitting-factor technique and the resulting behavior of these solutions can be seen in Figures 7.2 and 7.4. Also, for Example 7.3, it can be easily seen that the deviation in the graph (Figure 7.6) of the approximate solution using a fitting factor corresponding to two values of N (32 and 64) is much smaller than the corresponding deviation in the graph (Figure 7.5) of the approximate solution without a fitting factor corresponding to these values of N.

We have replaced [epsilon] by [[sigma].sub.j] in the normalized form and not in the original self-adjoint problem, with a(x) [not equal to] constant, because in that case [epsilon] is a multiple of both the second and first derivative terms, which will cause ill-conditioning in the tridiagonal system. In normalized form, [epsilon] multiplies the second derivative term only. Hence, the fitting-factor technique for the normalized form can be implemented easily. This shows the importance of reducing the original self-adjoint problem to normal form.

Acknowledgments. The authors are thankful to the referee for his valuable suggestions and comments, which improved this paper.

* Received April 3, 2008. Accepted for publication August 13, 2008. Published online on December 8, 2008. Recommended by K. Burrage.

REFERENCES

 M. H. PROTTER AND H. F. WEINBERGER, Maximum Principles in Differential Equations, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1967.

 A. H. SCHATZ AND L. B. WAHLBIN, On the finite element method for singularly perturbed reaction diffusion problems in two and one dimension, Math. Comp., 40 (1983), pp. 47-89.

 I. P. BOGLAEV, A variational difference scheme for a boundary value problems with a small parameter in highest derivative, USSR Comput. Math. Math. Phys., 21 (1981), pp. 71-81.

 J. J. H. MILLER, On the convergence uniformaly in [epsilon], of difference schemes for a two point boundary singular perturbation problem in Numerical Analysis of Singular Perturbation Problems, P. W. Hemker and J. J. H. Miller, eds., Academic Press, New York, 1979, pp. 467-474.

 K. NIIJIMA, On a three-point difference scheme for a singular perturbation problem without a first derivative term I, II, Mem. Numer. Math. 7 (1980), pp. 1-27.

 M. K. KADALBAJOO AND V. K. AGGARWAL, Fitted mesh B-spline collocation method for solving self-adjoint singularly perturbed boundary value problems, Appl. Math. Comput., 161 (2005), pp. 973-987.

 P. A. FARRELL, R. E. O'RIORDAN, AND G. I. SHISHKIN, A class of singularly perturbed semilinear differential equations with interior layers, Math. Comp., 74 (2005), pp. 1759-1776.

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

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

 H. G. ROOS, Layer adaptive grids for singular perturbation problems, Z. Angew. Math. Mech., 78 (1998), pp. 291-309.

 M. STYNES AND H. G. ROOS, The midpoint upwind scheme, Appl. Numer. Math., 23 (1997), pp. 361-374.

 H. G. ROOS, Global uniformly convergent schemes for a singularly perturbed boundary value problem using patched base spline functions, J. Comput. Appl. Math., 29 (1990), pp. 69-77.

 K. SURLA AND V. JERKOVIC, Solving singularly perturbed boundary value problems by spline in tension, J. Comput. Appl. Math., 24 (1988), pp. 355-363.

 R. E. O'RIORDAN AND M. STYNES, An analysis of a super convergence result for a singularly perturbed boundary value problem, Math. Comp., 46 (1986), pp. 81-92.

 R. E. O'RIORDAN, Singularly perturbed finite element methods, Numer. Math., 44 (1984), pp. 425-434.

 R. E. O'RIORDAN, AND M. STYNES, A uniformly accurate finite element method for a singularly perturbed one-dimensional reaction-diffusion problem, Math. Comp., 47 (1986), pp. 555-570.

 M. STYNES AND R. E. O'RIORDAN, Uniformly accurate finite element method for a singularly perturbed problem in conservative form, SIAM J. Numer. Anal., 23 (1986), pp. 369-375.

 M. STOJANOVIC, Spline collocation method for singular perturbation problem, Glas. Mat. Ser., 37 (2002), pp. 393-403.

 P. M. PRENTER, Spline and Variational Methods, John Wiley & Sons, New York, 1975.

 E. P. DOOLAN, J. J. H. MILLER AND W. H. A. SCHILDERS, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980.

MOHAN K. KADALBAJOO ([dagger]) AND DEVENDRA KUMAR ([dagger])

([dagger]) Department of Mathematics & Statistics, Indian Institute of Technology, Kanpur, Kanpur-208016, India (dpanwar@iitk.ac.in).
```TABLE 4.1

B-Spline basis values.

Nodal values

[x.sub.i-2] [x.sub.i-1] [x.sub.i]

[B.sub.i] (x) 0 1 4
[B'.sub.i] (x) 0 3/h 0
[B".sub.i](x) 0 6/[h.sup.2] -12/[h.sup.2]

Nodal values

[x.sub.i+1] [x.sub.i+2]

[B.sub.i] (x) 1 0
[B'.sub.i] (x) -3/h 0
[B".sub.i](x) 6/[h.sup.2] 0

TABLE 7.1 Maximum absolute error for Example 7.1,
without fitting factor.

Number of mesh points N
[epsilon]
= [2.sup.-k] 16 32 64 128

k = 4 7.10E-03 1.78E-03 4.45E-04 1.11E-04
8 1.64E-02 3.80E-03 9.34E-04 2.32E-04
12 1.52E-01 6.35E-02 1.69E-02 3.92E-03
16 2.57E-01 2.28E-01 1.52E-01 6.35E-02
20 2.67E-01 2.65E-01 2.57E-01 2.28E-01
24 2.68E-01 2.68E-01 2.67E-01 2.65E-01
[E.sup.N] 2.68E-01 2.68E-01 2.67E-01 2.65E-01

Number of mesh points N

[epsilon]
= [2.sup.-k] 256 512 1024 2048

k = 4 2.78E-05 6.96E-06 1.74E-06 4.35E-07
8 5.80E-05 1.45E-05 3.63E-06 9.07E-07
12 9.63E-04 2.40E-04 5.99E-05 1.50E-05
16 1.69E-02 3.92E-03 9.64E-04 2.40E-04
20 1.52E-01 6.35E-02 1.69E-02 3.92E-03
24 2.57E-01 2.28E-01 1.52E-01 6.35E-02
[E.sup.N] 2.57E-01 2.28E-01 1.52E-01 6.35E-02

TABLE 7.2 Maximum absolute error for Example 7.1, with fitting
factor.

Number of mesh points N
[epsilon]
= [2.sup.-k] 16 32 64 128

k = 4 8.10E-03 2.03E-03 5.08E-04 1.27E-04
8 6.69E-03 1.62E-03 4.03E-04 1.01E-04
12 9.41E-03 1.88E-03 4.21E-04 1.02E-04
16 1.24E-02 2.91E-03 5.93E-04 1.18E-04
20 1.27E-02 3.18E-03 7.84E-04 1.82E-04
24 1.27E-02 3.20E-03 8.01E-04 2.00E-04
[E.sup.N] 1.27E-02 3.20E-03 8.01E-04 2.00E-04

Number of mesh points N
[epsilon]
= [2.sup.-k] 256 512 1024 2048

k = 4 3.18E-05 7.94E-06 1.99E-06 4.96E-07
8 2.51E-05 6.28E-06 1.57E-06 3.92E-07
12 2.52E-05 6.28E-06 1.57E-06 3.92E-07
16 2.63E-05 6.35E-06 1.57E-06 3.92E-07
20 3.71E-05 7.36E-06 1.64E-06 3.97E-07
24 4.90E-05 1.14E-05 2.32E-06 4.60E-07
[E.sup.N] 4.90E-05 1.14E-05 2.32E-06 4.60E-07

TABLE 7.3 Maximum absolute error for Example 7.2, without
fitting factor.

Number of mesh points N
[epsilon]
= [2.sup.-k] 16 32 64 128

k = 4 4.56E-02 1.13E-02 2.85E-03 7.12E-04
8 2.47E-01 6.31E-02 1.44E-02 3.53E-03
12 8.57E-01 5.10E-01 2.04E-01 5.37E-02
16 1.01E+00 8.71E-01 7.25E-01 4.70E-01
20 1.02E+00 9.07E-01 8.46E-01 7.96E-01
24 1.02E+00 9.09E-01 8.55E-01 8.27E-01
[E.sup.N] 1.02E+00 9.09E-01 8.55E-01 8.27E-01

Number of mesh points N
[epsilon]
= [2.sup.-k] 256 512 1024 2048

k = 4 1.78E-04 4.45E-05 1.11E-05 2.78E-06
8 8.78E-04 2.19E-04 5.48E-05 1.37E-05
12 1.24E-02 3.05E-03 7.59E-04 1.90E-04
16 1.94E-01 5.14E-02 1.19E-02 2.93E-03
20 6.94E-01 4.60E-01 1.91E-01 5.08E-02
24 8.08E-01 7.78E-01 6.86E-01 4.58E-01
[E.sup.N] 8.08E-01 7.78E-01 6.86E-01 4.58E-01

TABLE 7.4 Maximum absolute error for Example 7.2, with
fitting factor.

Number of mesh points N
[epsilon] =
[2.sup.-k] 16 32 64 128

k = 4 3.96E-02 1.01E-02 2.53E-03 6.34E-04
8 2.18E-02 6.93E-03 1.67E-03 4.16E-04
12 4.47E-02 8.62E-03 6.34E-03 2.54E-03
16 5.04E-02 1.91E-02 5.03E-03 9.73E-04
20 5.07E-02 1.96E-02 5.69E-03 1.49E-03
24 5.07E-02 1.96E-02 5.73E-03 1.52E-03
[E.sup.N] 5.07E-02 1.96E-02 6.34E-03 2.54E-03

Number of mesh points N
[epsilon] =
[2.sup.-k] 256 512 1024 2048

k = 4 1.59E-04 3.96E-05 9.91E-06 2.48E-06
8 1.04E-04 2.60E-05 6.49E-06 1.62E-06
12 6.13E-04 1.51E-04 3.79E-05 9.46E-06
16 2.09E-03 7.23E-04 1.75E-04 4.31E-05
20 3.45E-04 4.06E-04 5.60E-04 1.87E-04
24 3.89E-04 9.68E-05 2.25E-05 1.15E-04
[E.sup.N] 2.09E-03 7.23E-04 5.60E-04 1.87E-04

TABLE 7.5 Comparisons of maximum absolute errors for
Example 7.2 with those in .

Results in 

[epsilon]= [epsilon]=
N [epsilon]=1 [(1/N).sup.0.5] [(1/N).sup.1.0]

8 4.2E-01 3.8E-01 3.3E-01
16 1.1E-01 9.5E-02 7.8E-02
32 2.7E-02 2.3E-02 1.8E-02
64 6.9E-03 5.6E-03 4.2E-03
128 1.7E-03 1.3E-03 1.0E-03
256 4.3E-04 3.1E-04 2.5E-04
512 1.1E-04 7.4E-05 6.3E-05

Our results

[epsilon]= [epsilon]=
N [epsilon]=1 [(1/N).sup.0.5] [(1/N).sup.1.0]

8 2.0E-01 1.9E-01 1.7E-01
16 5.4E-02 4.8E-02 4.0E-02
32 1.4E-02 1.2E-02 8.9E-03
64 3.5E-03 2.8E-03 1.9E-03
128 8.6E-04 6.7E-04 3.8E-04
256 2.2E-04 1.6E-04 1.0E-04
512 5.4E-05 3.7E-05 4.1E-05

TABLE 7.6 Maximum absolute error for Example 7.3,
without fitting factor.

Number of mesh points N
[epsilon]
= [2.sup.-k] 32 64 128

k = 2 1.39E-03 3.47E-04 8.67E-05
4 5.34E-03 1.34E-03 3.34E-04
6 1.86E-02 4.65E-03 1.16E-03
8 4.79E-02 1.20E-02 2.99E-03
10 1.20E-01 2.99E-02 7.46E-03
12 4.17E-01 1.04E-01 2.60E-02
14 1.68E+00 4.08E-01 1.02E-01
16 7.57E+00 1.70E+00 4.10E-01
18 2.39E+01 7.65E+00 1.71E+00
20 4.94E+01 2.42E+01 7.70E+00
22 6.91E+01 5.01E+01 2.43E+01
24 7.72E+01 7.01E+01 5.04E+01
26 7.95E+01 7.83E+01 7.06E+01
28 8.02E+01 8.08E+01 7.89E+01
[E.sup.N] 8.02E+01 8.08E+01 7.89E+01

Number of mesh points N
[epsilon]
= [2.sup.-k] 256 512 1024

k = 2 2.17E-05 5.42E-06 1.35E-06
4 8.35E-05 2.09E-05 5.22E-06
6 2.91E-04 7.27E-05 1.82E-05
8 7.48E-04 1.87E-04 4.68E-05
10 1.86E-03 4.66E-04 1.17E-04
12 6.49E-03 1.62E-03 4.06E-04
14 2.54E-02 6.34E-03 1.59E-03
16 1.02E-01 2.55E-02 6.38E-03
18 4.13E-01 1.03E-01 2.57E-02
20 1.72E+00 4.15E-01 1.04E-01
22 7.73E+00 1.73E+00 4.17E-01
24 2.44E+01 7.75E+00 1.73E+00
26 5.06E+01 2.45E+01 7.76E+00
28 7.08E+01 5.07E+01 2.45E+01
[E.sup.N] 7.08E+01 5.07E+01 2.45E+01

TABLE 7.7 Maximum absolute error for Example 7.3,
with fitting factor.

Number of mesh points N
[epsilon] =
[2.sup.-k] 32 64 128

k = 2 1.31E-03 3.28E-04 8.21E-05
4 4.93E-03 1.23E-03 3.08E-04
6 1.60E-02 4.00E-03 1.00E-03
8 3.71E-02 9.27E-03 2.32E-03
10 6.19E-02 1.54E-02 3.86E-03
12 9.39E-02 2.34E-02 5.83E-03
14 1.34E-01 3.29E-02 8.15E-03
16 1.90E-01 4.31E-02 1.05E-02
18 2.84E-01 7.89E-02 2.56E-02
20 4.08E-01 7.82E-02 5.26E-02
22 4.60E-01 1.06E-01 2.78E-02
24 4.62E-01 1.21E-01 2.78E-02
26 4.62E-01 1.21E-01 3.08E-02
28 4.62E-01 1.21E-01 3.10E-02
[E.sup.N] 4.62E-01 1.21E-01 5.26E-02

Number of mesh points N
[epsilon] =
[2.sup.-k] 256 512 1024

k = 2 2.05E-05 5.13E-06 1.28E-06
4 7.71E-05 1.93E-05 4.82E-06
6 2.50E-04 6.26E-05 1.56E-05
8 5.79E-04 1.45E-04 3.62E-05
10 9.65E-04 2.41E-04 6.03E-05
12 1.46E-03 3.64E-04 9.10E-05
14 2.03E-03 5.08E-04 1.27E-04
16 2.60E-03 6.50E-04 1.62E-04
18 6.20E-03 1.54E-03 3.84E-04
20 1.49E-02 3.61E-03 9.05E-04
22 2.97E-02 7.97E-03 1.94E-03
24 1.91E-02 1.57E-02 4.12E-03
26 7.09E-03 1.09E-02 8.09E-03
28 7.80E-03 1.79E-03 5.78E-03
[E.sup.N] 2.97E-02 1.57E-02 8.09E-03

TABLE 7.8 Rate of convergence for Example 7.3.

[epsilon] =
[2.sup.-k] r(0) r(1) r(2) r(3) r(4)

k = 2 2.00E+00 2.00E+00 2.00E+00 2.00E+00 2.00E+00
4 2.00E+00 2.00E+00 2.00E+00 2.00E+00 2.00E+00
6 2.00E+00 2.00E+00 2.00E+00 2.00E+00 2.00E+00
8 2.00E+00 2.00E+00 2.00E+00 2.00E+00 2.00E+00
10 2.01E+00 2.00E+00 2.00E+00 2.00E+00 2.00E+00
12 2.00E+00 2.00E+00 2.00E+00 2.00E+00 2.00E+00
14 2.03E+00 2.01E+00 2.01E+00 2.00E+00 2.00E+00
16 2.14E+00 2.04E+00 2.01E+00 2.00E+00 2.00E+00
```
COPYRIGHT 2008 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.