# Application of the Haar wavelet transform to solving integral and differential equations/Haari lainikute rakendamine integraalja diferentsiaalvorrandite lahendamiseks.

1. INTRODUCTION AND BACKGROUND

Let us first explain why we need wavelets at all. A signal x = x(t) is oftenanalysed with the aid of the Fourier transform

F([omega]) = [[integral].sup.+[infinity].sub.-[infinity]] x(t)[e.sup.i[omega]t] dt, (1)

where w denotes the frequency and F--the amplitude.

Two signals x(t) and their Fourier diagrams are plotted in Fig. 1. Although in both cases the Fourier diagrams are quite similar (both have two peaks), the signals are completely different and the time information is lost. This is a serious drawback: the Fourier method analyses the signal over the whole domain and is unable to characterize its behaviour in time.

[FIGURE 1 OMITTED]

Wavelet analysis allows us to represent a function in terms of a set of basic functions, called wavelets, which are localized both in space and time. Here a continuous function [psi](t), called the mother wavelet, is introduced. The corresponding wavelet family is formed by the dilation and translation of the function [psi](t):

[[psi].sub.j,k](t) = [2.sup.-j/2][psi]([2.sup.-j] t - k), (2)

where j and k are non-negative integers.

Depending upon the choice of the mother wavelet [psi](t), various wavelet families are obtained. The wavelets introduced by Daubechies [1] are quite frequently used for solving different problems. The Daubechies mother wavelet is plotted in Fig. 2. These wavelets are differentiable and have a minimum size support.

A shortcoming of the Daubechies wavelets is that they do not have an explicit expression and therefore analytical differentiation or integration is not possible. This complicates the solution of differential equations where integrals of the following type

[[integral].sup.b.sub.a] G (t, [[psi].sub.i,k], d[[psi].sub.i,k]/dt, [d.sup.2][[psi].sub.i,k]/d[t.sup.2], ...)dt (3)

must be computed (G is generally a nonlinear function). For calculating such integrals the conception of connection coefficients is introduced. Calculation of these coefficients is very complicated and must be carried out separately for different types of integrals (see, e.g., [2]). Besides, it can be done only for some simpler types of nonlinearities (mainly for quadratic nonlinearity). This remark holds also for other types of wavelets (such as Symlet, Coiflet, etc. wavelets).

[FIGURE 2 OMITTED]

The wavelet method was first applied to solving differential and integral equations in the 1990s. A survey of early results in this field can be found in [3]. Lately the number of respective papers has greatly increased and it is not possible to analyse them all here, but some are discussed in the following sections.

Due to the complexity of the wavelet solutions, some pessimistic estimates exist. So Strang and Nguyen write in their text-book [4] "... the competition with other methods is severe. We do not necessarily predict that wavelets will win" (p. 394). Jameson [5] writes "... nonlinearities etc., when treated in a wavelet subspace, are often unnecessarily complicated ... There appears to be no compelling reason to work with Galerkin-style coefficients in a wavelet method" (p. 1982).

Obviously attempts to simplify solutions based on the wavelet approach are wanted. One possibility is to make use of the Haar wavelet family.

In 1910 Alfred Haar introduced a function which presents a rectangular pulse pair. After that various generalizations were proposed (a state of the art about Haar transforms can be found in [6]). In the 1980s it turned out that the Haar function was in fact the Daubechies wavelet of order 1. It is the simplest orthonormal wavelet with compact support.

The Haar wavelet family is defined for t [member of] [0, 1] as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

Here [[xi].sub.1] = k/m, [[xi].sub.2] = (k + 0.5)/m, [[xi].sub.3] = (k + 1)m: The integer m = [2.sup.j] (j = 0, 1, ..., J) indicates the level of the wavelet; k = 0, 1, ..., m - 1 is the translation parameter. The maximal level of resolution is J. The index i in (4) is calculated according to the formula i = m + k + 1; in the case m = 1, k = 0 we have i = 2; the maximal value of i is i = 2M = [2.sup.J+1]. It is assumed that the value i = 1 corresponds to the scaling function for which [h.sub.1] = 1 for t [member of] [0, 1]. The first eight Haar functions are plotted in Fig. 3.

An essential shortcoming of the Haar wavelets is that they are not continuous. The derivatives do not exist in the points of discontinuity, therefore it is not possible to apply the Haar wavelets directly to solving differential equations.

There are at least two possibilities of ending this impasse. First, the piecewise constant Haar functions can be regularized with interpolation splines; this technique has been applied by Cattani [7,8]. This greatly complicates the solution process and the main advantage of the Haar wavelets--their simplicity--gets lost.

Another possibility was proposed by Chen and Hsiao [9,10]. They recommended to expand into the Haar series not the function itself, but its highest derivative appearing in the differential equation; the other derivatives (and the function) are obtained through integrations. All these ingredients are then incorporated into the whole system, discretized by the Galerkin or collocation method.

Chen and Hsiao demonstrated the possibilities of their method by solving linear systems of ordinary differential equations (ODEs) and partial differential equations (PDEs). In [10] an optimal control problem with the quadratic performance index is discussed.

In [11,12] Hsiao and Wang applied this method to solving singular bilinear and nonlinear systems. Nonlinear stiff systems were examined in [13].

In [14] Hsiao demonstrated that the Haar wavelet approach is effective also for solving variational problems. Of importance is also the paper [15] by Hsiao and Rathsfeld, in which wavelet collocation methods are applied to the boundary element solution of the first-kind integral equations arising in acoustic scattering.

[FIGURE 3 OMITTED]

The aim of the present paper is to acquaint the readers with the possibilities of the Haar wavelet method in solving integral and evolution equations. The illustrating examples are taken from papers published by the author.

2. INTEGRATION OF HAAR WAVELETS

If we want to integrate differential equations by the Hsiao-Chen method, we have to evaluate the integrals

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

Carrying out these integrations with the aid of (4), we find

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

It is necessary to evaluate these functions for a prescribed J only once and save them in the storage of the computer.

It is convenient to pass to the matrix formulation. For this purpose the interval t [member of] [0, 1] is divided into 2M parts of equal length [DELTA]t = 1/2M, where M = [2.sup.J]. The grid points are

t(l) = (l - 0.5)[DELTA]t, l = 1, 2, ..., 2M. (7)

The Haar coefficient matrix H is defined as H(i, l) = [h.sub.i](l). The integral matrices P[upsilon] have the elements P[upsilon](i, l) = [p.sub.i,[upsilon]](l).

Chen and Hsiao defined the integral matrix in a different way [5]. They introduced the row vector

[h.sub.([mu])](t) = [[h.sub.1](t), [h.sub.2](t), ..., [h.sub.[mu]](t)], (8)

where [mu] = 2M = [2.sup.J+1]. Now we have

[[integral].sup.t.sub.0] [h.sub.([mu])]([tau])d[tau] [approximately equal to] [P.sub.([mu]x[mu])][h.sub.([mu])](t). (9)

The square matrix [P([mu]x[mu])] is called the operational matrix of integration. Chen and Hsiao showed that the following recursive formula holds:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (10)

whereas [P.sub.(1x1)] = 0.5.

Formulae (8)-(10) can be applied to solving first-order differential equations. Higher-order equations have to be reduced to a system of first-order equations, but this increases the order of the Haar matrices and may cause complications of computational character.

3. INTEGRAL EQUATIONS

In this section application of wavelet methods to solution of integral equations is discussed. The following types of integral equations are considered: (i) Fredholm equation

u(x) - [[integral].sup.[beta].sub.[alpha]] K(x, t, u(t))dt = f(x), x, t [member of] [[alpha], [beta]], (11)

(ii) Volterra equation

u(x) - [[integral].sup.x.sub.[alpha]] K(x, t, u(t))dt = f(x), x, t, [member of] [[alpha], [beta]], (12)

(iii) integro-differential equation

au'(x) + bu(x) - [[integral].sup.x.sub.0] K(x, t, u(t), u'(t))dt = f(x),

for u'(0) = [u.sub.0], x, t [member of] [0, [beta]]. (13)

