Gregory Beylkin

James Keiser

On the Adaptive Numerical Solution of Nonlinear Partial Differential Equations in Wavelet Bases

requires the development of new algorithms, which are introduced in this paper.

This work develops fast and adaptive algorithms for numerically solving nonlinear partial differential equations of the form u t ϭ Any wavelet-expansion approach to solving differential L u ϩ N f (u), where L and N are linear differential operators and equations is essentially a projection method. In a projection f (u) is a nonlinear function. These equations are adaptively solved method the goal is to use the fewest number of expansion by projecting the solution u and the operators L and N into a coefficients to represent the solution since this leads to wavelet basis. Vanishing moments of the basis functions permit a sparse representation of the solution and operators. Using these efficient numerical computations. The number of coeffisparse representations fast and adaptive algorithms that apply opercients required to represent a function expanded in a Fouators to functions and evaluate nonlinear functions, are developed rier series (or similar expansions based on the eigenfuncfor solving evolution equations. For a wavelet representation of the tions of a differential operator) depends on the most solution u that contains N s significant coefficients, the algorithms singular behavior of the function. We are interested in update the solution using O(N s ) operations. The approach is applied to a number of examples and numerical results are given.

solutions of partial differential equations that have regions of smooth, nonoscillatory behavior interrupted by a number of well-defined localized shocks or shock-like structures. Therefore, expansions of these solutions, based upon

INTRODUCTION

the eigenfunctions of differential operators, require a large number of terms due to the singular regions. Alternately, This paper is concerned with the fast, adaptive numerical a localized representation of the solution, typified by frontsolution of nonlinear partial differential equations having tracking or adaptive grid methods, may be employed in solutions which exhibit both smooth and shock-like behavorder to distinguish between smooth and shock-like beior. The algorithms we describe take advantage of the havior. fact that wavelet expansions may be viewed as a localized

In this paper we use wavelet expansions in the develop-Fourier analysis with multiresolution structure that ''automent of adaptive numerical algorithms. Let the wavelet matically'' distinguishes between smooth and shock-like transform of a function consist of N s significant coefficients behavior. In smooth regions few wavelet coefficients are (those coefficients of size greater than some threshold needed and, in singular regions, large variations in the Ͼ 0) concentrated near shock-like structures. Our goal is function require more wavelet coefficients. The theoretical to design fully adaptive algorithms that perform numerical analysis of such functions by wavelet methods is well uncomputations in O(N s ) operations, using only the signifiderstood [17, [START_REF] Chui | An Introduction to Wavelets[END_REF][START_REF] Meyer | Wavelets and Operators[END_REF]. Additionally, there have been a numcant wavelet coefficients. In other words, we will look for ber of investigations into the use of wavelet expansions for a general ''spectral'' approach that has the desirable propnumerically computing solutions of differential equations erties of specialized adaptive algorithms. The resulting al- [START_REF] Schult | Using wavelets to solve the Burgers[END_REF][START_REF] Liandrat | Nonlinear Partial Differential Equations using the Wavelet Approach. Wavelets and Their Applications[END_REF][38]. However, these numerical approaches are simply gorithmic complexity of our approach is then proportional projection methods which do not address the important to the number of significant coefficients in the wavelet computational question of adaptively updating solutions expansions of functions and operators. of differential equations. Using wavelet expansions of func-We also note that in wavelet coordinates differential tions and operators for fast, adaptive numerical purposes operators may be preconditioned by a diagonal matrix, [22]. Moreover, a large class of operators, namely Caldero ´n-Zygmund and pseudo-differential operators, are sparse in wavelet bases. These observations make a good case for developing new numerical algorithms which

u t ϭ L u ϩ N f (u) (2.
1) take advantage of these properties.

We develop two new algorithms for computing solutions with the initial condition of partial differential equations, namely the adaptive application of operators to functions and the adaptive pointwise u(x, 0) ϭ u 0 (x), 0 Յ x Յ 1, (2.2) product of functions. These algorithms are necessary ingredients of any fast, adaptive numerical scheme for computand the periodic boundary condition ing solutions of partial differential equations. The algorithm for adaptively multiplying operators and functions u(0, t) ϭ u(1, t), 0 Յ t Յ T.

(2.3) is based on a vanishing-moment property associated with the nonstandard form representation of a class of opera-We explicitly separate the evolution Eq. (2.1) into a linear tors, which includes differential operators and Hilbert part, L u, and a nonlinear part, N f (u), where the operators transforms. We will use this property to develop a generic, L and N are constant-coefficient differential operators efficient, adaptive algorithm for applying differential operthat do not depend on time t. The function f (u) is typically ators to functions using only O(N s ) significant wavelet coefnonlinear, e.g., f (u) ϭ u p . ficients. We have also developed an adaptive algorithm for Examples of evolution Eq. (2.1) in 1 ϩ 1 dimensions computing the pointwise product of functions, again using include reaction-diffusion equations, e.g., only O(N s ) significant wavelet coefficients.

This paper is outlined as follows. In Section 2 we identify u t ϭ u xx ϩ u p , p Ͼ 1, Ͼ 0, (2.4) a class of partial differential equations for which we develop our methods. We use the semigroup method to reequations describing the buildup and propagation of place the differential equation by a nonlinear integral shocks, e.g., Burgers' equation equation and introduce a procedure for approximating the integral to any order of accuracy. In Section 3 we are u t ϩ uu x ϭ u xx , Ͼ 0 (2.5) concerned with the construction of and calculations with the operators appearing in the quadrature formulas derived in Section 2. Specifically, we describe a method for [START_REF] Burgers | A mathematical model illustrating the theory of turbuof matrices and fractals[END_REF], and equations having special soliton solutions, e.g., the constructing a wavelet representation of these operators, Korteweg-de Vries equation derive the vanishing-moment property of these operators and describe a fast, adaptive algorithm for applying these u t ϩ Ͱuu x ϩ ͱu xxx ϭ 0, (2.6) operators to functions expanded in a wavelet basis. In Section 3 we also provide a brief review of the notation and where Ͱ, ͱ are constant [35,1]. Finally, a simple example terminology associated with the wavelet representations of of Eq. (2.1) is the classical diffusion (or heat) equation functions and operators. In Section 4 we introduce a new adaptive algorithm for computing the pointwise product u t ϭ u xx , Ͼ 0.

(2.7) of functions expanded in a wavelet basis. In Section 5 we illustrate the use of these algorithms by providing the re-Remark. Although we do not address multidimensults of numerical experiments and comparing them with sional problems in this paper, we note that the Navierthe exact solutions. Finally, in Section 6 we draw a number Stokes equations may also be written in the form (2.1). of conclusions based on our results and indicate directions Consider of further investigation.

u t ϩ [u и ٌu ϩ ٌ(u и u)] ϭ ٌ 2 u Ϫ ٌp, (2.8)

SEMIGROUP APPROACH AND QUADRATURES

In this section we use the semigroup approach to recast where a partial differential equation as a nonlinear integral equation in time. We then approximate the integrals to arbitrary div u ϭ 0 (2.9) orders of accuracy by quadratures with operator valued coefficients. These operators have wavelet representations and p denotes the pressure. Applying divergence operator with a number of desirable properties described in Sections to both sides of (2.8) and using (2.9), we obtain 3.2 and 3.3.

⌬p ϭ f (u), (2.10) 2.1. The Model Equation

We consider the problem of computing numerical solu-where f (u) ϭ Ϫ ٌ[u и ٌu ϩ ٌ(u и u)] is a nonlinear function of u. Equation (2.1) is formally obtained by setting tions of Lu ϭ ٌ 2 u (2.11) disadvantage in that the linear and nonlinear contributions must be added in the same domain, either the physical space or the Fourier space. For solutions which exhibit and shock-like solutions such transformations between the domains is costly, whereas the transform into the wavelet Nu ϭ Ϫ [u и ٌu ϩ ٌ(u и u)] Ϫ ٌ(⌬ Ϫ1 f(u)). (2.12) domain is much less expensive due to the locality of the wavelet transform. This difficulty becomes significant when The term ⌬ Ϫ1 f (u) is an integral operator which introduces a long-range interaction and has a sparse representation one attempts to compute solutions of differential equations in multiple dimensions. We note that our wavelet approach in wavelet bases. A one-dimensional model that may be thought of as a prototype for the Navier-Stokes equation is is comparable to spectral methods in their accuracy and parallels general adaptive grid approaches in the automatic placement of significant wavelet coefficients in regions of u t ϭ H (u)u, (2.13) large gradients. where H (и) is the Hilbert transform (see [14]). The presence of the Hilbert transform in (2.13) introduces a long-2.2. The Semigroup Approach range interaction which models that found in the Navier-

The semigroup approach is a well-known analytical tool Stokes equations. Even though in this paper we develop which is used to convert partial differential equations to algorithms for one-dimensional problems, we take special nonlinear integral equations and to obtain estimates associcare that they generalize properly to several dimensions ated with the behavior of their solutions (see, e.g., [4,13]). so that we can address these problems in the future.

The solution of the initial value problem (2.1) is given by Several numerical techniques have been developed to compute approximate solutions of equations such as (2.1). These techniques include finite-difference, pseudo-spec-

u(x, t) ϭ e (tϪt 0 )L u 0 (x) ϩ ͵ t t 0 e (tϪ)L N f (u(x, )) d. (2.14)
tral, and adaptive grid methods (see, e.g., [START_REF] Dahlquist | Numerical Methods[END_REF]35]). An important step in solving Eq. (2.1) by any of these methods Expressing solutions of (2.1) in the form (2.14) is useful is the choice of time discretization. Explicit schemes (which for proving existence and uniqueness of solutions and comare easiest to implement) may require prohibitively small puting estimates of their magnitude, verifying dependence time steps (usually because of diffusion terms in the evoluon initial and boundary data, as well as performing asymption equation). On the other hand, implicit schemes allow totic analysis of the solution; see, e.g., [13]. for larger time steps but they require solving a system We are interested in using Eq. (2.14) as a starting point of equations at each time step and, for this reason, are for an efficient numerical algorithm. As far as we know, somewhat more difficult to implement in an efficient manthe semigroup approach has had limited use in numerical ner. In our approach we have used an implicit time intecalculations. A significant difficulty in designing numerical grator which is described below. We note that there are algorithms based directly on the solution (2.14) is that preconditioners available for the wavelet representation of the operators appearing in (2.14) are not sparse (i.e., the differential operators used in implicit numerical schemes, matrices representing these operators are dense). We show although we do not discuss this subject in this paper.

in Sections 3.2 and 3.3 that in the wavelet system of coordi-The main difficulty in computing solutions of equations nates these operators are sparse and have the desired proplike (2.1) is the resolution of shock-like structures. Straighterties for fast, adaptive numerical algorithms. forward refinement of a finite-difference scheme easily becomes computationally excessive. The specialized fronttracking or adaptive grid methods require some criteria to 2.3. Quadratures perform local grid refinement. Usually in such schemes As it follows from (2.14), we have to consider approxithese criteria are chosen in an ad hoc fashion (especially mating integrals of the form in multiple dimensions) and are generally based on the amplitudes or local gradients in the solution.

