Hölder and Sobolev regularity exponent of Daubechies’ wavelets as function of length 2 N of .
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.
- differential equations
- numerical analysis
- finite element method
- multiresolution analysis
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  and Meyer . Since then there have been a large number of papers ranging from very theoretical ones to application . In a multiresolution analysis, a space of signals (most commonly , 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.
is the domain where the functions of interest are defined. The boundary of is partitioned as follows: , , where or can be empty.
will denote a space of functions defined on .
is 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 such thatE1E2E3
In (1)–(3) is an operator that includes derivatives and is known. Boundary conditions are typically given as constraint about or 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, is 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 , where and the first coordinate represents time. Boundary conditions for are 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 and that satisfyE4
A wide class of eigenvalue problems is represented by Sturm-Liouville problems that can be written as
where the apostrophe denotes derivation, the unknowns are , and function is , while , , and , all belonging to , 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 .
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 and the derivatives are approximated with finite differences. By writing the differential equation for every , with the derivative approximated as finite differences, one obtains a system of equations in unknowns 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 will consider only few points around , 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 by approximating it with an -dimensional space , that is, unknown function is approximated as
where is a basis of .
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 and its basis .
One possibility is to choose functions that 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., -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 by 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 is constructed in a piecewise fashion by gluing, so to say, the approximations of over 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 . 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 and choosing as basis vectors the polynomials that satisfy the interpolation property
Remark 2.1 (generalized collocation method).
A generalization of this idea is to choose a set of functionals mapping functions defined over to and requiring
Eq. (8) gives back (7) if is defined as the functional that corresponds to evaluating the argument of the functional in . 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 . A typical meshless approach is to approximate with a discrete convolution with kernel as a chosen function , that is,
where are a set of points of and is a scale factor that allows to change the “resolution” of kernel function . Coefficient can be used to keep the energy constant as is changed.
2.1.3 Exploiting the discretization
After expressing as linear combination of , 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 is called the residual.
If we restrict to be a linear combination of , 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 and ask that the residual is zero on the chosen points, that is,
where, clearly, . Note that (12) is a linear system in unknowns .
With reference to Remark 2.1, one can generalize the collocation method by using a set of linear functionals . In this case one can obtain a generalized version of (12), namely,
Another approach is to solve in a least square sense, that is, to search for coefficients that minimize
Standard algebra allows to show that (14) is minimized when is orthogonal to for every , that is,
If is linear,
and we get
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 ; however, we do not know the solution, so we ask for the residual to be orthogonal to , that is,
Eq. (18) is equivalent to
which can be interpreted as the original differential equation in weak form. Form (19) is often exploited by integrating by parts the left-hand side scalar product, moving one differentiation from the unknown function to the test function . This is often useful when a piecewise linear approximation is employed and contains second-order differential operators (that cannot be applied on piecewise linear functions). Eq. (19) is verified for all if and only if it is verified for every vector in a basis of , that is, (19) is equivalent to
If is linear, from (20) one can easily derive the linear system in
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 used to approximate , we ask the residual to be orthogonal to a different -dimensional space where is clearly a basis of . One obtains
It is worth observing that from the weighted residual method, Galerkin and least square methods can be derived by a suitable choice of ; even collocation method can be derived if we allow to be a delta function (so that the scalar product needs to be interpreted as a distribution pairing). Moreover, for every since map is a functional, it is easy to recognize that every method can be considered like a generalized collocation method, as described in Remark 2.1.
The idea of multiresolution analysis is to approximate vectors of with 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 of subspaces of nested one inside the other, that is,
The idea is that if one approximates (in a least square sense) a function with vectors belonging to , the approximation error gets smaller as increases since every vector of also belongs to . Note, however, that (23) does not grant that we will be able to approximate with 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 ). Axiom (24) requires that every vector of is in the closure of the union in the left hand; this means that given any and , it is possible to find an element of the union whose distance from is less than . In other words, (24) means that whatever and whatever the chosen maximum approximation error allowed , one can find a space that approximates with the required precision.
An axiom dual to (24) is
that requires that there is only one “lowest resolution vector,” that is, the null vector.
In order to see that axiom (24) is not obvious, it is more convenient to work with Hilbert space . Recall that functions , , and and the constant are an orthogonal basis of .
Define as the set of all the even-numbered sines, and define as the space generated by , that is,
Now define spaces , by removing one vector at time from the basis of , and define spaces , by 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 is orthogonal to every . 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 vary twice as faster than the vectors in . In order to make this more precise, define operator as the rescaling operator . Note that because of the multiplication by , is unitary, that is, . The new axiom is
It follows that . With this position one can interpret (25) by saying that the slowest function is the constant (and the only constant in is the zero).
The last axiom puts a constraint on the structure of by asking that is generated by a function and its translations. In order to make this more precise, define the operator associated to a translation of as . Note that also is unitary and that the exponential notation is convenient since . Observe also the commutation relation . 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” . 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 , . Note also that set is an orthonormal basis of . It follows that one can write as linear combination of and its half-integer translations, that is,
for some sequence . Eq. (37) is known as two-scale equation and it is central to wavelet theory. Function is known as scaling function.
Note that from the orthonormality of follows
where the left-hand side scalar product is the usual scalar product in .
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 according to (36) and by 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 and apply repeatedly to a vector of in order to obtain . This is indeed possible, but the theoretical details are out of scope here; see .
Since one can consider the orthogonal complement of in ; call it , that is,
It is possible to find, starting from the two-scale Eq. (37), a function such that is an orthonormal basis of . This implies that it must be for all
This shows that sequences and are the impulse responses of a two-channel orthogonal filter bank. Moreover, if and are known, one can obtain by choosing according to (46) and computing according to (41). Function is known as wavelet, and it generates the whole with its translations and dilations.
This has also another interesting consequence. Suppose and that
are the coefficients of its projection on . Suppose we need the coefficients of the projection on . It is possible to exploit the two-scale equation
where is the time-reversed version of . Eq. (48) shows that it is possible to go to the space at lower resolution by means of a filtering by and a decimation by a factor of two. Similarly, by calling the coefficients relative to the projection of on , 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 exit from one branch, and the coefficients relative to the “missing details” space exit 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.”
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 ; if the function is continuous but not derivable, its Fourier transform decays as and 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 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 . 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 in the sense that if is a polynomial with degree less than , there exist coefficients such that
In other words, contains 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 satisfies a two-scale equation
with coefficients . Note that trivially . In order to create the corresponding wavelet, use prescription to 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 is the space of “low-pass” functions, that is, functions whose Fourier transform is zero outside interval . As well known, is generated by the Sinc function
and its translations, that is,
This suggests to use as scaling function. The fact that a two-scale equation is satisfied is easily checked in frequency since is the space of functions whose Fourier transform is zero outside ; therefore, every function of is contained in , as desired.
The corresponding wavelet is easily characterized in frequency as the function whose Fourier transform is
It is easy to verify that , , and .
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 ), and this introduces several practical issues. Moreover, sequences and are 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 . In this case is defined as the space of piecewise polynomial functions that are times differentiable (with continuous derivative), with the “breaking points” on the integer, more precisely
It is easy to see that is a similar space of piecewise polynomial functions but with the breaking points in half integers. It follows that every function in also belongs to , giving rise to a multiresolution analysis.
A generator for can easily be obtained as a suitable translation (necessary to align the breaking points) of
that is, the rect convolved with itself times. Function has 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.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 in the two-scale equation can be obtained as
According to (60) if has compact support, then has a finite number of coefficients that are different from zero. Since needs to be orthogonal to its even translations, its length (i.e., the number of nonzero coefficients) must be necessarily even .
Moreover, if the iteration of operator in (39) converges, it is easy to see that has a finite number of coefficients, and then the limit function has compact support. This suggests that it “suffices” to find a finite length sequence that is orthogonal to its own even translation and iterates operator to 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 and smoothness of the resulting ); see  for details.
Every member of the Daubechies family is identified by the length of the sequence (remember that the length of is necessarily even). It can be proven that the resulting scaling function has vanishing moments and its smoothness grows with ; see Table 1. See also Figure 2 that shows the results of the first three iterations of (first row), the final scaling function (second row), and the wavelet (third row) of three different Daubechies wavelets.
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.
Multiwavelets are a generalization of standard multiresolution analysis in the sense that now scaling functions and wavelets are vectors of functions. This means that is 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  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 . 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 is sometimes called the prediction step, and it is interpreted as a filter that predicts the odd samples from the even ones; filter is sometimes called update.
The advantage of this form is that, being similar to the Feistel structure used in cryptography , it is exactly invertible even if operations are implemented in fixed-point arithmetic. Actually, the invertibility does not depend on the detail of and ; 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,  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 . 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,  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 , and it is not guaranteed that will 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.
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 in (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.,
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 . 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  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  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 . 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  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  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  for the precise definition.
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.