If the kernel K is a linear function in regard of u(t) and u'(t), we have a linear integral equation, otherwise the equation is nonlinear.

In recent years interest in the solution of integral equations by the wavelet methods has greatly increased (we have counted more than 50 papers on this topic). The first paper in this field according to our information is the one by Alpert et al. [16], published in 1993. In the following papers mostly Fredholm or Volterra equations are considered (consult, e.g., [17-20]). The Abel equation is discussed in [18], the Hammerstein equation in [21,22], integro-differential equations are considered in [23]. Different types of wavelets are applied, such as Daubechies [17,20], Coifman [24], Cohen [25], and Legendre [18,26,27] wavelets, and wavelets of B-splines [28]. In most papers linear integral equations are solved; nonlinear equations are considered in [18,20,26,29].

The wavelet equations are discretized by the Galerkin or collocation method. In the case of nonlinear problems we get for the wavelet coefficients a nonlinear system of equations, which can be solved by the Newton method.

Now we pass to the solutions based on the Haar wavelets. Since the Haar wavelets are defined only for the interval [0,1], we must first normalize Eqs (11)-(13). This can be done by the change of the variables

[x.sub.*] = (x - [alpha])/([beta] - [alpha]); [t.sub.*] = (t - [alpha])/([beta] - [alpha]).

In [30] the Haar wavelet method is applied to solving different types of linear integral equations (Fredholm, Volterra, integro-differential, weakly singular integral equations), also the eigenvalue problem is solved. Nonlinear integral equations are considered in [31,32]. It is demonstrated that the recommended method suits well for solving boundary value problems of ODEs.

Maleknejad and Mirzaee [33] applied the Haar wavelet method to solving linear Fredholm integral equations of the second kind.

In the following the Haar wavelet approach is used to solve some types of integral equations. The efficiency of the method is demonstrated by some numerical examples; for getting the error estimates, problems for which the exact solution [u.sub.ex](x) is known are considered. The accuracy of the results is estimated by the error function

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

Computations carried out in [30] showed that the solutions obtained with the aid of the collocation method are simpler than those got by the Galerkin method, besides, their accuracy seems to be better. For this reason in the following the collocation method is used. All calculations were made with the aid of MATLAB programs.

(i) Linear Fredholm integral equation

Let us consider the equation

u(x) - [[integral].sup.1.sub.0] K(x, t)u(t)dt = f(x), x, t [member of] [0, 1] (15)

and seek the solution in the form

u(t) = [2M.summation over (i=1)] [a.sub.i][h.sub.i](t), (16)

where [a.sub.i] are the wavelet coefficients.

Putting this result into (15), we find

[2M.summation over (i=1)] [a.sub.i][h.sub.i](x) - [2M.summation over (i=1)] [a.sub.i][G.sub.i](x) = f(x), (17)

where

[G.sub.i](x) = [[integral].sup.1.sub.0] K(x, t)[h.sub.i](t)dt. (18)

Satisfying (17) in the collocation points (7), we get a linear system of equations for the coefficients [a.sub.i]:

[2M.summation over (i=1)] [a.sub.i][[h.sub.i]([x.sub.l]) - [G.sub.i]([x.sub.l])] = f([x.sub.l]), l = 1, 2, ..., 2M. (19)

More convenient is the matrix form of (16) and (19)

u = aH,

a(H - G) = F, (20)

where

u = {u([t.sub.l])}, F = {f([x.sub.l])},

G = {[G.sub.i]([x.sub.l])}.

Example 1. Let us take K = x + t, f(x) = [x.sup.2]. Equation (15) has the exact solution

[u.sub.ex] = [x.sup.2] - [5.sub.x] - 17/6.

In the present case

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Computations were carried out for different values of J. The errors are shown in Table 1.

(ii) Nonlinear Volterra equation

Let us consider Eq. (12) in the normalized form [alpha] = 0; [beta] = 1. Satisfying it in the collocation points

x(l) = (l - 0.5)[DELTA]t, l = 1, 2, ..., 2M,

we obtain

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

The solution is again sought in the form (16). The system of equations for evaluating the wavelet coefficients is now nonlinear and must be solved with the aid of some numerical method. We have applied for this purpose the Newton method, which brings us to the equation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (22)

where

[partial derivative]K/[partial derivative][a.sub.i] = [partial derivative]K/[partial derivative]u [h.sub.i](t). (23)

According to (4), in each subinterval [h.sub.i](t) is constant between two collocation points. This enables us to calculate the integrals in (22) analytically. First we carry out the computations for one subintegral ([t.sub.s], [t.sub.s+1]). Summing up all these results, we get the exact values of these integrals for the whole domain [0, x].

If some approximate solution [a.sup.([upsilon]).sub.i] is known, it can be corrected with the aid of

[a.sup.([upsilon]+1).sub.i] = [a.sup.([upsilon]).sub.i] + [DELTA][a.sup.([upsilon]).sub.i], i = 1, 2, ..., 2M.

For details consult the paper [32].

Example 2. Solve the equation

u(x) = 1/2 [[integral].sup.x.sub.0] u(t)u(x - t)dt + 1/2 sin x, x [member of] [0, 1]. (24)

It has the exact solution u(x) = [J.sub.i](x), where the symbol [J.sub.i](x) denotes the Bessel function of order 1.

In this case

K = 1/2 u(t)u(x - t),

[partial derivative]K/[partial derivative][a.sub.i] = 1/2[u(x - t)[h.sub.i](x - t) - u(t)[h.sub.i](t)]. (25)

Error estimates (14) for this case are shown in Table 2.

(iii) Boundary value problems of ODEs

The method of solution described in this section is also very effective for solving boundary value problems if we rewrite the equation in the form of an integro-differential equation.

Consider the equation

u'' = K(x, u, u'), x [member of] [0, 1] (26)

with the boundary conditions u(0) = A, u(1) = B.

By integrating (26) we get the integro-differential equation

u'(x) = [[integral].sup.x.sub.0] K[t, u(t), u'(t)dt + u'(0). (27)

We seek the solution in the form

u'(t) = [2M.summation over (i=1)] [a.sub.i][h.sub.i](t),

u(t) = [2M.summation over (i=1)] [a.sub.i][p.sub.i,1](t) + u(0), (28)

where [p.sub.i,1] denotes the integral (5).

Satisfying the boundary conditions, we get u(1) = [a.sub.1] + u(0) or

[a.sub.1] = B - A. (29)

Since [a.sub.1] is now prescribed, we have [DELTA][a.sub.1] = 0; it means that the matrix of system (22) must be replaced with a reduced matrix for which the first row and column are deleted. Besides, (23) obtains the form

[partial derivative]K/[partial derivative][a.sub.i] = [partial derivative]K/[partial derivative]u [p.sub.i,1](t) + [partial derivative]K/[partial derivative]u' [h.sub.i](t). (30)

For further details consult again [33].

Example 3. Consider the boundary value problem

u'' - uu' = 1/4 [square root of [(1 + x).sup.3]] - 1/2, u(0) = 1, u(1) = [square root of 2], (31)

which has the exact solution

u = [square root of 1 + x].

By integrating (31) we obtain

u'(x) = [[integral].sup.x.sub.0] uu'dt + u'(0) - 1/2 (1 + x) + 1/2[square root of 1 + x] (32)

with the initial condition u(0) = 1. It follows from boundary conditions that [a.sub.1] = u(1) - u(0) = [square root of 2] - 1.

According to (28),

u'(0) = [2M.summation over (i=1)] [a.sub.i][h.sub.i](0), (33)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now

[partial derivative]K/[partial derivative][a.sub.i] = u'(t)[p.sub.i,1](t) + u(t)[h.sub.i](t), i = 2, 3, ..., 2M.

Some numerical calculations were carried out. The errors of the obtained results are presented in Table 3.

4. DIFFERENTIAL EQUATIONS

Papers on solving ODEs and PDEs by the wavelet methods are quite numerous and it is not possible to cite all of them (some information can be found, e.g., in [34,35]). Let us restrict the research area. We shall mainly consider the solutions based on the Haar wavelets. Our analysis has shown [34] that in solving nonlinear ODEs the wavelet method has no preferences over the classical Runge-Kutta method. For this reason we shall concentrate our analysis on PDEs.

An evolution process u = u(x, t), where x is the space coordinate and t is time, is described by a PDE. The simplest equations are the diffusion (or heat) equation

[partial derivative]u/[partial derivative]t = A [[partial derivative].sup.2]u/[partial derivative][x.sup.2] (34)

and the Burgers equation

[partial derivative]u/[partial derivative]t + u [partial derivative]u/[partial derivative]x = [upsilon] [[partial derivative].sup.2]u/[partial derivative][x.sup.2]. (35)

Due to simplicity they have proved to be a touchstone for new numerical methods of solution. An outline of such papers can be found in [34,35].

There are several possibilities for solving PDEs by the wavelet method. First we can introduce the two-dimensional wavelets

[psi](x, t) = [psi](x) [direct sum] [psi](t). (36)

Here [psi](x), [psi](t) are 2M-dimensional vectors, composed of the wavelets (2); [psi](x, t) is a 2M x 2M matrix; the symbol [direct sum] denotes Kronecker multiplication.

This approach has been applied in [36] for harmonic wavelets and in [37] for the Haar wavelets. The shortcoming of it is that we have to deal with matrices of high dimension (e.g., if the resolution level is J = 5, then M = 16 and the wavelet coefficients matrix has 4096 elements!).

An alternative method is applied in [9] to compute an electrical transmission line and in [34] to solve the diffusion equation (34). The solution is sought in the form

[partial derivative]u/[partial derivative]t] = a(x)H(t), (37)

