Open access peer-reviewed chapter

Wavelets for Differential Equations and Numerical Operator Calculus

By Riccardo Bernardini

Submitted: September 7th 2018Reviewed: November 30th 2018Published: February 14th 2019

DOI: 10.5772/intechopen.82820

Downloaded: 278

Abstract

Differential equations are commonplace in engineering, and lots of research have been carried out in developing methods, both efficient and precise, for their numerical solution. Nowadays the numerical practitioner can rely on a wide range of tools for solving differential equations: finite difference methods, finite element methods, meshless, and so on. Wavelets, since their appearance in the early 1990s, have attracted attention for their multiresolution nature that allows them to act as a “mathematical zoom,” a characteristic that promises to describe efficiently the functions involved in the differential equation, especially in the presence of singularities. The objective of this chapter is to introduce the main concepts of wavelets and differential equation, allowing the reader to apply wavelets to the solution of differential equations and in numerical operator calculus.

Keywords

  • wavelets
  • differential equations
  • numerical analysis
  • finite element method
  • meshless
  • multiresolution analysis

1. Introduction

Partial differential equations (PDEs) are used commonplace in science and in engineering to model the behavior of physical systems. Because of their importance, many numerical techniques for their solutions have been developed: finite difference methods (FDMs), finite element methods (FEMs), spectral methods, Ritz/Galerkin approach, meshless approaches, and so on. The main characteristic of PDEs is that the “unknown” is a function, that is, an object with an infinite number of degrees of freedom. Because of this, it is usually impossible (even in principle) to get an exact solution by numerical means. The objective of every technique for PDE solving is to get a good approximation of the solution with limited computational resources (CPU time, memory, etc.).

The first step of every PDE solution algorithm is to discretize the PDE, that is, to approximate it with a finite-dimensional problem that can be solved by numerical means. A popular discretization technique is to discretize the space where the solution is searched by restricting the problem to a finite-dimensional vector space, that is, by writing the solution as the linear combination of several base functions. If the original PDE was linear, discretization will map it to a linear problem (typically a linear system or an eigenvalue problem). Even techniques such as FDM (that discretizes the domain) can be often reformulated as suitable discretization of the function space.

Intuitively, better approximations of the solution require finer discretization (the exact meaning of finer depends on the specific approach), especially if the solution has some regions of large variability. Since finer discretization implies larger problems (and, therefore, higher computational efforts), it is of interest to be able to change locally the discretization resolution to the solution variability, possibly in an adaptive way.

This need for different resolutions in different regions is the idea that links PDE with multiresolution analysis. The birth of multiresolution analysis goes back to 1990 with the works of Mallat [1] and Meyer [2]. Since then there have been a large number of papers ranging from very theoretical ones to application [3]. In a multiresolution analysis, a space of signals (most commonly L2IR, but not only) is represented as a nesting of spaces with different levels of “resolution.” This allows to write a signal as the sum of a “low-resolution” version plus some higher-resolution “details.”

Because of the ability of changing the resolution used to observe the signal (by adding or removing details), the multiresolution analysis is sometimes described as a mathematical zoom. This fact inspired many applications, including numerical solution of PDEs where they sound promising, especially for those problems that contains localized phenomena (e.g., shockwaves) or intra-scale interaction (e.g., turbulence).

The objective of this chapter is to introduce the reader to the application of wavelets to PDE solutions. This chapter can be ideally divided in three parts: in the first part, we recall briefly the main concepts about PDE and the main algorithms for solving PDE; successively we do a brief recall of multiresolution analysis and wavelets including also multiwavelets and second-generation wavelets that find often application in PDE solutions; and finally, we will illustrate few techniques that can be found in the literature.

1.1 Notation

ΩIRdis the domain where the functions of interest are defined. The boundary of Ωis partitioned as follows: ∂Ω=ΩDΩN, ΩDΩN=, where ΩDor ΩNcan be empty.

HΩwill denote a space of functions u:ΩIRdefined on ΩIRd.

IR0is the set of nonnegative reals.

2. Generalities on PDE

Most of the physical problems modeled by PDEs fall in one of the following three large classes: equilibrium, propagation, or eigenvalue problems.

  • In an equilibrium problem, we are interested in finding a function uHΩsuch that

    Du=fin ΩE1
    u=uDon ΩDE2
    un=uNon ΩNE3

