MATEMATIQKI VESNIK

59 (2007), 31–46

UDK 517.518.5 originalni nauqni rad

research paper

SPLINE-WAVELET SOLUTION OF SINGULARLY PERTURBED BOUNDARY PROBLEM

Desanka Radunović

Abstract. Boundary or interior layers typically appear in singularly perturbed boundary problems. Solution gradients are very sharp in these regions, and it seems natural to use wavelets to obtain numerical solution. As layers are usually positioned on boundaries, wavelets have to be modified appropriately. In this paper the collocation method based on spline-wavelets is derived and tested on simple linear one-dimensional singularly perturbed boundary problem.

1. Introduction

Singularly perturbed boundary problems appear in many areas of fluid and gas dynamics (transition from Navier-Stokes to Euler equations in fluid flow models, semiconductor device simulation, etc). Boundary layers may arise even in very simple flows modelled by one-dimensional linear singularly perturbed boundary problem

\begin{matrix}L u(x) ≡ ε u′′(x) + p(x) u′(x) & + q(x) u(x) = f(x), x ∈ [a, b] \\ u(a) = u_{0}, & u(b) = u_{1}, \end{matrix}  (1)

which will be analyzed. Coefficients p(x), q(x) and f(x) are smooth functions. When perturbation parameter ε, 0 < ε ≤ 1, is very small, close to zero, narrow regions of sharp gradients of the solution cause difficulties in numerical solving of such problems. Successful numerical techniques are those which exhibit so called ε-uniform convergence [3].

Definition 1. A numerical method is said to be ε-uniform of order p on the mesh Ω_{m} = {x_{i}, i = 0, …, m} if there exists a number m_{0}, independent of ε, such that for all m ≥ m_{0}

sup_{0<ε≤1} max_{Ω_{m}} |u(x) − u_{m}(x)| ≤ C m^{−p},  (2)

where u(x) is the solution of the differential equation, u_{m}(x) is the numerical approximation to u, and C and p are constants independent of ε and m.

AMS Subject Classification: 65L10, 65L60, 65T60

Keywords and phrases: Boundary layers, spline wavelets, collocation.

In the construction of ε-uniform methods two approaches have generally been taken [4]. The first approach is based on fitted finite difference operators applied on standard meshes. The second successful approach refers to fitted mesh methods, which apply standard difference operators on piecewise uniform fitted meshes [5]. This multiresolution idea leads naturally to wavelets, as wavelets are just oscillatory functions with compact support defined on fitted meshes.

2. Multiresolution

Let us clarify the last sentence.

Definition 2. Multiresolution analysis is the decomposition of the Hilbert space L_{2}(R) to the nested sequence of closed subspaces {V_{j}}_{j∈Z}, which satisfy the following relations [2]

⋯ ⊂ V_{2} ⊂ V_{1} ⊂ V_{0} ⊂ V_{−1} ⊂ V_{−2} ⊂ …

⋂_{j∈Z} V_{j} = \{0\}, \overline{⋃_{j∈Z} V_{j}} = L_{2}(R)

∀f ∈ L_{2}(R) and ∀j ∈ Z, f(x) ∈ V_{j} ⟺ f(2x) ∈ V_{j−1}
∀f ∈ L_{2}(R) and ∀k ∈ Z, f(x) ∈ V_{0} ⟺ f(x − k) ∈ V_{0}
∃φ ∈ V_{0} so that {φ(x − k)}_{k∈Z} is the Riesz basis in V_{0}.

Dilatations and translations of the function φ(x), so called scaling function, define bases of all subspaces V_{j},

{φ_{j,k}(x)}_{k∈Z}, where φ_{j,k}(x) = φ(2^{−j}x − k).

Let us denote by W_{j} the orthogonal complement of the subspace V_{j} in the space V_{j−1},

V_{j−1} = V_{j} ⊕ W_{j},  j ∈ Z.  (3)

For every k ≠ j it is W_{k} ⊥ W_{j}, which is the consequence of the relation between spaces V_{j} and definition (3).

Bases of all subspaces W_{j} are defined by dilatations and translations of one function ψ(x), so called wavelet,