where [partial derivative]u/[partial derivative]t and a are 2M-component vectors, H is the Haar matrix. Integrating (37) in regard of t and differentiating twice with respect to x, we get

u(x, t) = a(x)[P.sub.1](t) + u(x, 0)E,

[[partial derivative].sup.2]u/[partial derivative][x.sup.2] = [d.sup.2]a/d[x.sup.2] [P.sub.1](t) + [d.sup.2]u(x, 0)/d[x.sup.2] E, (38)

where [P.sub.1] is integral matrix and E is a 2M-dimensional unit vector. Inserting these results into (37), we get for a(x) a system of linear ODEs, which is solved analytically with the aid of the matrix exponential function expm. This method is quite effective, but applicable only in linear problems.

Now let us analyse the third method. Here the time interval t [member of] [0, T] is divided into N parts of equal length [DELTA]t and [t.sub.s] = s[DELTA]t, i = 0, 1, 2, ..., N. Following the suggestion of Chen and Hsiao [9], the highest derivative appearing in the evolution equation is developed into the Haar series. This term can be written in the matrix form

[a.sub.t](:)H(:, l), l = 1, 2, ..., 2M. (39)

It is assumed that [a.sub.t] = const for the subinterval t [member of] [[t.sub.s], [t.sub.s+1]]. This allows us to integrate (39) in regard of time t. Integration in regard of the space variable is carried out by the wavelet method. The evolution equation is satisfied by the collocation or Galerkin method. Transfer to the next instant [t.sub.s+1] is performed by using some finite-difference method (e.g. by the Euler or Crank-Nicholson method). This approach is applied in [38] for the Burgers equation and in [34] for the diffusion equation. It should be mentioned that the method has also a weakness. For calculating the wavelet coefficients we must invert some matrices, but they--especially for higher-order equations--are often nearly singular. This is due to the integral matrices [P[upsilon]], which are calculated according to (6). Even in the case of a small value J = 3 the determinants of these matrices are |[P.sub.1]| = 2.7E-20, [absolute value of [P.sub.2]] = 3.4E-49, [absolute value of [P.sub.3]] = 6.6E-81, [absolute value of [P.sub.4]]] = -2.2E-114. Here application of the theory of sparse matrices [38] may be of some help. Another possibility is to divide the space interval x [member of] [[x.sub.min], [x.sub.max]] into N segments and apply the Haar wavelet method to each segment (of course, some continuity conditions must be fulfilled on the boundaries of these segments). This approach was applied in [34,39].

Below we shall consider the third method. The application of this method is illustrated by solving the sine-Gordon equation.

5. SINE-GORDON EQUATION

Let us consider the sine-Gordon equation

1/[L.sup.2] [[partial derivative].sup.2]u/[partial derivative][x.sup.2] - [[partial derivative].sup.2]u/[partial derivative][t.sup.2] = x [member of] [0, 1], t [greater than or equal to] 0, (40)

with the boundary and initial conditions

u(x, 0) = f(x), [partial derivative]u/[partial derivative]t (x, 0) = g(x),

u(0, t) = [psi](x), [partial derivative]u/[partial derivative]x (0, t) = [psi](x).

Here f, g, [psi], [psi] are prescribed functions, L is a constant.

This equation has an analytical solitary wave solution

[u.sub.ex](x, t) = 4 arctan[exp(z)],