In (1)(3) Dis an operator that includes derivatives and fHΩis known. Boundary conditions are typically given as constraint about uor its derivatives on regions of ∂Ω. Typical examples of physical systems giving rise to this type of problem are systems in steady state (e.g., temperature distribution, potential distribution, steady flows, and so on). Typically equilibrium problems are elliptic, that is, Dis an elliptic operator (a generalization of the Laplacian).

  • In a propagation problem, we are interested in modeling the time evolution of a physical system. The PDE can still be written as in (1)(3), but the domain can typically be written as Ω=IR0×W, where WIRd1and the first coordinate represents time. Boundary conditions for t=0are known as initial conditions. Example physical problems are heat or wave propagation. Propagation problems are typically hyperbolic or parabolic.

  • Finally, in eigenvalue problems we are interested in finding uand λthat satisfy

    Du=λuE4

A wide class of eigenvalue problems is represented by Sturm-Liouville problems that can be written as

py+qy=λwyE5

where the apostrophe denotes derivation, the unknowns are λ, and function is yHab, while p, q, and w, all belonging to Hab, are known. Sturm-Liouville problems include Bessel differential equations (obtained by writing Laplace, Helmholtz, or Schrodinger equation in cylindrical coordinates) and Lagrange differential equation (obtained working in spherical coordinates).

2.1 Solution of differential equations

The field of numerical solution of differential equations is very wide, and many techniques have been developed. Nevertheless, a categorization in few large classes is possible. An important step in every solution algorithm is mapping the differential equation into a discrete version with only a finite number of degrees of freedom. A first distinction can be done between techniques that achieve this objective by discretizing the domain Ωor the function space HΩ.

2.1.1 Domain discretization

The most known technique based on a domain discretization is the FDM where the unknown function is sampled in a finite number of points p1,p2,,pNΩand the derivatives are approximated with finite differences. By writing the differential equation for every p, with the derivative approximated as finite differences, one obtains a system of Nequations in Nunknowns that can be solved with known techniques. If the original PDE was linear, the discretized system will be linear too.

FDM is maybe the simplest approach and the most intuitive, and it can work quite well for simple problems and geometries. Moreover, in the linear case, since any approximation of a derivative in pwill consider only few points around p, the matrix of the discretized linear system will be very sparse, allowing for a reduction in the computational effort. The application of FDM techniques becomes difficult, albeit possible, in the case of complex problems.

2.1.2 Function space discretization

Another class of techniques discretizes the function space HΩby approximating it with an n-dimensional space Hn, that is, unknown function uis approximated as

ui=1nαibi,αiIRE6

where bii=1nis a basis of Hn.

By exploiting approximation (6), one can transform PDE (1)(3) into a finite-dimensional problem. The different solution techniques differ in how (6) is used and in the way of choosing space Hnand its basis bii=1n.

One possibility is to choose functions bithat are infinitely differentiable and nonvanishing on the whole Ω. This gives rise to so-called spectral methods. Typical choices for basis functions can be complex exponential/sinusoidal functions (if the solution is expected to be periodic), Chebyshev polynomials (for separable domains, e.g., d-dimensional cubes), and spherical harmonics (for systems with spherical symmetry). Spectral methods can work very well if the solution is expected to be smooth; they can even converge exponentially fast. However, their spatial localization is not good, and if the functions involved are not smooth (e.g., they are discontinuous), they lose most of their interest.

Another approach, very popular, is FEM that chooses functions biby first partitioning the domain Ωinto a set of elements (triangles and their multidimensional counterpart are a popular choice) and assigning to every element a suitable finite-dimensional vector space. The final approximation of uis constructed in a piecewise fashion by gluing, so to say, the approximations of uover every single element.

In a typical implementation of FEM, all the elements are affine images of a single reference element. This simplifies the implementation since it suffices to choose only the vector space of the reference element T0. Another popular choice is to choose the space associated to the elements as spaces of polynomials. The basis is selected by choosing a set of control points in q1,q2,T0and choosing as basis vectors bithe polynomials that satisfy the interpolation property

biqj=δi,j=1ifi=j0ifijE7

Remark 2.1 (generalized collocation method).

A generalization of this idea is to choose a set of functionals σjmapping functions defined over T0to IRand requiring

σjbi=δi,jE8

Eq. (8) gives back (7) if σjis defined as the functional that corresponds to evaluating the argument of the functional in qj. Eq. (8) is, however, more general than (7) since it can be used, for example, to control the flow through a face of the element.

