Haar wavelet method for solving stiff differential equations.

1 Introduction

In the past two decades, stiff differential equations have been studied extensively and various methods for their solutions have been proposed. While the intuitive meaning of stiffness is clear to all specialists, a unique mathematical definition is lacking.

As a rule stiffness occurs in differential equations where there are two or more different time scales of the independent variables on which the dependent variables are changing. In the case of linear problems the stiffness is caused by eigenvalues of big negative values. To measure the degree of stiffness the following stiffness ratio SR = max [absolute value of [lambda]]/min [absolute value of [lambda]] is introduced: when SR < 20 the problem is not stiff, up to SR " 1000 the problem is classified as stiff, and when SR [greater than or equal to] 100000, the problem is very stiff. Stiff ODEs are called extremely stable if there is at least one eigenvalue with a large negative real part [2].

In the case of nonlinear problems the problem is more complicated since stiffness is a global problem and cannot be reduced to the solution structure in the neighbourhood of single points.

If we apply for the solution of stiff differential equations explicit solvers (e.q. MATLAB solvers ODE45, ODE23s) the solution may be oscillating or diverge, where oscillations of the solution are a result of the failure of the solution algorithm. For avoiding these oscillations implicit solvers or explicit solvers with a very small step size must be used. All this complicates the numerical work, besides small steps introduce too many round-off errors and cause numerical instabilities. Depending on the nature of stiffness different methods of solution should be used.

A detailed analysis of the methods for solving stiff ordinary differential equations can be found in the well known text-book [7] by Hairer and Wanner (here also an extensive list of publications is added). We would like to turn attention also to the paper [4] by Enright et al. in which different numerical methods have been tested for solving 25 systems of stiff equations. From the recent literature we refer here [1, 5, 6, 8, 12, 16, 18].

Of interest are also singular perturbation problems, where the differential equations contain a small parameter [epsilon] so that for [epsilon] [right arrow] 0 one of the equations loses its highest derivative [7]. In the limit case [epsilon] = 0 we get the so called differential algebraic system, which has been also investigated in many papers (see, e.g. [3, 5, 21]).

We are interested in solving stiff equations by the wavelet method. Hsiao [9, 10, 11] has proposed for solving linear and nonlinear differential equations the so called single-term Haar wavelet method (in fact it is a method of piecewise constant approximation). The method is very simple and effective, but its exactness in the region of rapid changes may turn out to be insufficient (see, Fig. 1 in [9]).

From recently published papers let us cite [19, 20]. In the these papers the numerical solutions of stiff differential equations, especially the Robertson problem, by using wavelet-based methods are investigated.

The main goal of the paper is to develop further Hsiao's results and demonstrate that the multi-term Haar wavelets are a powerful tool for solving stiff problems. The paper is organized as follows. In Section 2 the Haar wavelet method is shortly referred. In Section 3 this method is applied for solving linear stiff ODEs. Nonlinear problems are discussed in Section 4. Section 5 is devoted to solution of the Robertson's problem. Singular perturbation problems are considered in Section 6. In Section 7 some recommendations for further research are given.

2 Haar Wavelet Method

Let us consider the interval x [member of] [A, B], where A and B are given constants. We shall define the quantity M = [2.sup.J], where J is the maximal level of resolution. The interval [A, B] is divided into 2M subintervals of equal length; the length of each subinterval is [DELTA]x = (B - A)/(2M). Next two parameters are introduced: the dilatation parameter j = 0,1, ..., J and the translation parameter k = 0,1, ..., m - 1 (here the notation m = [2.sup.j] is introduced). The wavelet number i is identified as i = m + k + 1.

The i-th Haar wavelet is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.1)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The case i = 1 corresponds to the scaling function: [h.sub.1](x) = 1 for x [member of] [A, B] and [h.sub.1](x) = 0 elsewhere.

If we want to solve a n-th order ODE we need the integrals

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The case v = 0 corresponds to function [h.sub.i](t). Taking into account (2.1), these integrals can be calculated analytically; by doing it we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.2)

These formulas hold for i > 1. In the case i = 1 we have [[xi].sub.1] = A, [[xi].sub.2] = [[xi].sub.3] = B and 1

[P.sub.[alpha]!] = 1/[alpha]![(x - A).sup.[alpha]].