z = [alpha](x - [beta]t), [alpha] = L/[square root of (1 - [L.sup.2][[beta].sup.2]]. (41)

Solution of the sine-Gordon equation is discussed in many papers. According to Science Direct, 368 papers have been published on this subject since 1996. Among these we found only one paper [40] in which the wavelet method (Gaussian wavelets) was applied.

Now let us present the solution based on the Haar wavelet method [35]. According to Chen and Hsiao [9], the solution of (40) is sought in the form

[??]''(x, t) = [2M.summation over (i=1)] [a.sub.s][h.sub.i](x), t [member of] [[t.sub.s], [t.sub.s+1]], x [member of] [0, 1]. (42)

This equation is integrated twice in regard of x in the limits [0, x] and in regard of t in the limits [[t.sub.s], [t.sub.s+1]]. By doing this we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

These results are discretized by replacing x [right arrow] [x.sub.l], t [right arrow] [t.sub.s+1]. To simplify the writing of the formulae, the matrix formulation is used and the notations [DELTA]t = [t.sub.s+1] - [t.sub.s], u(l, s) = u([x.sub.l], [t.sub.s]) etc. are introduced. Taking into account the initial conditions (40), Eqs (43) can be put into the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (44)

Inserting these results into (40), we get a linear matrix equation for calculating the wavelet coefficients [a.sub.s](:):

[a.sub.s](:)[P.sub.2](:, l) = 1/[L.sub.2] u''(s, l) + sin u(s, l) - [x.sub.l] [??](s + 1) - [??](s + 1). (45)

Example 4. Computer simulation was carried out for t [member of] [10, 30], L = 20, [beta] = 0.025. If we want to get the classical solitary wave solution, we must take

f(x) = 4 arctan[exp([alpha]x)],

g(x) = [alpha]V ([alpha]x),

[phi](t) = 4 arctan[exp(-[alpha][beta]t)],

[psi](t) = -[alpha][beta]V (-[alpha][beta]t), (46)

where

V (z) = 4[e.sup.z]/1 + [e.sup.2z].

The accuracy of our approach is estimated by the error function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (47)

The calculations show that the function v(t) increases monotonically, therefore [v]([t.sub.max]) is taken as the error estimate. Some computer results are presented in Table 4 and Fig. 4.

[FIGURE 4 OMITTED]

It follows from Table 4 that already in the cases J = 4 or J = 5 we get the results which visually coincide with the exact solution. The time step was taken as [DELTA]t = 0.005 or [DELTA]t = 0.001; further decrease in [DELTA]t gave no considerable effect. If the accuracy of these results is insufficient, more precise results could be obtained by the segmentation method proposed in [34].

We would like to draw attention to Fig. 4b. It follows from here that most of the wavelet coefficients are zero (or near to zero) and the number of significant coefficients is quite small. This is one of the reasons for rapid convergence of the Haar wavelet series.

6. CONCLUSIONS

The main goal of this paper was to demonstrate that the Haar wavelet method is a powerful tool for solving different types of integral equations and partial differential equations. The method with far less degrees of freedom and with smaller CPU time provides better solutions than classical ones.

The main advantage of this method is its simplicity and small computation costs, resulting from the sparsity of the transform matrices and the small number of significant wavelet coefficients. The method is also very convenient for solving the boundary value problems, since the boundary conditions are taken care of automatically.

ACKNOWLEDGEMENT

Financial support from the Estonian Science Foundation (grant No. 6697) is gratefully acknowledged.

Received 13 November 2006

REFERENCES

[1.] Daubechies, I. Orthonormal bases of compactly supported wavelets. Comm. Pure Appl. Math., 1988, 41, 909-996.

[2.] Chen, M.-Q., Hwang, C. and Shin, Y.-P. The computation of wavelet-Galerkin approximation on a bounded interval. Int. J. Numer. Methods Eng., 1996, 39, 2921-2944.

[3.] Dahmen, W., Kurdila, A. J. and Oswald, P. (eds). Multiscale Wavelet Methods for Partial Differential Equations. Academic Press, San Diego, 1997.

[4.] Strang, G. and Nguyen, T. Wavelets and Filter Banks. Wellesley-Cambridge Press, Wellesley (MA), 1997.

[5.] Jameson, L. A wavelet-optimized, very much order adaptive grid and order numerical method. SIAM J. Sci. Comput., 1998, 19, 1980-2013.

[6.] Stankovic, R. S. and Falkowski, B. J. The Haar wavelet transform: its status and achievements. Comput. Electr. Eng., 2003, 29, 25-44.

[7.] Cattani, C. Haar wavelet splines. J. Interdisciplinary Math., 2001, 4, 35-47.

[8.] Cattani, C. Haar wavelets based technique in evolution problems. Proc. Estonian Acad. Sci. Phys. Math., 2004, 53, 45-65.

[9.] Chen, C. F. and Hsiao, C.-H. Haar wavelet method for solving lumped and distributed-parameter systems. IEE Proc. Control Theory Appl., 1997, 144, 87-94.

[10.] Chen, C. F. and Hsiao, C.-H.Wavelet approach to optimising dynamic systems. IEE Proc. Control Theory Appl., 1997, 146, 213-219.

[11.] Hsiao, C.-H. and Wang, W.-J. State analysis of time-varying singular bilinear systems via Haar wavelets. Math. Comput. Simul., 2000, 52, 11-20.

[12.] Hsiao, C.-H. and Wang, W.-J. State analysis of time-varying singular nonlinear systems via Haar wavelets. Math. Comput. Simul., 1999, 51, 91-100.

[13.] Hsiao, C.-H. and Wang, W.-J. Haar wavelet approach to nonlinear stiff systems. Math. Comput. Simul., 2001, 57, 347-353.

[14.] Hsiao, C.-H. Haar wavelet direct method for solving variational problems. Math. Comput. Simul., 2004, 64, 569-585.

[15.] Hsiao, G.-C. and Rathsfeld, A. Wavelet collocation methods for a first kind boundary integral equation in acoustic scattering. Adv. Comput. Math., 2002, 17, 281-308.

[16.] Alpert, B., Beylkin, G., Coifman, R. and Rokhlin, V. Wavelet-like bases for the fast solution of second-kind integral equations. SIAM J. Sci. Comput., 1993, 14, 159-184.

[17.] Vainikko, G., Kivinukk, A. and Lippus, J. Fast solvers of integral equations of the second kind: wavelet methods. J. Complexity, 2005, 21, 243-273.

[18.] Yousefi, S. and Razzaghi, M. Legendre wavelets method for the nonlinear Volterra-Fredholm integral equations. Math. Comput. Simul., 2005, 70, 1-8.

[19.] Maleknejad, K., Aghazadeh, N. and Molapourasl, F. Numerical solution of Fredholm integral equation of the first kind with collocation method and estimation of error bound. Appl. Math. Comput., 2006, 179, 352-359.

[20.] Xiao, J.-Y., Wen, L.-H. and Zhang, D. Solving second kind Fredholm integral equation by periodic wavelet Galerkin method. Appl. Math. Comput., 2006, 175, 508-518.

[21.] Kaneko, H., Noren, R. D. and Novaprateep, B. Wavelet applications to the Petrov-Galerkin method for Hammerstein equations. Appl. Numer. Math., 2003, 45, 255-273.

[22.] Maleknejad, K. and Derili, H. The collocation method for Hammerstein equations by Daubechies wavelets. Appl. Math. Comput., 2006, 172, 846-864.

[23.] Avudainayagam, A. and Vano, C. Wavelet-Galerkin method for integro-differential equations. Appl. Numer. Math., 2000, 32, 247-254.

[24.] Fedorov, M. V. and Chyev, G. N. Wavelet method for solving integral equations of simple liquids. J. Mol. Liq., 2005, 120, 159-162.

[25.] Liang, X.-Z., Liu, M.-C. and Che, X.-J. Solving second kind integral equations by Galerkin methods with continuous orthogonal wavelets. J. Comput. Appl. Math., 2001, 136, 149-161.

[26.] Mahmoudi, Y. Wavelet Galerkin method for numerical solution of nonlinear integral equation. Appl. Math. Comput., 2005, 167, 1119-1129.

[27.] Khellat, F. and Yousefi, S. A. The linear Legendre mother wavelets operational matrix of integration and its application. J. Franklin Inst., 2006, 143, 181-190.

[28.] Maleknejad, K. and Lotfi, T. Expansion method for linear integral equations by cardinal B-spline wavelet and Shannon wavelet bases to obtain Galerkin system. Appl. Math. Comput., 2006, 175, 347-355.

[29.] Maleknejad, K. and Karami, M. Numerical solution of non-linear Fredholm integral equations by using multiwavelets in the Petrov-Galerkin method. Appl. Math. Comput., 2005, 168, 102-110.

[30.] Lepik, U. and Tamme, E. Application of the Haar wavelets for solution of linear integral equations. In 5-10 July 2004, Antalya, Turkey--Dynamical Systems and Applications, Proceedings. 2005, 395-407.

[31.] Lepik, U. and Tamme, E. Solution of nonlinear integral equations via the Haar wavelet method. Proc. Estonian Acad. Sci. Phys. Math., 2007, 56, 17-27.

[32.] Lepik, U. Haar wavelet method for non-linear integro-differential equations. Appl. Math. Comput., 2006, 176, 324-333.

[33.] Maleknejad, K. and Mirzaee, F. Using rationalized Haar wavelet for solving linear integral equations. Appl. Math. Comput., 2005, 160, 579-587.

[34.] Lepik, U. Numerical solution of differential equations using Haar wavelets. Math. Comput. Simul., 2003, 68, 127-143.

[35.] Lepik, U. Numerical solution of evolution equations by the Haar wavelet method. Appl. Math. Comput. (accepted).

[36.] Newland, D. E. An Introduction to Random Vibrations, Spectral and Wavelet Analysis. Longman Scientific and Technical, New York, 1993.

[37.] Kim, B. H., Kim, H. and Park, T. Nondestructive damage evaluation of plates using the multiresolution analysis of two-dimensional Haar wavelet. J. Sound Vibration, 2006, 292, 82-104.

[38.] Kumar , B. V. R. and Mehra, M. Wavelet based preconditioners for sparse linear systems. Appl. Math. Comput., 2005, 171, 203-224.

[39.] Cattani, C. Wavelet analysis of dynamical systems. Electron. Commun. (Kiev), 2002, 17, 115-124.

[40.] Forinash, K. and Willis, C. R. Nonlinear response of the sine-Gordon breather to an a.c. driver. Physica D, 2001, 149, 95-106.

Ulo Lepik

Institute of Applied Mathematics, University of Tartu, Liivi 2, 50409 Tartu, Estonia; ulo.lepik@ut.ee

Let us first explain why we need wavelets at all. A signal x = x(t) is oftenanalysed with the aid of the Fourier transform

F([omega]) = [[integral].sup.+[infinity].sub.-[infinity]] x(t)[e.sup.i[omega]t] dt, (1)

where w denotes the frequency and F--the amplitude.

Two signals x(t) and their Fourier diagrams are plotted in Fig. 1. Although in both cases the Fourier diagrams are quite similar (both have two peaks), the signals are completely different and the time information is lost. This is a serious drawback: the Fourier method analyses the signal over the whole domain and is unable to characterize its behaviour in time.

[FIGURE 1 OMITTED]

Wavelet analysis allows us to represent a function in terms of a set of basic functions, called wavelets, which are localized both in space and time. Here a continuous function [psi](t), called the mother wavelet, is introduced. The corresponding wavelet family is formed by the dilation and translation of the function [psi](t):

[[psi].sub.j,k](t) = [2.sup.-j/2][psi]([2.sup.-j] t - k), (2)

where j and k are non-negative integers.

Depending upon the choice of the mother wavelet [psi](t), various wavelet families are obtained. The wavelets introduced by Daubechies [1] are quite frequently used for solving different problems. The Daubechies mother wavelet is plotted in Fig. 2. These wavelets are differentiable and have a minimum size support.

A shortcoming of the Daubechies wavelets is that they do not have an explicit expression and therefore analytical differentiation or integration is not possible. This complicates the solution of differential equations where integrals of the following type

[[integral].sup.b.sub.a] G (t, [[psi].sub.i,k], d[[psi].sub.i,k]/dt, [d.sup.2][[psi].sub.i,k]/d[t.sup.2], ...)dt (3)

must be computed (G is generally a nonlinear function). For calculating such integrals the conception of connection coefficients is introduced. Calculation of these coefficients is very complicated and must be carried out separately for different types of integrals (see, e.g., [2]). Besides, it can be done only for some simpler types of nonlinearities (mainly for quadratic nonlinearity). This remark holds also for other types of wavelets (such as Symlet, Coiflet, etc. wavelets).

[FIGURE 2 OMITTED]

The wavelet method was first applied to solving differential and integral equations in the 1990s. A survey of early results in this field can be found in [3]. Lately the number of respective papers has greatly increased and it is not possible to analyse them all here, but some are discussed in the following sections.

Due to the complexity of the wavelet solutions, some pessimistic estimates exist. So Strang and Nguyen write in their text-book [4] "... the competition with other methods is severe. We do not necessarily predict that wavelets will win" (p. 394). Jameson [5] writes "... nonlinearities etc., when treated in a wavelet subspace, are often unnecessarily complicated ... There appears to be no compelling reason to work with Galerkin-style coefficients in a wavelet method" (p. 1982).

Obviously attempts to simplify solutions based on the wavelet approach are wanted. One possibility is to make use of the Haar wavelet family.

In 1910 Alfred Haar introduced a function which presents a rectangular pulse pair. After that various generalizations were proposed (a state of the art about Haar transforms can be found in [6]). In the 1980s it turned out that the Haar function was in fact the Daubechies wavelet of order 1. It is the simplest orthonormal wavelet with compact support.

The Haar wavelet family is defined for t [member of] [0, 1] as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

Here [[xi].sub.1] = k/m, [[xi].sub.2] = (k + 0.5)/m, [[xi].sub.3] = (k + 1)m: The integer m = [2.sup.j] (j = 0, 1, ..., J) indicates the level of the wavelet; k = 0, 1, ..., m - 1 is the translation parameter. The maximal level of resolution is J. The index i in (4) is calculated according to the formula i = m + k + 1; in the case m = 1, k = 0 we have i = 2; the maximal value of i is i = 2M = [2.sup.J+1]. It is assumed that the value i = 1 corresponds to the scaling function for which [h.sub.1] = 1 for t [member of] [0, 1]. The first eight Haar functions are plotted in Fig. 3.

An essential shortcoming of the Haar wavelets is that they are not continuous. The derivatives do not exist in the points of discontinuity, therefore it is not possible to apply the Haar wavelets directly to solving differential equations.

There are at least two possibilities of ending this impasse. First, the piecewise constant Haar functions can be regularized with interpolation splines; this technique has been applied by Cattani [7,8]. This greatly complicates the solution process and the main advantage of the Haar wavelets--their simplicity--gets lost.

Another possibility was proposed by Chen and Hsiao [9,10]. They recommended to expand into the Haar series not the function itself, but its highest derivative appearing in the differential equation; the other derivatives (and the function) are obtained through integrations. All these ingredients are then incorporated into the whole system, discretized by the Galerkin or collocation method.

Chen and Hsiao demonstrated the possibilities of their method by solving linear systems of ordinary differential equations (ODEs) and partial differential equations (PDEs). In [10] an optimal control problem with the quadratic performance index is discussed.

In [11,12] Hsiao and Wang applied this method to solving singular bilinear and nonlinear systems. Nonlinear stiff systems were examined in [13].

In [14] Hsiao demonstrated that the Haar wavelet approach is effective also for solving variational problems. Of importance is also the paper [15] by Hsiao and Rathsfeld, in which wavelet collocation methods are applied to the boundary element solution of the first-kind integral equations arising in acoustic scattering.

[FIGURE 3 OMITTED]

The aim of the present paper is to acquaint the readers with the possibilities of the Haar wavelet method in solving integral and evolution equations. The illustrating examples are taken from papers published by the author.

2. INTEGRATION OF HAAR WAVELETS

If we want to integrate differential equations by the Hsiao-Chen method, we have to evaluate the integrals

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

Carrying out these integrations with the aid of (4), we find

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

It is necessary to evaluate these functions for a prescribed J only once and save them in the storage of the computer.

It is convenient to pass to the matrix formulation. For this purpose the interval t [member of] [0, 1] is divided into 2M parts of equal length [DELTA]t = 1/2M, where M = [2.sup.J]. The grid points are

t(l) = (l - 0.5)[DELTA]t, l = 1, 2, ..., 2M. (7)

The Haar coefficient matrix H is defined as H(i, l) = [h.sub.i](l). The integral matrices P[upsilon] have the elements P[upsilon](i, l) = [p.sub.i,[upsilon]](l).

Chen and Hsiao defined the integral matrix in a different way [5]. They introduced the row vector

[h.sub.([mu])](t) = [[h.sub.1](t), [h.sub.2](t), ..., [h.sub.[mu]](t)], (8)

where [mu] = 2M = [2.sup.J+1]. Now we have

[[integral].sup.t.sub.0] [h.sub.([mu])]([tau])d[tau] [approximately equal to] [P.sub.([mu]x[mu])][h.sub.([mu])](t). (9)

The square matrix [P([mu]x[mu])] is called the operational matrix of integration. Chen and Hsiao showed that the following recursive formula holds:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (10)

whereas [P.sub.(1x1)] = 0.5.

Formulae (8)-(10) can be applied to solving first-order differential equations. Higher-order equations have to be reduced to a system of first-order equations, but this increases the order of the Haar matrices and may cause complications of computational character.

3. INTEGRAL EQUATIONS

In this section application of wavelet methods to solution of integral equations is discussed. The following types of integral equations are considered: (i) Fredholm equation

u(x) - [[integral].sup.[beta].sub.[alpha]] K(x, t, u(t))dt = f(x), x, t [member of] [[alpha], [beta]], (11)

(ii) Volterra equation

u(x) - [[integral].sup.x.sub.[alpha]] K(x, t, u(t))dt = f(x), x, t, [member of] [[alpha], [beta]], (12)

(iii) integro-differential equation

au'(x) + bu(x) - [[integral].sup.x.sub.0] K(x, t, u(t), u'(t))dt = f(x),