{ψ_{j,k}(x)}_{k∈Z}, where ψ_{j,k}(x) = ψ(2^{−j}x − k).

It follows, from (3) for j = 0, that the scaling function φ(x) (“father wavelet”) is the solution of the dilatation equation

φ(x) = ∑_{k=0}^{N} c(k) φ(2x − k),  (4)

and that the wavelet ψ(x) (“mother wavelet”) can be expressed by the wavelet equation

ψ(x) = ∑_{k} d(k) φ(2x − k).  (5)

By recursion of the relation (3), we obtain following orthogonal decomposition of the space V_{j−1}

V_{j−1} = W_{j} ⊕ ⋯ ⊕ W_{J−1} ⊕ W_{J} ⊕ V_{J},  j < J.  (6)

So, the approximation f_{j−1}(x) ∈ V_{j−1} of a function f(x) ∈ L_{2}(R) can be expressed in a multiresolution form

f_{j−1}(x) = ∑_{k∈Z} a_{J,k} φ_{J,k}(x) + ∑_{l=j}^{J} ∑_{k∈Z} b_{l,k} ψ_{l,k}(x),  (7)

i.e. as a sum of the coarse approximation (first sum) and details on various levels of resolution (double sum). When j → −∞ the expression (7) is equal to f(x) in L_{2} sense, as, according to the definition 2, V_{j−1} → L_{2} in (6).

3. Spline-wavelets

If coefficients c(k) in the dilatation equation (4) are chosen to be binomial coefficients of the order N, the scaling function φ^{(N)}(x) is B-spline of the order (N − 1), defined on the integer division Δ = {x_{k} = k, k ∈ Z}. The compact support of B-spline φ^{(N)}(x) is the interval [0, N].

Theorem 1. [6] The following recursion is valid