An issue with FEM is that creating the grid of elements can be expensive. This is especially true in those problems where the geometry is not fixed but needs to be updated. An example of this type of system is free-surface fluid flows, where the interface between air and fluid changes with time, requiring a continuous update of the mesh. In order to solve this problem, meshless methods have been developed [4]. A typical meshless approach is to approximate uwith a discrete convolution with kernel as a chosen function φ, that is,

ux=CρIαIϕxxIρE9

where XIare a set of points of Ωand ρis a scale factor that allows to change the “resolution” of kernel function ϕ. Coefficient Cρcan be used to keep the energy constant as ρis changed.

2.1.3 Exploiting the discretization

After expressing uas linear combination of bi, we are left with the problem of determining the coefficients of the linear combination. Several approaches are possible; the easiest way to briefly present them is by rewriting the differential equation as

RuDuf=0E10

where operator R:HΩHΩis called the residual.

If we restrict uto be a linear combination of bi, most probably we will not be able to make residual (10) exactly zero; therefore, we will aim to make it as small as possible. Since the result of the residual operator is a function, there are many possible approaches in minimizing it.

With the collocation approach, we choose a number of points of the domain p1,p2,,pnΩand ask that the residual is zero on the chosen points, that is,

0=Rupj=Dupjfpjj=1,,nE11

Eq. (11) represents a system of nequations having as unknown the coefficients αi, i=1,,n. For example, if Dis linear, (11) becomes

fpj=Di=1nαibipj=i=1nαiDbipj=i=1nαiAj,ij=1,,nE12

where, clearly, Aj,i=Dbipj. Note that (12) is a linear system in unknowns αi.

Remark 2.2.

With reference to Remark 2.1, one can generalize the collocation method by using a set of linear functionals σj:HΩIR. In this case one can obtain a generalized version of (12), namely,

σjf=i=1nαiσjDbiAj,ij=1,,nE13

Another approach is to solve Ru=0in a least square sense, that is, to search for coefficients αithat minimize

Ru2=RuRuE14

Standard algebra allows to show that (14) is minimized when Ruis orthogonal to Ru/αifor every i, that is,

RuRuαi=0i=1,,nE15

If Dis linear,

Ruαj=αjDi=1nαibif=DbjE16

and we get

0=RuRuαj=Di=1nαibifDbj=i=1nαiDbiDbjfDbjE17

which is still a linear system.

The Galerkin method is inspired on the idea that in a least square approximation, the error is orthogonal to the space where the approximating function lives. We would like to approximate the solution of the PDE with a vector of Hn; however, we do not know the solution, so we ask for the residual to be orthogonal to Hn, that is,

Ruv=0vHnE18

Eq. (18) is equivalent to

Duv=fvvHnE19

which can be interpreted as the original differential equation Du=fin weak form. Form (19) is often exploited by integrating by parts the left-hand side scalar product, moving one differentiation from the unknown function uto the test function v. This is often useful when a piecewise linear approximation is employed and Dcontains second-order differential operators (that cannot be applied on piecewise linear functions). Eq. (19) is verified for all vHnif and only if it is verified for every vector in a basis of Hn, that is, (19) is equivalent to

Dubj=fbjj=1,,nE20

If Dis linear, from (20) one can easily derive the linear system in αi

fbj=i=1nαiDbibjj=1,,nE21

Finally, it is worth citing the method of weighted residuals that can be seen as a generalization of the Galerkin PDE method. The idea is that instead of asking the residual being orthogonal to the space Hnused to approximate u, we ask the residual to be orthogonal to a different n-dimensional space Kn=spanβ1βnwhere βii=1nis clearly a basis of Kn. One obtains

fβj=i=1nαiDbiβjj=1,,nE22

Remark 2.3.

It is worth observing that from the weighted residual method, Galerkin and least square methods can be derived by a suitable choice of βi; even collocation method can be derived if we allow βito be a delta function (so that the scalar product needs to be interpreted as a distribution pairing). Moreover, for every vsince map xxvis a functional, it is easy to recognize that every method can be considered like a generalized collocation method, as described in Remark 2.1.

3. Wavelets

The idea of multiresolution analysis is to approximate vectors of L2IRwith variable degrees of resolution. This is achieved through a multiresolution analysis scheme defined by means of some axioms. The first axiom is the existence of a sequence VnnZof subspaces of L2IRnested one inside the other, that is,

V2V1V0V1V2E23