The pseudo-spectral method usually splits the evolution

I(x, t) ϭ ͵ t t 0 e (tϪ)L N f (u(x, )) d.
(2.15) equation into linear and nonlinear parts and updates the solution by adding the linear contribution, calculated in the Fourier space, and the nonlinear contribution, calculated in As mentioned earlier, the differential operator N is assumed to be independent of t and the function f (u) is the physical space [34,35]. Pseudo-spectral schemes have the advantages that they are easy to understand analyti-nonlinear. For example, in the case of Burgers' equation N ϭ Ѩ/Ѩx and f (u) ϭ u 2 , so that N f (u) ϭ uu x appears cally, spectrally accurate, and relatively straightforward to implement. However, pseudo-spectral schemes have a as products of u and its derivative. In the case of quadratic nonlinearity we seek approximations to integrals of the For m ϭ 1, we approximate (2.19) by form 2 ), (2.20)

I(t) ϭ O L,1 (u(t 0 )u x (t 0 ) ϩ u(t 1 )u x (t 1 )) ϩ O((⌬t)
I(t) ϭ ͵ t t 0 e (tϪ)L u()v() d,
(2.16) or 2 ), (2.21) where we have suppressed the explicit x-dependence of u(x, t). In order to derive an approximation to this integral, where we partition the interval of integration [t 0 , t] into m equal subintervals with grid points at t i ϭ t 0 ϩ i ⌬t, for i ϭ 0, O L,m ϭ (e m⌬tL Ϫ I)L Ϫ1 (2.22) 1, ..., m, and we denote u(t i ) and v(t i ) by u i and v i , respectively.

I(t) ϭ O L,1 (u(t 0 )u x (t 1 ) ϩ u(t 1 )u x (t 0 )) ϩ O((⌬t)
and where I is the identity operator. Note that (2.20) is Remark. We do not address adaptive time integration equivalent to the standard trapezoidal rule. For m ϭ 2 our in this paper, but we note that it can be accommodated procedure yields an analogue of Simpson's rule, by our algorithms.

We seek an approximation to (2.16) of the form 3 ), (2.23)

