Open access peer-reviewed chapter

Wavelets for Differential Equations and Numerical Operator Calculus

Written By

Riccardo Bernardini

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

DOI: 10.5772/intechopen.82820

Chapter metrics overview

978 Chapter Downloads

View Full Metrics


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.


  • 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 discretizethe 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 finerdiscretization (the exact meaning of finerdepends 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 problemswe are interested in finding uand λthat satisfy


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


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


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 spectralmethods. 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


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


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, meshlessmethods have been developed [4]. A typical meshless approach is to approximate uwith a discrete convolution with kernel as a chosen function φ, that is,


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


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 collocationapproach, 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,


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


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,


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


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


If Dis linear,


and we get


which is still a linear system.

The Galerkin methodis 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 residualto be orthogonal to Hn, that is,


Eq. (18) is equivalent to


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


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


Finally, it is worth citing the method of weighted residualsthat 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


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 schemedefined by means of some axioms. The first axiom is the existence of a sequence VnnZof subspaces of L2IRnested one inside the other, that is,


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


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


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,


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




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


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


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


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


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,


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

Remark 3.2.

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


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


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,


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


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


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


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


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


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


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


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


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


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


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


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


and its translations, that is,


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


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


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


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


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.


Table 1.

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

Figure 2.

Daubechies’ wavelets. First three iterations ofO(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 wavelethas 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 predictionstep, 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 productor Kronecker product) approach that create a multidimensional function by the product of several one-dimensional ones, e.g.,


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 waveletwhere the dilation is replaced by a diffusion operatorthat 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 alland that every practitioner needs to find the solution specific for the problem at hand using knowledge in both fields and some ingenuity.


  1. 1. Mallat S. Multiresolution approximations and wavelet orthonormal bases ofL2(R). Transactions of the American Mathematical Society. 1989;315:69-87
  2. 2. Meyer Y. Ondelettes, Vol. 1 of Ondelettes et Opérateurs. Paris: Hermann; 1990
  3. 3. Vetterli M, Kovačević J. Wavelets and Subband Coding. Signal Processing. Englewood Cliffs, NJ: Prentice-Hall; 1995
  4. 4. Huerta A, Belytschko T, Fernández-Méndez S, Rabczuk T, Zhuang X, Arroyo M. Meshfree Methods. In: Stein E, Borst R, Hughes TJ, editors. Encyclopedia of Computational Mechanics. 2nd ed. 2017. DOI: 10.1002/9781119176817.ecm2005
  5. 5. Daubechies I. Ten Lectures on Wavelets. Philadelphia: SIAM; 1992
  6. 6. Mallat S, Hwang WL. Singularity detection and processing with wavelets. IEEE Transactions on Information Theory. 1992;38:617-643
  7. 7. Keinert F. Studies in advanced mathematics. In: Wavelets and Multiwavelets. Chapman and Hall/CRC; 2003
  8. 8. Schneier B. Applied Cryptography: Protocols, Algorithms, and Source Code in C. Hoboken, New Jersey: John Wiley and Sons; 1994
  9. 9. Goyal K, Mehra M. An adaptive meshfree diffusion wavelet method for partial differential equations on the sphere. Journal of Computational Physics. 2014;272:747-771
  10. 10. Li B, Chen X. Wavelet-based numerical analysis: A review and classification. Finite Elements in Analysis and Design. 2014;81:14-31
  11. 11. Hussain NEA. 2D wavelets Galerkin method for the computation of EM field on seafloor excited by a point source. International Journal of Applied Electromagnetics and Mechanics. 2017;53(4):631-644
  12. 12. Qian S, Weiss J. Wavelets and the numerical solution of partial differential equations. Journal of Computational Physics. 1993:155-175
  13. 13. Amaratunga K, Williams JR. Wavelet based green’s function approach to 2d pdes. Engineering Computations. 1993;10(4):349-367
  14. 14. Ho S, Yang S. Wavelet-galerkin method for solving parabolic equations in finite domains. Finite Elements in Analysis and Design. 2001;37(12):1023-1037
  15. 15. AL-Qassab M, Nair S. Wavelet-galerkin method for free vibrations of elastic cable. Journal of Engineering Mechanics. 2003;129(3):350-357
  16. 16. Monasse P, Perrier V. Orthonormal wavelet bases adapted for partial differential equations with boundary conditions. SIAM Journal on Mathematical Analysis. 1998;29(4):1040-1065
  17. 17. Nastos CV, Theodosiou TC, Rekatsinas CS, Saravanos DA. A 2d daubechies finite wavelet domain method for transient wave response analysis in shear deformable laminated composite plates. Computational Mechanics. 2018;62:1187-1198
  18. 18. Dahmen W, Kunoth A, Urban K. A wavelet galerkin method for the stokes equations. Computing. 1996;56:259-301
  19. 19. Yang S, Ni G, Ho SL, Machado JM, Rahman MA, Wong HC. Wavelet-galerkin method for computations of electromagnetic fields-computation of connection coefficients. IEEE Transactions on Magnetics. 2000;36:644-648
  20. 20. Yongping S, Pu L, Rufu H. A wavelet-galerkin method for the simulation of torsion mems devices under the effect of squeeze film damping. In: 2011 International Conference on Electric Information and Control Engineering; 2011. pp. 5671-5674
  21. 21. Yang S, Ni G, Cardoso JR, Ho SL, Machado JM. A combined wavelet-element free galerkin method for numerical calculations of electromagnetic fields. IEEE Transactions on Magnetics. 2003;39:1413-1416
  22. 22. Liu Y, Liu Y, Ding K. A new coupling technique for the combination of wavelet-galerkin method with finite element method in solids and structures. International Journal for Numerical Methods in Engineering.112(10):1295-1322
  23. 23. Venini P, Morana P. An adaptive wavelet-galerkin method for an elastic-plastic-damage constitutive model: 1d problem. Computer Methods in Applied Mechanics and Engineering. 2001;190(42):5619-5638
  24. 24. Sannomaru S, Tanaka S, ichiro Yoshida K, Bui TQ, Okazawa S, Hagihara S. Treatment of dirichlet-type boundary conditions in the spline-based wavelet galerkin method employing multiple point constraints. Applied Mathematical Modelling. 2017;43:592-610
  25. 25. Yang Z, Liao S. On the generalized wavelet-galerkin method. Journal of Computational and Applied Mathematics. 2018;331:178-195
  26. 26. Bu H, Wang D, Zhou P, Zhu H. An improved wavelet–galerkin method for dynamic response reconstruction and parameter identification of shear-type frames. Journal of Sound and Vibration. 2018;419:140-157
  27. 27. Choudhury A, Deka R. Wavelet–galerkin solutions of one dimensional elliptic problems. Applied Mathematical Modelling. 2010;34(7):1939-1951
  28. 28. Zhang T, Tian Y-C, Tadé MO. Wavelet-based collocation method for stiff systems in process engineering. Journal of Mathematical Chemistry. 2008;44:501-513
  29. 29. Lepik Ü. Haar wavelet method for solving stiff differential equations. In: Mathematical Modelling and Analysis. Vilnius, Lithuania: Vilnius Gediminas Technical University; 2009;14:467-481
  30. 30. Vasilyev OV, Bowman C. Second-generation wavelet collocation method for the solution of partial differential equations. Journal of Computational Physics. 2000;165(2):660-693
  31. 31. Vasilyev OV, Kevlahan NK-R. An adaptive multilevel wavelet collocation method for elliptic problems. Journal of Computational Physics. 2005;206(2):412-431
  32. 32. Adigun BJ, Buchan AG, Adam A, Dargaville S, Goffin MA, Pain CC. A haar wavelet method for angularly discretising the boltzmann transport equation. Progress in Nuclear Energy. 2018;108:295-309
  33. 33. Patton RD, Marks PC. One-dimensional finite elements based on the daubechies family of wavelets. AIAA Journal. 1996;34:1696-1698
  34. 34. Ma J, Xue J, Yang S, He Z. A study of the construction and application of a daubechies wavelet-based beam element. Finite Elements in Analysis and Design. 2003;39(10):965-975
  35. 35. Li B, Hongrui C, He Z. The construction of one-dimensional Daubechies wavelet-based finite elements for structural response analysis. Journal of Vibroengineering. 2011;13:729-738
  36. 36. He W-Y, Wang Y, Zhu S. Adaptive reconstruction of a dynamic force using multiscale wavelet shape functions. Shock and Vibration. 2018;2018:8213105
  37. 37. Mitra M, Gopalakrishnan S. Wave propagation analysis in anisotropic plate using wavelet spectral element approach. Journal of Applied Mechanics. 2008;75(1):014504
  38. 38. Han J-G, Ren W-X, Huang Y. A spline wavelet finite-element method in structural mechanics. International Journal for Numerical Methods in Engineering.66(1):166-190
  39. 39. Chen WH, Wu CW. A spline wavelets element method for frame structures vibration. Computational Mechanics. 1995;16:11-21
  40. 40. Pian THH, Chen D-P. Alternative ways for formulation of hybrid stress elements. International Journal for Numerical Methods in Engineering. 1982;18(11):1679-1684
  41. 41. Jia R-Q, Liu S-T. Wavelet bases of hermite cubic splines on the interval. Advances in Computational Mathematics. 2006;25:23-39
  42. 42. Xiang J-w, Chen X-f, Li X-k. Numerical solution of poisson equation with wavelet bases of hermite cubic splines on the interval. Applied Mathematics and Mechanics. 2009;30:1325
  43. 43. Chen X, Xiang J. Solving diffusion equation using wavelet method. Applied Mathematics and Computation. 2011;217(13):6426-6432
  44. 44. Xiang J, Wang Y, Jiang Z, Long J, Ma G. Numerical simulation of plane crack using hermite cubic spline wavelet. CMES: Computer Modeling in Engineering and Sciences. 2012;88:1-16
  45. 45. He Y, Chen X, Xiang J, He Z. Adaptive multiresolution finite element method based on second generation wavelets. Finite Elements in Analysis and Design. 2007;43(6):566-579
  46. 46. He Y, Chen X, Xiang J, He Z. Multiresolution analysis for finite element method using interpolating wavelet and lifting scheme. Communications in Numerical Methods in Engineering;24(11):1045-1066
  47. 47. Sudarshan R, D’Heedene S, Amaratunga K. A multiresolution finite element method using second generation hermite multiwavelets. In: Bathe K, editor. Computational Fluid and Solid Mechanics 2003. Oxford: Elsevier Science Ltd; 2003. pp. 2135-2140
  48. 48. Liu Y, Liu Y, Cen Z. Daubechies wavelet meshless method for 2-d elastic problems. Tsinghua Science and Technology. 2008;13:605-608
  49. 49. Yousefi MR, Jafari R, Moghaddam HA. A combined wavelet based mesh free method for solving the forward problem in electrical impedance tomography. In: 2012 IEEE International Symposium on Medical Measurements and Applications Proceedings; 2012. pp. 1-4

Written By

Riccardo Bernardini

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