The idea is that if one approximates (in a least square sense) a function fwith vectors belonging to Vn, the approximation error gets smaller as nincreases since every vector of Vnalso belongs to Vn+1. Note, however, that (23) does not grant that we will be able to approximate fwith an error as small as desired; in order to grant this, we need another axiom

nZVn¯=L2IRE24

where the overline denotes set closure (in the topology induced by the norm on L2IR). Axiom (24) requires that every vector of L2IRis in the closure of the union in the left hand; this means that given any ϵ>0and fL2IR, it is possible to find an element of the union whose distance from fis less than ϵ. In other words, (24) means that whatever fL2IRand whatever the chosen maximum approximation error allowed ϵ, one can find a space Vnthat approximates fwith the required precision.

An axiom dual to (24) is

nZVn=0E25

that requires that there is only one “lowest resolution vector,” that is, the null vector.

Remark 3.1.

In order to see that axiom (24) is not obvious, it is more convenient to work with Hilbert space L201. Recall that functions xcos2πnx, xsin2πnx, and nINand the constant 1are an orthogonal basis of L201.

Define S0=sin2π2ktkINas the set of all the even-numbered sines, and define V0as the space generated by S0, that is,

V0spanS0E26

Now define spaces Vn, n<0by removing one vector at time from the basis of V0, and define spaces Vn, n<0by adding one odd harmonic at time. More precisely, define

Vn=spanSnE27

where

Sn=Sn+1\sin2π2ntifn<0Sn1sin2π2n1tifn>0E28

It is clear that the sequence of spaces defined in this way satisfies axiom (23), but not (24), since, for example, function cos2πtis orthogonal to every Vn. Note that this construction can be repeated for any Hilbert space using an orthonormal basis of the space instead of sines and cosines. Another axiom makes more precise the idea of “increasing resolution” by asking that vectors in Vnvary twice as faster than the vectors in Vn1. In order to make this more precise, define operator S:L2IRL2IRas the rescaling operator Sfx=2f2x. Note that because of the multiplication by 2, Sis unitary, that is, Sf=f. The new axiom is

fVnSfVn+1E29

It follows that Vn=SnV0. With this position one can interpret (25) by saying that the slowest function is the constant (and the only constant in L2IRis the zero).

The last axiom puts a constraint on the structure of V0by asking that is generated by a function ϕand its translations. In order to make this more precise, define the operator τt:L2IRL2IRassociated to a translation of tas τtfx=fxt. Note that also τtis unitary and that the exponential notation is convenient since τaτb=τa+b. Observe also the commutation relation Sτt=τt/2S. The last axiom can be written as

ϕL2IR:V0=spanτiϕiZE30

Often as part of axioms, it is required that ϕis orthogonal to its translations, that is,

τiϕτjϕ=δi,jE31

However, it is not necessary to include (31) explicitly in the axioms since, given a ϕthat satisfies (30), it is possible to orthonormalize it, so that it satisfies (31), with a well-known “Fourier trick” [3]. Therefore, we will suppose (31) satisfied.

It is worth to summarize here the axioms

V2V1V0V1V2E32
nZVn¯=L2IRE33
nZVn=0E34
Vn+1=SVnE35
V0=spanτiϕiZϕL2IRE36

The axioms above allow us to determine a property of ϕ. Note that since V1V0, ϕV1. Note also that set Sτiϕ=τi/2SϕiZis an orthonormal basis of V1. It follows that one can write ϕas linear combination of Sϕand its half-integer translations, that is,

ϕ=iZgiτi/2SϕE37

for some sequence gi:ZIR. Eq. (37) is known as two-scale equation and it is central to wavelet theory. Function ϕis known as scaling function.

Remark 3.2.

Note that from the orthonormality of τi/2SϕiZfollows

giτ2kgi=φτkφ=δk,0E38

where the left-hand side scalar product is the usual scalar product in 2Z.

Remark 3.3.

Note that starting from a ϕthat satisfies a two-scale equation like (37), it is possible to recover a full multiresolution analysis. Indeed, one defines V0according to (36) and Vnby repeated applications of (35). Two-scale Eq. (37) grants that the nesting axiom (32) is satisfied.

Note also that (37) shows that ϕis the fixed point of operator

OiZgiτi/2SE39

This suggests that maybe one could start from a sequence giand apply repeatedly Oto a vector of L2IRin order to obtain ϕ. This is indeed possible, but the theoretical details are out of scope here; see [5].