I(t) ϭ 2 iϭ0 c i,i u(t i )u x (t i ) ϩ O((⌬t)
I(t) ϭ I ˆ(t) ϩ O(⌬t mϩ1 ), (2.17) where where c 0,0 ϭ O L,2 Ϫ L , (2.24) c 1,1 ϭ O L,2 , (2.25)
I ˆ(t) ϭ m i, jϭ0 c i, j u i v j (2.18) c 2,2 ϭ O L,2 ϩ L . (2.26)
Higher order quadratures are accommodated by this proand where the coefficients c c, j are time-independent, operacedure by considering m Ͼ 2. Detailed derivation and tor-valued functions of the operator L . Observe that we stability analysis of these schemes is outside the scope of have included in (2.18) cross terms of the form u i v j , i ϶ j;

this paper and we refer to [43] for details. typically, quadrature approximations only involve products u i v i , e.g., the trapezoidal rule. We would like to use

WAVELET REPRESENTATIONS OF

the fewest number of terms of the form u i v j in (2.18). We

OPERATOR FUNCTIONS

reduce the number of such terms (from (m ϩ 1) 2 to m ϩ 1) by reducing the number of nonzero operator coef-In this section we are concerned with the construction ficients c i, j . The coefficients c i, j are determined by comparof and calculations with the nonstandard form (NS-form) ing (2.17) and (2.18) with a scheme of known order of of the operator functions (see, e.g., (2.16)). We begin by accuracy. One such comparison scheme is constructed ussetting our notation and refer to Appendix A for details. ing Lagrange polynomial approximations in t of the func-

We then show how to compute the NS-form of the operator tions u(t) and v(t). The coefficients c i, j in (2.18) are then functions and establish the vanishing-moment property of determined by straightforward expansion techniques. This the wavelet representation of these operators. Finally, we leads to a system of equations for determining the operator describe a fast, adaptive algorithm for applying operators coefficients c i, j that, in general, has more than one solution.

to functions in the wavelet system of coordinates. We then choose a solution of this system of equations that consists of m ϩ 1 nonzero coefficients c i, j . Substituting 3.1. Notation these c i, j into Eq. (2.17) yields an approximation to (2.16) which is O((⌬t) mϩ1 ) accurate and involves m ϩ 1 terms of We begin by setting our notation associated with multithe form u i v j . resolution analysis and expansions of functions and opera-Applying this procedure to Burgers' equation (2.5), we tors into a wavelet basis (see also Appendix A and [18, are led to approximate [START_REF] Chui | An Introduction to Wavelets[END_REF][START_REF] Wickerhauser | On the fast algorithm for multiplication of functions ware[END_REF]). We consider a multiresolution analysis (MRA) of L 2 )ޒ( as an infinite chain of subspaces

I(t) ϭ ͵ t t 0 e (tϪ)L u()u x () d. (2.19) и и и ʚ V 2 ʚ V 1 ʚ V 0 ʚ V Ϫ1 ʚ V Ϫ2 ʚ и и и . (3.1)
As usual, we define an associated sequence of subspaces J Յ n. Once again we refer the reader to Appendix A and [18, [START_REF] Chui | An Introduction to Wavelets[END_REF][START_REF] Wickerhauser | On the fast algorithm for multiplication of functions ware[END_REF] for additional introductory and background ma-W j as the orthogonal complements of V j in V jϪ1 , terial. V jϪ1 ϭ V j W j .

(3.2) 3.2. The Nonstandard Form of Operator Functions We denote by (и) the scaling function and (и) the wave-In this section we construct the NS-forms of analytic let. The family of functions ͕ j,k (x) ϭ 2 Ϫj/2 (2 Ϫj x Ϫ k)͖ kʦޚ functions of the differential operator Ѩ x . Following [25, 21] forms an orthonormal basis of V j and ͕ j,k (x) ϭ we introduce two approaches for computing the NS-forms 2 Ϫj/2 (2 Ϫj x Ϫ k)͖ kʦޚ , an orthonormal basis of W j . of operator functions: (i) compute the projection of the We consider representations of operators in the NSoperator function on V 0 , form [25,21]. The NS-form of an operator T is obtained by expanding T in the ''telescopic'' series,

P 0 f (Ѩ x )P 0 , (3.7)
T ϭ jʦޚ (Q j TQ j ϩ Q j TP j ϩ P j TQ j ), (3.3) or, (ii) compute the function of the projection of the operator, where P j and Q j are projection operators on subspaces V j and W j , respectively. The NS-form of T is, thus, repref (P 0 Ѩ x P 0 ). (3.8) sented by the set of operators T ϭ ͕A j , B j , ⌫ j ͖ jʦޚ , (3.4) The difference between these two approaches depends on how well ͉ ˆ()͉ 2 acts as a cutoff function, where (x) is the where the operators A j , B j , and ⌫ j act on subspaces V j and scaling function associated with a wavelet basis. It might W j as be convenient to use either (3.7) or (3.8) in applications.

The operator functions we are interested in are those

A j ϭ Q j TQ j : W j Ǟ W j ,
appearing in solutions of the partial differential Eq. (2.1). For example, using (2.14) with (2.21), solutions of Burgers' B j ϭ Q j TP j : V j Ǟ W j , (3.5) equation can be approximated to order (⌬t) 2 by ⌫ j ϭ P j TQ j : W j Ǟ V j .

u(x, t ϩ ⌬t) For numerical purposes we define a ''finest'' scale, ϭ e ⌬tL u(x, t) (3.9) j ϭ 0, and a ''coarsest'' scale, j ϭ J, such that the infinite chain (3.1) is restricted to

Ϫ O L,1 [u(x, t)Ѩ x u(x, t ϩ ⌬t) ϩ u(x, t ϩ ⌬t)Ѩ x u(x, t)], V J ʚ V JϪ1 ʚ и и и ʚ V 0 . (3.6)
where L ϭ Ѩ 2 x and O L,1 is given by (2.22). Therefore, we are interested in constructing the NS-forms of the opera-We also consider a periodized version of the multiresolutor functions, tion analysis that is obtained if we consider periodic functions. This periodization is the simplest (but not the most e ⌬tL (3.10) efficient or elegant) way to consider the multiresolution analysis of a function on an interval. The problem with periodization is that we might introduce an artificial singuand larity at the boundary. A more elegant approach would use wavelets on the interval, [17], or multiwavelets, [19].

O L,1 ϭ (e ⌬t L Ϫ I)L Ϫ1 , (3.11) We choose to consider the periodization described here since it is the easiest way to describe our adaptive algorithms and our approach does not change substantially if for example. In computing solutions of (2.1) (via, e.g., (3.9)), we precompute the NS-forms of the operator func-we use other bases. We will therefore consider functions having projections on V 0 which are periodic of period tions and apply them as necessary.

We note that if the operator function is homogeneous N ϭ 2 n , where N is the dimension of V 0 . With a slight abuse of notation we will denote these periodized subspaces also of degree m (e.g., m ϭ 1 and 2 for the first and second derivative operators), then the coefficients appearing in by V j and W j . We can then view the space V 0 as consisting of 2 n ''samples'' or lattice points and each space V j and the NS-form at different scales are simply related (see (A.14) and (A. 19)), W j as consisting of 2 nϪj lattice points, for j ϭ 1, 2, ...,

Ͱ j l ϭ 2 Ϫmj Ͱ 0 l , ͉ ˆ()͉ 2 Ͻ
for ͉͉ Ͼ for some Ͼ 0. Therefore, Eq. (3.17) is approximated to within by ͱ j l ϭ 2 Ϫmj ͱ 0 l , (3.12)

Ͳ j l ϭ 2 Ϫmj Ͳ 0 l , g ˜() ϭ K kϭϪK f (Ϫi2 Ϫj ( ϩ 2ȏk))͉ ˆ( ϩ 2ȏk)͉ 2 (3.20) s j l ϭ 2 Ϫmj s 0 l ,
for some K. Using (3.20) in place of g() in (3.18), we On the other hand, if the operator function is not homogeobtain an approximation to the coefficients s j l , neous then we compute s 0 k,kЈ via (A.14) and compute the coefficients Ͱ j k,kЈ , ͱ j k,kЈ , and Ͳ j k,kЈ via equations (A.20) for each scale, j ϭ 1, 2, ..., J Յ n. We note that if the operator

s ˜j l ϭ 1 N NϪ1 nϭ0 g ˜( n )e i n l .
(3.21) function is a convolution then the formulas for s 0 kϪkЈ are considerably simplified (see [21]).

We first describe computing the NS-form of an operator The coefficients s ˜j l are computed by applying the FFT to function f (Ѩ x ) by projecting the operator function into the the sequence ͕g ˜( n )͖, computed via (3.20).

wavelet basis via (3.7). To compute the coefficients

In order to compute the NS-form of an operator function via (3.8), we use the DFT to diagonalize the differential operator Ѩ x and apply the spectral theorem to compute (3.13) the operator functions. Starting with the wavelet representation of Ѩ x on V 0 (see Section A.2 or [21]) of the discretizalet us consider tion of Ѩ x , we write the eigenvalues explicitly as

s j k,kЈ ϭ 2 Ϫj ͵ ϩȍ Ϫȍ (2 Ϫj x Ϫ k) f (Ѩ x )(2 Ϫj x Ϫ kЈ) dx,
f (Ѩ x )(2 Ϫj x Ϫ kЈ) (3.14) k ϭ s 0 ϩ L lϭ1 (s l e 2ȏi(kl/N) ϩ s Ϫl e Ϫ2ȏi(kl/N) ), (3.22) ϭ 1 ͙2ȏ ͵ ȍ Ϫȍ f (Ϫi 2 Ϫj ) ˆ()e ikЈ e i2 Ϫj x d,
where the wavelet coefficients of the derivative, s l ϭ s 0 l , where ˆ() is the Fourier transform of (x), are defined by (A.14). Since

f ( A ) ϭ F f (⌳) F Ϫ1 , (3.23) ˆ() ϭ 1 ͙2ȏ ͵ ȍ Ϫȍ (x)e ix dx. (3.15)
where ⌳ is a diagonal matrix and F is the Fourier transform Substituting (3.14) into (3.13) and noting that s j k,kЈ ϭ s j kϪkЈ , we compute f ( k ) and apply the inverse Fourier transform we arrive at to the sequence f ( k ),

s j l ϭ ͵ ϩȍ Ϫȍ f (Ϫi 2 Ϫj )͉ ˆ()͉ 2 e il d, (3.16) s 0 l ϭ N kϭ1 f ( k )e 2ȏi((kϪ1)(lϪ1)/N) , (3.24)
We evaluate (3.16) by setting to arrive at the wavelet coefficients s 0 l . The remaining elements of the NS-form are then recursively computed using (3.17) Eqs. (A.19). or

s j l ϭ ͵ 2ȏ 0 kʦޚ f (Ϫi2 Ϫj ( ϩ 2ȏk))͉ ˆ( ϩ 2ȏk)͉ 2 e il d,

Vanishing Moments of the B-Blocks

We now establish the vanishing-moment property of the s j l ϭ ͵ 2ȏ 0 g()e il d, (3.18) B-blocks of the NS-form representation of functions of a differential operator and the Hilbert transform. We note that a similar result also holds for the B-blocks of some where classes of pseudo-differential operators; see, e.g., [START_REF] Meyer | Wavelets and numerical split-step method: A global adaptive scheme[END_REF]. Additionally, we note that the results of this section do not

g() ϭ kʦޚ f (Ϫi2 Ϫj ( ϩ 2ȏk))͉ ˆ( ϩ 2ȏk)͉ 2 . (3.19)
require compactly supported wavelets. These results are used to design an adaptive algorithm for multiplying the NS-form of an operator and the wavelet expansion of a We now observe that for a given accuracy the function ͉ ˆ()͉ 2 acts as a cutoff function in the Fourier domain, i.e., function. for m ϭ 0, 1, 2, ..., M Ϫ 1 and j ϭ 1, 2, ..., J.

Proof. See Appendix A.3.

FIG. 1.

For the operators considered in Section 3.3 the vanishingmoment property of the rows of the B-block yields a sparse result (up LEMMA 2. Under the conditions of Lemma 1, the to a given accuracy ) when applied to a smooth and dense vector ͕s j ͖. Since the averages are ''smoothed'' versions of the function for 0 Յ m Յ M Ϫ 1 and j ϭ 1, 2, ..., J.

B-blocks of the NS-form of the Hilbert transform

( H f )(x) ϭ 1 ȏ p.v. ͵ ȍ Ϫȍ f (s) s Ϫ x ds (3.26) d ˆJ k ϭ l A J kϩl d J kϩl ϩ l B J kϩl s J kϩl (3.30) s ˆJ k ϭ l ⌫ J kϩl d J kϩl ϩ l T J kϩl s J kϩl ( 3
itself, these vectors are not necessarily sparse and may consist of 2 nϪj significant coefficients on scale j. Our algo-Proof. See Appendix A.3.

rithm uses the fact that for the operator functions considered in Section 3.2, the rows of the B-blocks have M 3.4. Adaptive Calculations with the Nonstandard Form vanishing moments. This means that when the row of a In [25] it was shown that Caldero ´n-Zygmund and B-block is applied to the ''smooth'' averages ͕s j ͖ the repseudo-differential operators can be applied to functions sulting vector is sparse (for a given accuracy ), as is illusin O(ϪN log ) operations, where N ϭ 2 n is the dimension trated in Fig. 1. of the finest subspace V 0 and is the desired accuracy. In Since each row of the B-block has the same number of this section we describe an algorithm for applying operavanishing moments as the filter G, we can use the ͕d j ͖ tors to functions with sublinear complexity, O(CN s ), where coefficients of the wavelet expansion to predict significant N s is the number of significant coefficients in the wavelet contributions to (3.28). In this way we can replace the representation of the function.

calculations with a dense vector ͕s͖ in (3.28) by calculations

We are interested in applying operators to functions with a sparse vector ͕s ˜͖, that are solutions of partial differential equations having regions of smooth, nonoscillatory behavior interrupted by

d ˜j k ϭ l A j kϩl d j kϩl ϩ l B j kϩl s ˜j kϩl , (3.32
) a number of well-defined localized shocks or shock-like structures. The wavelet expansion of such functions (see, e.g., (A.10)) then consists of differences ͕d j ͖ that are sparse for j ϭ 1, 2, ..., J Ϫ 1 and k ʦ ކ 2 nϪj . In what follows we and averages ͕s j ͖ that may be dense. Adaptively applying describe a method for determining the indices of ͕s ˜j͖ using the NS-form representation of an operator to a function the indices of the significant wavelet coefficients ͕d j ͖. expanded in a wavelet basis requires rapid evaluation of

The formal description of the procedure is as follows. For the functions under consideration the magnitude of

d ˆj k ϭ l A j kϩl d j kϩl ϩ l B j kϩl s j kϩl
(3.28) many wavelet coefficients ͕d j ͖ are below a given threshold of accuracy . The representation of f on V 0 , (A.10), using

s ˆj k ϭ l ⌫ j kϩl d j kϩl
(3.29) only coefficients above the threshold is

(P 0 f ) (x) ϭ J jϭ1 ͕k:͉d j k ͉Ͼ͖ d j k j,k (x) ϩ kʦކ 2 nϪJ s J k J,k (x), (3.33) for j ϭ 1, 2, ..., J Ϫ 1 and k ʦ ކ 2 nϪj ϭ ͕0, 1, 2, ..., 2 nϪJ Ϫ 1͖
and on the final, coarse scale whereas for the error we have If u(x) is expanded in a basis,

ʈ(P 0 f ) (x) Ϫ (P 0 f )(x)ʈ 2 (3.34) u(x) ϭ N iϭ1 u i b i (x), (4.1) ϭ ͩ J jϭ1 ͕k:͉d j k ͉Յ͖ ͉d j k ͉ 2 ͪ 1/2 Ͻ N 1/2 r ,
where u i are the coefficients and b i (x) are the basis functions, then in general where N r is the number of coefficients below the threshold. The number of significant wavelet coefficients is defined as N s ϭ N Ϫ N r , where N is the dimension of the space V 0 .

f (u(x)) ϶ N iϭ1 f (u i )b i (x).
(4.2) We define the -accurate subspace for f, denoted D f ʚ V 0 , as the subspace spanned by only those basis functions present in (3.33),

Clearly, this is the case for Fourier expansions.

Let us now assume that u and f (u) are both elements of V 0 . Then D

f ϭ V J ʜ ͕span ͕ j,k (x)͖ : ͉d j k ͉ Ͼ ͖ (3.35) u(x) ϭ k s 0 k (x Ϫ k), (4.3) for 1 Յ j Յ J and k ʦ ކ 2 nϪj . Associated with D f are sub- spaces S
f, j determined using the two-scale difference relation, e.g., Eq. (A.2). Namely, for each j ϭ 0, 1, ..., J Ϫ 1, where s k 0 are defined by s 0

k ϭ ͐ ȍ Ϫȍ u(x)(x Ϫ k) dx.
In addition, let us assume that the scaling function is interpolating,

S f, j ϭ ͕span ͕ j,2kϩ1 (x)͖ : jϩ1,k (x) ʦ D f ͖. (3.36) so that s 0 k ϭ u(k)
. Thus, we obtain For j ϭ J we define the space S f, j as

f (u) ϭ k f (s 0 k )(x Ϫ k); (4.4) S f, J ϭ V J . (3.37)
i.e., f (u) is evaluated by computing the function of the In terms of the coefficients d jϩ1 k the space S f, j may be deexpansion coefficients f (s 0 k ). Below we will describe how fined by to relax the requirement that the scaling function be interpolating and still have property (4.4) as a quantifiable ap-S f, j ϭ ͕span ͕ j,2kϩ1 (x)͖ : ͉d jϩ1 k ͉ Ͼ ͖.

(3.38) proximation.

We point out that typically f (u) is not in the same sub-In this way we can use D f to ''mask'' V 0 forming S f, j ; in space as u. In what follows we describe an adaptive algopractice all we do is manipulate indices. The subset of rithm for computing the pointwise square of a function, coefficients ͕s ˜j͖ that contribute to the sum (3.32) may now f (u) ϭ u 2 , where we split f (u) into projections on different be identified by indices of the coefficients corresponding subspaces. Working with ''pieces'' of the wavelet expanto basis functions in S f, j . In appendix A.4 we show that we sion of u we calculate contributions to f (u) using an approxmay indeed use the coefficients of the ͕d j ͖ to determine imation to (4.4). This is in direct contrast with calculating the ͕s ˜j͖ that contribute to (3.32). f (u) in a basis where the entire expansion must first be projected into a ''physical'' space, e.g., pseudo-spectral

EVALUATING FUNCTIONS IN WAVELET BASES

methods. In Section 4.2 we briefly discuss an algorithm for adaptively evaluating an arbitrary function f (u). In this section we describe our adaptive algorithm for evaluating the pointwise product of functions represented 4.1. Adaptive Calculation of u 2 in wavelet bases. More generally, our results may be applied to computing functions f (u), where f is an analytic Since the product of two functions can be expressed as a difference of squares, it is sufficient to explain an algo-function and u is expanded in a wavelet basis. We start by noting that since pointwise multiplication is a diagonal rithm for evaluating u 2 . The algorithm we describe is an improvement over that found in [40,23]. operator in the ''physical'' domain, computing the pointwise product in any other domain appears to be less effi-In order to compute u 2 in a wavelet basis, we first recall that the projections of u on subspaces V j and W j are given cient. In other words, a successful and efficient algorithm should at some point compute f (u) in the physical domain by P j u ʦ V j and Q j u ʦ W j for j ϭ 0, 1, 2, ..., J Յ n, respectively (see the discussion in Appendix A). Let j f , using values of u and not the expansion coefficients of u.

u 2 ϭ jʦޚ 2(P j u)(Q j u) ϩ (Q j u) 2 , (4.8)
which is essentially the paraproduct; see [27].

Evaluating (4.7) requires computing (Q j u) 2 and (P j u)(Q j u), where Q j u and P j u are elements of subspaces on the same scale and, thus, have basis functions with the same size support. In addition, we need to compute (P J u) 2 which involves only the coarsest scale and is not computationally expensive. The difficulty in evaluating (4.7) is that the terms (Q j u) 2 and (P j u) do not necessarily belong to the same subspace as the multiplicands. However, since

V j W j ϭ V jϪ1 ʚ V jϪ2 ʚ и и и ʚ V jϪj 0 ʚ ʚ и и и , (4.9)
we may think of both P j u ʦ V j and Q j u ʦ W j as elements of a finer subspace, that we denote V jϪj 0 , for some j 0 Ն 1. We compute the coefficients of P j u and Q j u in V jϪj 0 using the reconstruction algorithm, e.g., (A.10), and on V jϪj 0 we can calculate contributions to (4.7) using (4.4). The key observation is that, in order to apply (4.4), we may always choose j 0 in such a way that, to within a given accuracy , FIG. 2. The adaptive pseudo-wavelet algorithm. Averages on V j are (Q j u) 2 and (P j u)(Q j u) belong to V jϪj 0 . It is sufficient to ''masked'' by corresponding differences on W j . These coefficients are demonstrate this fact for j ϭ 0, which we do in Appenthen projected onto a finer subspace V jϪj 0 , Eq. (4.10) is evaluated, and the result is projected into the wavelet basis. dix A.5.

Remark. In practice j 0 must be small, and in our numerical experiments j 0 ϭ 3. We note that for the case of multi-1 Յ j f Յ J (see, e.g., Fig. 2, where j f ϭ 5 and J ϭ 8) be wavelets [20, 24] the proof using the Fourier domain does the finest scale having significant wavelet coefficients that not work since basis functions may be discontinuous. Howcontribute to the -accurate approximation of u; i.e., the ever, one can directly use the piecewise polynomial repreprojection of u can be expressed as sentation of the basis functions instead. For spline wavelets both approaches are available.

(P 0 u) (x) ϭ J jϭj f ͕k:͉d j k ͉Ͼ͖ d j k j,k (x) ϩ kʦކ 2 nϪJ s J k J,k (x). (4.5)
To describe the algorithm for computing the pointwise product, let us denote by R j j 0 (и) the operator to reconstruct (represent) a vector on subspace V j or W j in the subspace Let us first consider the case where u and u 2 ʦ V 0 , so that V jϪj 0 . On V jϪj 0 we can then use the coefficients R j j 0 (P j u) we can expand (P 0 u) 2 in a ''telescopic'' series, and R j j 0 (Q j u) to calculate contributions to the product (4.7) using ordinary multiplication as in (4.4). To this end, the contributions to (4.7) for j ϭ j f , j f ϩ 1, ..., J Ϫ 1 are

(P 0 u) 2 Ϫ (P J u) 2 ϭ J jϭj f (P jϪ1 u) 2 Ϫ (P j u) 2 .
(4.6) computed as

P jϪj 0 (u 2 ) ϭ 2( R j j 0 (P j u))( R j j 0 (Q j u)) (4.10) Decoupling scale interactions in (4.6) using P jϪ1 ϭ Q j ϩ P j , we arrive at ϩ ( R j j 0 (Q j u)) 2 , (P 0 u) 2 ϭ (P J u) 2 ϩ J jϭj f 2(P j u)(Q j u) ϩ (Q j u) 2 . (4.7)
where P j f (u) is the contribution to f (u) on subspace V j (see (4.7). On the final coarse scale J, we compute Later we will remove the condition that u and u 2 ʦ V 0 .

P JϪj 0 (u 2 ) ϭ ( R j j 0 (P J u)) 2 ϩ 2( R j j 0 (P J u))( R j j 0 (Q J u)) (4.11) Remark. Equation (4.7
) is written in terms of a finite number of scales. If j ranges over ,ޚ then (4.7) can be ϩ ( R j j 0 (Q J u)) 2 . written as

We then project the representation on subspaces V jϪj 0 , for rapidly converging Taylor series expansions, e.g., f (u) ϭ sin(u) for ͉u͉ sufficiently small. In this case, for a given j ϭ j f , ..., J into the wavelet basis. This procedure is completely equivalent to the decomposition one has to perform accuracy we fix an N so that ͉E j,N ( f, u)͉ Ͻ . We note that the partial differential Eq. (2.1) typically involves func-after applying the NS-form. The algorithm for computing the projection of u 2 in a wavelet basis is illustrated in Fig. tions f (и) that are not only analytic but in many cases are p-degree polynomials in u. If this is the case then for each 2. In analogy with ''pseudo-spectral'' schemes, as in, e.g., [34,35], we refer to this as an adaptive pseudo-wavelet fixed j the series in (4.13) is of degree p and E j,N ( f, u) ϭ 0 for N Ͼ p. In any event we are led to evaluate the algorithm.

To demonstrate that the algorithm is adaptive, we recall double sum in (4.14), which can be done using the adaptive pseudo-wavelet algorithm described in Section 4.1. that u has a sparse representation in the wavelet basis. Thus, evaluating (Q j u) 2 for j ϭ 1, 2, ..., J requires manipu-If the function f is not analytic, e.g., f (u) ϭ ͉u͉, then the primary concern is how to quantify an appropriate value lating only sparse vectors. Evaluating the square of the final coarse scale averages (P J u) 2 is inexpensive. The difficulty of j 0 , i.e., how fine a reconstruction (or how much ''oversampling'') is needed to take advantage of the inter-in evaluating (4.10) lies in evaluating the products R j j 0 (P j u)( R j j 0 Q j u) since the vectors P j u are typically dense. polating property s k 0 ϭ u(k). On the other hand, determining j 0 may become a significant problem even if f is analytic. The adaptivity of the algorithm comes from an observation that, in the products appearing in (4.10), we may use the For example if the Taylor series expansion of f (u) does not converge rapidly, as in the case of f (u) ϭ e u for large coefficients Q j u as a ''mask'' of the P j u (this is similar to the algorithm for adaptively applying operators to func-u, we have to consider alternate approaches.

For example, expanding e u in the ''telescopic'' series tions). In this way contributions to (4.10) are calculated, based on the presence of significant wavelet coefficients Q j u and, therefore, significant products R j j 0 (P j u)( R j j 0 Q j u). The complexity of our algorithm is automatically adaptable e P 0 u Ϫ e P J u ϭ J jϭ1 e P jϪ1 u Ϫ e P j u , (4.15) to the complexity of the wavelet representation of u.

Remarks on the Adaptive Calculation of

and using P jϪ1 ϭ Q j ϩ P j to decouple scale interactions, General f (u)

we arrive at This section consists of a number of observations regarding the evaluation of functions other than f (u) ϭ u 2 in wavelet bases. For analytic f (u) we can apply the same e P 0 u ϭ e P J u ϩ J jϭ1 e P j u (e Q j u Ϫ 1).

(4.16) approach as in Section 4.1, wherein we assume f (P 0 u) ʦ V 0 and expand the projection f (P 0 u) in the ''telescopic'' series Since the wavelet coefficients Q j u are sparse, the multipli-

f (P 0 u) Ϫ f (P J u) ϭ J jϭ1
f (P jϪ1 u) Ϫ f (P j u). (4.12) cand e Q j u Ϫ 1 is significant only where Q j u is significant.

Therefore, we can evaluate (4.16) using the adaptive pseudo-wavelet algorithm described in Section 4.1, where Using P jϪ1 ϭ Q j ϩ P j to decouple scale interactions in in this case the mask is determined by significant values (4.12) and assuming f (и) to be analytic, we substitute the of e Q j u Ϫ 1. The applicability of such an approach depends Taylor series on the relative size (or dynamic range) of the variable u.

For example, if u(x) ϭ Ͱ sin(2ȏx) on 0 Յ x Յ 1 then e ϪͰ Յ f (u) Յ e Ͱ .
It is clear that even for relatively moderate

f (Q j u ϩ P j u) ϭ N nϭ0 f (n) (P j u) n! (Q j u) n ϩ E j,N ( f, u) (4.13)
values of Ͱ the function e u may range over several orders of magnitude.

In order to take the dynamic range into account, we to arrive at apply a scaling and squaring method. Instead of computing e u directly, one calculates e u2 Ϫk and repeatedly squares the

f (P 0 u) ϭ f (P J u) ϩ J jϭ1 N nϭ1 f (n) (P j u) n! (Q j u) n ϩ E j,N ( f, u) result k times.
The constant k depends on the magnitude of u and is chosen so that the variable u is scaled as Ϫ1 Յ (4.14) 2 Ϫk u Յ 1, for example. In this interval, calculating e u2 Ϫk can be accomplished as described by Eq. (4.16) and the adaptive pseudo-wavelet algorithm of Section 4.1. One then For f (u) ϭ u 2 , j f ϭ 1, and N ϭ 2 we note that (4.14) and (4.7) are identical.

repeatedly applies the algorithm for squaring e u2 Ϫk to arrive at the wavelet expansion of e u . This approach can be used for functions f (u) that have

RESULTS OF NUMERICAL EXPERIMENTS

x dependence. In (5.4) E(и) is the explicit part of the approximation to (2.14) and I(и) is the implicit part. In this section we present the results of numerical experi-One can use specialized techniques for solving (5.4), ments in which we compute approximations to the solue.g., accelerating the convergence of the iteration by using tions of the heat equation, Burgers' equation, and two preconditioners (which may be readily reconstructed in a generalized Burgers' equations. In each of the examples wavelet basis; see, e.g., [22]). However, in our experiments we replace the initial value problem (2.1) with (2.2) and we use a straightforward fixed-point method to compute (2.3) by a suitable approximation, e.g., (2.17). The wavelet U(t jϩ1 ). We begin by setting representation of the operators appearing in this approximation are computed via (3.8). In order to illustrate the use U 0 (t jϩ1 ) ϭ E(U(t j )) ϩ I(U(t j ), U(t j )), (5.5) of our adaptive algorithm for computing f (u) developed in Section 4, we choose the basis having a scaling function and repeatedly evaluate with M shifted vanishing moments (see (A.4)) the so-called ''coiflets.'' This allows us to use the approximate interpolat-

U kϩ1 (t jϩ1 ) ϭ E(U(t j )) ϩ I(U(t j ), U k (t jϩ1 ))
(5.6) ing property; see, e.g., (5.2), below.

In each experiment we use a cutoff of ϭ 10 Ϫ6 , roughly for k ϭ 0, 1, 2, .... We terminate the iteration when corresponding to single precision accuracy. The number of vanishing moments is then chosen to be M ϭ 6 and the ʈU kϩ1 (t jϩ1 ) Ϫ U k (t jϩ1 )ʈ Ͻ , (5.7) corresponding length of the quadrature mirror filters H ϭ ͕h k ͖ L f kϭ1 and G ϭ ͕g k ͖ L f kϭ1 for ''coiflets'' satisfies L f ϭ 3M where (see, e.g., [18]). The number of scales n in the numerical realization of the multiresolution analysis depends on the most singular behaviour of the solution u(x, t). The specific

ʈU kϩ1 (t jϩ1 ) Ϫ U k (t jϩ1 )ʈ ϭ ͩ 2 Ϫn 2 n iϭ1 (U kϩ1 (x i , t jϩ1 ) Ϫ U k (x i , t jϩ1 )) 2 ͪ 1/2 .
(5.8) value of n used in our experiments is given with each example. We fix J, the depth of the wavelet decomposition, satisfying 2 nϪJ Ͼ L f , so that there is no ''wrap-around'' of the filters H and G on the coarsest scale.

Once (5.7) is satisfied, we update the solution and set Each of our experiments begins by projecting the initial condition (2.2) on V 0 , which amounts to evaluating U(t jϩ1 ) ϭ U kϩ1 (t jϩ1 ).

(5.9)

s 0 l ϭ ͵ ȍ ȍ u 0 (x)(x Ϫ l) dx.
(5.1) Again we note that one can use a more sophisticated iterative scheme and different stopping criteria for evaluating (5.4) (e.g., simply compute (5.6) for a fixed number of itera-For smooth initial conditions we approximate the integral tions). (5.1) (using the shifted vanishing moments of the scaling function (и)) to within via 5.1. The Heat Equation

s 0 l Ȃ u(l Ϫ Ͱ) (5.2)
We begin with this simple linear example in order to illustrate several points and provide a bridge to the nonlinear problems discussed below. In particular we show (see the discussion in Section 4.1). We note that in this that, in the wavelet system of coordinates, higher order case the discretization of the initial condition is similar to schemes do not necessarily require more operations than traditional discretizations, where one sets lower order schemes. We consider the heat equation on the unit interval, U(x i , t 0 ) ϭ u 0 (i ⌬x)

(5.3)

u t ϭ u xx , 0 Յ x Յ 1, 0 Յ t Յ 1, ( 5
.10) for i ϭ 0, 1, 2, ..., 2 n Ϫ 1, where ⌬x ϭ 2 Ϫn , and where U(x i , t) is the numerical approximation of the solution at grid point x i ϭ i ⌬x and time t.

for Ͼ 0, with the initial condition Since approximations to the integral in (2.14) are implicit in time, we solve an equation of the form u(x, 0) ϭ u 0 (x), 0 Յ x Յ 1, (5.11)

U(t jϩ1 ) ϭ E(U(t j )) ϩ I(U(t j ), U(t jϩ1 ))
(5.4) and the periodic boundary condition u(0, t) ϭ u(1, t). There are several well-known approaches for solving (5.10) and more general equations of this type having variable coeffi-for U(t jϩ1 ) by iteration, where we have dropped the explicit cients. Equation (5.10) can be viewed as a simple representative of this class of equations and we emphasize that the following remarks are applicable to the variable coefficient case, ϭ (x) (see also [START_REF] Engquist | Fast wavelet based algorithms[END_REF]).

For diffusion-type equations, explicit finite difference schemes are conditionally stable with the stability condition ⌬t/(⌬x) 2 Ͻ 1 (see, e.g., [START_REF] Stoer | les e ´quations aux de ´rive ´es partielles non-line ´aires[END_REF][START_REF] Dahlquist | Numerical Methods[END_REF]), where ⌬t ϭ 1/N t , ⌬x ϭ 1/N, and N t is the number of time steps. This condition tends to require prohibitively small time steps. An alternate, implicit approach is the Crank-Nicholson scheme [START_REF] Stoer | les e ´quations aux de ´rive ´es partielles non-line ´aires[END_REF][START_REF] Dahlquist | Numerical Methods[END_REF], which is unconditionally stable and accurate to O((⌬t) 2 ϩ (⌬x) 2 ). At each time step, the Crank-Nicholson scheme requires solving a system of equations, (5.12) method (5.12) with ⌬t ϭ ⌬x ϭ 2 Ϫ9 and ϭ 1.0. Note the slowly decaying peak in the solution that is due to the eigenvalue N ϭ ϭ Ϫ 0.99902344.

for j ϭ 0, 1, 2, ..., N t Ϫ 1, where we have suppressed the dependence of U(x, t) on x. The matrices A and B are given by A ϭ diag(ϪͰ/2, 1 ϩ Ͱ, ϪͰ/2) and B ϭ u 0 (x) ϭ ͭ x, 0ՅxՅ , 1Ϫx, ՅxՅ1, (5.15) diag(Ͱ/2, 1 Ϫ Ͱ, Ͱ/2), where Ͱ ϭ (⌬t/(⌬x) 2 ).

Alternatively, we can write the solution of (5.10) as that has a discontinuous derivative at x ϭ . Figure 3 illusu(x, t) ϭ e tL u 0 (x), (5.13) trates the evolution of (5.15) via (5.12) with ⌬t ϭ ⌬x and ϭ 1 and the slow decay of high frequency components of the initial condition. We have implemented Eq. (5.14) where L ϭ Ѩ xx , compute (5.13) by discretizing the time and display the result in Fig. 4 for the case where ϭ 1, interval [0, 1] into N t subintervals of length ⌬t ϭ 1/N t , and ⌬t ϭ ⌬x ϭ 2 Ϫn ϭ 1/N, and n ϭ 9. We note that there is a by repeatedly applying the NS-form of the operator e ⌬tL via proper decay of the sharp peak in the initial condition.

In order to illustrate the difference between the results U(t jϩ1 ) ϭ e ⌬tL U(t j )

(5.14) of our wavelet based approach and those of the Crank-Nicholson scheme, we construct the NS-form of the operator A Ϫ1 B and compare it with that of e ⌬tL . The NS-form of for j ϭ 0, 1, 2, ..., N t Ϫ 1, where U(t 0 ) ϭ U(0). The numerical an operator explicitly separates blocks of the operator that method described by (5.14) is explicit and unconditionally act on the high frequency components of u. These finer stable since the eigenvalues of e ⌬tѨ 2

x are less than one. The fact that the Crank-Nicholson scheme is unconditionally stable allows one to choose ⌬t independently of ⌬x; in particular one can choose ⌬t to be proportional to ⌬x. In order to emphasize our point we set ⌬x ϭ ⌬t and ϭ 1. Although the Crank-Nicholson scheme is secondorder accurate and such choices of the parameters ⌬x, ⌬t, and appear to be reasonable, by analyzing the scheme in the Fourier domain, we find that high frequency components in an initial condition decay very slowly. By diagonalizing matrices A and B in (5.12), it is easy to find the largest eigenvalue of A Ϫ1 B, N ϭ (1 Ϫ 2Ͱ)/(1 ϩ 2Ͱ). For the choice of parameters ϭ 1 and ⌬t ϭ ⌬x, we see that as Ͱ becomes large, the eigenvalue N tends to Ϫ1. We note that there are various ad hoc remedies (e.g., smoothing) used in conjunction with the Crank-Nicholson scheme to remove these slowly decaying high frequency components.

For example, let us consider the following initial con- dition Let us conclude by reiterating that the wavelet based scheme via (5.13) is explicit and unconditionally stable. The accuracy in the spatial variable of our scheme is O((⌬x) 2M ), where M is the number of vanishing moments, ⌬x ϭ 2 Ϫn , and n is the number of scales in the multiresolution analysis. Additionally, our scheme is spectrally accurate in time. Also it is adaptive simply by virtue of using a sparse data structure to represent the operator e ⌬tѨ xx , the adaptive algorithm developed in Section 3.4 and the sparsity of the solution in the wavelet basis. Finally, we note that if we were to consider (5.10) with variable coefficients, e.g., in, e.g., [START_REF] Beylkin | Wavelets in numerical[END_REF] (see also [43]).

u t ϭ (x)u xx , ( 5
Crank-Nicolson scheme (5.12). Entries of absolute value greater than

u t ϩ uu x ϭ u xx , 0 Յ x Յ 1, t Ն 0,
(5.17) scale or high frequency blocks are located in the upper for Ͼ 0, together with an initial condition, left corner of the NS-form. Therefore, the blocks of the NS-form of the operator A Ϫ1 B that are responsible for the u(x, 0) ϭ u 0 (x), 0 Յ x Յ 1, (5.18) high frequency components in the solution are located in the upper left portion of Fig. 5. One can compare Fig. 5 and periodic boundary conditions u(0, t) ϭ u(1, t). Burgers' with Fig. 6, illustrating the NS-form of the exponential equation is the simplest example of a nonlinear partial operator used in (5.14). Although the Crank-Nicholson differential equation incorporating both linear diffusion scheme is not typically used for this regime of parameters and nonlinear advection. Solutions of Burgers' equation (i.e., ϭ 1 and ⌬t ϭ ⌬x), a similar phenomena will be consist of stationary or moving shocks and capturing such observed for any low-order method. Namely, for a given behavior is an important simple test of a new numerical cutoff, the NS-form representation of the matrix for the method (see, e.g., [START_REF] Schult | Using wavelets to solve the Burgers[END_REF][START_REF] Liandrat | Nonlinear Partial Differential Equations using the Wavelet Approach. Wavelets and Their Applications[END_REF]42]). low-order scheme will have more entries than that of the Burgers' equation may be solved analytically by the corresponding exponential operator in the wavelet basis. Cole-Hopf transformation [START_REF] Hopf | The partial differential equation u t ϩ uu x ϭ Ȑu xx[END_REF][START_REF] Cole | On a quasilinear parabolic equation occurring in aerody-37[END_REF], wherein it is observed Referring to Fig. 5 and6 it is clear that the NS-form of that a solution of (5.17) may be expressed as the operator e ⌬tL in our high order scheme is sparser than the NS-form for the operator A Ϫ1 B in the second-order Crank-Nicholson scheme. The matrix in Fig. 5 has approxiu(x, t) ϭ Ϫ2 x , (5.19) mately 3.5 times as many entries as the matrix in Fig. 6.

where ϭ (x, t) is a solution of the heat equation with initial condition (x, 0) ϭ e Ϫ(1/4 ȏ) ͐ u(x,0)) dx .

(5.20)

Remark. We note that if is small, e.g., ϭ 10 Ϫ3 ; then using (5.19) as the starting point for a numerical method dition u(x, 0) ϭ sin(2ȏx)). Consequently, the finite arithme-tic involved in a numerical scheme leads to a loss of accu-increases during the formation of the shock, yet it never exceeded 10 over the entire simulation. The compression racy in directly calculating u(x, t) via (5.19), most notably within the vicinity of the shock. ratios of the NS-form representation of the first derivative, exponential, and nonlinear operator O L,m are 442.2, 3708.5, Our numerical scheme for computing approximations to the solution of (5.17) consists of evaluating and 1364.9, respectively, where the compression ratio is defined as N 2 /N s , where N is the dimension of the finest subspace V 0 and N s is the number of significant entries. U(t iϩ1 ) ϭ e ⌬tL U(t i )

(5.21) EXAMPLE 2. In this example we illustrate the wavelet

Ϫ O L,1 [U(t i )Ѩ x U(t iϩ1 ) ϩ U(t iϩ1 )Ѩ x U(t i )],
analogue of the Gibbs phenomena encountered when one does not use a sufficiently resolved basis expansion of the subject to the stopping criterion (5.7). Since the solution solution. In this example n ϭ 10, J ϭ 4, ⌬t ϭ 0.001, ϭ is expressed as the sum (5.21) and the linear part is equiva-0.001, and ϭ 10 Ϫ6 , and we refer to Figs. 9 and 10. Using lent to the operator used in the solution of the heat equan ϭ 10 scales to represent the solution in the wavelet basis tion, the linear diffusion in (5.17) is accounted for in an is insufficient to represent the high frequency components essentially exact way. Thus, we may attribute all numerical present in the solution. Figure 9 illustrates the projection artifacts in the solution to the nonlinear advection term of the solution on V 0 beyond the point in time where the in (5.17).

solution is well represented by n ϭ 10 scales. We see that For each of the following examples, we illustrate the high frequency oscillations have appeared in the projection accuracy of our approach by comparing the approximate which may be viewed as a local analogue of the Gibbs solution U w with the exact solution U e using phenomenon. Figure 10 illustrates the number of significant coefficients and the number of iterations per time step required to satisfy the stopping criterion (5.7). The

ʈU w Ϫ U e ʈ ϭ ͩ 2 Ϫn 2 n Ϫ1 iϭ0 (U w (x i , t) Ϫ U e (x i , t)) 2 ͪ 1/2
G(; x, t) ϭ ͵ 0 F(Ј) dЈ ϩ (x Ϫ ) 2 2t
(5.24) basis per time step. The number of operations per time step used to update the solution is proportional to the number of significant coefficients in the wavelet represen-and F() ϭ u 0 () is the initial condition (5.18) (see, e.g., [START_REF] Whitham | Global regularity of solutions; II. Local regularity, infinite products[END_REF]). The initial conditions have been chosen so that (5.24) tation of the solution. may be evaluated analytically and we compute the integrals 5.3. Generalized Burgers' Equation in (5.23) using a high order quadrature approximation.

In this section we consider the numerical solution of the EXAMPLE 1. In this example we set n ϭ 15, J ϭ 9, generalized Burgers' equation ⌬t ϭ 0.001, ϭ 0.001, and ϭ 10 Ϫ6 . The subspace V 0 may be viewed as a discretization of the unit interval into 2 15 u t ϩ u ͱ u x ϩ u Ͱ ϭ u xx , 0 Յ x Յ 1, t Ն 0, (5.26) grid points with the step size ⌬x ϭ 2 Ϫ15 . We refer to Figs. 7 and8. Figure 7 illustrates the projection of the solution on V 0 , and Fig. 8 illustrates the error (5.22) and the number for constants Ͱ, ͱ, Ͼ 0 and real , together with an initial condition u(x, 0), and periodic boundary conditions of significant coefficients per time step. The number of operations needed to update the solution is proportional u(0, t) ϭ u(1, t). This equation is thoroughly studied in [START_REF] Sachdev | Nonlinear Diffusive Waves[END_REF] and we illustrate the results of a number of experi-to the number of significant coefficients. The number of iterations required to satisfy the stopping criterion (5.7) ments which may be compared with [START_REF] Sachdev | Nonlinear Diffusive Waves[END_REF]. EXAMPLE 4. In this example we set ͱ ϭ Ͱ ϭ 1 and is computed via ϭ Ϫ1, and consider the evolution of a gaussian initial condition centered on the interval 0 Յ x Յ 1, e.g., U(t iϩ1 ) ϭ e ⌬t/(Ѩ 2

x ϩI) U(t i ) (5.28) u(x, 0) ϭ u 0 e Ϫ((xϪ1/2)) 2 . On the interval, the decay of u(x, 0) is suffciently fast that we can consider the initial