for u'(0) = [u.sub.0], x, t [member of] [0, [beta]]. (13)

If the kernel K is a linear function in regard of u(t) and u'(t), we have a linear integral equation, otherwise the equation is nonlinear.

In recent years interest in the solution of integral equations by the wavelet methods has greatly increased (we have counted more than 50 papers on this topic). The first paper in this field according to our information is the one by Alpert et al. [16], published in 1993. In the following papers mostly Fredholm or Volterra equations are considered (consult, e.g., [17-20]). The Abel equation is discussed in [18], the Hammerstein equation in [21,22], integro-differential equations are considered in [23]. Different types of wavelets are applied, such as Daubechies [17,20], Coifman [24], Cohen [25], and Legendre [18,26,27] wavelets, and wavelets of B-splines [28]. In most papers linear integral equations are solved; nonlinear equations are considered in [18,20,26,29].

The wavelet equations are discretized by the Galerkin or collocation method. In the case of nonlinear problems we get for the wavelet coefficients a nonlinear system of equations, which can be solved by the Newton method.

Now we pass to the solutions based on the Haar wavelets. Since the Haar wavelets are defined only for the interval [0,1], we must first normalize Eqs (11)-(13). This can be done by the change of the variables

[x.sub.*] = (x - [alpha])/([beta] - [alpha]); [t.sub.*] = (t - [alpha])/([beta] - [alpha]).