Since Vn+1Vnone can consider the orthogonal complement of Vnin Vn+1; call it Wn, that is,

Vn+1=VnWnE40

It is possible to find, starting from the two-scale Eq. (37), a function ψsuch that τiψiZis an orthonormal basis of W0. This implies that it must be for all

ψ=iZhiτi/2SϕψV1E41
ψτiψ=δiorthonormalbasisE42
τjψτiϕ=0i,jZW0orthogonaltoV0E43

By using (41) and the orthonormality of τi/2Sϕ, it is possible to rewrite (42) and (43) as conditions on hi, namely,

τ2khihi=δk,0E44
τ2khigi=0E45

It is easy to verify that, given gi, a possible hithat satisfies (44) and (45) is

hi=1ign+1E46

This shows that sequences hiand giare the impulse responses of a two-channel orthogonal filter bank. Moreover, if giand ϕare known, one can obtain ψby choosing hiaccording to (46) and computing ψaccording to (41). Function ψis known as wavelet, and it generates the whole L2IRwith its translations and dilations.

This has also another interesting consequence. Suppose fL2IRand that

γi1=τi/2Sϕf=SτiϕfE47

are the coefficients of its projection on V1. Suppose we need the coefficients γk0=τkϕfof the projection on V0. It is possible to exploit the two-scale equation

γk0=τkϕf=τkiZgiSτiϕf=iZgiSτi+2kϕf=iZgiγi+2k1=gγi12kE48

where gis the time-reversed version of gi. Eq. (48) shows that it is possible to go to the space at lower resolution by means of a filtering by gand a decimation by a factor of two. Similarly, by calling ηk0=τkψfthe coefficients relative to the projection of fon W0, one can obtain

ηk0=iZhiγi+2k1=hγi12kE49

Figure 1a shows this idea: the sequence of high-resolution coefficients are processed with a two-channel orthogonal filter bank, and the coefficients relative to the lower resolution space Vnexit from one branch, and the coefficients relative to the “missing details” space Wnexit from the other. The idea can be iterated several times; see Figure 1b. This is the basis of the well-known fast algorithm to compute wavelet coefficients and also the origin of the minor, and very common, misnomer in calling Figure 1b a “discrete-time wavelet transform.”

Figure 1.

(a) Splitting coefficient sequence into a low-resolution and a high-resolution one using a two-channel filter bank. (b) Iteration of structure (a) makes a fast algorithm for computing the wavelet coefficients.

An interesting characteristic of wavelets is that they can be used to detect the local regularity of a function. This is similar to what happened with Fourier transform where a function that is discontinuous has a Fourier transform that decays as 1/ω; if the function is continuous but not derivable, its Fourier transform decays as 1/ω2and so on. With the wavelet transform happens something similar, with the scale playing the role of frequency. The interesting difference is that while a Fourier transform that decays as 1/ωtells us that there is at least one discontinuity, but not where, with the wavelet transform the slow decay with the scale is localized around the discontinuity. The precise claim of this property requires the introduction of the concept of Lipschitz regularity and would take us too far; see [6]. This suggests that when approximating the unknown function in a PDE, we can keep high-resolution coefficients only in the neighborhood of singularities, saving on computational effort.

We will say that wavelet ψhas vanishing moments if

IRxkψxdx=0k=0,1,,1E50

An interesting property of compactly supported wavelets with vanishing moments is that the corresponding scaling function (not the wavelet itself) can reproduce polynomials of degree at most 1in the sense that if Pxis a polynomial with degree less than , there exist coefficients cisuch that

Px=iZciτiϕxE51

In other words, V0contains all the polynomials of degree less than .

Example 3.1. (Haar wavelet).

The simplest example of wavelet is the Haar wavelet whose scaling function is

ϕHx=1ifx010elseE52

It is immediate to verify that ϕHsatisfies a two-scale equation

ϕH=12SϕH+τ1/2SϕHE53

with coefficients g0=g1=1/2. Note that trivially τ2kϕHϕH=δk. In order to create the corresponding wavelet, use prescription hi=1igi+1to get

ϕH=12SϕHτ1/2SϕHE54

Note that the Haar wavelet is compactly supported, but it is discontinuous. This makes it not well suited to approximate smooth functions.

Example 3.2 (Sinc wavelet).

An example in some sense opposite to the Haar wavelet is the Sinc wavelet. In this case V0is the space of “low-pass” functions, that is, functions whose Fourier transform is zero outside interval ππ. As well known, V0is generated by the Sinc function

