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 eigenvalues
statistical term meaning latent root. 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 eigenvalue
In mathematical analysis, one of a set of discrete values of a parameter, k, in an equation of the form Lx = kx. Such characteristic equations are particularly useful in solving differential equations, integral equations, and systems of with a large negative real part .
In the case of nonlinear A system in which the output is not a uniform relationship to the input.
nonlinear - (Scientific computation) A property of a system whose output is not proportional to its input. 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 (MATrix LABoratory) A programming language for technical computing from The MathWorks, Natick, MA (www.mathworks.com). Used for a wide variety of scientific and engineering calculations, especially for automatic control and signal processing, MATLAB runs on Windows, Mac and solvers ODE ode, elaborate and stately lyric poem of some length. The ode dates back to the Greek choral songs that were sung and danced at public events and celebrations. 45, ODE23s) the solution may be oscillating os·cil·late
intr.v. os·cil·lat·ed, os·cil·lat·ing, os·cil·lates
1. To swing back and forth with a steady, uninterrupted rhythm.
2. or diverge diverge - If a series of approximations to some value get progressively further from it then the series is said to diverge.
The reduction of some term under some evaluation strategy diverges if it does not reach a normal form after a finite number of reductions. , where oscillations oscillations See Cortical 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  by Hairer and Wanner (here also an extensive list of publications is added). We would like to turn attention also to the paper  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 In mathematics, more precisely in perturbation theory, a singular perturbation problem is a problem containing a small parameter that cannot be approximated by setting the parameter value to zero. 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 . In the limit case [epsilon] = 0 we get the so called differential algebraic 1. (language) ALGEBRAIC - An early system on MIT's Whirlwind.
[CACM 2(5):16 (May 1959)].
2. (theory) algebraic - In domain theory, a complete partial order is algebraic if every element is the least upper bound of some chain of compact elements. 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 (mathematics) wavelet - A waveform that is bounded in both frequency and duration. Wavelet tranforms provide an alternative to more traditional Fourier transforms used for analysing waveforms, e.g. sound. 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 approximation /ap·prox·i·ma·tion/ (ah-prok?si-ma´shun)
1. the act or process of bringing into proximity or apposition.
2. a numerical value of limited accuracy. ). 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 ).
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 See wavelet compression.
The elementary building blocks in a mathematical tool for analyzing functions. The functions can be very diverse; examples are solutions of a differential equation, and one- and two-dimensional signals. 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 max·i·mal
1. Of, relating to, or consisting of a maximum.
2. Being the greatest or highest possible. 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 dilatation /dil·a·ta·tion/ (dil?ah-ta´shun)
1. the condition, as of an orifice or tubular structure, of being dilated or stretched beyond normal dimensions.
2. the act of dilating or stretching. parameter j = 0,1, ..., J and the translation parameter k = 0,1, ..., m - 1 (here the notation notation: see arithmetic and musical notation.
How a system of numbers, phrases, words or quantities is written or expressed. Positional notation is the location and value of digits in a numbering system, such as the decimal or binary system. 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 ASCII or American Standard Code for Information Interchange, a set of codes used to represent letters, numbers, a few symbols, and control characters. Originally designed for teletype operations, it has found wide application in computers. ] (2.1)
[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 In mathematics, a collocation method is a method for the numerical solution of ordinary differential equation and partial differential equations and integral equations. The idea to choose a finite-dimensional space of candidate solutions (usually, polynomials up to a certain for solving the ODEs is applied. The collocation collocation - co-location 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 summation n. the final argument of an attorney at the close of a trial in which he/she attempts to convince the judge and/or jury of the virtues of the client's case. (See: closing argument) 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.  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  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.  solved by the Adomian decomposition method The Adomian decomposition method (ADM) is a non-numerical method for solving nonlinear differential equations, both ordinary and partial. The general direction of this work is towards a unified theory for Partial Differential Equations (PDE). 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  who applied the power series method. Hojjati et al.  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 Noun 1. nonlinear system - a system whose performance cannot be described by equations of the first degree
system, scheme - a group of independent but interrelated elements comprising a unified whole; "a vast system of production and distribution and consumption
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 coefficient /co·ef·fi·cient/ (ko?ah-fish´int)
1. an expression of the change or effect produced by variation in certain factors, or of the ratio between two different quantities.
2. , which is selected to guarantee the decrease of [parallel]F[[parallel].sub.[infinity infinity, in mathematics, that which is not finite. A sequence of numbers, a1, a2, a3, … , is said to "approach infinity" if the numbers eventually become arbitrarily large, i.e. ]] (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 .
[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  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 gradient
In mathematics, a differential operator applied to a three-dimensional vector-valued function to yield a vector whose three components are the partial derivatives of the function with respect to its three variables. The symbol for gradient is ∇. 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 : [[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 . It is shown that explicit methods give oscillating solutions (see Fig. 1.3 in ). Stability was guaranteed in the case of implicit methods (in  the Runge-Kutta codes DOPRI 5 and RADAU 5 were applied). In  the system (5.1) is integrated by the single term Haar wavelet method and in  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 Qualitative Analysis
Securities analysis that uses subjective judgment based on nonquantifiable information, such as management expertise, industry cycles, strength of research and development, and labor relations. 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 ex·pe·di·ent
1. Appropriate to a purpose.
a. Serving to promote one's interest: was merciful only when mercy was expedient.
b. 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 according to
1. As stated or indicated by; on the authority of: according to historians.
2. In keeping with: according to instructions.
3. 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 One repetition of a sequence of instructions or events. For example, in a program loop, one iteration is once through the instructions in the loop. See iterative development.
(programming) iteration - Repetition of a sequence of instructions. 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 
[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 differential equation
Mathematical statement that contains one or more derivatives. It states a relationship involving the rates of change of continuously changing quantities modeled by functions. is combined with an algebraic equation algebraic equation
Mathematical statement of equality between algebraic expressions. An expression is algebraic if it involves a finite combination of numbers and variables and algebraic operations (addition, subtraction, multiplication, division, raising to a power, and . 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 
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6.2)
Boundary condition boundary condition
The set of conditions specified for behavior of the solution to a set of differential equations at the boundary of its domain. 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).
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, , 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 ad·vis·a·ble
Worthy of being recommended or suggested; prudent.
ad·visa·bil to put into series the higher derivative [y.sup.(n+1)]; this idea is implemented in .
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  in which the Burgers equation was solved.
The Haar wavelet transform and PCA (single term Haar wavelet series In mathematics, a wavelet series is a representation of a square-integrable (real- or complex-valued) function by a certain orthonormal series generated by a wavelet. This article provides a formal, mathematical definition of an orthonormal wavelet and of the ) 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 Computation is a general term for any type of information processing that can be represented mathematically. This includes phenomena ranging from simple calculations to human thinking. costs, resulting from the sparsity sparse
adj. spars·er, spars·est
Occurring, growing, or settled at widely spaced intervals; not thick or dense.
[Latin sparsus, past participle of spargere, to scatter. 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.
Financial support from The Estonian Science Foundation under Grant ETF ETF
See Exchange Traded Fund.
See exchange-traded fund (ETF). 6697 is gratefully acknowledged.
Received March 19, 2009; revised September 22, 2009; published online November 10, 2009
 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)
 N.M. Bujurke, C.S. Salimath and S.C. Shirashetti. Numerical solution of stiff systems from nonlinear dynamics nonlinear dynamics, study of systems governed by equations in which a small change in one variable can induce a large systematic change; the discipline is more popularly known as chaos (see chaos theory). using single-term Haar wavelet series. Nonlinear Dynamics, 51:595-605, 2008. (doi:10.1007/s11071-007-9248-8)
 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)
 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)
 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)
 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)
 E. Hairer and G. Wanner. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. Springer-Verlag, 1991.
 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)
 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)
 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)
 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)
 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)
 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)
 U. Lepik. Solving differential and integral equations by the Haar wavelet method; revisited. Int. J. Math. Comput., 1:43-52, 2008.
 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)
 A.S. Mahmood, L. Casasus and W. Al-Hayani. The decomposition method In constraint satisfaction, a decomposition method translates a constraint satisfaction problem into another constraint satisfaction problem that is binary and acyclic. Decomposition methods work by grouping variables into sets, and solving a subproblem for each set. for stiff systems of ordinary differential equations. Appl. Math. Comp., 167:964- 975, 2005. (doi:10.1016/j.amc.2004.06.134)
 W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling. Numerical Recipes. Cambridge University Press Cambridge University Press (known colloquially as CUP) is a publisher given a Royal Charter by Henry VIII in 1534, and one of the two privileged presses (the other being Oxford University Press). , 1987.
 M. Shacham and N. Brauner. Preventing oscillatory oscillatory
characterized by oscillation.
see pendular nystagmus. behavior in error control for ODEs. Comput. Chem. Engng, 32:409-419, 2008. (doi:10.1016/j.amc.2004.06.134)
 M. Tade and T. Zhang. Wavelet approach incorporated with optimization optimization
Field of applied mathematics whose principles and methods are used to solve quantitative problems in disciplines including physics, biology, engineering, and economics. for solving stiff systems. J. Math. Chem., 43:1553-1548, 2008.
 Y.-C. Tian Tian
In indigenous Chinese religion, the supreme power reigning over humans and lesser gods. The term refers to a deity, to impersonal nature, or to both. 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)
 P.K. Vijalapura, J. Strain and S. Covindjee. Fractional fractional
size expressed as a relative part of a unit.
fractional catabolic rate
the percentage of an available pool of body component, e.g. protein, iron, which is replaced, transferred or lost per unit of time. step methods for index--1 differential-algebraic equations. J. Comput. Phys., 203:305-320, 2005. (doi:10.1016/j.jcp.2004.08.015)
Institute of Mathematics, Tartu University, Estonia
Liivi 2, 50409 Tartu, Estonia