In [30] the Haar wavelet method is applied to solving different types of linear integral equations (Fredholm, Volterra, integro-differential, weakly singular integral equations), also the eigenvalue problem is solved. Nonlinear integral equations are considered in [31,32]. It is demonstrated that the recommended method suits well for solving boundary value problems of ODEs.

Maleknejad and Mirzaee [33] applied the Haar wavelet method to solving linear Fredholm integral equations of the second kind.

In the following the Haar wavelet approach is used to solve some types of integral equations. The efficiency of the method is demonstrated by some numerical examples; for getting the error estimates, problems for which the exact solution [u.sub.ex](x) is known are considered. The accuracy of the results is estimated by the error function

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

Computations carried out in [30] showed that the solutions obtained with the aid of the collocation method are simpler than those got by the Galerkin method, besides, their accuracy seems to be better. For this reason in the following the collocation method is used. All calculations were made with the aid of MATLAB programs.

(i) Linear Fredholm integral equation

Let us consider the equation

u(x) - [[integral].sup.1.sub.0] K(x, t)u(t)dt = f(x), x, t [member of] [0, 1] (15)

and seek the solution in the form

u(t) = [2M.summation over (i=1)] [a.sub.i][h.sub.i](t), (16)

where [a.sub.i] are the wavelet coefficients.

Putting this result into (15), we find

[2M.summation over (i=1)] [a.sub.i][h.sub.i](x) - [2M.summation over (i=1)] [a.sub.i][G.sub.i](x) = f(x), (17)

where

[G.sub.i](x) = [[integral].sup.1.sub.0] K(x, t)[h.sub.i](t)dt. (18)

Satisfying (17) in the collocation points (7), we get a linear system of equations for the coefficients [a.sub.i]:

[2M.summation over (i=1)] [a.sub.i][[h.sub.i]([x.sub.l]) - [G.sub.i]([x.sub.l])] = f([x.sub.l]), l = 1, 2, ..., 2M. (19)

More convenient is the matrix form of (16) and (19)

u = aH,

a(H - G) = F, (20)

where

u = {u([t.sub.l])}, F = {f([x.sub.l])},

G = {[G.sub.i]([x.sub.l])}.

Example 1. Let us take K = x + t, f(x) = [x.sup.2]. Equation (15) has the exact solution

[u.sub.ex] = [x.sup.2] - [5.sub.x] - 17/6.

In the present case

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Computations were carried out for different values of J. The errors are shown in Table 1.

(ii) Nonlinear Volterra equation

Let us consider Eq. (12) in the normalized form [alpha] = 0; [beta] = 1. Satisfying it in the collocation points

x(l) = (l - 0.5)[DELTA]t, l = 1, 2, ..., 2M,

we obtain

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

The solution is again sought in the form (16). The system of equations for evaluating the wavelet coefficients is now nonlinear and must be solved with the aid of some numerical method. We have applied for this purpose the Newton method, which brings us to the equation

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (22)

where

[partial derivative]K/[partial derivative][a.sub.i] = [partial derivative]K/[partial derivative]u [h.sub.i](t). (23)

According to (4), in each subinterval [h.sub.i](t) is constant between two collocation points. This enables us to calculate the integrals in (22) analytically. First we carry out the computations for one subintegral ([t.sub.s], [t.sub.s+1]). Summing up all these results, we get the exact values of these integrals for the whole domain [0, x].

If some approximate solution [a.sup.([upsilon]).sub.i] is known, it can be corrected with the aid of

[a.sup.([upsilon]+1).sub.i] = [a.sup.([upsilon]).sub.i] + [DELTA][a.sup.([upsilon]).sub.i], i = 1, 2, ..., 2M.

For details consult the paper [32].

Example 2. Solve the equation

u(x) = 1/2 [[integral].sup.x.sub.0] u(t)u(x - t)dt + 1/2 sin x, x [member of] [0, 1]. (24)

It has the exact solution u(x) = [J.sub.i](x), where the symbol [J.sub.i](x) denotes the Bessel function of order 1.

In this case

K = 1/2 u(t)u(x - t),

[partial derivative]K/[partial derivative][a.sub.i] = 1/2[u(x - t)[h.sub.i](x - t) - u(t)[h.sub.i](t)]. (25)

Error estimates (14) for this case are shown in Table 2.

(iii) Boundary value problems of ODEs

The method of solution described in this section is also very effective for solving boundary value problems if we rewrite the equation in the form of an integro-differential equation.

Consider the equation

u'' = K(x, u, u'), x [member of] [0, 1] (26)

with the boundary conditions u(0) = A, u(1) = B.

By integrating (26) we get the integro-differential equation

u'(x) = [[integral].sup.x.sub.0] K[t, u(t), u'(t)dt + u'(0). (27)

We seek the solution in the form

u'(t) = [2M.summation over (i=1)] [a.sub.i][h.sub.i](t),

u(t) = [2M.summation over (i=1)] [a.sub.i][p.sub.i,1](t) + u(0), (28)

where [p.sub.i,1] denotes the integral (5).

Satisfying the boundary conditions, we get u(1) = [a.sub.1] + u(0) or

[a.sub.1] = B - A. (29)

Since [a.sub.1] is now prescribed, we have [DELTA][a.sub.1] = 0; it means that the matrix of system (22) must be replaced with a reduced matrix for which the first row and column are deleted. Besides, (23) obtains the form

[partial derivative]K/[partial derivative][a.sub.i] = [partial derivative]K/[partial derivative]u [p.sub.i,1](t) + [partial derivative]K/[partial derivative]u' [h.sub.i](t). (30)

For further details consult again [33].

Example 3. Consider the boundary value problem

u'' - uu' = 1/4 [square root of [(1 + x).sup.3]] - 1/2, u(0) = 1, u(1) = [square root of 2], (31)

which has the exact solution

u = [square root of 1 + x].

By integrating (31) we obtain

u'(x) = [[integral].sup.x.sub.0] uu'dt + u'(0) - 1/2 (1 + x) + 1/2[square root of 1 + x] (32)

with the initial condition u(0) = 1. It follows from boundary conditions that [a.sub.1] = u(1) - u(0) = [square root of 2] - 1.

According to (28),

u'(0) = [2M.summation over (i=1)] [a.sub.i][h.sub.i](0), (33)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now

[partial derivative]K/[partial derivative][a.sub.i] = u'(t)[p.sub.i,1](t) + u(t)[h.sub.i](t), i = 2, 3, ..., 2M.

Some numerical calculations were carried out. The errors of the obtained results are presented in Table 3.

4. DIFFERENTIAL EQUATIONS

Papers on solving ODEs and PDEs by the wavelet methods are quite numerous and it is not possible to cite all of them (some information can be found, e.g., in [34,35]). Let us restrict the research area. We shall mainly consider the solutions based on the Haar wavelets. Our analysis has shown [34] that in solving nonlinear ODEs the wavelet method has no preferences over the classical Runge-Kutta method. For this reason we shall concentrate our analysis on PDEs.

An evolution process u = u(x, t), where x is the space coordinate and t is time, is described by a PDE. The simplest equations are the diffusion (or heat) equation

[partial derivative]u/[partial derivative]t = A [[partial derivative].sup.2]u/[partial derivative][x.sup.2] (34)

and the Burgers equation

[partial derivative]u/[partial derivative]t + u [partial derivative]u/[partial derivative]x = [upsilon] [[partial derivative].sup.2]u/[partial derivative][x.sup.2]. (35)

Due to simplicity they have proved to be a touchstone for new numerical methods of solution. An outline of such papers can be found in [34,35].

There are several possibilities for solving PDEs by the wavelet method. First we can introduce the two-dimensional wavelets

[psi](x, t) = [psi](x) [direct sum] [psi](t). (36)

Here [psi](x), [psi](t) are 2M-dimensional vectors, composed of the wavelets (2); [psi](x, t) is a 2M x 2M matrix; the symbol [direct sum] denotes Kronecker multiplication.

This approach has been applied in [36] for harmonic wavelets and in [37] for the Haar wavelets. The shortcoming of it is that we have to deal with matrices of high dimension (e.g., if the resolution level is J = 5, then M = 16 and the wavelet coefficients matrix has 4096 elements!).

An alternative method is applied in [9] to compute an electrical transmission line and in [34] to solve the diffusion equation (34). The solution is sought in the form

[partial derivative]u/[partial derivative]t] = a(x)H(t), (37)