φ_{k}^{(N)}(x) = \frac{x − k}{N − 1} φ_{k}^{(N−1)}(x) + \frac{k + N − x}{N − 1} φ_{k+1}^{(N−1)}(x), N = 2, 3, …,
φ_{k}^{(1)}(x) = ℵ_{[k, k+1)}(x) = \{ \begin{matrix}1, & k ≤ x < k + 1 \\ 0, & otherwise, \end{matrix}

(8)

where φ_{k}^{(N)}(x) ≡ φ^{(N)}(x − k).

Wavelet bases are defined on an infinite domain. To derive spline-wavelet approximation of the solution of the problem (1), which is defined on the interval, we have to modify scaling functions and wavelets on the boundaries of the interval. We shall call a boundary spline every spline whose compact support do not whole belong to the observed interval. Without lost of generality, we shall define interior and boundary splines on the integer division {x_{k} = k, k = 0, …, M} of the interval [0, M].

Theorem 2. [1] Quadratic boundary splines, which form spline basis on [0, M] together with interior quadratic B-splines ((8), for N = 3), are

\begin{aligned}φ_{−2}^{(3)}(x) & = (1 − x)^{2} ℵ_{[0, 1)}(x), \\ φ_{−1}^{(3)}(x) & = (2x − \frac{3}{2}x^{2}) ℵ_{[0, 1)}(x) + \frac{1}{2}(2 − x)^{2} ℵ_{[1, 2)}(x) \\ φ_{M−2}^{(3)}(x) & = φ_{−1}^{(3)}(M − x),  φ_{M−1}^{(3)}(x) = φ_{−2}^{(3)}(M − x). \end{aligned}  (9)

Fig. 1. Quadratic splines

Suppose that the length of the interval is even integer, M = 2m. For shortness, let us denote the scaling function on the finer resolution level by u_{k}(x) ≡ φ(2x−k+N), and the scaling function on the coarser resolution level by v_{k}(x) ≡ φ(x − k + N). If we express interior quadratic splines (8) and boundary splines (9) in the form of the dilatation equation (4),

v_{k}(x) = ∑_{l=1}^{2m+N−1} p(k, l) u_{l}(x), k = 1, …, m + N − 1,  (10)

for N = 3, coefficients p(k, l) are given, up to the scaling factor, by nonzero elements of row vectors

p(1, 1 : 2) = ( 4 2 )
p(2, 2 : 4) = ( 2 3 1 )
p(k, 2k − 3 : 2k) = ( 1 3 3 1 ) ,
p(m + 1, 2m − 1 : 2m + 1) = ( 1 3 2 )
p(m + 2, 2m + 1 : 2m + 2) = ( 2 4 )

Transformation matrix P = {p(k, l)} is a rectangular matrix whose dimension is (m + 2) × (2m + 2).

Analogous results are derived for cubic splines.

Theorem 3. [1] Cubic boundary splines, which form spline basis on [0, M] together with interior cubic B-splines ((8), for N = 4), are

\begin{aligned}φ_{−3}^{(4)}(x) & = (1 − x)^{3} ℵ_{[0, 1)}(x), \\ φ_{−2}^{(4)}(x) & = (3x − \frac{9}{2}x^{2} + \frac{7}{4}x^{3}) ℵ_{[0, 1)}(x) + \frac{1}{4}(2 − x)^{3} ℵ_{[1, 2)}(x) \\ φ_{−1}^{(4)}(x) & = (\frac{3}{2}x^{2} − \frac{11}{12}x^{3}) ℵ_{[0, 1)}(x) + (\frac{7}{12} + \frac{1}{4}(x − 1) − \frac{5}{4}(x − 1)^{2} \\ & + \frac{7}{12}(x − 1)^{3}) ℵ_{[1, 2)}(x) + \frac{1}{6}(3 − x)^{3} ℵ_{[2, 3)}(x) \\ φ_{M−3}^{(4)}(x) & = φ_{−1}^{(4)}(M − x), \\ φ_{M−2}^{(4)}(x) & = φ_{−2}^{(4)}(M − x),  φ_{M−1}^{(4)}(x) = φ_{−3}^{(4)}(M − x). \end{aligned}  (11)

Fig. 2. Qubic splines

Dimension of the transformation matrix P, defined by coefficients p(k, l) in (10) for N = 4, is (m + 3) × (2m + 3). Nonzero entries of its row vectors, up to the scaling factor, are

p(1, 1 : 2) = ( 8 4 )
p(2, 2 : 4) = ( 4 6 \frac{3}{2} )
p(3, 3 : 6) = ( 2 \frac{11}{2} 4 1 )
p(k, 2k − 4 : 2k) = ( 1 4 6 4 1 ) ,

p(m + 1, 2m − 2 : 2m + 1) = ( 1 4 \frac{11}{2} 2 )
p(m + 2, 2m : 2m + 2) = ( \frac{3}{2} 6 4 )
p(m + 3, 2m + 2 : 2m + 3) = ( 4 8 )

Now, we shall construct semiorthogonal wavelets ψ_{jk}^{(N)}(x) for B-spline scaling function φ^{(N)}(x). Semiorthogonal means that wavelets on the same scale are not orthogonal, but they are orthogonal to the scaling functions, and also wavelets on various scales (for various j) are orthogonal between themselves.

Theorem 4. [1] Compactly supported semiorthogonal spline-wavelet ψ^{(N)}(x) is given by the wavelet equation (5), where scaling function is B-spline φ^{(N)}(x) ≡ φ_{0}^{(N)}(x) given in (8), and coefficients d(k) are equal to

d(k) = \frac{(−1)^{k}}{2^{N−1}} ∑_{l=0}^{N} \binom{N}{l}φ^{(2N)}(k − l + 1), k = 0, …, 3N − 2.  (12)

If we denote, as before, wavelets on the coarser resolution level by w_{k}(x) ≡ ψ(x − k), the wavelet equation has the form similar to (10),

w_{k}(x) = ∑_{l=1}^{2m+N−1} q(k, l) u_{l}(x), k = 1, …, m.  (13)

The coefficients q(k, l) for boundary wavelets (k = 1, …, N − 1, and k = m − N + 2, …, m,) are calculated from the request that wavelets have to be orthogonal to the scaling functions on the same resolution level, (w_{k}(x), v_{l}(x)) = 0 ((·, ·) denotes scalar product). These coefficients for interior wavelets are calculated from (12).

Theorem 5. Quadratic spline-wavelets, which form wavelet basis on [0, M] together with quadratic B-spline scaling functions (8) and (9), are given by wavelet equation (13) for N = 3. Wavelet coefficients q(k, l) are given, up to the scaling factor, by nonzero elements of row vectors

q(1, 1 : 6) = ( \frac{4576}{15} \frac{−5564}{15} \frac{9934}{45} \frac{−4166}{45} \frac{812}{45} \frac{−28}{45} ) q(2, 2 : 8) = ( \frac{−338}{5} \frac{29222}{173} \frac{−29867}{99} \frac{21853}{75} \frac{−7432}{53} \frac{4147}{150} \frac{−143}{150} ) q(k, 2k − 3 : 2k + 4) = ( 1 −29 147 −303 303 −147 29 −1 ) q(m − 1, 2m − 5 : 2m + 1) = ( \frac{143}{150} \frac{−4147}{150} \frac{7432}{53} \frac{−21853}{75} \frac{29867}{99} \frac{−29222}{173} \frac{338}{5} ) q(m, 2m − 3 : 2m + 2) = ( \frac{28}{45} \frac{−812}{45} \frac{4166}{45} \frac{−9934}{45} \frac{5564}{15} \frac{−4576}{15} )

Q = {q(k, l)} is the rectangular matrix with dimension m × (2m + 2).

Theorem 6. Cubic spline-wavelets, which form wavelet basis on [0, M] together with cubic B-spline scaling functions (8) and (11), are given by wavelet equation (13) for N = 4. Wavelet coefficients q(k, l) are given, up to the scaling factor, by nonzero entries of row vectors

q(1, 1 : 8) = ( 22642 \frac{−260258}{9} \frac{91274}{5} \frac{−547933}{84} \frac{66042}{19} \frac{−45387}{59} \frac{4387}{77} \frac{−1485}{3232} )

Fig. 3. Quadratic spline-wavelets

q(2, 2 : 10) =

( \frac{−64591}{17} \frac{113671}{10} \frac{−185221}{11} \frac{45219}{2} \frac{−68345}{4} \frac{167457}{23} \frac{−12349}{8} \frac{16093}{141} \frac{−833}{905} )

q(3, 3 : 12) =

( \frac{4881}{427} \frac{−7289}{77} \frac{81806}{67} −15124 \frac{417353}{15} \frac{−471341}{20} \frac{166557}{16} \frac{−11114}{5} \frac{23509}{143} \frac{−1876}{1415} )

q(k, 2k − 4 : 2k + 6) =

(1 − 124 1677 − 7904 18482 − 24264 18482 − 7904 1677 − 124 1) q(m − 2, 2m − 8 : 2m + 1) =

( \frac{1876}{1415} \frac{−23509}{143} \frac{11114}{5} \frac{−166557}{16} \frac{471341}{20} \frac{−417353}{15} 15124 \frac{−81806}{67} \frac{7289}{77} \frac{−4881}{427} )

q(m − 1, 2m − 6 : 2m + 2) =

( \frac{833}{905} \frac{−16093}{141} \frac{12349}{8} \frac{−167457}{23} \frac{68345}{4} \frac{−45219}{2} \frac{185221}{11} \frac{−113671}{10} \frac{64591}{17} )

q(m, 2m − 4 : 2m + 3) =

( \frac{1485}{3232} \frac{−4387}{77} \frac{45387}{59} \frac{−66042}{19} \frac{547933}{84} \frac{−91274}{5} \frac{260258}{9} −22642 )

Dimension of the transformation matrix Q = {q_{k,l}} is m × (2m + 3).

Previous results can be expressed as discrete spline-wavelet transformation by matrix relation

( \begin{matrix}v \\ w \end{matrix} ) = ( \begin{matrix}P \\ Q \end{matrix} ) u.  (14)

Vector u = {u_{k}(x)} is (2m + N − 1)-dimensional vector of approximations on the finer level, vector v = {v_{k}(x)} is (m+N −1)-dimensional vector of approximations

Fig. 4. Qubic spline-wavelets

on the coarser level and w = {w_{k}(x)} is m-dimensional vector of details on the coarser level. P is low-frequency matrix with dimension (m+N −1)×(2m+N −1), and Q is high-frequency matrix with dimension m × (2m + N − 1).

4. Collocation based on spline-wavelets

The linear one-dimensional singularly perturbed boundary problem (1), with boundary or interior layers, is solved by the collocation method based on B-spline or spline-wavelet bases.

B-splines, quadratic or cubic, are defined on the uniform dyadic division of the interval [a, b], i.e. mesh step is h = (b − a)/2^{m}. By linear transformation of the interval [a, b] to the interval [0, 2^{m}], all calculations are done on integer divisions, with the finest resolution level division

Δ = {x_{k} = k, k = 0, …, 2^{m}}.

The approximate solution of the transformed problem (1) is searched in a spline form (only first sum in (7)),

u_{s}(x) = ∑_{k=1−N}^{L_{0}−1} a_{0,k} φ_{0,k}(x)  (15)

or in a spline-wavelet form (expression (7)),

u_{w}(x) = ∑_{k=1−N}^{L_{J} −1} a_{J,k} φ_{J,k}(x) + ∑_{j=1}^{J} ∑_{k=1−N}^{L_{j}−N} b_{j,k} ψ_{j,k}(x),  (16)

where L_{j} = 2^{m−j}, j = 0, …, J. The number of levels J is limited by the request that on every level it must exist at least one interior spline or wavelet. The total number L of free parameters a_{J,k}, b_{j,k} in both options is equal. It depends only on the number m in division Δ and the order of the spline N, and does not depend on the number of resolution levels J,

\begin{aligned}L & = (L_{J} + N − 1) + ∑_{j=1}^{J} L_{j} = (L_{J} + N − 1) + L_{J} ∑_{j=0}^{J−1} 2^{j} \\ & = (L_{J} + N − 1) + L_{J}(2^{J} − 1) = 2^{m} + N − 1. \end{aligned}  (17)

Free parameters are determined by the collocation method. The error function, in the case of spline approximation (15), is equal to

R_{s}(x; a_{0}) = −f_{T} (x)
+ ∑_{k=1−N}^{L_{0}−1} a_{0,k} (ε c_{0}^{2} φ′′_{0,k}(x) + p_{T} (x) c_{0} φ′_{0,k}(x) + q_{T} (x) φ_{0,k}(x)) ,(18)

and, in the case of spline-wavelet approximation (16), it is equal to

\begin{aligned}R_{w}(x; b_{1} & , …, b_{J}, a_{J}) = −f_{T} (x) \\ & + ∑_{k=1−N}^{L_{J} −1} a_{J,k}(ε c_{J}^{2} φ′′_{J,k}(x) + p_{T} (x) c_{J} φ′_{J,k}(x) + q_{T} (x) φ_{J,k}(x)) \\ & + ∑_{j=1}^{J} ∑_{k=1−N}^{L_{j}−N} b_{j,k}(ε c_{j}^{2} ψ′′_{j,k}(x) + p_{T} (x) c_{j} ψ′_{j,k}(x) + q_{T} (x) ψ_{j,k}(x)), \end{aligned} (19)

where c_{j} = L_{j}/(b − a), b_{j} = {b_{j,k}}_{k}, a_{J} = {a_{J,k}}_{k}, and the subscript T denotes transformed coefficients of the problem (1). Derivatives of splines are calculated by use of well known relation for the l-th derivative of the B-spline of order (N − 1) defined on the integer division [6],

(φ^{(N)}(x))^{(l)} = ∑_{k=0}^{l}(−1)^{k}\binom{l}{k}φ^{(N−l)}(x − k).

The L-dimensional system of linear equations per unknown vector a_{0} in (18), or vectors b_{j}, j = 1, …, J, and a_{J} in (19), is obtained from requirements that the solution satisfies two boundary conditions (1), and that the error function is zero in collocation points x_{l}, l = 1, …, L − 2. For the spline problem, relations (15) and (18) give the system

\begin{aligned}a_{0,1−N} & = u_{0} \\ ∑_{k=1−N}^{L_{0}−1} a_{0,k}(ε c_{0}^{2} φ′′_{0,k}(x_{l}) & + p_{T} (x_{l}) c_{0} φ′_{0,k}(x_{l}) + q_{T} (x_{l}) φ_{0,k}(x_{l})) \\ & = f_{T} (x_{l}),  l = 1, …, L − 2, \\ a_{0,L_{0}−1} & = u_{1} \end{aligned}  (20)

and for the spline-wavelet problem, from (16) and (19) we obtain the system

\begin{aligned}a_{J,1−N} + ∑_{j=1}^{J} b_{j,1−N} & = u_{0} \\ ∑_{k=1−N}^{L_{J} −1} a_{J,k}(ε c_{J}^{2} φ′′_{J,k}(x_{l}) & + p_{T} (x_{l}) c_{J} φ′_{J,k}(x_{l}) + q_{T} (x_{l}) φ_{J,k}(x_{l})) \\ + ∑_{j=1}^{J} ∑_{k=1−N}^{L_{j}−N} b_{j,k}(ε c_{j}^{2} ψ′′_{j,k}(x_{l}) & + p_{T} (x_{l}) c_{j} ψ′_{j,k}(x_{l}) + q_{T} (x_{l}) ψ_{j,k}(x_{l})) \\ & = f_{T} (x_{l}),  l = 1, …, L − 2. \\ a_{J,L_{J} −1} + ∑_{j=1}^{J} b_{j,L_{j}−N} & = u_{1} \end{aligned}  (21)

Form of the first and the last equations in systems (20) and (21) follows from the fact that only the first and the last splines and wavelets are not equal to zero on the boundaries. We suppose that these functions are normalized, so that they are equal to one in relevant boundary point.

The number of collocation points L − 2 = 2^{m} + N − 3 depends on the order of the spline. For the quadratic spline (N = 3) this number is equal to the number of intervals, and for the cubic spline (N = 4) it is equal to the number of nodes of the division Δ.

5. Numerical results

In this section we shall present a few examples, whose approximate solutions are computed both in spline (15) and spline-wavelet form (16) (for N = 3 and N = 4), for various numbers of resolution levels and for different distributions of collocation points. Examples are chosen to be simple, with known solutions, and, for small values of the parameter ε, they represent convection-diffusion and reaction-diffusion problems with boundary and interior layers.

Example 1. The convection-diffusion boundary problem

ε u′′(x) + u′(x) = 0,  u(0) = 0, u(1) = 1,  (22)

has the solution

u(x) = \frac{1 − e^{−x/ε}}{1 − e^{−1/ε}} ,

with the boundary layer at the left boundary, as

lim_{ε→0} lim_{x→0} u(x) ≠ lim_{x→0} lim_{ε→0} u(x).

Example 2. The convection-diffusion boundary problem

ε u′′(x) + 2x u′(x) = 0,  u(−1) = −1, u(1) = 2,  (23)

has the interior layer in the turning point x = 0, as the solution tends to

lim_{ε→0} u(x) = \{ \begin{matrix}−1, & x < 0 \\ 2, & x > 0 \end{matrix} .

Example 3. The reaction-diffusion boundary problem

ε u′′(x) − u(x) = 0,  u(0) = u(1) = 1,  (24)

has the solution with boundary layers at both boundaries,

u(x) = \frac{e^{−x/\sqrt{ε}} + e^{−(1−x)/\sqrt{ε}}}{1 + e^{−1/\sqrt{ε}}} .

As it was expected, numerical solutions obtained by splines and by splinewavelets of the same type, for the same division of the interval, are identical. In both cases calculations are numerically stable, and condition numbers of linear systems (20) and (21) (Tab. 1) do not grow with decreasing of the small parameter or increasing of the dimension of the system for the optimal choice of collocation points (which will be elaborated later).

ε\h  2^{−2}  2^{−3}  2^{−4}  2^{−5}  2^{−6}^{2^{−6} 1.06e + 002 2.47e + 002 7.59e + 002 2.61e + 003 1.02e + 004}^{2^{−7} 3.54e + 002 2.54e + 002 6.21e + 002 1.84e + 003 5.85e + 003}^{2^{−8} 1.54e + 003 9.33e + 002 6.46e + 002 1.65e + 003 4.79e + 003}^{2^{−9} 6.29e + 003 4.31e + 003 2.44e + 003 1.71e + 003 4.52e + 003}

Table 1. Condition numbers

The optimal choice of collocation points are division points in cubic case and middle points of division intervals in quadratic case. In both cases the matrix of the system (20) is a three-diagonal matrix (Fig. 5(a)). The matrix of the system (21) is also sparse with multi-diagonal blocks, where the number of blocks is equal to the number of resolution levels (Fig. 5(b)). The last block corresponds to spline coefficients on the coarsest level, and the others correspond to wavelet coefficients on all resolution levels. Fig. 5(c) and Fig. 5(d) represent one resolution level matrices of the system (21) in the quadratic and cubic case, respectively.

Although the system (21) is more complex, the advantage of this approach is the possibility of data compression, as a lot of components of the solution vector are negligible. Fig. 6(a) shows magnitude of coefficients in the wavelet representation (16) for the convection problem (22), while Fig. 6(b) shows magnitude of coefficients in the spline representation (15) for the same example.

A few words about the optimal choice of the uniform distribution of collocation points, which gives the best accuracy. In the cubic case, the best result is obtained when collocation points are the division points, so that the ends of the interval are

Fig. 5 Matrices assigned to spline and spline-wavelet bases

also collocation points. Every other choice takes with less weight the first boundary basis function, the only one which is different from zero in the relevant end point. Figure 7 represents the numerical solution of the example 3 for ε = 2^{−10}, obtained by use of cubic spline-wavelets and 65 collocation points, which (a) include ends,

(b) are far from ends at least for half of the division interval.

In the quadratic case, on the contrary, the best choice is when collocation points are middle points of division intervals. In this case, the ends are not collocation points. The instability grows when we approach with collocation points to division points, as the second derivative of the quadratic spline is discontinuous in these points. The worst option is when ends are collocation points, when the system of linear equations (20) or (21) becomes very ill-conditioned, and the calculation error grows.

Reaction-diffusion problems manifest similar behavior, but they are not so sensitive to the position of collocation points as convection-diffusion problems are. The reason is that the width of the boundary layer is of the order \sqrt{ε} for these problems, while it is of the order ε in the case of convection-diffusion problems.

Fig. 6. Magnitude of wavelet and spline coefficients

Fig. 7. The effect of the collocation points distribution in cubic case

Fig. 8 shows the influence of the collocation points distribution to the numerical solution of the reaction problem (24) for ε = 10^{−4}. Quadratic spline-wavelets and 32 collocation points are used, when (a) ends are not collocation points, and (b) ends are collocation points.

Even when we make optimal choice of collocation points, the numerical solution will be poor, specially in the layer, if the number of collocation points is not adequate. It means, the less is parameter ε more collocations points, i.e. more basis functions, we need.

Fig. 9 illustrates a behavior of the numerical solution of the problem (22), obtained by cubic spline-wavelets. First two figures, (a) and (b), show the solutions for ε = 10^{−2} and ε = 10^{−3} for the same mesh parameter h = 2^{−6}. Other two figures,

Fig. 8. The effect of the collocation points distribution in quadratic case

Fig. 9. The influence of small parameter and mesh size

(c) and (d), illustrate the dependance of the solution from the grid parameter, (h = 2^{−7} and h = 2^{−8}), for the same small value of parameter ε = 10^{−3}. Reactiondiffusion problems are not so sensitive to the values of small parameter, i.e. the same accuracy can be obtained with the greater mesh parameter h for the same value of the small parameter ε.

Numerical estimates give that the order of accuracy is about two for both kinds of splines and both kinds of problems, although numerical tests show that uniform errors are less when we use cubic spline-wavelets. Tables 2 and 3 show that the method is not ε-uniform convergent, in the sense of Definition 1.

\begin{matrix}ε\h & 2^{−3} & 2^{−4} & 2^{−5} & 2^{−6} & 2^{−7} & 2^{−8} \\ 1 & 2e − 4 & 4e − 5 & 1e − 5 & 2e − 6 & 6e − 7 & 2e − 7 \\ 0.1 & 6e − 2 & 1e − 2 & 3e − 3 & 7e − 4 & 2e − 4 & 5e − 5 \\ 0.01 & 9e − 1 & 1e + 0 & 2e + 0 & 1e − 1 & 2e − 2 & 5e − 3 \\ 0.001 & 8e + 0 & 2e + 0 & 9e − 1 & 6e − 1 & 6e − 1 & 2e − 1 \\ 0.0001 & 8e + 1 & 2e + 1 & 5e + 0 & 1e + 0 & 7e − 1 & 6e − 1 \\ 0.00001 & 8e + 2 & 2e + 2 & 5e + 1 & 1e + 1 & 3e + 0 & 1e + 0 \end{matrix}

Table 2. Uniform errors for problem (22) in cubic case

\begin{matrix}ε\h & 2^{−3} & 2^{−4} & 2^{−5} & 2^{−6} & 2^{−7} & 2^{−8} \\ 1 & 6e − 4 & 2e − 4 & 4e − 5 & 1e − 5 & 2e − 6 & 5e − 7 \\ 0.1 & 1e − 1 & 3e − 2 & 8e − 3 & 2e − 3 & 5e − 4 & 1e − 4 \\ 0.01 & 3e + 0 & 1e + 0 & 4e − 1 & 1e − 1 & 3e − 2 & 7e − 3 \\ 0.001 & 1e + 2 & 2e + 1 & 7e + 0 & 3e + 0 & 5e − 1 & 4e − 2 \\ 0.0001 & 1e + 4 & 2e + 3 & 2e + 2 & 4e + 1 & 2e + 1 & 7e + 0 \\ 0.00001 & 1e + 6 & 1e + 5 & 2e + 4 & 3e + 3 & 4e + 2 & 1e + 2 \end{matrix}

Table 3. Uniform errors for problem (22) in quadratic case

6. Conclusion

We can see that, used for solving singularly perturbed boundary problems, spline and spline-wavelet approximations of the same type give identical results. The advantage of splines is a simpler system of linear equations (three-diagonal), and the advantage of spline-wavelets is a significant possibility of data compression.

The accuracy for fixed mesh depends on the type of the problem (convective or reactive), on the order of spline (quadratic or cubic) and on the distribution of collocation points. Convective problems are more sensitive to the value of small parameter, so we need finer mesh to get the same accuracy as in the case of reactive problems. Better results are obtained by cubic splines for both types of problems. Best results in the case of convective problems are obtained when collocation points are mesh points, while, in the case of reactive problems, best results are obtained when collocation points are middle interval points.

Numerical order of convergence for both types of splines is about two, but scheme is not ε-uniform (Definition 1). A future work will go towards construction of ε-uniform schemes by use of different types of wavelets and their various modifications on boundaries.

REFERENCES

[1] Chui C., Wavelets: A Mathematical Tool for Signal Analysis, SIAM, Philadelphia, 1997.

[2] Daubechies I., Ten Lectures on Wavelets, SIAM, Philadelphia, 1992.

[3] Farrell P., Miller J., O’Riordan E., Shishkin G., A uniformly convergent finite difference scheme for a singularly perturbed semilinear equation, SIAM J. Numer. Anal., 33, 1135–1140, 1996.

[4] Roos H.-G., Stynes M., Tobiska L., Mumerical Methods for Singularly Perturbed Differential Equations, Springer, Berlin, 1996.

[5] Shishkin G.I., Discrete approximation of singularly perturbed elliptic and parabolic equations, Russian Academy of Sciences, Ural Section, Ekaterinburg, 1992.

[6] Zavjalov Y., Kvasov I., Miroshnicenko L., Spline Function Methods, Moscow, 1980 (in Russian).

(received 19.02.2006)

Faculty of Mathematics, University of Belgrade, Studentski trg 16, Belgrade, Serbia

E-mail: dradun@matf.bg.ac.yu