sincx=sinπxπxE55

and its translations, that is,

V0=spanτksinckZE56

This suggests to use ϕS=sincas scaling function. The fact that a two-scale equation is satisfied is easily checked in frequency since V1is the space of functions whose Fourier transform is zero outside 2π2π; therefore, every function of V0is contained in V1, as desired.

The corresponding wavelet ψSis easily characterized in frequency as the function whose Fourier transform is

ΨSω=1ifπ<ω<2π0otherwiseE57

It is easy to verify that ψSV1, ψSV1, and τ2kψSψS=δk.

As said above, the Sinc example is somehow the opposite of Haar wavelet: it is arbitrarily differentiable, but it has infinite support; actually, it decays very slowly (as O1/x), and this introduces several practical issues. Moreover, sequences giand hiare of infinite length, and they decay slowly too (they do not even have a z-transform), making it difficult to implement it.

Example 3.3 (spline wavelet).

An example intermediate between Haar and Sinc wavelet is represented by spline spaces of degree d. In this case V0dis defined as the space of piecewise polynomial functions that are dtimes differentiable (with continuous derivative), with the “breaking points” on the integer, more precisely

V0d={fL2IRCd1fkk+1=polynomialdegreedkZ}E58

It is easy to see that V1d=SV0dis a similar space of piecewise polynomial functions but with the breaking points in half integers. It follows that every function in V0dalso belongs to V1d, giving rise to a multiresolution analysis.

A generator for V0dcan easily be obtained as a suitable translation (necessary to align the breaking points) of

rectdx=rectx=1ifx<1/20otherwiseifd=0rectd1rectxifd>0,E59

that is, the rect convolved with itself dtimes. Function rectdhas compact support but it is not, however, orthogonal to its own translations. It can be orthogonalized with the Fourier trick, but the result has no compact support. For more details about this case, see [3].

3.1 Compactly supported wavelets: Daubechies’ wavelets

In the examples above, we found multiresolution analysis whose scaling function had at most two out of the following three desirable characteristics: orthogonality, smoothness, and compact support. Is it possible to find a wavelet that has all these three characteristics? The answer was given by Daubechies. It turns out that imposing all the three characteristics is very demanding and only a small family of wavelets exists.

An easy observation is that if ϕis orthogonal to its own translations, the coefficients giin the two-scale equation can be obtained as

gi=ϕτi/2SϕE60

According to (60) if ϕhas compact support, then gihas a finite number of coefficients that are different from zero. Since gineeds to be orthogonal to its even translations, its length (i.e., the number of nonzero coefficients) must be necessarily even [3].

Moreover, if the iteration of operator Oin (39) converges, it is easy to see that gihas a finite number of coefficients, and then the limit function has compact support. This suggests that it “suffices” to find a finite length sequence githat is orthogonal to its own even translation and iterates operator Oto obtain the desired scaling functions. It actually turns out that this can be done, although there are lots of technical details to be taken care of (e.g., about convergence of Okand smoothness of the resulting ϕ); see [5] for details.

Every member of the Daubechies family is identified by the length 2Nof the sequence gi(remember that the length of giis necessarily even). It can be proven that the resulting scaling function has Nvanishing moments and its smoothness grows with N; see Table 1. See also Figure 2 that shows the results of the first three iterations of O(first row), the final scaling function (second row), and the wavelet (third row) of three different Daubechies wavelets.

N2345678910
βS11.4151.7752.0962.3882.6582.9143.1613.402
βH0.5500.9151.2751.5961.8882.1582.4152.6612.902

Table 1.

Hölder βHand Sobolev βSregularity exponent of Daubechies’ wavelets as function of length 2 N of gi.

Figure 2.

Daubechies’ wavelets. First three iterations of O (first row), the final scaling function (second row), and the wavelet (third row) of three different Daubechies’ wavelets.

3.2 Extensions

The construction given above is the original idea of multiresolution analysis. Since the early 1990s, many researchers worked in this field, and many variations and extensions have been introduced. Here we briefly recall those that have more interest in the field of differential equation solutions.

3.2.1 Multiwavelets

Multiwavelets are a generalization of standard multiresolution analysis in the sense that now scaling functions and wavelets are vectors of functions. This means that Vnis not generated by the translations of a single function but from the translations of many functions. Every idea of standard multiresolution analysis can be reformulated without much difficulty in this case, with the most notable difference that the two-scale equation now has vector function and coefficients that are matrices. Multiwavelets can accommodate scaling factors different from two and there is a larger choice for compact support wavelets. See [7] for more details.