where [partial derivative]u/[partial derivative]t and a are 2M-component vectors, H is the Haar matrix. Integrating (37) in regard of t and differentiating twice with respect to x, we get

u(x, t) = a(x)[P.sub.1](t) + u(x, 0)E,

[[partial derivative].sup.2]u/[partial derivative][x.sup.2] = [d.sup.2]a/d[x.sup.2] [P.sub.1](t) + [d.sup.2]u(x, 0)/d[x.sup.2] E, (38)

where [P.sub.1] is integral matrix and E is a 2M-dimensional unit vector. Inserting these results into (37), we get for a(x) a system of linear ODEs, which is solved analytically with the aid of the matrix exponential function expm. This method is quite effective, but applicable only in linear problems.

Now let us analyse the third method. Here the time interval t [member of] [0, T] is divided into N parts of equal length [DELTA]t and [t.sub.s] = s[DELTA]t, i = 0, 1, 2, ..., N. Following the suggestion of Chen and Hsiao [9], the highest derivative appearing in the evolution equation is developed into the Haar series. This term can be written in the matrix form

[a.sub.t](:)H(:, l), l = 1, 2, ..., 2M. (39)

It is assumed that [a.sub.t] = const for the subinterval t [member of] [[t.sub.s], [t.sub.s+1]]. This allows us to integrate (39) in regard of time t. Integration in regard of the space variable is carried out by the wavelet method. The evolution equation is satisfied by the collocation or Galerkin method. Transfer to the next instant [t.sub.s+1] is performed by using some finite-difference method (e.g. by the Euler or Crank-Nicholson method). This approach is applied in [38] for the Burgers equation and in [34] for the diffusion equation. It should be mentioned that the method has also a weakness. For calculating the wavelet coefficients we must invert some matrices, but they--especially for higher-order equations--are often nearly singular. This is due to the integral matrices [P[upsilon]], which are calculated according to (6). Even in the case of a small value J = 3 the determinants of these matrices are |[P.sub.1]| = 2.7E-20, [absolute value of [P.sub.2]] = 3.4E-49, [absolute value of [P.sub.3]] = 6.6E-81, [absolute value of [P.sub.4]]] = -2.2E-114. Here application of the theory of sparse matrices [38] may be of some help. Another possibility is to divide the space interval x [member of] [[x.sub.min], [x.sub.max]] into N segments and apply the Haar wavelet method to each segment (of course, some continuity conditions must be fulfilled on the boundaries of these segments). This approach was applied in [34,39].

Below we shall consider the third method. The application of this method is illustrated by solving the sine-Gordon equation.

5. SINE-GORDON EQUATION

Let us consider the sine-Gordon equation

1/[L.sup.2] [[partial derivative].sup.2]u/[partial derivative][x.sup.2] - [[partial derivative].sup.2]u/[partial derivative][t.sup.2] = x [member of] [0, 1], t [greater than or equal to] 0, (40)

with the boundary and initial conditions

u(x, 0) = f(x), [partial derivative]u/[partial derivative]t (x, 0) = g(x),

u(0, t) = [psi](x), [partial derivative]u/[partial derivative]x (0, t) = [psi](x).

Here f, g, [psi], [psi] are prescribed functions, L is a constant.

This equation has an analytical solitary wave solution

[u.sub.ex](x, t) = 4 arctan[exp(z)],