Ϫ O ˜Ѩ2 x ,1 [U(t i )Ѩ x U(t iϩ1 ) ϩ U(t iϩ1 )Ѩ x U(t i )],
condition to be periodic. We set n ϭ 15, J ϭ 4, ⌬t ϭ 0.001, and ϭ 10 Ϫ6 . For easy comparison with the results of [START_REF] Sachdev | Nonlinear Diffusive Waves[END_REF], we choose ϭ 0.0005. The approximation to the where solution of

u t ϩ uu x Ϫ u ϭ u xx , 0 Յ x Յ 1, t Ն 0, (5.27) O ˜Ѩ2
x ,1 ϭ

e ⌬t(Ѩ 2 x ϩI) Ϫ I Ѩ 2 x ϩ I , ( 5.29)
and I is the identity operator. We have chosen to use the operator L in the form L ϭ Ѩ 2 x ϩ I (see the development in, e.g., Section 2). We note that the NS-forms of the operators e ⌬t(Ѩ 2

x ϩI) and (5.29) are computed as described in Section 3.

Due to the negative damping in (5.27), the operator Ѩ 2 x ϩ I is no longer negative definite. Therefore, if the nonlinear term were not present, the solution would grow without bound as t increased. The solution of the nonlinear Eq. (5.27) evolves to form a single shock which grows as it moves to the right. Figure 13 illustrates the evolution of the projection of the solution and Fig. 14 illustrates the number of significant wavelet coefficients needed to represent the solution over the course of the experiment. On the other hand, the presence of the nonlinearity may affect the growth of the solution, depending on the size of the coefficient . We have increased the diffusion coefficient to ϭ 0.005; Fig. 15 illustrates the evolution of the projection of the solution and Fig. 16 illustrates the number of significant wavelet coefficients. We point out that the number of operations required to update the solution is FIG. 8. The error (5.22) per sample (Fig. 7) and the number of significant wavelet coefficients per time step in the approximation (5.21). proportional to the number of significant coefficients. x ,1 is given by (2.22). The only difference in (5.31), as compared with the approximation to Burgers' equations, mations to the solution of the so-called cubic Burgers' equation (5.21), is the presence of the cubic nonlinearity. We have computed approximations to the solution using our algorithms with n ϭ 13, J ϭ 6, ⌬t ϭ 0.001, ϭ 0.001, and ϭ x U(t i ) (5.31) a moving shock, and the sinusoidal initial condition evolves