In the present paper the collocation method for solving the ODEs is applied. The collocation points are [x.sub.l] = 0.5[[[??].sub.l-1] + [[??].sub.l]], l = 1, 2, ..., 2M: the symbol [[??].sub.l] denotes the l-th grid point [[??].sub.l] = A + l[DELTA]x, / = 1, 2, ..., 2M. Eqs. (2.1), (2.2) are discretized by replacing x [right arrow] [[??].sub.l]. It is convenient to introduce the Haar matrices H(i,l) = [h.sub.i]([[??].sub.l]), [P.sub.v](i,l) = [p.sub.v,I]([[??].sub.l]). In the following sections computer simulations are carried out with the aid of the Matlab programs for which the matrix representation is effective.

3 Linear Problems

If we have to solve a n-th order linear ODE we seek the solution in the form

[y.sup.(n)]([[??].sub.l]) = [2M.summation over (i=1)] [a.sub.i][h.sub.i]([[??].sub.l]) = aH. (3.1)

Lower order derivatives (and the function y(x)) are obtained through integrations of (3.1). All these ingredients are incorporated into (3.1) which is discretized by the collocation method. In this way we get a system of equations from which the wavelet coefficients [a.sub.i] are calculated.

We investigate the accuracy of such algorithms. First, let us consider the case where exact solution [y.sub.ex](x) of the problem is known. We define in the collocation points the error function u(l) = y([??].sub.l]) = [y.sub.ex]([[??].sub.l]). Next the norm

[parallel]u[[parallel].sub.[] [([2M.summation over (l=1)] [[absolute value of u(l)].sup.p]).sup.1/p]

is introduced. The following two error estimates are applied (J denotes the resolution level):

(i) local estimate [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

(ii) global estimate [[sigma].sub.J] = [parallel]v[[parallel].sub.2]/(2M).

Error estimates in the case when the exact solution is not known is considered further on in Sections 4-6.

Example 1. Press et al. [17] investigated the system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.2)

We reduce (3.2) to the second order ODE

y" + 1001y' + 1000y = 0, x [member of] [0,1], (3.3) y(0) = 1, y'(0) = 0.

The exact solution of (3.3) is

[y.sub.ex] = 1/999 (1000[e.sup.-x] - [e.sup.-1000x]). (3-4)

The eigenvalues of (3.3) are [[lambda].sub.1] = -1, [[lambda].sub.2] = -1000; since [absolute value of [[lambda].sub.2]] >> [absolute value of [[lambda].sub.1]] the system (3.2) is stiff. It is demonstrated in [17] that explicit methods of integration bring to oscillations and the solution is unstable.

Now let us solve (3.3) by the Haar wavelet method. The solution is sought in the form

y" = aH, y' = a[P.sub.1] + [y'.sub.0]E, y = a[P.sub.1] + [y'.sub.0]x + [y.sub.0]E. (3.5)

Here E = (1,1, ..., 1) and all quantities are calculated in the collocation points. Substituting (3.5) into (3.3) we get the matrix equation

d(H + 1001[P.sub.1] + 1000[P.sub.2]) = -1000E

from which the wavelet coefficients [a.sub.j] are calculated. The solution of the problem y = y(x) is found from [(3.5).sub.3].

Computer simulations gave the following error estimates

[[delta].sub.3] = 4.5e - 4, [[sigma].sub.3] = 2.8e - 5, [[delta].sub.4] = 1.7e - 4, [[sigma].sub.4] = 5.3e - 6, [[delta].sub.5] = 4.1e - 5, [[sigma].sub.5] = 6.5e - 7.

It follows from here that high precision of the results is guaranteed already for a small number of collocation points (e.g. if J = 3 we have only 16 points).

Example 2. Mahmood et al. [16] solved by the Adomian decomposition method the linear initial value problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.6)

which has an exact solution

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (3.7)

The solution (3.7) is plotted in Fig. 1.

[FIGURE 1 OMITTED]

This problem was as well solved by Guzel and Bayram [6] who applied the power series method. Hojjati et al. [8] solved (3.4) with the aid of a predictor-corrector method, based on backward differentiation.

The wavelet solution is sought in the form