z = [alpha](x - [beta]t), [alpha] = L/[square root of (1 - [L.sup.2][[beta].sup.2]]. (41)

Solution of the sine-Gordon equation is discussed in many papers. According to Science Direct, 368 papers have been published on this subject since 1996. Among these we found only one paper [40] in which the wavelet method (Gaussian wavelets) was applied.

Now let us present the solution based on the Haar wavelet method [35]. According to Chen and Hsiao [9], the solution of (40) is sought in the form

[??]''(x, t) = [2M.summation over (i=1)] [a.sub.s][h.sub.i](x), t [member of] [[t.sub.s], [t.sub.s+1]], x [member of] [0, 1]. (42)

This equation is integrated twice in regard of x in the limits [0, x] and in regard of t in the limits [[t.sub.s], [t.sub.s+1]]. By doing this we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

These results are discretized by replacing x [right arrow] [x.sub.l], t [right arrow] [t.sub.s+1]. To simplify the writing of the formulae, the matrix formulation is used and the notations [DELTA]t = [t.sub.s+1] - [t.sub.s], u(l, s) = u([x.sub.l], [t.sub.s]) etc. are introduced. Taking into account the initial conditions (40), Eqs (43) can be put into the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (44)

Inserting these results into (40), we get a linear matrix equation for calculating the wavelet coefficients [a.sub.s](:):

[a.sub.s](:)[P.sub.2](:, l) = 1/[L.sub.2] u''(s, l) + sin u(s, l) - [x.sub.l] [??](s + 1) - [??](s + 1). (45)

Example 4. Computer simulation was carried out for t [member of] [10, 30], L = 20, [beta] = 0.025. If we want to get the classical solitary wave solution, we must take

f(x) = 4 arctan[exp([alpha]x)],

g(x) = [alpha]V ([alpha]x),

[phi](t) = 4 arctan[exp(-[alpha][beta]t)],

[psi](t) = -[alpha][beta]V (-[alpha][beta]t), (46)

where

V (z) = 4[e.sup.z]/1 + [e.sup.2z].

The accuracy of our approach is estimated by the error function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (47)

The calculations show that the function v(t) increases monotonically, therefore [v]([t.sub.max]) is taken as the error estimate. Some computer results are presented in Table 4 and Fig. 4.

[FIGURE 4 OMITTED]

It follows from Table 4 that already in the cases J = 4 or J = 5 we get the results which visually coincide with the exact solution. The time step was taken as [DELTA]t = 0.005 or [DELTA]t = 0.001; further decrease in [DELTA]t gave no considerable effect. If the accuracy of these results is insufficient, more precise results could be obtained by the segmentation method proposed in [34].

We would like to draw attention to Fig. 4b. It follows from here that most of the wavelet coefficients are zero (or near to zero) and the number of significant coefficients is quite small. This is one of the reasons for rapid convergence of the Haar wavelet series.

6. CONCLUSIONS

The main goal of this paper was to demonstrate that the Haar wavelet method is a powerful tool for solving different types of integral equations and partial differential equations. The method with far less degrees of freedom and with smaller CPU time provides better solutions than classical ones.

The main advantage of this method is its simplicity and small computation costs, resulting from the sparsity of the transform matrices and the small number of significant wavelet coefficients. The method is also very convenient for solving the boundary value problems, since the boundary conditions are taken care of automatically.

ACKNOWLEDGEMENT

Financial support from the Estonian Science Foundation (grant No. 6697) is gratefully acknowledged.

Received 13 November 2006

REFERENCES

[1.] Daubechies, I. Orthonormal bases of compactly supported wavelets. Comm. Pure Appl. Math., 1988, 41, 909-996.

[2.] Chen, M.-Q., Hwang, C. and Shin, Y.-P. The computation of wavelet-Galerkin approximation on a bounded interval. Int. J. Numer. Methods Eng., 1996, 39, 2921-2944.

[3.] Dahmen, W., Kurdila, A. J. and Oswald, P. (eds). Multiscale Wavelet Methods for Partial Differential Equations. Academic Press, San Diego, 1997.

[4.] Strang, G. and Nguyen, T. Wavelets and Filter Banks. Wellesley-Cambridge Press, Wellesley (MA), 1997.

[5.] Jameson, L. A wavelet-optimized, very much order adaptive grid and order numerical method. SIAM J. Sci. Comput., 1998, 19, 1980-2013.

[6.] Stankovic, R. S. and Falkowski, B. J. The Haar wavelet transform: its status and achievements. Comput. Electr. Eng., 2003, 29, 25-44.

[7.] Cattani, C. Haar wavelet splines. J. Interdisciplinary Math., 2001, 4, 35-47.

[8.] Cattani, C. Haar wavelets based technique in evolution problems. Proc. Estonian Acad. Sci. Phys. Math., 2004, 53, 45-65.

[9.] Chen, C. F. and Hsiao, C.-H. Haar wavelet method for solving lumped and distributed-parameter systems. IEE Proc. Control Theory Appl., 1997, 144, 87-94.

[10.] Chen, C. F. and Hsiao, C.-H.Wavelet approach to optimising dynamic systems. IEE Proc. Control Theory Appl., 1997, 146, 213-219.

[11.] Hsiao, C.-H. and Wang, W.-J. State analysis of time-varying singular bilinear systems via Haar wavelets. Math. Comput. Simul., 2000, 52, 11-20.

[12.] Hsiao, C.-H. and Wang, W.-J. State analysis of time-varying singular nonlinear systems via Haar wavelets. Math. Comput. Simul., 1999, 51, 91-100.

[13.] Hsiao, C.-H. and Wang, W.-J. Haar wavelet approach to nonlinear stiff systems. Math. Comput. Simul., 2001, 57, 347-353.

[14.] Hsiao, C.-H. Haar wavelet direct method for solving variational problems. Math. Comput. Simul., 2004, 64, 569-585.

[15.] Hsiao, G.-C. and Rathsfeld, A. Wavelet collocation methods for a first kind boundary integral equation in acoustic scattering. Adv. Comput. Math., 2002, 17, 281-308.

[16.] Alpert, B., Beylkin, G., Coifman, R. and Rokhlin, V. Wavelet-like bases for the fast solution of second-kind integral equations. SIAM J. Sci. Comput., 1993, 14, 159-184.

[17.] Vainikko, G., Kivinukk, A. and Lippus, J. Fast solvers of integral equations of the second kind: wavelet methods. J. Complexity, 2005, 21, 243-273.

[18.] Yousefi, S. and Razzaghi, M. Legendre wavelets method for the nonlinear Volterra-Fredholm integral equations. Math. Comput. Simul., 2005, 70, 1-8.

[19.] Maleknejad, K., Aghazadeh, N. and Molapourasl, F. Numerical solution of Fredholm integral equation of the first kind with collocation method and estimation of error bound. Appl. Math. Comput., 2006, 179, 352-359.

[20.] Xiao, J.-Y., Wen, L.-H. and Zhang, D. Solving second kind Fredholm integral equation by periodic wavelet Galerkin method. Appl. Math. Comput., 2006, 175, 508-518.

[21.] Kaneko, H., Noren, R. D. and Novaprateep, B. Wavelet applications to the Petrov-Galerkin method for Hammerstein equations. Appl. Numer. Math., 2003, 45, 255-273.

[22.] Maleknejad, K. and Derili, H. The collocation method for Hammerstein equations by Daubechies wavelets. Appl. Math. Comput., 2006, 172, 846-864.

[23.] Avudainayagam, A. and Vano, C. Wavelet-Galerkin method for integro-differential equations. Appl. Numer. Math., 2000, 32, 247-254.

[24.] Fedorov, M. V. and Chyev, G. N. Wavelet method for solving integral equations of simple liquids. J. Mol. Liq., 2005, 120, 159-162.

[25.] Liang, X.-Z., Liu, M.-C. and Che, X.-J. Solving second kind integral equations by Galerkin methods with continuous orthogonal wavelets. J. Comput. Appl. Math., 2001, 136, 149-161.

[26.] Mahmoudi, Y. Wavelet Galerkin method for numerical solution of nonlinear integral equation. Appl. Math. Comput., 2005, 167, 1119-1129.

[27.] Khellat, F. and Yousefi, S. A. The linear Legendre mother wavelets operational matrix of integration and its application. J. Franklin Inst., 2006, 143, 181-190.

[28.] Maleknejad, K. and Lotfi, T. Expansion method for linear integral equations by cardinal B-spline wavelet and Shannon wavelet bases to obtain Galerkin system. Appl. Math. Comput., 2006, 175, 347-355.

[29.] Maleknejad, K. and Karami, M. Numerical solution of non-linear Fredholm integral equations by using multiwavelets in the Petrov-Galerkin method. Appl. Math. Comput., 2005, 168, 102-110.

[30.] Lepik, U. and Tamme, E. Application of the Haar wavelets for solution of linear integral equations. In 5-10 July 2004, Antalya, Turkey--Dynamical Systems and Applications, Proceedings. 2005, 395-407.

[31.] Lepik, U. and Tamme, E. Solution of nonlinear integral equations via the Haar wavelet method. Proc. Estonian Acad. Sci. Phys. Math., 2007, 56, 17-27.

[32.] Lepik, U. Haar wavelet method for non-linear integro-differential equations. Appl. Math. Comput., 2006, 176, 324-333.

[33.] Maleknejad, K. and Mirzaee, F. Using rationalized Haar wavelet for solving linear integral equations. Appl. Math. Comput., 2005, 160, 579-587.

[34.] Lepik, U. Numerical solution of differential equations using Haar wavelets. Math. Comput. Simul., 2003, 68, 127-143.

[35.] Lepik, U. Numerical solution of evolution equations by the Haar wavelet method. Appl. Math. Comput. (accepted).

[36.] Newland, D. E. An Introduction to Random Vibrations, Spectral and Wavelet Analysis. Longman Scientific and Technical, New York, 1993.

[37.] Kim, B. H., Kim, H. and Park, T. Nondestructive damage evaluation of plates using the multiresolution analysis of two-dimensional Haar wavelet. J. Sound Vibration, 2006, 292, 82-104.

[38.] Kumar , B. V. R. and Mehra, M. Wavelet based preconditioners for sparse linear systems. Appl. Math. Comput., 2005, 171, 203-224.

[39.] Cattani, C. Wavelet analysis of dynamical systems. Electron. Commun. (Kiev), 2002, 17, 115-124.

[40.] Forinash, K. and Willis, C. R. Nonlinear response of the sine-Gordon breather to an a.c. driver. Physica D, 2001, 149, 95-106.

Ulo Lepik

Institute of Applied Mathematics, University of Tartu, Liivi 2, 50409 Tartu, Estonia; ulo.lepik@ut.ee

Table 1. Error function [e.sub.J] for Eq. (15) J 2M [e.sub.J] 2 8 7.2E-2 3 16 1.7E-2 4 32 4.3E-3 5 64 1.3E-3 Table 2. Error estimates for Eq. (24) J 2M [e.sub.J] 2 8 1.2E-3 3 16 1.3E-3 4 32 7.9E-4 5 64 4.3E-4 6 128 2.2E-4 Table 3. Error estimates for Eq. (30) J 2M [e.sub.J] 2 8 1.1E-3 3 16 3.2E-4 4 32 8.6E-5 5 64 2.2E-5 6 128 5.6E-6 Table 4. Error estimate v([t.sub.max]) for Example 4 J 2M v v([t.sub.max]) [DELTA]t = 0.005 [DELTA]t = 0.001 4 32 0.051 0.038 5 64 0.018 0.009 6 128 0.0096 0.0036