u t ϩ u 2 u x ϭ u xx , 0 Յ x Յ 1, t Ն 0, (5.30
Ϫ O Ѩ 2 x ,1 [U 2 (t i )Ѩ x U(t iϩ1 ) ϩ U 2 (t iϩ1 )Ѩ x U(t i )],
into two right-moving shocks. We note that, although the number of grid points in a uniform discretization of such an initial value problem is, in this case, N ϭ 2 13 , we are using only a few hundred significant wavelet coefficients to update the solution.

CONCLUSIONS

In this paper we have synthesized the elements of numerical wavelet analysis into an overall approach for solving nonlinear partial differential equations. We have demonstrated an approach which combines the desirable features of finite difference approaches, spectral methods, and front-tracking or adaptive grid approaches usually applied to such problems. Specifically, we have considered the construction of and adaptive calculations with operator functions in wavelet bases, and we have developed an algorithm for the adaptive calculation of nonlinear functions, e.g., f (u) ϭ u 2 .

We used the semigroup method to replace the nonlinear partial differential equation (2.1) by a nonlinear integral equation (2.14), and outlined our approach for approximating such integrals. These approximations are expressed in terms of functions of differential operators, and we have (NS-form) representation. We then presented a fast, adap-In order to verify our approach, we have included the results of a number of numerical experiments, including tive algorithm for multiplying operators in the NS-form and functions expanded in wavelet bases. Additionally, the approximation to the solutions of the heat equation, Burgers' equation, and the generalized Burgers' equation. we have introduced an adaptive algorithm for computing functions f (u), in particular the pointwise product, where The heat equation was included to illustrate a number of simple observations made available by our approach. u is expanded in a wavelet basis. Both of these algorithms have an operation count which is proportional to the num-Burgers' equation and its generalization were included to illustrate the adaptivity inherent in wavelet-based ap-ber of significant wavelet coefficients in the expansion of u, and we note that both of these algorithms are necessary proaches, namely the ''automatic'' identification of sharp gradients inherent in the solutions of such equations. Since ingredients in any basis-expansion approach to numerically solving PDEs.

Burgers' equation is the simplest nonlinear example incorporating both diffusion and advection, it is typically a first example researchers investigate when introducing a new numerical method. There are several directions for this course of work which we have left for the future. One may consider nonperiodic boundary conditions instead of the periodic boundary condition (2.3). This may be accomplished by simply using a wavelet (or multiwavelet) basis on an interval rather than a periodized wavelet basis. Also, we note that variable coefficients in the linear terms of the evolution equation (2.1) (see, e.g., (5.16)) may be accommodated by computing the NS-form of the corresponding operators as outlined in, e.g., [START_REF] Beylkin | Wavelets in numerical[END_REF]. Another direction has to do with the choice of the wavelet basis. One of the conclusions which we have drawn from this study is that there seem to be a number of advantages to using basis functions which are piecewise polynomial. In particular the spline family of bases appears to be attractive as well as multiwavelets (see, e.g., [19]). In both cases there are also disadvantages and an additional study would help to understand such a trade-off. Yet another extension, which of course is the ultimate goal, is to consider multidimensional problems, e.g., the Navier-Stokes equations. Finally, although we did not address in this paper the puting the nonlinear contribution to the solution can be directly applied to this problem, the NS-form representa-In this appendix we provide a brief review of notions associated with multiresolution analysis (MRA); see, e.g., tion of the operator functions associated with this problem, e.g., e ͱ⌬tѨ 3

x , may be dense, even for rather small values of [18, [START_REF] Chui | An Introduction to Wavelets[END_REF][START_REF] Wickerhauser | On the fast algorithm for multiplication of functions ware[END_REF] for more details. Introducing the MRA as in Let P j denote the projection operator onto subspace V j The set of coefficients ͕s j k ͖ kʦޚ , which we refer to as ''averages,'' is computed via the inner product and let Q j ϭ P jϪ1 Ϫ P j be the projection operator onto subspace W j . The projection of a function f (x) onto subspace V j is given by

s j k ϭ ͵ ϩȍ Ϫȍ f (x) j,k (x) dx, (A .8)
(P j f )(x) ϭ kʦޚ s j k j,k (x). (A.6) and the set of coefficients ͕d j k ͖ kʦޚ , which we refer to as ''differences,'' is computed via the inner product Alternatively, it follows from (3.2) and (A.6) that we can also write (P j f )(x) as a sum of projections of f (x) onto The operators A j , B j , ⌫ j , and T J appearing in the NSform are represented by matrices Ͱ j , ͱ j , Ͳ j , and s j with entries defined by Ͱ j k,kЈ ϭ ͵ ͵ K(x, y) j,k (x) j,kЈ (y) dx dy, FIG. 22. Illustration of the application of the nonstandard form to a vector. ͱ j k,kЈ ϭ ͵ ͵ K(x, y) j,k (x) j,kЈ (y) dx dy, (A.14) Ͳ j k,kЈ ϭ ͵ ͵ K(x, y) j,k (x) j,kЈ (y) dx dy, 22 we see that the NS-form is applied to both averages and differences of the wavelet expansion of a function. In s j k,kЈ ϭ ͵ ͵ K(x, y) j,k (x) j,kЈ (y) dx dy.

this case we can view the multiplication of the NS-form and a vector as an embedding of matrix-vector multiplica-The operators in (3.5) are organized as blocks of a matrix tion into a space of dimension as shown in Fig. 21.

The price of uncoupling the scale interactions in (3.3) M ϭ 2 nϪJ (2 Jϩ1 Ϫ 1), (A.15) is the need for an additional projection into the wavelet basis of the product of the NS-form and a vector. The term where n is the number of scales in the wavelet expansion nonstandard form comes from the fact that the vector to and J Յ n is the depth of the expansion. This result must which the NS-form is applied is not a representation of then be projected back into the original space of dimension the original vector in the wavelet basis. Referring to Fig. N ϭ 2 n . We note that in general M Ͼ N, and for J ϭ n we have M ϭ 2N Ϫ 1.

It follows from (3.3) that after applying the NS-form to a vector we arrive at the representation

(T 0 f 0 )(x) ϭ J jϭ1 kʦކ 2 nϪj d ˆj k j,k (x) ϩ J jϭ1 kʦކ 2 nϪj s ˆj k j,k (x). (A.16)
The representation (A.16) consists of both averages and differences on all scales which can either be projected into the wavelet basis or reconstructed to space V 0 . In order to project (A.16) into the wavelet basis we form the representation We note that if we were to use any other finite-difference representation as coefficients on V 0 , the coefficients on V j scale j ϭ 2, 3, ..., J Ϫ 1, we decompose ͕s j ͖ ϭ ͕s ˆj ϩ s ˜j͖ into would not be related by scaling and would require individ-͕s ˜jϩ1 ͖ and ͕d ˜jϩ1 ͖ and form the sums ͕s jϩ1 ͖ ϭ ͕s ˆjϩ1 ϩ s ˜jϩ1 ͖ ual calculations for each j.

(T 0 f 0 )(x) ϭ J jϭ1 kʦކ 2 nϪj d j k j,k (x) ϩ kʦކ 2 nϪJ
and ͕d jϩ1 ͖ ϭ ͕d ˆjϩ1 ϩ d ˜jϩ1 ͖. The sets ͕s J ͖ and ͕d j ͖ J jϭ1 are the Using the two-scale difference equations (A.1) and coefficients of the wavelet expansion of (T 0 f 0 )(x), i.e., the (A.2), we are led to coefficients appearing in (A.17). This procedure is illustrated in Fig. 23.

Ͱ j l ϭ 2 LϪ1 kϭ0 LϪ1 kЈϭ0 g k g kЈ s jϪ1
2iϩkϪkЈ , Remark. An alternative to projecting the representation (A.16) into the wavelet basis is to reconstruct (A.16) to space V 0 , i.e., form the representation (A.6)

ͱ j l ϭ 2 LϪ1 kϭ0 LϪ1 kЈϭ0 g k h kЈ s jϪ1 2iϩkϪkЈ , (A.20) (P 0 f )(x) ϭ kʦޚ s 0 k 0,k (x), (A.18) Ͳ j l ϭ 2 LϪ1 kϭ0 LϪ1 kЈϭ0 h k g kЈ s jϪ1 2iϩkϪkЈ .
using the reconstruction algorithm described in Section A Therefore, the representation of Ѩ p x is completely deteras follows. Given the coefficients ͕s ˆj͖ J jϭ1 and ͕d ˆj͖ J jϭ1 , we mined by s 0 l in (A.14), or in other words, by the representareconstruct ͕d ˆJ͖ and ͕s ˆJ͖ into ͕s ˜JϪ1 ͖ and form the sum tion of Ѩ p x projected on the subspace V 0 . ͕s JϪ1 ͖ ϭ ͕s ˆJϪ1 ϩ s ˜JϪ1 ͖. Then on each scale j ϭ J Ϫ 1, J Ϫ To compute the coefficients s 0 l corresponding to the pro-2, ..., 1 we reconstruct ͕s ˆj͖ and ͕d ˆj͖ into ͕s ˜jϪ1 ͖ and form the jection of Ѩ p x on V 0 , it is sufficient to solve the system of sum ͕s jϪ1 ͖ ϭ ͕s ˆjϪ1 ϩ s ˜jϪ1 ͖. The final reconstruction (of ͕d 1 ͖ linear algebraic equations and ͕s 1 ͖) forms the coefficients ͕s 0 ͖ appearing in (A.18). This procedure is illustrated in Fig. 24.

s 0 l ϭ 2 p ͫ s 0 2l ϩ 1 2 L/2 kϭ1 a 2kϪ1 (s 0 2lϪ2kϩ1 ϩ s 0 2lϩ2kϪ1 ) ͬ (A.21) A.2

. The Nonstandard Form of Differential Operators

In this appendix we recall the wavelet representation of and differential operators Ѩ p x in the NS-form. The rows of the NS-form of differential operators may be viewed as finitel l p s 0 l ϭ (Ϫ1) p p!, (A.22) difference approximations on subspace V 0 of order 2M Ϫ 1, where M is the number of vanishing moments of the wavelet (x). This material is a review of material where a 2kϪ1 are the autocorrelation coefficients of H de- found in [21].

fined by The NS-form of the operator Ѩ p x consists of matrices A j , B j , ⌫ j for j ϭ 0, 1, ..., J and a ''coarse scale'' approximation T J . We denote the elements of these matrices by Ͱ j i,l ,

a n ϭ 2 LϪ1Ϫn iϭ0 h i h iϩn , n ϭ 1, ..., L Ϫ 1. (A.23) ͱ j i,l
, and Ͳ j i,l for j ϭ 0, 1, ..., J, and s J i,l . Since the operator Ѩ p

x is homogeneous of degree p, it is sufficient to compute We note that the autocorrelation coefficients a n with even the coefficients on scale j ϭ 0 and use indices are zero, a 2k ϭ 0, k ϭ 1, ..., L/2 Ϫ 1. (A.24)

The resulting coefficients s 0 l corresponding to the projection of the operator Ѩ p

x on V 0 may be viewed as a finitedifference approximation of order 2M Ϫ 1. Further details In order to show that (H)(x) has M vanishing moments, where P m (x) is a polynomial of degree m for 0 Յ m Յ we recall that in the Fourier domain vanishing moments M Ϫ 1; see [START_REF] Meyer | Wavelets and Operators[END_REF].

are characterized by Since the function f (и) is an analytic function of Ѩ x , we can expand f in terms of its Taylor series. The series for d m d m ˆ()͉ ϭ0 ϭ 0 for m ϭ 0, 1, ..., M Ϫ 1, (A.33) f (Ѩ x )P m (x) is finite and yields a polynomial of degree less than or equal to m, where ˆ() is the Fourier transform of (x). Setting g ˆ() ϭ f (Ѩ x )P m (x) ϭ P ˜mЈ (x), (A.27) ˆ() in (A.32), the sum on the right-hand side of (A.32) is zero. We also observe that the integrand on the righthand side of (A. where z ϭ z(x, j, l) lies between x and x ϩ 2 j l, we compute which again is of the same order as d jϩ1 kЈ . Therefore, if ͉d jϩ1 kЈ ͉ Ͻ for kЈ ʦ ކ 2 JϪ( jϩ1) , then for some constant C, d j k ϭ ͚ l ͱ j kϩl s j kϩl using (A.36) and obtain ͉d j k ͉ Ͻ C for k ʦ ކ 2 JϪj . choose j 0 in such a way that, to within a given accuracy , (Q j u) 2 and (P j u)(Q j u) belong to V jϪj 0 . It is sufficient to demonstrate this fact for j ϭ 0. In order to show that such (A.37) j 0 Ն 1 exists, we begin by assuming u ʦ V 0 ʚ V Ϫj 0 . This assumption implies that, in the Fourier domain, the support Due to the vanishing-moment property of the B-block of ˆ(2 Ϫj 0 ) ''overlaps'' the support of u ˆ(). Then, for scaling (Lemmas 1 and 2), the first term in (A. 2nϩl are of the same order, the differences satisfy 18]). We then write ͉d j 2n ͉ Ͻ C for some constant C. On the other hand, if k ϭ 2n ϩ 1 for n ϭʦ ކ 2 JϪ( jϩ1) , we find

LEMMA 1 .

 1 If the wavelet basis has M vanishing moments, then the B-blocks of the NS-form of the analytic operator function f (Ѩ x ), described in Section 3.2, satisfy ϩȍ lϭϪȍ l m ͱ j l ϭ 0 (3.25)

  .31) (where p.v. indicates the principle value) satisfy for k ʦ ކ 2 nϪJ . The difficulty in adaptively applying the NS-ϩȍ lϭϪȍ l m ͱ j l ϭ 0, (3.27) form of an operator to such functions is the need to apply the B-blocks of the operator to the averages ͕s j ͖ in (3.28).

FIG. 3 .

 3 FIG. 3. Solution of the heat equation using the Crank-Nicholson

FIG. 4 .

 4 FIG. 4. Solution of the heat equation using the NS-form of the exponential with ⌬t ϭ ⌬x ϭ 2 Ϫ9 and ϭ 1.0, i.e., Eq. (5.14).

  FIG. 5. NS-form representation of the operatorA Ϫ1 B used in the

10

  Ϫ8 are shown in black. The wavelet basis is Daubechies with M ϭ 6 5.2. Burgers' Equation vanishing moments (L f ϭ 18), the number of scales is n ϭ 9 and J ϭ 7. We have set ϭ 1.0 and ⌬t ϭ ⌬x ϭ 2 Ϫ9 . Note that the top left portion Our next example is the numerical calculation of soluof the figure contains nonzero entries which indicate high frequency tions of Burgers' equation components present in the operator A Ϫ1 B.

FIG. 6 .

 6 FIG. 6. NS-form representation of the operator e ⌬tL used in (5.14). turns out to be a poor approach. This is due to the large Entries of absolute value greater than 10 Ϫ8 are shown in black. The dynamic range of the transformed initial condition (5.20) wavelet basis is Daubechies with M ϭ 6 vanishing moments (L f ϭ 18), (approximately 70 orders of magnitude for the initial conthe number of scales is n ϭ 9 and J ϭ 7. We have set ϭ 1.0 and ⌬t ϭ ⌬x ϭ 2 Ϫ9 .

FIG. 7 .

 7 FIG. 7. The projection on V 0 of the solution of Burgers' equation at various time steps computed via the iteration (5.21). In this experiment n ϭ 15, J ϭ 9, ⌬t ϭ 0.001, ϭ 0.001, and ϭ 10 Ϫ6 . This figure corresponds to Example 1 of the text.

FIG. 9 .

 9 FIG. 9. The projection on V 0 of the solution of Burgers' equation at various time steps computed via the iteration (5.21). In this experiment n ϭ 10, J ϭ 4, ⌬t ϭ 0.001, ϭ 0.001, and ϭ 10 Ϫ6 . An analogue of the Gibbs phenomenon begins because the shock cannot be accurately represented by n ϭ 10 scales. Observe that the scheme remains stable in spite of the oscillations. This figure corresponds to Example 2 of the text.

  ) 10 Ϫ6 . Figures 17 and 18 illustrate the evolution of the solution for a gaussian initial condition, and Figs. 19 and 20 via illustrate the evolution of the solution for a sinusoidal initial condition. The gaussian initial condition evolves to U(t iϩ1 ) ϭ e ⌬tѨ 2

FIG. 10 .

 10 FIG.10. The total number of significant wavelet coefficients and the number of iterations needed to satisfy the stopping criterion (5.7) per shown how to expand these operator functions into a wave-

FIG. 12 .

 12 FIG.12. The error (5.22) per sample (Fig.11) and the number of significant wavelet coefficients per time step in the approximation (5.21).

FIG. 13 .

 13 FIG. 13. The projection on V 0 of the solution of (5.27) at various time steps. In this experiment n ϭ 15, J ϭ 4, ⌬t ϭ 0.001, ϭ 10 Ϫ6 , and ϭ 0.005. This figure corresponds to Example 4 of the text.

FIG. 14 .

 14 FIG. 14. The total number of significant wavelet coefficients per time step. This figure corresponds to Example 4 of the text.

FIG. 15 .

 15 FIG. 15. The projection on V 0 of the solution of (5.27) at various time steps. In this experiment n ϭ 15, J ϭ 4, ⌬t ϭ 0.001, ϭ 10 Ϫ6 , and ϭ 0.005. This figure corresponds to Example 4 of the text.

( 3 .

 3 FIG. 16. The total number of significant wavelet coefficients per time step. This figure corresponds to Example 4 of the text.

FIG. 17 .

 17 FIG. 17.The projection on V 0 of the solution of cubic Burgers' equation (5.30) at various time steps, computed using a gaussian initial condition. In this experiment n ϭ 13, J ϭ 6, ⌬t ϭ 0.001, ϭ 0.001, and ϭ 10 Ϫ6 . This figure corresponds to Example 5 of the text.

  FIG. 18. The total number of significant wavelet coefficients per time step. This figure corresponds to Example 5 of the text.

FIG. 19 .FIG. 20 .

 1920 FIG. 19.The projection on V 0 of the solution of cubic Burgers' equation (5.30) at various time steps, computed using a sinusoidal initial condition. In this experiment n ϭ 13, J ϭ 6, ⌬t ϭ 0.001, ϭ 0.001, and ϭ 10 Ϫ6 . This figure corresponds to Example 5 of the text.

ͱ j l ϭ 2 FIG. 23 .j l ϭ 2

 2232 FIG. 21. Organization of the nonstandard form of a matrix, A j , B j , and ⌫ j , j ϭ 1, 2, 3, and T 3 are the only nonzero blocks. sums ͕s 2 ͖ ϭ ͕s ˆ2 ϩ s ˜2͖ and ͕d 2 ͖ ϭ ͕d ˆ2 ϩ d ˜2͖. Then on each

FIG. 24 .A. 3 .

 243 FIG. 24. Reconstruction of the product of the NS-form and a function to space V 0 . are found in [25].

  x Ϫ l) ϭ P m (x), (A.[START_REF] Beylkin | Wavelets in numerical[END_REF]

  [START_REF] Meyer | Wavelets and numerical split-step method: A global adaptive scheme[END_REF], i.e., sign() ˆ(m) ()f ˆ(), is continuous where mЈ Յ m. Due to the M Ͼ m vanishing moments of at ϭ 0, once again because (x) has M vanishing mo-(x), the integrals (A.25) are zero and (3.25) is verified.ments. We can then define functions W ˆ(m) () for m ϭ 0, Proof of Lemma 2. The ͱ l elements of the NS-form of 1, ..., M Ϫ 1 as the Hilbert transform are given byͱ l ϭ ͵ ϩȍ Ϫȍ (x Ϫ l)(H )(x) dx, (A.28) W ˆ(m) () ϭ Ά Ϫi ˆ(m) (),as in Lemma 1, we find such that W ˆ(m) () coincides with the mth derivative of the generalized function (A.31) on the test functions f ʦ C ȍ 0 .)ޒ( Since W ˆ(m) () are continuous functions for m ϭ 0, 1, ..., M Ϫ 1, we obtain, instead of (A.30), ϩȍ lϭϪȍ l m ͱ l ϭ ϩȍ lϭϪȍ l m ͵ ϩȍ Ϫȍ (x Ϫ l)(H )(x) dx (x)x m e ix dx ϭ W ˆ(m) (). (A.35) Since W ˆ(m) ()͉ ϭ0 ϭ 0 the integrals (A.29) are zero and where, once again, we have used (A.26). (3.27) is established. To show that the integrals in (A.29) are zero, we estab-A.4. Masking ͕s j ͖ Coefficients Using ͕d j ͖ Coefficients lish that (H )(x) has at least M vanishing moments. Let us consider the generalized function Let us now show that we may indeed use significant wavelet coefficients ͕d jϩ1 ͖ to find coefficients ͕s ˜j͖ that contribute to (3.28). Expanding f (x ϩ 2 j l ) into the Taylor ͵ ȍ Ϫȍ (H )(x)x m e ix d ϭ i Ϫm Ѩ m (H )(). (A.30) series, In the Fourier domain the Hilbert transform of the function f (x ϩ 2 j l) ϭ MϪ1 mϭ0 f (m) (x) m! 2 jm l m ϩ f (M) (z) M! (z Ϫ x) M , (A.36) g defined by

ͪ dx A. 5 . 2 2

 522 Estimating the Amount of OversamplingIn this appendix we demonstrate that we may alwaysϩ ϪjϪj x Ϫ k)f (M) (z)(z Ϫ x) M dx.

2

 2 [START_REF] Liandrat | Nonlinear Partial Differential Equations using the Wavelet Approach. Wavelets and Their Applications[END_REF]) is zero and functions with a sufficient number of vanishing moments, the coefficients s Ϫj 0 l and the values u(x l ) for some x l may d f (M) (z)(z Ϫ 2 j (x ϩ k)) M dx be made to be within of each other. In this way we may then apply (4.4).The coefficients s Ϫj 0 l of the projection of u on V Ϫj 0 are(A.38) given by for k ʦ ކ 2 JϪj . To compute the differences d jϩ1 kЈ ϭ ͚ l g l s j 2kЈϩl , we uses Ϫj 0 l ϭ 2 j 0 /2 ͵ ȍ Ϫȍ u(x)(2 j 0 x Ϫ l) dx, (A.42) the averages which can be written in terms of u ˆ() as s j 2kЈϩl ϭ 2 Ϫj/2 ͵ ȍ Ϫȍ (2 Ϫj x Ϫ 2kЈ)f (x ϩ 2 j l) dx.(A.39)s Ϫj 0 l ϭ 2 j 0 /2 ͵ ȍ Ϫȍ u ˆ(2 j 0 ) ˆ()e Ϫil d.(A.43) Substituting (A.36) into (A.39), we obtain Replacing the integral in (A.43) by that over [Ϫȏ, ȏ], wed jϩ1 kЈ ϭ 2 Ϫj/2 ͵ ȍ Ϫȍ (2 Ϫj x Ϫ 2kЈ) Ϫj x Ϫ 2kЈ)f (M) (z)(z Ϫ x) M dx. s Ϫj 0 l ϭ 2 j 0 /2 kʦޚ ͵ ȏ Ϫȏ u ˆ(2 j 0 ( ϩ 2ȏk)) ˆ( ϩ 2ȏk)e Ϫil d.and, using the vanishing moments of the filter G ϭ ͕g l ͖, (A.44)Since u ʦ V 0 for any Ͼ 0, there is a j 0 such that thed f (M) (z)(z Ϫ 2 j (x ϩ 2kЈ)) M dxinfinite sum in (A.44) may be approximated to within by the first term, (A.40)for kЈ ʦ ކ 2 JϪ( jϩ1) . s Ϫj 0 l ϭ 2 j 0 /2 ͵ ȏ Ϫȏ u ˆ(2 j 0 ) ˆ()e Ϫil d. (A.45) To show that ͉d jϩ1 kЈ ͉ Ͻ implies ͉d j k ͉ Ͻ C, we consider two cases. First, if ͉d jϩ1 kЈ ͉ Ͻ and k is even, i.e., k ϭ 2n for In order to evaluate (A.45), we consider scaling functions n ʦ ކ 2 JϪ( jϩ1) , then we see that d j 2n and d jϩ1 kЈ given by (A.40) (x) having M-shifted vanishing moments, i.e., ͐ ȍ Ϫȍ (x Ϫ only differ in the coefficients g l and ͱ j 2nϩl . Since g l and Ͱ) m (x) dx ϭ 0, where Ͱ ϭ ͐ ȍ Ϫȍ x(x) dx (see, e.g., [25, ͱ j

  ͵ ȍ Ϫȍ (x Ϫ Ͱ) m (x) dx 1)f (M) (z) (z Ϫ 2 j (x ϩ 2n)) M dx,

Since u ʦ V 0 , the support of u ˆ(2 j 0 ( ϩ 2ȏk)) occupies a smaller portion of the support of ˆ( ϩ 2ȏk) as j 0 increases, 17. I. Daubechies, Orthonormal bases of compactly supported wavelets, Commun. Pure Appl. Math. 41, 909 (1988). and there exists a sufficiently large j 0 such that the coeffi- is the error term that is controlled by choosing j 0 sufficiently large.