[y'.sub.1] = aH, [y.sub.1] = a[P.sub.1] + E [y'.sub.2] = bH, [y.sub.2] = b[P.sub.1] [y'.sub.3] = cH, [y.sub.3] = c[P.sub.1] - E. (3.8)

Substituting (3.8) into the system (3.6) we obtain the matrix system

a(H + 20[P.sub.1]) + 0.25b[P.sub.1] + 19.75c[P.sub.1] = -0.25E, -20a[P.sub.1] + b(H + 20.25[P.sub.1]) - 0.25c[P.sub.1] = 19.75E, -20a[P.sub.1] + 19.75b[P.sub.1] + c(H + 0.25[P.sub.1]) = 20.25E. (3.9)

Solving (3.9) we calculate the vectors o, b, c, the solution of the problem [y.sub.1](x), [y.sub.2](x), [y.sub.3](x) is obtained from (3.8).

Computations for J = 5 (with 64 collocation points) were carried out and the results were compared with the exact solution (3.5). The following error estimates were obtained:

[[delta].sub.5] ([y.sub.1]) = 2.0e - 5, [[delta].sub.5]([y.sub.1]) = 3.1e - 7, [[delta].sub.5] ([y.sub.2]) = 4.4e - 5, [[delta].sub.5]([y.sub.2]) = 6.9e - 7, [[delta].sub.5] ([y.sub.3]) = 1.6e - 3, [[delta].sub.5]([y.sub.3]) = 2.5e - 5.

4 Nonlinear Problems

Consider the equation

[y.sup.(n)] = f (x,y,y', ..., [y.sup.(n-1)]), x [member of] [A, B] (4.1)

with initial conditions [y.sup.(i)](A) = [y.sup.i.sub.0], i = 0,1, ..., 2M - 1. Here f is a nonlinear function.

The wavelet solution of the problem is sought in the form (3.1). Calculating again lower order derivatives by integration, replacing these results into (4.1) and satisfying these equations in the collocation points we get a nonlinear system 

F ([x.sub.l], [a.sub.l], ..., [a.sub.2M]) = 0, l = 1, 2, ..., 2M, (4.2)

from which the wavelet coefficients oj are calculated.

The system (4.2) is solved by the Newton method. If some initial guess for these coefficients [[??].sub.j] is known, then the corrected values of the wavelet coefficients are computed as

o = [??] + [lambda][DELTA]a,

here [DELTA]a = -- F/[partial derivative]F/[partial derivative]a and [lambda] is a coefficient, which is selected to guarantee the decrease of [parallel]F[[parallel].sub.[infinity]] (in the case of the exact solution we have [parallel]F(a)[[parallel].sub.[infinity]] = 0). This Newton step is repeated until [parallel]F[[parallel].sub.[infinity]] is sufficiently close to zero.

It is well known that the Newton method converges only when the initial guess is sufficiently good. To find such an initial estimate may be rather difficult. We have applied the following approach. We begin with a small number of collocation points (1, 2 or 4 points) for which the solution of (4.2) is more simple. Let us assume that we have found the solutions [a.sup.(0).sub.i] for some level of resolution J. The wavelet coefficients at the next level J + 1 are estimated as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Such an estimation is motivated by the fact that higher order coefficients of the sequence [[alpha].sub.i] are usually small (compare e.g. Fig. 2). These estimates are corrected by the Newton method and after that we can step over to the next level of resolution. Other possibilities for the initial guess of the wavelet coefficients are considered in Sections 5-6. For further details about this approach consult [13].

[FIGURE 2 OMITTED]

Example 3. Consider the initial value problem

y' = -[y.sup.2]/1 + x, y(0) = 1, x [member of][0, 200], (4.3)

which has the exact solution y = 1/(1 + ln(1 + x)). This problem is classified by Hsiao [10] as a problem of "strong instability". Hsiao's solution is based on the single-term Haar wavelet method.

Let us solve (4.3) with the aid of multi-term Haar wavelets. At the beginning, the solution y(x) rapidly decreases and then calms down. Therefore it is reasonable to solve the problem in two stages: first consider the phase of rapid decrease x [member of] [0, [delta]] and after that find solution for x [member of] [[delta], 200]. In the following calculations the value [delta] = 20 was chosen. In both regions we take

y' = aH, y = a[P.sub.1] + [y.sub.in]E. (4.4)

The function (4.2) and its gradient get the form

F(x) = aH + y/(E + x), [partial derivative]F/[partial derivative]a = H + 2y x [P.sub.1]/(E + x).

Here x = ([x.sub.1], [x.sub.2], ..., [x.sub.2]M) and the point denotes element- by-element multiplication:

[(y.[sup.2]).sub.l] = y[(l).sup.2], [(y x [P.sub.1]).sub.i.l] = y(l)[P.sub.1](i,l).

In the first phase x [member of] [0, 20] we have [y.sub.in] = 1. For the initial solution we take [a.sup.0] [equivalent to] 0. This value is corrected by the Newton method. By increasing step by step the resolution level J the solution is made more and more exact. Computer simulation for J = 5 gave the error estimates [[delta].sub.5] = 6.9e - 3, [[delta].sub.5] = 6.0e - 5. The value y(20) = 0.2475 was calculated, it is taken for the initial value [y.sub.in] for the second stage x [member of] [20, 200]. It is assumed [a.sup.0] [equivalent to] 0 and this value is again corrected by the Newton method. Here the number of collocation points may be smaller since already [[delta].sub.3] = 8. 5e - 3, [[delta].sub.3] = 1.6e - 4. The wavelet solution y = y(x) and the wavelet coefficients [a.sub.i], i = 1, 2, ..., 2M are plotted in Fig. 2. It follows from this figure that the number of significant wavelet coefficients is rather small.

5 Robertson's Problem

In 1966 Robertson investigated a chemical system containing fast and slow motions at the same time. By modelling this system he got the following mathematical model

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5.1)

The eigenvalues of (5.1) are given in [1]: [[lambda].sub.1] = 0, [[lambda].sub.2] = - 0.405, [[lambda].sub.3] = -2189.6 and the third one produces the stiffness.

The system (5.1) has been investigated in several papers; it is often taken as a benchmark for numerical methods. A detailed analysis about it can be found in the text-book [7]. It is shown that explicit methods give oscillating solutions (see Fig. 1.3 in [7]). Stability was guaranteed in the case of implicit methods (in [7] the Runge-Kutta codes DOPRI 5 and RADAU 5 were applied). In [11] the system (5.1) is integrated by the single term Haar wavelet method and in [14] the Adomian decomposition method was used. Solution obtained by the ODE45 code is plotted in Fig. 3. Small oscillations of the curve y2 indicate instability of the solution.

[FIGURE 3 OMITTED]

Qualitative analysis of (5.1) indicates that the component [y.sub.2] rapidly reaches its maximal value for which [y'.sub.2] =0 and 3 x [10.sup.7][y.sup.2.sub.2] [approximately equal to] 0.04, consequently max([y.sub.2]) [approximately equal to] 3. 65e - 5; after that the function [y.sub.2] decreases very slowly. Due to this fact it is again expedient to solve (5.1) separately for the regions x [member of] [0, [delta]] and x [member of] [[delta], [x.sub.max]]. The wavelet solution is sought in the form

[y'.sub.1] = aH, [y'.sub.2] = bH, [y'.sub.3] = cH, (5.2) [y.sub.1] = [y.sub.1](0)E + a[P.sub.1], [y.sub.2] = [y.sub.2](0)E + b[P.sub.1], [y.sub.3] = [y.sub.3](0)E + c[P.sub.2].

Next the vectors [F.sub.1], [F.sub.2], [F.sub.3] are introduced

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

On account of (5.2) these vectors are functions of the wavelet coefficients o, b, c. We have to find such values o, b, c for which [F.sub.1] = [F.sub.2] = [F.sub.3] = 0. This can be realized according to the following procedure.

(i) We assume that some initial estimates [[??].sub.1], [[??].sub.2], [[??].sub.3] are known. From (5.2) we calculate the estimates [??], [??], [??] and also [[??]'.sub.1], [??]'.sub.2], [??]'.sub.3].

(ii) Evaluate [F.sub.2].

(iii) Since a rapid growth takes place for the function [y.sub.2], we shall vary only the coefficients b (o and c are fixed) and calculate the gradient

[partial derivative][F.sub.2] = H + [10.sup.4][y.sub.3] x [P.sub.1] + 6 x [10.sup.7][y.sub.2] x [P.sub.1]. (5.3)

(iv) Using the Newton method next approximation for b is calculated:

b = [??] [lambda][F.sub.2] [([partial derivative][F.sub.2]/[partial derivative]b).sup.-1].

For the coefficient [lambda] the value which minimizes [parallel][F.sub.2][parallel] is taken.

(v) Calculate from [(5.1).sub.3] and (5.2) the corrected values for c, [y.sub.3], [y'.sub.3].

(vi) Calculate from [(5.1).sub.1] and (5.2) the corrected values for o, [y.sub.1], [y'.sub.1].

(vii) Evaluate [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; if these norms are not sufficiently near to zero repeat the steps (ii)-(vii).

Computer simulations were carried out for [delta] = 0.005, [x.sub.max] = 0.3. In the first phase x [member of] [0, [delta]] the argument x is very small and we can develop the function [y.sub.1], [y.sub.2], [y.sub.3] into power series. Taking into account only the first terms of this series we find

[[??].sub.1] = E - 0.04x, [[??].sub.2] = 0.04E - 1.6 x [10.sup.4] x. [sup.3], [y.sub.3] = 1.6 x [10.sup.4] x. [sup.3].

According to (5.2) we obtain

[??] = -0.04E/H, [??] = -4.8 x [10.sup.4] x. [sup.2]/H, [??] = -[??].

In order to correct these values the Newton iteration steps are done. Restricting ourselves to two iterations we get for J = 5:

[y.sub.1]([delta]) = 0.9998, [y.sub.2]([delta]) = 3.648e - 5, [y.sub.3]([delta]) = 1.626e - 5.

These values were taken as initial values for the second phase x [member of] [0. 005, 0. 3]. In this phase the functions [y.sub.1], [y.sub.3] decrease very slowly, therefore we take

[[??].sub.1] = [y.sub.1] ([delta])E, [[??].sub.2] = [y.sub.2]([delta])E.

Calculating [y.sub.3] from [(5.1).sub.3] we find

[[??].sub.3] = [y.sub.3]([delta])E + 3 x [10.sup.7][y.sub.2][([delta]).sup.2]x.

Correcting once again the values by the Newton method we find

[y.sub.1] (0.3) = 0.9888, [y.sub.2](0.3) = 3.446e - 5, [y.sub.3](0.3) = 1.138e - 4.

Carrying out the calculations by the ODE45 code we get

[y.sub.1] (0.3) = 0.9987, [y.sub.2](0.3) = 3.444e - 5, [y.sub.3](0.3) = 1.129e - 4.

These results are in accordance with our results.

It is interesting to note that in the second phase only the first wavelet coefficients are significant. To demonstrate this, the first three wavelet coefficients are written down:

a: -0.0376, -0.0011, -0.0006, ... 105b: -0.6851, -0.0183, -0.0092, ... c: 0.0377, 0.0011, 0.0006, ...

All the other coefficients are considerably smaller.

6 Singular Perturbation Problems

A singular perturbation problem of Index 1 is defined as [7]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6.1)

If [epsilon] is small, then the second equation of (5.3) is stiff. In the case [epsilon] = 0 we get a differential-algebraic equation, when a differential equation is combined with an algebraic equation. Eqs. (6.1) for [epsilon] = 0 are called also the reduced system. It is usually much easier to analyze the reduced system than the primary system.

Consider the phase space (y, z). The slope of the phase curves is

dz/dy = g(y, z)/[epsilon]f(y, z).

We assume that the initial point (y0, z0) is not placed on the curve g(y, z) = 0 and [epsilon] << 1; in this case the phase curve is very steep. We are mainly interested in the case where the phase point (y, z) moves towards the curve g(y, z) = 0 and reaches it at some instant x = #. Now z'([delta]) = 0 and the following motion proceeds along the curve g(y, z) = 0. Numerical difficulties can appear in the first phase of motion x [less than or equal to] [delta], since for [epsilon] << 1 the second equation of (6.1) is stiff.

Example 4. Consider the van der Pol equation

[epsilon]z" + ([z.sup.2] - 1) z' + z = 0,

which can be rewritten in the form [7]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6.2)

Boundary condition for z' can be calculated from [(6. 2).sub.2]:

ez'([x.sub.0]) = ([y.sub.0] - [z.sup.3.sub.0]/3 + [z.sub.0]). (6.3)

The reduced problem [epsilon] = 0 can be easily solved:

ln [absolute value of z] - 0.5[z.sup.2] = x + c.

Computer simulations were carried out for [epsilon] = 0.001, [y.sub.0] = 0.5, [z.sub.0] = 1. According to (6.3) we have [z'.sub.0] = 7000/6. Since z' > 0 the phase point (y, z) moves toward the phase curve

y = [z.sup.3]/3 - z (6.4)

and reaches it at some instant x = [delta].

Numerical results obtained by the MATLAB program ODE23s and with the aid of the piecewise constant approximation method (PCA) are plotted in Fig. 4. It follows from this figure that at this number of calculation points the PCA solution is unstable.

[FIGURE 4 OMITTED]

Now let us present the Haar wavelet solution. Again it is suitable to consider separately two phases of motion: (i) the phase of rapid changes x [member of] [0, [delta]] and (ii) motion along the curve (6.4) for x [member of] [[delta], [x.sub.fin]]. In the following calculations it was taken [delta] = 0.003, [x.sub.fin] = 0.02.

In the first phase the solution is sought in the form

z" = aH, z' = a[P.sub.1] + [z'.sub.0]E, z = a[P.sub.2] + [z'.sub.0]x + [z.sub.0]E.

Next the 2M-dimensional vector F is defined:

F(a) = [epsilon]z" + (z. [sup.2] - E). z' + z;

its gradient is given by

[partial derivative]F/[partial derivative]a = [epsilon] H + [z. [sup.2] - E). [P.sub.1] + (E + 2z. z'). [P.sub.2].

The corrected estimate for wavelet coefficients is given as

a = [??] - [lambda]F [([partial derivative]F/[partial derivative]a).sup.-1]. do /

Computer simulations were carried out for [epsilon] = 4, for the initial approximation the value [??] [equivalent to] 0 was taken and three Newton iteration steps were made. For the final result we got [parallel]F[[parallel].sub.[infinity]] = 0.0026, y([delta]) = 0.4938, z([delta]) = 1.940, z'([delta]) = 0.9224.

In the second phase x [member of] [[delta], [x.sub.fin]] we assume

z' = bH, z = b[P.sub.1] + z([delta])E.

Now the function

G(b) = (z. [sup.2] - E). z' + z (6.5)

is to be minimized. The gradient of (6.5) is given by

[partial derivative]G/[partial derivative]b = ([z. [sup.2] - E).H + (E + 2z.z'). [P.sub.1].

Again we assume b [equivalent to] 0; after two Newton steps we find [parallel]G[[parallel].sub.[infinity]] = 8.5e - 5, y(0. 02) = 0. 4629, z(0.02) = 1 . 928. This is in a good accordance with the ODE23s solution which gives y(0.02) = 0.4618, z(0. 02) = 1 . 928.

As to the wavelet coefficients then the first coefficient is [b.sub.1] = - 0.7058; all other coefficients are very small (less than 0.27% from b1). This is due to the quasi-linearity of the function z(x).

7 Conclusions

The examples of this paper demonstrate that in solving stiff systems the Haar wavelet method can successfully compete with the other efficient methods. The main benefits of the Haar approach are simplicity (already a small number of grid points guarantees the necessary accuracy) and universality (the same approach is applicable for a wide class of differential and integral equations). The subprograms, developed for calculation of integrals [p.sub.[alpha]i](x) can be used without changes for solving different problems.

In the case of stiff problems the Haar wavelet method seems to be more stable than many other methods. There are some complementary possibilities to raise the stability of the solution. First, we could use the wavelets with a variable step size (see, [15], especially the example in Section 5). If we develop into the Haar series the highest derivative [y.sup.(n)], then it is not continuous: this fact may also cause some instabilities. In such cases it is advisable to put into series the higher derivative [y.sup.(n+1)]; this idea is implemented in [14].

In this paper only ordinary differential equations were considered, but the same approach is applicable also for stiff partial differential equations. We recommend to consult the paper [15] in which the Burgers equation was solved.

The Haar wavelet transform and PCA (single term Haar wavelet series) both consist of piecewise constant functions and therefore have the same convergence rate O([M.sup.-2]). Nevertheless the multi-Haar wavelet method has by our mind some advantages:

(i) Very high accuracy fast transformation exists and therefore a possibility of implementation of fast algorithms compared with other methods;

(ii) The simplicity and small computation costs, resulting from the sparsity of the transform matrices and small number of significant wavelet coefficients;

(iii) The method is also very convenient for solving boundary value problems, since the boundary conditions are taken into account automatically.

Acknowledgement

Financial support from The Estonian Science Foundation under Grant ETF 6697 is gratefully acknowledged.

Received March 19, 2009; revised September 22, 2009; published online November 10, 2009

References

[1] S. Abelman and K.C. Patidar. Comparison of some recent numerical methods for initial-value problems for stiff ordinary differential equations. Comput. Math. with Applications, 55:733-744, 2008. (doi:10.1016/j.camwa.2007.05.012)

[2] N.M. Bujurke, C.S. Salimath and S.C. Shirashetti. Numerical solution of stiff systems from nonlinear dynamics using single-term Haar wavelet series. Nonlinear Dynamics, 51:595-605, 2008. (doi:10.1007/s11071-007-9248-8)

[3] E. Celik. On the numerical solution of differential-algebraic equations with index--2. Appl. Math. Comp., 156:541-550, 2004. (doi:10.1016/j.amc.2003.08.018)

[4] W.H. Enright, T.E. Hull and B. Lindberg. Comparing numerical methods for stiff systems of ODEs. BIT, 15(1):10-48, 1975. (doi:10.1007/BF01932994)

[5] J. Gao and Y.-L. Jiang. An adaptive wavelet method for nonlinear differential-algebraic equations. Appl. Math. Comp., 189:208-220, 2007. (doi:10.1016/j.amc.2006.11.102)

[6] N. Guzel and M. Bayram. On the numerical solution of stiff systems. Appl. Math. Comp., 170:230-236, 2005. (doi:10.1016/j.amc.2004.11.035)

[7] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer-Verlag, 1991.

[8] G. Hojjati, M.Y. Rahimi Ardabili and S.M. Hosseini. A-EBDF: an adaptive method for numerical solution of stiff systems of ODEs. Math. Computers in Simulation, 66:33-41, 2004. (doi:10.1016/j.matcom.2004.02.019)

[9] C.-H. Hsiao. Haar wavelet approach to linear stiff systems. Math. Computers in Simulation, 64:561-567, 2004. (doi:10.1016/j.matcom.2003.11.011)

[10] C.-H. Hsiao. Numerical solution of stiff differential equations via Haar wavelets. Int. J. Comp. Math., 82(9):1117-1123, 2005. (doi:10.1080/00207160512331323308)

[11] C.-H. Hsiao and W.-J. Wang. Haar wavelet approach to nonlinear stiff systems. Math. Computers in Simulation, 57:347-353, 2001. (doi:10.1016/S0378-4754(01)00275-0)

[12] U. Lepik. Haar wavelet method for nonlinear integro-differential equations. Appl. Math. Comput., 176:324-333, 2006. (doi:10.1016/j.amc.2005.09.021)

[13] U. Lepik. Numerical solution of evolution equations by the Haar wavelet method. Appl. Math. Comput., 185(1):695-704, 2007. (doi:10.1016/j.amc.2006.07.077)

[14] U. Lepik. Solving differential and integral equations by the Haar wavelet method; revisited. Int. J. Math. Comput., 1:43-52, 2008.

[15] U. Lepik. Solving integral and differential equations by the aid of nonuniform Haar wavelets. Appl. Math. Comput., 198:326-332, 2008. (doi:10.1016/j.amc.2007.08.036)

[16] A.S. Mahmood, L. Casasus and W. Al-Hayani. The decomposition method for stiff systems of ordinary differential equations. Appl. Math. Comp., 167:964- 975, 2005. (doi:10.1016/j.amc.2004.06.134)

[17] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling. Numerical Recipes. Cambridge University Press, 1987.

[18] M. Shacham and N. Brauner. Preventing oscillatory behavior in error control for ODEs. Comput. Chem. Engng, 32:409-419, 2008. (doi:10.1016/j.amc.2004.06.134)

[19] M. Tade and T. Zhang. Wavelet approach incorporated with optimization for solving stiff systems. J. Math. Chem., 43:1553-1548, 2008.

[20] Y.-C. Tian and T. Zhang. Wavelet-based collocation method for stiff systems in process engineering. J. Math. Chem., 44:501-513, 2008. (doi:10.1007/s10910-007-9324-9)

[21] P.K. Vijalapura, J. Strain and S. Covindjee. Fractional step methods for index--1 differential-algebraic equations. J. Comput. Phys., 203:305-320, 2005. (doi:10.1016/j.jcp.2004.08.015)

Doi: 10.3846/1392-6292.2009.14

U. Lepik

Institute of Mathematics, Tartu University, Estonia

Liivi 2, 50409 Tartu, Estonia

E-mail: [e-mail]