3.2.2 Second-generation wavelets

Two-scale Eq. (37) and the resulting filter bank-based procedure work well when the data are sampled on a regular grid and/or the functions of interest are defined on IRd. Since there are many applications that do not satisfy this requirement (e.g., differential equations on general manifolds), the idea of second-generation wavelet has been introduced.

The starting point is the so-called lifting form of filter bank (Figure 1). It is possible to show that any two-channel filter bank (Figure 1) can be implemented as shown in Figure 3. In the lifting approach, the input signal is split into odd and even samples by a serial-to-parallel converter. The first branch is filtered, and the result combined with the other branch; the result of this operation is filtered again and combined with the first branch, and this iterated as long as necessary. Filter Pis sometimes called the prediction step, and it is interpreted as a filter that predicts the odd samples from the even ones; filter Uis sometimes called update.

Figure 3.

Lifting implementation of the two-channel filter banks associated with a wavelet analysis.

The advantage of this form is that, being similar to the Feistel structure used in cryptography [8], it is exactly invertible even if operations are implemented in fixed-point arithmetic. Actually, the invertibility does not depend on the detail of Piand Ui; that can be anything, even nonlinear. Another interesting advantage of this predict/update idea is that it does not require a regular domain, allowing to bring the wavelet concept to more general contexts. For example, [9] uses this idea to solve differential equations on the sphere.

4. Some examples of application of wavelets to PDE

By now it should be clear how multiresolution analysis can be applied to differential equation solution: by using scaling functions and/or wavelets as basis functions in approximation (6). All the approaches described in Section 2 can be used with wavelet: collocation, Galerkin PDE method, weighted residual method, meshless methods, etc. Before describing some details of few approaches described in the literature, it is worth to do some general remarks.

What makes wavelet interesting is their multiresolution property and the fact that a wisely chosen wavelet (smooth and/or with many vanishing gradients) has interesting “singularity sensing” properties: in the neighborhood of a singularity (discontinuity, nondifferentiability, etc.), the coefficients decay as a function of scale with a speed that depends on the singularity involved (similar to what Fourier transform does, only on a local level), but away from the singularity, they decade fast [3]. This implies that good approximations can be obtained with few coefficients, using high-resolution decomposition only where it is necessary, reducing the size of the matrices involved in the solution of the PDE. A similar effect can be obtained, for example, in FEM by using a finer mesh around points of large variation. However, using this approach in an adaptive way would require to adjust at running time the mesh, a potentially heavy operation. Wavelets have the potential of employing an adaptive resolution in an easier way. See, for example, [10] for few examples of adaptive techniques employing wavelets.

While orthogonality is considered an important feature in many theoretically wavelet papers, in the context of differential equation solution, it plays a smaller role. The reason is that basis functions enter in the scalar products associated with the various methods via the differential operator D, and it is not guaranteed that Dwill preserve orthogonality (that would give rise to many zero entries, that is, sparser matrices).

Actually, orthogonality is preserved if the two basis functions have disjoint support in space (since differential operators do not extend the support) or in frequency (since differential operators are translation-invariant and in frequency they become a product). This suggests that in the context of differential equations, compact support and well-localization in frequency are more important than just orthogonality. In a sense, they represent a “robust” orthogonality condition.

Remark 4.1.

It is true that true compact support in frequency is less common than compact support in space. With the exception of few very special and theoretical cases (e.g., Sinc), the best we can get is a rapid decay in frequency. This means that the scalar product of two basis functions separated in frequency will be maybe very small, but not zero. Nevertheless, even this kind of “almost sparseness” can be exploited.

A general issue with wavelets is that it can be difficult to impose boundary conditions since they have no natural interpolation property that would make boundary condition handling simpler. Another problem that can arise is that many wavelets have no closed-form description, for example, Daubechies wavelets that are described as the result of the iteration of operator Oin (39). This can make their application to PDE more difficult, for example, when computing scalar products involved in weighted residual and other methods.

Finally, another common issue is that most of the known wavelets are defined on a one-dimensional domain, while many physical systems are on a multidimensional domain. The easiest way to create a multidimensional multiresolution analysis by a one-dimensional one is the separable (or tensor product or Kronecker product) approach that create a multidimensional function by the product of several one-dimensional ones, e.g.,

ϕ3Dxyz=ϕxϕyϕzE61

This kind of approach, however, produces “cube-like” wavelets, and their application to FEM schemes based on triangular elements can be difficult.

4.1 Some schemes from the literature

In this section we briefly summarize some interesting wavelet-based schemes that can be found in the literature. As said above, wavelets and scaling functions can be used as basis in the approximation used in collocation, weighted residuals, and other methods.

Wavelets in Galerkin and weighted residual methods bring the advantage of their multiresolution and localization properties while, however, suffering from difficulties in handling complex boundary conditions. Moreover, nonlinear equations can turn out to be difficult to handle. Nevertheless, there have been many successful examples in the application to elliptic, hyperbolic, and parabolic PDE [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. Wavelet-based collocation methods, where wavelet functions are used as shape functions, also registered some success. The advantage of collocation methods is that they are more easily applicable in nonlinear cases [28, 29] and irregular boundary conditions [30]. A collocation method based on second-generation wavelet and lifting is applied to a nonlinear vibration problem in [30, 31].

Much more popular seems to be the application of wavelets to FEM techniques. In this case wavelets or scaling functions are used as shape functions instead of the more traditional polynomials. Daubechies wavelets are particularly popular most probably because of their compact support property. Also of interest is the fact that Daubechies’ wavelets can have any number of null moments, making possible the perfect interpolation of polynomials. Some examples of successful application Daubechies wavelets to PDE (mostly mechanical problems) are [32, 33, 34, 35, 36]. Of special interest is the proposal of Mitra [37] where wavelet-based FEM is used to transform a wave propagation problem into ordinary differential equations that are successively solved.

Another popular solution for wavelet-based FEM is the wavelets based on spline spaces. Although spline bases cannot have both compact support and orthogonality, in differential equations, as explained above, we gladly give up on orthogonality if we can get compact support and smoothness. Another important advantage of splines is that a simple closed-form expression is known. Examples of spline applications can be found in [38, 39, 40]. Of special interest is the application of Hermite cubic splines (HCS), a kind of multiwavelet [41] that shows promise in handling in a numerically robust way boundary conditions. The HCS is a multiwavelet with four smooth (twice differentiable) components defined on interval 01. Some examples of application can be found in [42, 43, 44]. A problem with the application of wavelets to FEM techniques is the difficulty to adapt the wavelet construction to complex meshes. In this case the use of second-generation wavelet based on an extension of the lifting idea has attracted some attention [45, 46, 47].

Wavelets have attracted some interest also in the context of meshless methods [9, 48, 49]. Of some interest for the special problem is [9] that uses a wavelet approach to implement a meshless solver for differential equations defined on the sphere. A problem with applying wavelets to generic manifolds like a sphere is that it is not clear what a “rescaling by 2” should mean for a manifold that is not a Euclidean space. The idea used in [9] is to use a so-called diffusion wavelet where the dilation is replaced by a diffusion operator that looks like a kind of “low-pass filtering” that smear out the details; see [9] for the precise definition.

5. Conclusions

This chapter introduced the reader to the field of applying wavelets to the numerical solution of differential equations. Both wavelets and differential equations are research fields with many applications, contributions, and results. Their combination gives rise to wide varieties of methods, each one suited for specific applications. By looking at the literature, we can see that wavelets can be a very powerful tool for solving PDE especially because of their multiresolution nature that allows to optimize the level of detail where it is needed. Wavelets, however, are not a silver bullet for all problems either, since they can have some characteristics (multidimensional construction via tensor product, nonexistence of a closed-form expression, difficulty in handling some boundary conditions, etc.) that can make their application not trivial in some cases. We can say that this is a field where, more than ever, no single solution fits all and that every practitioner needs to find the solution specific for the problem at hand using knowledge in both fields and some ingenuity.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Riccardo Bernardini (February 14th 2019). Wavelets for Differential Equations and Numerical Operator Calculus, Wavelet Transform and Complexity, Dumitru Baleanu, IntechOpen, DOI: 10.5772/intechopen.82820. Available from:

chapter statistics

278total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

DWT-Based Data Hiding Technique for Videos Ownership Protection

By Farhan Al-Enizi and Awad Al-Asmari

Related Book

Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology

Edited by Dumitru Baleanu

First chapter

Real-Time DSP-Based License Plate Character Segmentation Algorithm Using 2D Haar Wavelet Transform

By Zoe Jeffrey, Soodamani Ramalingam and Nico Bekooy

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us