Wavelet-based Methods for Numerical Solutions of Differential Equations
Bin Han Michelle Michelle Yau Shu Wong


Wavelet Methods for Differential Equations (ICCM 2019)]Wavelet-based Methods for Numerical Solutions of Differential Equations



Research supported in part by
Natural Sciences and Engineering Research Council (NSERC) of Canada and Alberta Innovates

Department of Mathematical and Statistical Sciences,
University of Alberta, Edmonton,  Alberta, Canada T6G 2G1.
  bhan@ualberta.ca  mmichell@ualberta.ca   yauwong@ualberta.ca


addtoresetequationsection 
Wavelet theory has been well studied in recent decades. Due to their appealing features such as sparse multiscale representation and fast algorithms, wavelets have enjoyed many tremendous successes in the areas of signal/image processing and computational mathematics. This paper primarily intends to shed some light on the advantages of using wavelets in the context of numerical differential equations. We shall identify a few prominent problems in this field and recapitulate some important results along these directions.
Wavelet-based methods for numerical differential equations offer the advantages of sparse matrices with uniformly bounded small condition numbers.
We shall demonstrate wavelets' ability in solving some one-dimensional differential equations: the biharmonic equation and the Helmholtz equation with high wave numbers (of magnitude (10^4) or larger).





[
    Bin Han, Michelle Michelle, and Yau Shu Wong
    December 30, 2023
================================================



arabic



§ INTRODUCTION AND MOTIVATIONS

Being unconditional bases in many function spaces,
wavelets are sparse multiscale representation systems, which serve as excellent approximation tools for various types of functions and signals. The theory of wavelet analysis can be traced back as early as Haar's discovery in 1909, which is now known as the Haar orthogonal wavelet system. It is the simplest example of a wavelet, since it takes the form of an indicator function supported on the unit interval. However, wavelet analysis only experienced a surge of interest and rapid development after the emergence of
Morlet <cit.>, Meyer <cit.>, and Daubechies <cit.> wavelets as well as the notion of multiresolution analysis <cit.>. At the present time, the theory of wavelet analysis is considerably well-established (e.g., see <cit.>) and its applications are far-reaching. In the context of computational mathematics, wavelets have demonstrated their ability in tackling problems related to signal/image processing and numerical differential equations (DEs). The main focus of this paper is to explore the upper hand we gain from utilizing wavelets in the context of numerical DEs. We refer interested readers to the recent book <cit.> and references therein for a similar discussion in the context of signal and image processing.

To elucidate the concept of wavelets, we recall some basic facts and present some examples. Define :=∪{0} and let m ∈. Recall that the Sobolev space H^m() contains all functions f on  such that f, f', …, f^(m-1) are absolutely continuous on  and f, …, f^(m)∈ L_2(). When m=0, we naturally have H^0() = L_2(). Let ϕ=(ϕ_1,…,ϕ_r)^ and ψ=(ψ_1,…,ψ_s)^ be vectors of tempered distributions or functions on . For J∈, a wavelet affine system for H^m() is defined as


    _J^m(ϕ;ψ):=
       {2^J(1/2-m)ϕ_ℓ(2^J ·-k)  k∈, 1≤ℓ≤ r}
         ∪{2^j(1/2-m)ψ_ℓ(2^j·-k)  j≥ J, k∈, 1≤ℓ≤ s}.


We define _J(ϕ;ψ):=_J^0(ϕ;ψ) for m=0. We often refer j and k in (<ref>) as the scale/resolution level and the integer shift, respectively. Hence, a wavelet affine system is simply a set that contains dilated and translated versions of functions in all the entries of ϕ and ψ. We shall see later that the wavelet function ψ is also generated from a linear combination of dilated and translated versions of the refinable (vector) function ϕ. Note that ϕ in Figure <ref> is the hat function or B-spline of order 2. In fact, B-splines and Hermite splines are two popular examples of refinable functions, ϕ, due to their analytic expressions and wide usage in computational mathematics.
The wavelet affine system _J^m(ϕ;ψ) is a Riesz basis for m if the following two conditions are satisfied: (1) the linear span of _J^m(ϕ;ψ) is dense in m, and (2) there exist positive constants C_1 and C_2 such that


    C_1 (∑_ℓ=1^r ∑_k∈ |v_ℓ,k|^2
       +∑_j=J^∞∑_ℓ=1^s ∑_k∈
    |w_ℓ,j;k|^2) ≤∑_ℓ=1^r ∑_k∈ v_ℓ,k 2^J(1/2-m)ϕ_ℓ(2^J ·-k)
       + ∑_j=J^∞∑_ℓ=1^s ∑_k∈ w_ℓ,j;k 2^j(1/2-m)ψ_ℓ(2^j·-k)^2_m≤ C_2 ( ∑_ℓ=1^r ∑_k∈ |v_ℓ,k|^2
    + ∑_j=J^∞∑_ℓ=1^s ∑_k∈ |w_ℓ,j;k|^2)



for all finitely supported sequences { v_ℓ,k}_1≤ℓ≤ r,k∈ and {w_ℓ,j;k}_1≤ℓ≤ s,j≥ J, k∈. The best possible constants C_1 and C_2 are respectively called the lower and upper Riesz bounds of _J^m(ϕ;ψ).
It is known in <cit.> that _J^m(ϕ;ψ) is a Riesz basis for H^m() if and only if it is a Riesz basis for H^m() for all J∈.
Hence, we call {ϕ;ψ} a Riesz wavelet in the Sobolev space m if _0^m(ϕ;ψ) is a Riesz basis for H^m().
Moreover, the ratio C_2/C_1 of its Riesz bounds is called the condition number of the Riesz wavelet {ϕ;ψ} (or of the Riesz wavelet basis
_0^m(ϕ;ψ)) in the Sobolev space m.
See Figure <ref> for sample elements taken from such Riesz wavelet bases.
If r=1, we often call {ϕ;ψ} a scalar Riesz wavelet. Meanwhile, if r>1, we often call {ϕ;ψ} a Riesz multiwavelet due to the fact that ϕ is a vector function. For the sake of convenience, we simply use wavelets or Riesz wavelets to refer to both scalar Riesz wavelets and Riesz multiwavelets.



We now discuss why it is favourable to employ wavelets in solving numerical DEs. Wavelets are commonly used in the finite element/Galerkin method. Firstly, we have a well-established theory to guarantee that they completely span the space of interest. Thus, we can confidently conclude that the approximated solution approaches the true solution as the relative error decreases. Secondly, Riesz wavelets are well-conditioned bases, which give rise to coefficient matrices with uniformly bounded (possibly) small condition numbers. This property is essential for fast convergence in numerical schemes. Moreover, one standard technique to handle multidimensional DEs is to take the tensor product of a univariate Riesz wavelet with a small condition number. This will prevent the condition number from growing too rapidly with respect to the dimension. In the finite element method, we typically use the hat function as our basis. We shall demonstrate (see Table <ref>) that the condition number of a wavelet basis is far smaller than that of the hat functions (i.e., the standard finite element method), even though both of them generate the same space. In fact, the later condition number seems to go unbounded as the scale level increases. Thirdly, there is a lot of flexibility in designing wavelets. The speed of convergence in the continuous Galerkin method hinges on the polynomial reproduction order of the employed basis. Hence, we can control the speed by simply picking a (primal) refinable function with a suitable polynomial reproduction order. Given this choice, we can construct a (dual) refinable function having a prescribed polynomial reproduction order, which in turn affects the wavelet's vanishing moments. The latter plays an indispensable role in enhancing the sparsity of the coefficient matrix. Typically, we may even end up with some freedom to construct our wavelet basis so that ratio of the Riesz bounds (i.e., the condition number of the Riesz wavelet basis) is as small as possible. Note that high polynomial reproduction order often comes at the cost of longer supports. However, we have the ability to strike the right balance depending on the problem in hand. Once a wavelet basis is aptly chosen, we can use a suitable construction procedure (e.g., <cit.>) to ensure the boundary conditions (e.g., Dirichlet, Neumann, etc.) are satisfied. In some instances, wavelets can be constructed so that in order to obtain the numerical solution, no linear system needs to be solved. The fourth advantage that wavelets offer comes in the assembly of coefficient matrices in the continuous Galerkin method. In some cases, there is effectively no need to use quadrature to calculate the inner products of the basis functions. Such information can be efficiently obtained from the eigenvector of the transition operator and exploiting the refinability structure of our wavelet basis (i.e., via the fast wavelet transform). Furthermore, we can design an efficient quadrature to compute the inner products of our wavelet basis and the source term as in <cit.>.

There are two fundamental problems in wavelet-based methods for numerical DEs, which act as the focal points of our discussion. The first one is in the construction of suitable wavelets to handle the problem in hand. Having a sparse coefficient matrix is very much desired in solving numerical DEs. Motivated by the model problem

    u^(2m)(x)+ α u(x)=f(x),    x ∈ℐ,

where α∈ and ℐ is a bounded interval on , <cit.> provides the necessary and sufficient conditions for the construction of wavelets whose m-th order derivatives are orthogonal. Prior to <cit.>, these conditions had been unknown. However, examples of such wavelets have long been known in the literature <cit.>. Recall that a Riesz wavelet {ϕ;ψ} in the Sobolev space H^m() is called an m-th order derivative orthogonal Riesz wavelet in H^m() if it satisfies the following two properties:


    ψ^(m), ϕ^(m)(·-k)=0,    ∀  k∈,


and


    ψ^(m)(2^j·-k), ψ^(m)(2^j'·-k')=0,    ∀  k,k'∈, j,j'∈    j j'.


It is not hard to see that the weak formulation of (<ref>) is

    ⟨ u^(2m),v ⟩
    +α u,v
    = ∑_k=0^m-1 (-1)^k.(u^(2m-k-1)v^(k))|_ℐ + (-1)^m⟨ u^(m),v^(m)⟩ + α⟨ u,v ⟩ = ⟨ f,v ⟩

for all v ∈ H^m(ℐ). Indeed, employing an m-th order derivative-orthogonal Riesz wavelets gives rise to a stiffness matrix with a nice block diagonal matrix structure. On the other hand, the corresponding mass matrix has a sparse finger-like structure. See Figure <ref>. Let κ be the condition number of a matrix; i.e., κ is the ratio of the largest and smallest singular values. For a symmetric positive definite matrix, κ is simply the ratio of the largest and smallest eigenvalues. The contribution of κ((-1)^m⟨ u^(m),v^(m)⟩) in the condition number of the coefficient matrix, κ((-1)^m⟨ u^(m),v^(m)⟩ + α⟨ u,v ⟩), often overpowers that of κ(⟨ u,v ⟩). Hence, κ((-1)^m⟨ u^(m),v^(m)⟩ + α⟨ u,v ⟩) is typically small (i.e., well-conditioned). We shall revisit this result in Section 2.

In the same vein, <cit.> impose an orthogonality condition on the first derivative of the wavelets (generated from Hermite cubic splines) to achieve better sparsity. However, their methods differ from <cit.> in that they substitute the wavelet functions at even integers with a pair that have a longer support.

Just like in many other applications, numerical DEs are typically defined on a bounded domain. For rectangular bounded domains, one strategy is to take the tensor product of wavelets constructed on a bounded interval. Hence, one critical question that we need to resolve is how to optimally construct wavelets on a bounded interval. This serves as the second fundamental problem. Our starting point is a Riesz wavelet on the real line. Ideally, we want to have a tractable algorithm at our disposal that gives us boundary elements with simple structure and allows us to retain all desirable properties that our starting wavelet has. Additionally, we want to ensure that the wavelet formed on the interval of interest indeed spans the whole space. This is in fact a long-standing problem in wavelet analysis. The fact that the analysis of multiwavelets is far more intricate than their scalar counterparts just adds to the difficulty. Nevertheless, any construction we propose needs to be applicable to multiwavelets, since they have a key advantage over scalar wavelets in that they generally have higher smoothness and vanishing moments for a given support. In other words, multiwavelets may attain the desired smoothness and vanishing moments with a much shorter support. Having a short support is indeed necessary to produce boundary elements with simple structure. Tremendous research efforts have gone into developing various construction procedures, but virtually all of them suffer from some serious shortcomings in one form or another. It is not until recently that <cit.> presents a general construction that gives rise to all possible square integrable compactly supported biorthogonal wavelets on a semi-infinite or bounded interval satisfying maximum vanishing moments and/or boundary conditions. Many existing constructions in the literature serve as special cases of the construction proposed in <cit.>. We shall elaborate on this further in Section 3.

Next, we mention other equally significant studies in the topic of wavelet-based methods for numerical DEs. First off, the book <cit.> provides a comprehensive treatment of the subject. Some key references in relation to the construction of wavelets on general bounded domains for numerical DEs are <cit.>. In addition to the two important research directions described above, a substantial amount of noteworthy work has gone towards developing adaptive wavelet algorithms. Such a scheme is useful in handling problems with singularities and nonlinear problems. There is a subtle distinction between adaptive finite element and adaptive wavelet methods. The first one directly relies on iterative mesh refinements driven by a posteriori local errors; meanwhile the latter involves iterative trial space refinements by proper inclusions of extra basis functions <cit.>. Two of the seminal papers in this field are <cit.>. Further developments of adaptive wavelet algorithms, to name a few, are reported in <cit.> for hyperbolic problems, <cit.> for high dimensional elliptic problems, and <cit.> for parabolic evolution problems. The use of wavelets in efficiently dealing with integral operators can be traced back to the influential papers <cit.> and more recent results can be found in <cit.>. Multilevel preconditioning techniques for linear systems stemming from Galerkin methods were studied in <cit.>.

Towards the end of this paper, we shall consider the Helmholtz equation, which is a widely used model in wave propagation. This equation appears in seismic imaging (e.g., full waveform inversion), acoustics, and electromagnetism. One predominant challenge in solving the Helmholtz equation is caused by the infamous pollution effect; i.e., small enough hk (here h and k refer to the mesh size and wave number respectively) provides no guarantee that the relative error has a bound independent of k <cit.>. To make matters worse, the linear system associated with the discretization process is ill-conditioned. There is an abundance of numerical schemes in the literature that deal with the Helmholtz equation with high wave numbers. For example, see <cit.> and references therein. To this day, developing a numerical scheme for the Helmholtz equation continues to be an active research area. In this paper, we aim to examine a wavelet-based method's potential in solving the one-dimensional Helmholtz equation with wave numbers of magnitude (10^4) or larger.

The paper is organized as follows. In Section 2, we shall review the main result of <cit.>. In Section 3, we shall revisit the main result of <cit.> and provide a new example of a biorthogonal wavelet formed on the unit interval with homogeneous Dirichlet boundary conditions and maximum vanishing moments. In Section 4, we shall use this example to solve the Helmholtz equation with high wave numbers and consider the biharmonic equation. Finally, we outline some future work in Section 5.



§ DERIVATIVE-ORTHOGONAL RIESZ WAVELETS IN H^M()

In this section, we review the main result presented in <cit.>. For complete details and proofs, we refer interested readers to the foregoing paper.

To facilitate our discussion, we introduce some definitions. The Fourier transform of f∈1 is defined as f(ξ):=∫_ f(x) e^-ixξ dx for all ξ∈. We then extend the definition to tempered distributions by means of duality. Let τ∈. The Sobolev space H^τ() contains all tempered distributions on  such that

    f ^2_H^τ() := 1/2 π∫_ |f(ξ)|^2(1+|ξ|^2)^τ dξ < ∞.

Let (l_0())^r × s be the space of all finitely supported matrix-valued filters a: →^r × s. For a ∈ (l_0())^r × s, we define a(ξ):=∑_k ∈ a(k) e^-ikξ, which is an r× s matrix of 2π-periodic trigonometric polynomials. A wavelet affine system and a Riesz wavelet in the Sobolev space H^τ() are defined the same way as in Section 1; we only need to replace m ∈∪{0} with τ∈. A compactly supported Riesz wavelet {ϕ;ψ} is typically derived from a compactly supported refinable vector function via the following refinability structure:

    ϕ(2ξ) = a(ξ) ϕ(ξ)   and  ψ(2ξ) = b(ξ) ϕ(ξ),   a.e. ξ∈

for some filters a ∈ (l_0())^r× r and b ∈ (l_0())^s× r.
For J∈, we say that (_J^-τ(ϕ̃;ψ̃), _J^τ(ϕ;ψ)) is a biorthogonal wavelet in (H^-τ(),H^τ()) if (1) ^τ_J(ϕ;ψ) is a Riesz basis for H^τ() and ^-τ_J(ϕ̃;ψ̃) is a Riesz basis for H^-τ(), as well as (2) ^τ_J(ϕ;ψ) and ^-τ_J(ϕ̃;ψ̃) are biorthogonal to each other.
It is known in <cit.> that
(_J^-τ(ϕ̃;ψ̃), _J^τ(ϕ;ψ)) is a biorthogonal wavelet in (H^-τ(),H^τ()) for some J∈ if and only if (_J^-τ(ϕ̃;ψ̃), _J^τ(ϕ;ψ)) is a biorthogonal wavelet in (H^-τ(),H^τ()) for all J∈. Consequently,
we say that ({ϕ̃;ψ̃},{ϕ;ψ}) is a biorthogonal wavelet in (H^-τ(),H^τ()) if (_0^-τ(ϕ̃;ψ̃), _0^τ(ϕ;ψ)) is a biorthogonal wavelet in (H^-τ(),H^τ()).

The bracket product is a key tool in the analysis of shift-invariant spaces.
Let f ∈ (H^τ())^r × t and g ∈ (H^-τ())^s × t be vector functions. For τ∈, we define the bracket product by

    [f,g]_τ(ξ):= ∑_k ∈f(ξ + 2π k) g(ξ + 2π k)^(1+|ξ|^2)^τ,   ξ∈.

When τ =0, [f,g](ξ):=[f,g]_0(ξ). We say that the integer shifts of a vector function ϕ∈ (H^τ())^r× 1 are stable in H^τ() if there exists a positive constant C such that

    C^-1 I_r ≤ [ϕ,ϕ]_τ(ξ)≤ C I_r,      a.e. ξ∈.

If ϕ∈ H^τ() has compact support, then the integer shifts of ϕ are stable in H^τ() if and only if
span{ϕ(ξ + 2π k)  k∈}=^r for all ξ∈ (see <cit.>).
Next, we introduce the definitions of sum rules and vanishing moments. The former is closely related to polynomial reproduction orders, while the latter is linked to polynomial annihilation. For a given filter a ∈ (l_0())^r × r, we say that a satisfies order m sum rules if there is a matching filter v ∈  (l_0())^1 × r such that v(0) ≠ 0,

    v(2ξ) a(ξ) = v(ξ) + (|ξ|^m),   v(2ξ)a(ξ + π)=(|ξ|^m),   ξ→ 0.

We also use the notation (a) to denote the highest order of sum rules satisfied by a filter a. Here, f(ξ) = g(ξ) + (|ξ|^m) as ξ→ 0 is equivalent to saying that f^(j)(0) = g^(j)(0) for all j=0,…,m-1. On the other hand, given a compactly supported function ψ∈ (L_2())^s × 1, we say that ψ has order n vanishing moments if ⟨ψ,x^j ⟩ = ∫_ψ(x) x^j dx =0 for all j=0,…,n-1; or in other words, ψ(0)=…=ψ^(n-1)(0)=0. Notation-wise, we write (ψ) to indicate the highest vanishing moments attained by ψ.

Let ϕ be a tempered distribution on . Recall that the smoothness/regularity of ϕ is measured by its smoothness exponent:

    (ϕ):=sup{τ∈ϕ∈ H^τ()}.

By convention, we set (ϕ):=-∞ if {τ∈ϕ∈ H^τ()} = ∅.

The theorem immediately below provides a full characterization of Riesz wavelets in the Sobolev space H^τ(), where τ∈. Riesz wavelets and dual wavelet frames in the Sobolev space H^τ(), where τ∈, were initially studied in <cit.>.

 
	(a special case of <cit.>)
	Let a, b, ã, b̃∈ (l_0())^r × r.
Assume that


  (i)
1 is a simple eigenvalue of a(0) and (2^jI_r-a(0))≠ 0 for all j∈;

  (ii) 1 is a simple eigenvalue of ã(0) and (2^jI_r-ã(0))≠ 0 for all j∈.

Let ϕ, ϕ̃ be r × 1 vectors of compactly supported distributions satisfying

    ϕ(2ξ) = a(ξ)ϕ(ξ),   ϕ̃(2ξ) = ã(ξ)ϕ̃(ξ).

Define ψ and ψ̃ by
	
    ψ(ξ) = b(ξ/2)ψ(ξ/2),   ψ̃(ξ) = b̃(ξ/2)ψ̃(ξ/2).

	For τ∈, the pair ({ϕ̃;ψ̃},{ϕ;ψ}) is a biorthogonal wavelet in (H^-τ(),H^τ()) if and only if


  (1) ({ã;b̃},{a;b}) is a biorthogonal wavelet filter bank; that is,
	
    [ ã(ξ)    ã(ξ+π)
    b̃(ξ)    b̃(ξ+π)
    	]
    	[ a(ξ)^   b(ξ)^
    a(ξ+π)^   b(ξ+π)^]=
    	I_2r,


  (2) ϕ(0)^ϕ̃(0)=1, ϕ∈ (H^τ())^r× 1 and ϕ̃∈ (H^-τ())^r× 1.

  (3) The integer shifts of ϕ and ϕ̃ are biorthogonal to each other:


ϕ̃, ϕ(·-k):=
∫_ ϕ̃(x) ϕ(x-k)^dx=(k) I_r,    ∀  k∈,

where (0)=1 and (k)=0 for all k 0.

  (4) ψ(ξ) = o(|ξ|^-τ) as ξ→ 0 (i.e., (ψ) > -τ) if τ≤ 0, and ψ̃(ξ) = o(|ξ|^τ) as ξ→ 0 (i.e., (ψ̃) > τ) if τ≥ 0.



The conditions in items (i) and (ii) of Theorem <ref> are not assumed in <cit.>. It is well known (e.g., see <cit.>)
that the condition in item (i) guarantees the existence and uniqueness (up to a multiplicative constant) of a compactly supported refinable vector function ϕ satisfying ϕ(2ξ)=a(ξ)ϕ(ξ).
Item (4) puts conditions on vanishing moments. For τ>0, no vanishing moments are required for ψ at all, while for τ<0,
no vanishing moments are required for
ψ̃.
If ϕ∈ (H^τ())^r× 1 satisfies ϕ(2ξ)=a(ξ)ϕ(ξ) with a∈ (l_0())^r × r and τ≥ 0, then item (i) must hold if the integer shifts of ϕ are stable in H^τ().
Note that the biorthogonality condition in (<ref>) is equivalent to saying that [ϕ̃,ϕ](ξ)=I_r for almost every ξ∈.
Since ϕ∈ (H^τ())^r× 1 and ϕ̃∈ (H^-τ())^r× 1 have compact support, item (3) implies that the integer shifts of ϕ are stable in H^τ(), while the integer shifts of ϕ̃ are stable in H^-τ().
Moreover, the biorthogonality condition (<ref>) in item (3) can be fully characterized by (a)>τ and (ã)>-τ, where the smoothness exponent (a) is defined in <cit.>.
For more details, see <cit.> and <cit.>.

Built on Theorem <ref>,
the following result
characterizes derivative-orthogonal Riesz wavelets.



(<cit.>)
Let ϕ=(ϕ_1,…,ϕ_r)^ be a compactly supported refinable vector function in m with m∈ such that ϕ(2ξ)=a(ξ)ϕ(ξ) for some a∈0rr. Then


  (i) there exists a finitely supported high-pass filter b∈0rr such that {ϕ;ψ} with ψ(ξ):=b(ξ/2)ϕ(ξ/2) is an mth-order derivative-orthogonal Riesz wavelet in the Sobolev space m satisfying (<ref>) and (<ref>) if and only if the integer shifts of ϕ are stable and the filter a has at least order 2m sum rules (i.e., (a)≥ 2m).

  (ii) Under the condition that the integer shifts of ϕ are stable, for any b∈0rr, {ϕ;ψ} with ψ(ξ):=b(ξ/2)ϕ(ξ/2) is an mth-order derivative-orthogonal Riesz wavelet in the Sobolev space m  satisfying (<ref>) and (<ref>) if and only if
 

    b(ξ) [ϕ^(m),ϕ^(m)](ξ)
    a(ξ)^+b(ξ+π) [ϕ^(m),ϕ^(m)](ξ+π)
    a(ξ+π)^=0


and


    ({a;b})(ξ):=([ a(ξ)    a(ξ+π)
    b(ξ)    b(ξ+π)
    ]) 0,     ∀ ξ∈.


are satisfied. Moreover, for every J∈, ^τ_J(ϕ;ψ) is a Riesz basis in the Sobolev space τ for all τ in the nonempty open interval (2m-(ϕ), (ϕ)) with (ψ)=(a)-2m.



For scalar filters, there is in fact an explicit formula to find their corresponding m-th order derivative orthogonal Riesz wavelets. See <cit.> for details. To make the presentation of this paper self-contained, we reproduce the second-order derivative-orthogonal Riesz wavelet generated from Hermite cubic splines, since it will be used in Section 4.2.

 
	Let ϕ:=(ϕ_1,ϕ_2) be the well-known Hermite cubic splines with ϕ_1 and ϕ_2 given below
	
    ϕ_1(x)=
    	(1-x)^2(1+2x),    x∈ [0,1],
    
    	(1+x)^2(1-2x),    x∈ [-1,0),
    
    	0,    otherwise,
    	    ϕ_2(x)=
    	(1-x)^2x,    x∈ [0,1],
    
    	(1+x)^2x,    x∈ [-1,0),
    
    	0,    otherwise.

	The integer shifts of ϕ are indeed stable and (ϕ)=2.5. Moreover, (a)=4, where the filter a is defined below. Note that ϕ, ψ satisfy equations ϕ(2ξ) = a(ξ)ϕ(ξ) and ψ(ξ) = b(ξ/2)ϕ(ξ/2) respectively with the filters a,b ∈ (l_0())^2 × 2 defined as follows:
	
    a={[   14   38; -116 -116 ],  [ 12  0;  0 14 ],   [   14  -38;  116 -116 ]}_[-1,1],   
    	b={[ 12  0;  0 12 ]}_[1,1].

	Then, {ϕ;ψ} is a second-order derivative-orthogonal Riesz wavelet in H^2(). Also, the wavelet affine system _0^τ(ϕ;ψ) is a Riesz basis in H^τ() for all τ∈ (3/2,5/2). See Figure <ref> for plots of ϕ and ψ.
	
	


Other examples of zeroth, first, and second-order derivative-orthogonal Riesz wavelets generated from B-splines and Hermite splines are presented in <cit.>.



§ CONSTRUCTION OF BIORTHOGONAL WAVELETS IN L_2([0,N]) WITH MAXIMUM VANISHING MOMENTS AND BOUNDARY CONDITIONS


In this section, we shift our attention to a general construction of compactly supported biorthogonal wavelets in L_2([0,N]) with N∈. For complete details and proofs, we refer readers to <cit.>.

First and foremost, we review relevant studies that precede ours. A pioneering work in the construction of Daubechies orthogonal wavelets in L_2([0,1]) can be attributed to Meyer <cit.>, which was followed up by <cit.> and other references. A more general construction of orthogonal wavelets in L_2([0,N]) is available in <cit.>. The use of orthogonal wavelets in L_2([0,N]) within the framework of numerical DEs is not preferable for two reasons. Firstly, they often do not have any analytic expressions with a few exceptions like the Haar wavelet. Even then, this particular example has low smoothness. Secondly, their orthogonality property often conflicts with the boundary conditions. I.e., orthogonal wavelets often fail to satisfy the boundary conditions. If orthogonal wavelets in L_2([0,N]) would satisfy a given boundary condition, then they cannot have high orders of vanishing moments, which are the key property for sparse representations.
Consequently, implementing orthogonal wavelets in a numerical scheme becomes a cumbersome task. In contrast, biorthogonal wavelets are preferable mainly because they do not suffer from the previous disadvantages. Within the scope of biorthogonal wavelets, some existing constructions can solely be applied to a particular refinable function ϕ; e.g., Hermite cubic splines in <cit.>. Scalar biorthogonal wavelets generated from B-splines can be formed on the unit interval by utilizing the methods in <cit.>. Additionally, the method explained in <cit.> yields spline wavelets on the unit interval with homogeneous boundary conditions of arbitrary order. On a similar note, <cit.> provides a means to construct biorthogonal wavelets satisfying complementary boundary conditions. The constructions proposed in <cit.> rely on an operator introduced in <cit.>. They can be applied to a relatively large class of wavelets and framelets (i.e., a generalization of wavelets with redundancy) satisfying some symmetry property. However, these constructions sacrifice the vanishing moments near the boundary in order to fulfill the boundary conditions. Limited applicability as well as inability to concurrently retain maximum vanishing moments and meet the boundary conditions are two common deficiencies in the above constructions. This is where our contribution comes in.

To aid the discussion of our main result, we first introduce a notation and an important result that describes the structure of compactly supported Riesz wavelets on [0,∞). Suppose f is a compactly supported vector function. Then, (f) is defined to be the shortest interval with integer endpoints such that f vanishes outside (f). Also, denote f_j;k:=2^j/2f(2^j · -k).

The theorem below essentially guarantees the existence of compactly supported vector functions ϕ̃^L, ψ̃^L given compactly supported vector functions ϕ^L, ψ^L: a critical fact for the construction procedure under consideration.


(<cit.>)
	Let ({ϕ̃;ψ̃},{ϕ;ψ}) be a compactly supported biorthogonal wavelet in 2 with ϕ,ψ,ϕ̃,ψ̃∈ (2)^r. Then Theorem <ref> holds with τ =0.
	Define
	
	
	[l_ϕ, h_ϕ]:=(ϕ),  [l_ψ, h_ψ]:=(ψ),  [l_ϕ̃,
	 h_ϕ̃]:=(ϕ̃),  [l_ψ̃, h_ψ̃]:=(ψ̃)
	
	and
	
	
	[l_a, h_a]:=(a),  [l_b, h_b]:=(b),  [l_ã, h_ã]:=(ã),  [l_b̃, h_b̃]:=(b̃).
	
	
	Let ϕ^L and ψ^L be vectors of compactly supported functions in L_2([0,∞)). Define
	
	Φ:= {ϕ^L} ∪{ϕ(·-k): k ≥n_ϕ},   Ψ:= {ψ^L} ∪{ψ(·-k): k ≥n_ψ}
	with n_ϕ≥max(-l_ϕ,-l_a) and n_ψ≥max(-l_ψ, n_ϕ-l_b/2).
Define

    _J(Φ;Ψ)_[0,∞):={2^J/2φ(2^J·) φ∈Φ}∪{ 2^j/2η(2^j·)  j≥ J, η∈Ψ},     J∈.

	Suppose that _0(Φ;Ψ)_[0,∞) is a Riesz basis of L_2([0,∞))
	
	and satisfies
	
	
    ϕ^L=2A_Lϕ^L(2·)+2∑_k=n_ϕ^m_ϕ A(k) ϕ(2·-k),
    	
       ψ^L=2B_Lϕ^L(2·)+2∑_k=n_ϕ^m_ψ B(k) ϕ(2·-k),

	
	for some matrices A_L,B_L and finitely supported sequences A,B of matrices.
	Then
	
		
  (1) there must exist compactly supported vector functions ϕ̃^L, ψ̃^L in L_2([0,∞)) and integers n_ϕ̃≥max(-l_ϕ̃,-l_ã,n_ϕ) and n_ψ̃≥max(-l_ψ̃,n_ϕ̃-l_b̃/2,n_ψ) such that
		 _0(Φ̃;Ψ̃)_[0,∞) is the dual Riesz basis of _0(Φ;Ψ)_[0,∞) in L_2([0,∞)),
		where _0(Φ̃;Ψ̃)_[0,∞):=Φ̃∪{ 2^j/2η̃(2^j·)  j∈, η̃∈Ψ̃} with :=∪{0} and
		
		
		Φ̃:={ ϕ̃^L}∪{ ϕ̃(·-k) k≥n_ϕ̃},  Ψ̃:={ ψ̃^L}∪{ ψ̃(·-k) k≥n_ψ̃};
		
		
		
  (2) there exist matrices Ã_L,B̃_L and finitely supported sequences Ã,B̃ of matrices such that
		
		
    ϕ̃^L=2Ã_Lϕ̃^L(2·)+2∑_k=n_ϕ̃^m_ϕ̃Ã(k) ϕ̃(2·-k), 
       ψ̃^L=2B̃_Lϕ̃^L(2·)+2∑_k=n_ϕ̃^m_ψ̃B̃(k) ϕ̃(2·-k),

		
		and
		
		
    ϕ̃(·-k_0)=2∑_k=n_ϕ̃^∞ã(k-2k_0) ϕ̃(2·-k),    ∀  k_0≥ n_ϕ̃,
       ψ̃(·-k_0)=2∑_k=n_ϕ̃^∞b̃(k-2k_0) ϕ̃(2·-k),    ∀  k_0≥ n_ψ̃.

		
	


Before we present the complete algorithms, we shall outline several key steps to better illustrate how the construction operates. For the following, let ({ϕ̃;ψ̃},{ϕ;ψ}) be a compactly supported biorthogonal wavelet in L_2().

	
  (1) We commence by forming a compactly supported L_2([0,∞)) biorthogonal wavelet.
	
		
  (1.1) (Algorithm <ref>) Construct the space Φ as in (<ref>) satisfying a prescribed polynomial reproduction property (cannot be higher than that of ϕ) and the refinability condition.
		
  (1.2) (Algorithm <ref>) Compute the inner products between the shifted versions of ϕ and ϕ̃ restricted to [0,1]. This will help us to recover any inner products between the shifts of ϕ and ϕ̃.
		
  (1.3) (Algorithm <ref>) Construct the space Φ̃ as in (<ref>) satisfying a prescribed polynomial reproduction (cannot be higher than that of ϕ̃), the refinability condition, and the biorthogonality condition with respect to Φ.
		
  (1.4) (Algorithm <ref>) Construct the space Ψ as in (<ref>) such that all elements in Ψ are perpendicular to all elements in Φ̃. Additionally, we need to ensure that ψ^L and all interior primal wavelets, {ψ(· -k): k ≥ n_ψ}, are linearly independent. Construct the space Ψ̃ as in (<ref>) such that the biorthogonality condition with respect to Ψ is satisfied and all elements in Ψ̃ are perpendicular to all elements in Φ. Additionally, we need to ensure that ψ̃^L and all interior dual wavelets, {ψ̃(· -k): k ≥ n_ψ̃}, are linearly independent.
	
	
  (2) Repeat Step (1) for the reflected compactly supported L_2() biorthogonal wavelet ({ϕ̃,ψ̃},{ϕ,ψ}), where ϕ:=ϕ(-·), ψ:=ψ(-·), ϕ̃:=ϕ̃(-·), and  ψ̃:=ψ̃(-·), in order to obtain compactly supported Φ, Ψ, Φ̃, and Ψ̃.
	
  (3) Set J_0 to be the smallest nonnegative integer such that for all j ≥ J_0,
	
    max(m_ϕ+n_ϕ, m_ψ+n_ϕ, m_ϕ+n_ϕ, m_ψ+n_ϕ)≤ 2^j+1N

	and each element in Φ(2^j), Ψ(2^j), Φ̃(2^j), Ψ̃(2^j) does not essentially touch both endpoints 0 and N simultaneously. I.e., neither 0 nor N is an interior point of (h) for all h ∈Φ(2^j ·), Ψ(2^j ·), Φ̃(2^j ·), Ψ̃(2^j ·). Similarly, set J̃_0 to be the smallest nonnegative integer such that for all j ≥J̃_0
	
    max(m_ϕ̃+n_ϕ̃,
    	 m_ψ̃+n_ϕ̃,
    	 m_ϕ̃+n_ϕ̃,
    	 m_ψ̃+n_ϕ̃)
    	≤ 2^j+1N,

	each element in Φ(2^j), Ψ(2^j), Φ̃(2^j), Ψ̃(2^j) does not essentially touch both endpoints 0 and N simultaneously, as well as
	
    ψ^L_j;0, ϕ̃^R_j;2^jN-N=0,  ψ^L_j;0, ψ̃^R_j;2^jN-N=0,
    		  ψ̃^L_j;0, ϕ^R_j;2^jN-N=0,  ψ̃^L_j;0, ψ^R_j;2^jN-N=0,
    ϕ̃^L_j;0, ψ^R_j;2^jN-N=0,  ψ̃^L_j;0, ψ^R_j;2^jN-N=0,
    		  ϕ^L_j;0, ψ̃^R_j;2^jN-N=0,  ψ^L_j;0, ψ̃^R_j;2^jN-N=0,

	where
	
    ϕ^R:=ϕ^L(N-·),  ψ^R:=ψ^L(N-·),  ϕ̃^R:=ϕ̃^L(N-·),  ψ̃^R:=ψ̃^L(N-·).


By <cit.>, the above Steps (1) and (2) ensure that _J(Φ̃,Ψ̃)|_[0,∞) and _J(Φ,Ψ)|_[0,∞) form a pair of compactly supported Riesz bases in L_2([0,∞)). Similarly, _J(Φ̃,Ψ̃)|_[0,∞) and _J(Φ,Ψ)|_[0,∞) form a pair of compactly supported Riesz bases in L_2([0,∞)). Without loss of generality, suppose that J̃_0≥ J_0. Define

    Φ_j:={ϕ^L_j;0}∪{ϕ^R_j;2^jN-N}∪{ϕ_j;k n_ϕ≤ k≤ 2^jN-n_ϕ}, 
       Ψ_j:={ψ^L_j;0}∪{ψ^R_j;2^jN-N}∪{ψ_j;k n_ψ≤ k≤ 2^jN-n_ψ}, 
       Φ̃_j:={ϕ̃^L_j;0}∪{ϕ̃^R_j;2^jN-N}∪{ϕ̃_j;k n_ϕ̃≤ k≤ 2^jN-n_ϕ̃}, 
       Ψ̃_j:= {ψ̃^L_j;0}∪{ψ̃^R_j;2^jN-N}∪{ψ̃_j;k n_ψ̃≤ k≤ 2^jN-n_ψ̃}.

Furthermore, define _J:=Φ_J ∪{Ψ_j   j≥ J}, and _J:=Φ̃_J ∪{Ψ̃_j   j≥ J}. Following the above three key steps, we have by <cit.> that (_J,_J) constitutes a pair of biorthogonal Riesz bases of L_2([0,N]) for all J ≥J̃_0. Moreover, there exist matrices A_j, B_j, Ã_j, B̃_j such that the following refinable structures hold

    Φ_j = A_jΦ_j+1,   Ψ_j = B_jΦ_j+1,   Φ̃_j = Ã_jΦ̃_j+1,   Ψ̃_j = B̃_jΦ̃_j+1,

and [A_j^,B_j^] is an invertible square matrix with

    [ Ã_j
    B̃_j]=[A_j^, B_j^]^-1

for j ≥ J. Note that orthogonal wavelets can also be constructed by the foregoing procedure. We now present Algorithms <ref> to <ref> in their full form.









Various examples of compactly supported biorthogonal wavelets in L_2([0,1]) are available in <cit.>. We next turn to a compactly supported biorthogonal scalar wavelet on L_2([0,1]) satisfying homogeneous Dirichlet boundary conditions and maximum vanishing moments. Such an example will be used later in our numerical experiments for solving Helmholtz equations.

 
	Consider the scalar biorthogonal wavelet ({ϕ̃;ψ̃},{ϕ;ψ}) in <cit.> satisfying (<ref>) and (<ref>) with ϕ(0)=ϕ̃(0)=1 and an associated
	biorthogonal wavelet filter bank ({ã;b̃},{a;b}) given by
	
	
    a=   {14,12,14}_[-1,1],    b={-18,-14,34,-14,-18}_[-1,3],
    ã=   {-18, 14, 34, 14, -18}_[-2,2],   b̃={-14, 12, -14}_[0,2].

	Note that ϕ is the hat function or the centered B-spline of order 2, B_2, which is defined as follows
	
	ϕ(x)=B_2(x)=
	1+x,    x∈[-1,0),

	1-x,    x∈[0,1],

	0,    otherwise.
	
	By calculation, we have (ϕ)=(a)=1.5, (ϕ̃)=(ã) ≈ 0.440765 and (a)=(ã)=2, where
(a) is defined in <cit.>. Using Algorithm <ref> with (x)=x, we have n_ϕ=3 and
	
	ϕ^L=A_c ϕ^c      ϕ^c=[   ϕ(·-2);   ϕ(·-1); ϕχ_[0,∞) ],
	  A_c=[ 1 0 0; 0 1 0;   ].
	Setting (x)=x is what helps us to fulfill the homogeneous Dirichlet boundary condition. Similarly, using Algorithm <ref> with (x)=x, we have n_ϕ=3 and (<ref>) with ϕ being replaced by ϕ. Using Algorithm <ref> with m=(ã)=2 and n_ϕ̃=3, we have
	
	ϕ̃^L=Ã_c ϕ̃^c      ϕ̃^c=
	[ 
	ϕ̃(·-2)

	ϕ̃(·-1) χ_[0,∞)

	ϕ̃ χ_[0,∞)

	ϕ̃(·+1) χ_[0,∞)
	 ],   Ã_c =
	[ [  1  0 -1 -2;  0  1  2  3 ] ].
	Setting m=(ã)=2 is what helps us to achieve the maximum vanishing moments. Similarly, using Algorithm 3 with m=(ã)=2 and n_ϕ̃=3, we have (<ref>) with ϕ̃^L, ϕ̃^c, Ã_c, ϕ̃ being replaced by ϕ̃^L, ϕ̃^c, Ã_c, ϕ̃ respectively. Moreover, ϕ^L, ψ^L, ϕ̃^L, and ψ̃^L satisfy
	
    ϕ^L =     2[   0   0; 1/2 1/4 ]ϕ^L(2·)
    	+ 2[ 1/4; 1/4 ]ϕ(2·-3)
    	+ 2 [ 12;  0 ]ϕ(2·-4)
    	+ 2[ 1/4;   0 ]ϕ(2·-5), 
    ϕ̃^L=     2 [ [    0 -1/4;  1/2  3/4 ]] ϕ̃^L(2·)
    	+ 2 [ [ 1/4; 1/4 ]] ϕ̃(2·-3)
    	+ 2 [ [      3/4; -
    			1/8 ]] ϕ̃(2·-4)
    	+ 2 [ [ 1/4;   0 ]] ϕ̃(2·-5) 
       + 2 [ [ -1/8;    0 ]] ϕ̃(2·-6), 
    ψ^L =     2[   b(0)  b(-1); -27/64 37/128 ]ϕ^L(2·)
    	+ 2 [ b(1); 1/64 ]ϕ(2·-3)
    	+ 2[ b(2); 5/64 ]ϕ(2·-4)
    	+ 2 [  b(3); 5/128 ]ϕ(2·-5), 
    ψ̃^L=     2 [ [ -5/32 -3/16;  -1/2     1 ]] ϕ̃^L(2·)
    	+ 2 [ 1/2;   0 ]ϕ̃(2·-3)
    	+ 2 [ -1/4;    0 ]ϕ̃(2·-4).

	Additionally, ϕ^L, ϕ̃^L respectively satisfy (<ref>) and (<ref>) with ϕ, ϕ̃, ϕ^L, ϕ̃^L being properly replaced by ϕ, ϕ̃, ϕ^L, ϕ̃^L. On the other hand, ψ^L and ψ̃^L satisfy
	
    ψ^L =     2[   2/7     0; -9/28   1/4 ]ϕ^L(2·)
    	+ 2 [  -1/2; -1/16 ]ϕ(2·-3)
    	+ 2[  1/7; 5/56 ]ϕ(2·-4)
    	+ 2 [  1/14; 5/112 ]ϕ(2·-5),  
    ψ̃^L=     2 [ [ 5/16 3/32; -1/2  5/4 ]] ϕ̃^L(2·)
    	+ 2 [ -23/32;   -1/4 ]ϕ̃(2·-3) + [ 23/64;   1/8 ]ϕ̃(2·-4).

	
	Let _J=Φ_J ∪{Ψ_j  j≥ J} for all J≥ 2, where Φ_j and Ψ_j are defined in (<ref>) and (<ref>) respectively, with n_ϕ=n_ϕ=n_ψ=3, n_ψ=2, ϕ^R=ϕ^L(1-·), and ψ^R=ψ^L(1-·). Let _J=Φ̃_J ∪{Ψ̃_j  j≥ J} for all J ≥ 3, where Φ̃_j and Ψ̃_j are defined in (<ref>) and (<ref>) respectively, with n_ϕ̃=n_ϕ̃=n_ψ̃=3, n_ψ̃=2, ϕ̃^R=ϕ̃^L(1-·), and ψ̃^R=ψ̃^L(1-·). According to <cit.> with N=1, (_J,_J) form a biorthogonal Riesz basis of L_2([0,1]) for every J ≥ 3. Note that (ψ^L)=(ψ^R)=(ψ)=2=(ã). However, (ψ̃^L)=(ψ̃^R)=0 and (ψ̃)=2=(a). See Figure <ref> for the generators of (_J,_J) for J ≥ 3.
	
	As an implication of <cit.>, (_2,_2) in Example <ref> is also a biorthogonal Riesz basis of L_2([0,1]). Though _2 is unchanged and is still given by (<ref>) and (<ref>), we point out that the boundary dual functions in _2 may be different from (<ref>), (<ref>), and (<ref>).
	
	It is indeed possible to replace the second component of ψ^L with one that has a shorter support and symmetry. The same thing applies to ψ^R=ψ^L(1-·). However, our calculation suggests that such a choice leads to a much larger condition number, assuming the first components of ψ^L, ψ^L are the same as our present choice.






§ NUMERICAL EXAMPLES USING RIESZ WAVELETS IN L_2([0,1])


In this section, we consider two different model problems. One of which is the Helmholtz equation and the other is the biharmonic equation. Without loss of generality, we consider ℐ=[0,1] as our domain of interest. 

The performance is measured by the relative L_2 error given by

    e_N_L_2/ u_L_2:= u_N -u_L_2/ u_L_2         e_N:=u_N-u,

where u is the true solution and u_N is the numerically computed approximated solution with N being its corresponding finest scale level for computing u_N. The meaning of N will be further clarified below.
Because both u and u_N in our Helmholtz examples are highly oscillating functions, for the purpose of high accuracy and fair comparison, we symbolically compute both u_N-u_L_2 and u_L_2 with extremely high precision.

Let us first exemplify some advantages of wavelet-based methods over the standard finite element method by comparing the condition numbers of the following system _2,N truncated at the finest scale levels N of the wavelet basis _2 in L_2([0,1]) constructed in Example <ref> of Section 3:


_2,N:=Φ_2 ∪{Ψ_j 2≤j≤N-1},    N≥2

to the condition numbers of the standard finite element method using only the shifts of the refinable function ϕ=B_2 defined in (<ref>) as follows:


    _N:={B_2(2^N·-k)  k=1,…,2^N-1},     N≥ 2.


Note that these two systems _2,N and _N
satisfy the homogeneous Dirichlet boundary conditions on [0,1] and generate the same finite-dimensional space for all N≥ 2. See Table <ref> for details.



	
		Scale N    Size
    Mass (FEM)     Stiffness (FEM)     Mass (Wavelet)     Stiffness (Wavelet) 

		11    2047     3.0000     1.6999 × 10^6     18.4336     16.9644 

		12    4095     3.0000     6.7929 × 10^6      19.2825     17.2715 

		13    8191     3.0000     2.7198 × 10^7      20.0209     17.5118 

		14    16383     3.0000     1.0879 × 10^8      20.6658     17.7025 

		
	tableCondition numbers of coefficient matrices using the finite element basis _N in (<ref>) and the wavelet basis _2,N in (<ref>) of Example <ref> with Dirichlet homogeneous boundary conditions. The mass and stiffness matrices have different normalizations. For the mass matrix, we normalize each element of  _N and _2,N such that it has norm equal to 1. For the stiffness matrix, we normalize each element of _N and _2,N such that its derivative has norm equal to 1.









 §.§ Helmholtz equation

Let N≥ 2 be the finest scale level for computation.
In this section, we shall modify the truncated system _2,N in (<ref>) of the wavelet basis _2 for L_2([0,1]) constructed in Example <ref> of Section 3 to solve the Helmholtz equation.

Consider the following 1D model problem


    -u”-k^2u=f   on   (0,1), 
        u(0)=0,    u'(1)-ik u(1)=0,


where f ∈ L_2([0,1]) and κ >0.
Since the generators ϕ,ϕ̃ for the biorthogonal wavelets (_2,_2) in Example <ref> has smoothness (ϕ)=(a)=1.5 and (ϕ̃)=(ã)≈ 0.440765, according to Theorem <ref>, _2^τ(ϕ;ψ) must be a Riesz basis in the Sobolev space H^τ() for all -0.440765<τ<1.5.
Hence, after renormalizing each element in _2,N in H^1
and modifying the boundary wavelets to satisfy (<ref>) (more details below),
then we can use the modified _2,N satisfying the boundary conditions in (<ref>) to solve Helmholtz equation in (<ref>). The modification is admittedly heuristic by nature, but the modified _2,N still appears to be a Riesz basis in H^1([0,1]).

To capture the highly oscillating waves in the solution of the Helmholtz equation,
we additionally supplement our wavelet basis on [0,1] with special waves. Suppose we have a non-overlapping partition of the unit interval; i.e., [0,1] = ∪_l=1^M [a_l,b_l], where M denotes the number of partitions. Then, these special waves take the following form

    s^+_l(x):= (e^ikx-(λ^+_1,l+λ^+_2,lx))|_[a_l,b_l],   
    s^-_l(x):= (e^-ikx-(λ^-_1,l+λ^-_2,lx))|_[a_l,b_l],
      
    l=1,…,M.

Furthermore, the parameters λ^+_1,l, λ^+_2,l, λ^-_1,l, λ^-_2,l are chosen so that

    s^+_l(a_l)=s^+_l(b_l)=s^-_l(a_l)=s^-_l(b_l)=0,    l=1,…,M-1, 
        s^+_M(a_M)=s^-_M(a_M)=0,    (s^+_M)'(1)-iks^+_M(1)=(s^-_M)'(1)-iks^-_M(1)=0.


Let us disregard the normalization constants for a moment, we note that our original right boundary elements are

    ℛ_N:= {ϕ(2^2·-3)}∪{ψ^R(2^j·-(2^j-1))  2 ≤ j ≤ N-1}.

Denote ψ^R,n to be the n-th component of ψ^R. To ensure all right boundary elements satisfy the radiation boundary condition, we introduce a heuristic modification to a subset of these elements. Define

    𝒮_N:={ϕ(2^2·-3)}∪{ψ^R,2(2^j·-(2^j-1))  2 ≤ j ≤ N-1},

where ψ^R=ψ^L(1-·) and ψ^L is defined in (<ref>). For each g ∈𝒮_N, let supp(g) be the support of g, supp(g):=[l_g,1], and define g^new(x):=g(x)-(λ^g_1+λ^g_2 x), where λ^g_1, λ^g_2∈. Instead of ℛ_N, we shall use the following modified right boundary elements

    ℛ_N^new:={ψ^R,1(2^j·-(2^j-1))  j ≥ 2}∪{g^new g^new(l_g)=0, (g^new)'(1)-ikg^new(1)=0, g ∈𝒮_N}.

To further reduce the condition number, we may first find a pair of linear combinations of the original right boundary wavelet functions at a given scale level, perform the modification, and find another pair of linear combinations of the modified right boundary wavelet functions at the same scale level. The same pairs of linear combinations are then applied to all scale levels.

Next, we normalize our modified wavelet basis on [0,1] and special waves so that each of their first derivatives has norm equal to 1. Let {g_n}_n=1^2^N+2M-1 be the enumerated basis elements starting from the wavelet basis (ordered first by scales then by shifts) and ending with the special waves. The Galerkin formulation for (<ref>) and (<ref>) is

    ∑_n=1^2^N+2M-1( ⟨ g'_l,g'_n⟩-k^2⟨ g_l,g_n⟩-g_l(1)g_n'(1)) c_n = ⟨ g_l,f ⟩,    l=1,…,2^N+2M-1.


Consider the linear system of equations Ax=b induced by the above Galerkin formulation with

    A=[ A_1 A_2; A_3 A_4 ],   
    x=[ x_1; x_2 ],   
    b=[ b_1; b_2 ],

such that A_1∈^n_1× n_1, A_2∈^n_1× n_2, A_3∈^n_2× n_1, A_4∈^n_2× n_2, x_1,b_1∈^n_1, and x_2,b_2∈^n_2. Then, (A_1-A_2A_4^-1A_3) x_1 = b_1 - A_2A_4^-1b_2 and A_4 x_2 = b_2-A_3x_1. For the two following examples, we shall use two different conditions numbers κ and κ^*. The former corresponds to the condition number of the entire coefficient matrix A; meanwhile, the latter corresponds to the condition number of A_1-A_2A_4^-1A_3. The submatrix A_4 coincides with the part of the coefficient matrix, where the special waves reside. The reason for having such a quantity is because generally the number of special waves is very small. Consequently, the matrix A_4 is often well-conditioned and can be inverted without any problem. Therefore, we can also treat the condition number of (A_1-A_2A_4^-1A_3) as a good indicator of how well-conditioned the entire coefficient matrix A is and we may use it to reduce the size of A with improved condition numbers.

To assess the performance of our enriched wavelet basis (the modified _2,N in (<ref>) plus the special waves), we compare it with the standard wavelet-based Galerkin method (only the modified _2,N in (<ref>)) and the (pollution free) finite difference method proposed in <cit.>. Additionally, we utilize two distinct error measurements for the finite difference method: discrete and interpolation errors. Let h be the grid size and H_FD=1/h. Define the discrete (relative) error as

    Dis. Err.:= . (∑_n=1^H_FD |u(n h)-U_H_FD(n)|^2)^1/2/ (∑_n=1^H_FD |u(n h)|^2)^1/2, .

where u is the exact solution and U_H_FD is the solution obtained from the finite difference method. The reason for introducing another measure of errors is because in a finite difference scheme, we do not know the values in between the discrete approximated solution. To partially reconcile the L_2 relative error of the wavelet-based Galerkin method and the discrete error of the finite difference method, we simply interpolate the discrete points obtained in the finite difference method. To calculate the interpolation error, we first interpolate the discrete points U_H_FD with the hat function

    u_N_FD^int:=∑_n=1^H_FD U_H_FD(n) B_2(2^N_FD·-n)     with     N_FD:=log_2(H_FD),

and then calculate the relative L_2 error

    e_N_FD^int_L_2/ u_L_2:= u_N_FD^int -u_L_2/ u_L_2         e_N_FD^int:=u_N_FD^int-u.


 
	Suppose f:=2√(2)× 10^8× (χ_[3/16,5/16] + χ_[11/16,13/16]) and k=20000 in (<ref>). The unit interval is partitioned into [0,3/16] ∪ [3/16,5/16] ∪ [5/16,11/16] ∪ [11/16,13/16] ∪ [13/16,1]. Table <ref> lists sizes of coefficient matrices, condition numbers, and errors. We see that the relative error corresponding to the wavelet-based Galerkin method with special waves hits 0% at the scale level N=4, which implies that the solution is recovered exactly. See Figure <ref> for plots of the source term as well as real and imaginary approximated solutions.






 
	Suppose f:=3000(2824/5χ_[0,1/3] + (-1743+5729/5x) χ_[1/3,3/5] + (-3525x^2+700x+949) χ_[3/5,5/7] + (37000x^3 - 9370211/175x^2+3803172/245x+2397)χ_[5/7,1]) and k=200000 in (<ref>). The unit interval is partitioned into [0,3/8] ∪ [3/8,5/8] ∪ [5/8,1]. Table <ref> lists sizes of coefficient matrices, condition numbers, and errors. See Figure <ref> for plots of the source term as well as the real and imaginary true solutions.






The condition numbers of coefficient matrices coming from the enriched wavelet basis with special waves are large. Our wavelet basis is definitely not the source of this problem, as its condition numbers are small (see κ in Tables <ref> and <ref> under `Wavelet only'), but rather the large condition numbers are due to the enrichment by special waves. This phenomenon has been reported in the literature of numerical methods for the Helmholtz equation, particularly those that use the homogeneous solution of the underlying DE in the trial and test spaces <cit.>.



 §.§ Biharmonic equation

In this subsection, we shall apply Example <ref> of Section 2 to solve the following biharmonic equation

    u^(4)=f   on   (0,1) 
        u(0)=u(1)=u'(0)=u'(1)=0,

where f ∈ L_2([0,1]). To justify our construction on [0,1], we observe {ϕ;ψ} is an orthogonal wavelet in L_2(), where ψ:=diag(√(3)/24,1/8) ψ” and ϕ:=(χ_[0,1], (√(3)(2x-1)χ_[0,1])^ such that {ϕ(·-k)}_k∈ and
{ϕ”(·-k)}_k∈ generate the same shift-invariant space. Using the construction in Section 3, we obtain an orthogonal wavelet in L_2([0,1]) derived from {ϕ;ψ}. As a consequence, we can deduce a Riesz wavelet in H^2([0,1]) derived from {ϕ;ψ}. Unlike in Section 3.1 with the Helmholtz equation, we choose our coarsest scale level to be J=1. That is, define

    Φ_1:={2^-1/2ϕ_1(2·-1),2^-1/2ϕ_2(2·-1)}, 
       Ψ_j:={2^-j/2ψ_1(2^j·-k),2^-j/2ψ_2(2^j·-k)  0≤ k ≤ 2^j-1},    j ≥ 1, 
       _1:= Φ_1∪{Ψ_j j ≥ 1}.

After an appropriate renormalization in H^2, the system _1 forms a Riesz basis in H^2([0,1]). As before, assuming that N is our finest scale level, we have the truncated system:

    _1,N:= Φ_1∪{Ψ_j: 1 ≤ j ≤ N-1},    N ≥ 1,

where Φ_1, Ψ_j are defined in (<ref>) and (<ref>), respectively. Let {g_n}_n=1^2^N+2-2 be the ordered normalized wavelet basis _1,N (first by scales then by shifts). See <cit.> for the explicit ordering and normalization. The Galerkin formulation for (<ref>) and (<ref>) is

    ∑_n=1^2^N+2-2⟨ g_l”,g_n”⟩ c_n = ⟨ g_l,f⟩,    l=1,…,2^N+2-2.


 
	Suppose f:=-6250000π^2((-3/625+(x^2-x)π^2)sin(50π x) -4/25πcos(50π x)(x-1/2)) in (<ref>). The true solution is u(x)=sin(50π x)(-x^2+x). Table <ref> lists sizes of coefficient matrices, condition numbers, and errors. As can be seen from the table, the condition number of the coefficient matrix is identically equal to 1. I.e., the coefficient matrix is an identity matrix. Therefore, we do not need to solve any linear systems.




We conclude this section with one remark. Notice that we have deliberately used a derivative-orthogonal Riesz wavelet only for the biharmonic equation. There certainly exists a first-order derivative-orthogonal Riesz wavelet generated from the hat function, whose stiffness matrix has a condition number identically equal to 1. See <cit.>. However, such a wavelet is a Riesz basis only in H^1(), but not in L_2(). On the contrary, the wavelet studied in Example <ref> of Section 3 is a Riesz basis in both L_2() and H^1(). To ensure that the condition numbers of mass and stiffness matrices are uniformly bounded, the wavelet studied in Example <ref> of Section 3, which is a Riesz basis in both H^1() and L_2() after renormalization,
is a more suitable choice for solving the Helmholtz equation with very large wave numbers.



§ CONCLUSION

In summary we have reviewed some important research directions in the context of wavelet-based methods for numerical DEs. Two main results revisited were the construction of m-th order derivative orthogonal Riesz wavelets in the Sobolev space H^m() as well as a general construction of wavelets in L_2([0,N]) satisfying maximum vanishing moments and given boundary conditions. Furthermore, we pointed out that adaptive wavelet algorithms have received a lot of interest over the years, even though we did not delve into this topic in this paper. Earlier on, we also discussed the advantages that wavelets bring in solving numerical DEs and demonstrated their ability in solving some model problems. One thing worth emphasizing is that we do not claim wavelets can efficiently solve all problems in numerical DEs. However, for selected problems, wavelets may serve as an effective approach.

Going forward, there are a few interesting problems to consider. Developing wavelets with simple structures that are capable of effectively tackling high dimensional problems remains as a critical research direction. We restrict the domain of interest to [0,1]^d, since the construction of wavelets on a general bounded domain poses a great challenge. Currently, the existence of derivative-orthogonal Riesz wavelet in ^d, where d ≥ 2, is unknown. If such wavelets existed, then we would encounter a desired situation like Example <ref> in a multidimensional setting. I.e., we do not need to solve any linear systems for basic Poisson and biharmonic problems. Otherwise, we need to explore the possibility of finding Riesz wavelets in ^d with near orthogonality in their derivatives. The Helmholtz equation serves as a good motivation to study the construction of wavelets on a bounded interval satisfying mixed (Robin) boundary conditions and maximum vanishing moments. Unfortunately, we do not have any theoretical justification for the boundary elements modification we adopted in Section 4.1. One difficulty that arises from mixed boundary conditions is that we now have to deal with nonstationary wavelet systems, since the boundary elements are no longer scale invariant. We are curious to see how we can extend our construction in Section 3 to accommodate mixed boundary conditions. Still on the topic of the Helmholtz equation, we earlier observe that wavelets combined with special waves yield encouraging results. However, the corresponding condition numbers are still large and they will continue to grow in high dimensions. Naturally, we want to know if we can design a better wavelet-based method that captures both low and high frequencies more efficiently.

99
	abcr93
	B. Alpert, G. Beylkin, R. Coifman, and V. Rokhlin, Wavelet-like bases for the fast solution of second-kind integral equations. SIAM J. Sci. Comput. 14 (1993), 159–184.
	
	ak12
	A. Altürk and F. Keinert, Regularity of boundary wavelets. Appl. Comput. Harmon. Anal. 32 (2012), no. 1, 65–85.
	
	ahjp94
	L. Andersson, N. Hall, B. Jawerth, and G. Peters, Wavelets on closed subsets of the real line. Recent advances in wavelet analysis, 1–61, Wavelet Anal. Appl., 3, Academic Press, Boston, MA, 1994.
	
ahl17
E. Ashpazzadeh, B. Han, and M. Lakestani, Biorthogonal multiwavelets on the interval for numerical solutions of Burgers' equation. J. Comput. Appl. Math. 317 (2017), 510–534.

	a93
	P. Auscher, Ondelletes á support compact et conditions aux limites. J. Funct. Anal. 111 (1993), 29-43.
	
	bcr91
	G. Beylkin, R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithm I. Comm. Pure Appl. Math. 44 (1991), 141–183.
	
	ctu99
	C. Canuto, A. Tabacco, and K. Urban, The wavelet element method. Appl. Comp. Harmon. Anal. 6 (1999), 1–52.
	
	cr00
	M. A. Christon and D. W. Roach, The numerical performance of wavelets for PDEs: the multi-scale finite element. Comput. Mech. 25 (2000), no. 2, 230–244.
	
	cohenbook
	A. Cohen, Wavelet methods in numerical analysis. Handbook of numerical analysis, Vol. VII, Handb. Numer. Anal., VII, North-Holland, Amsterdam, 2000.
	
	cdd00
	A. Cohen, W. Dahmen, and R. DeVore, Adaptive wavelet methods for elliptic operator equations: convergence rates. Math. Comp. 70 (2000), no. 233, 27–75.
	
	cdd02
	A. Cohen, w. Dahmen, and R. DeVore, Adaptive wavelet methods II–Beyond the elliptic case. Found. Comput. Math. 2 (2002), 203–245.
	
	cdf92
	A. Cohen, I. Daubechies, and J. C. Feauveau, Biorthogonal bases of compactly supported wavelets. Comm. Pure Appl. Math.
	45 (1992), 485–560.
	
	cdv93
	A. Cohen, I. Daubechies, and P. Vial, Wavelets on the interval and fast wavelet transforms. Appl. Comput. Harmon. Anal. 1 (1993), no. 1, 54–81.
	
	cf16
	D. Černá, and V. Finěk, On a sparse representation of an n-Dimensional Laplacian in wavelet coordinates. Results. Math. 69 (2016), no. 1, 225–243.
	
	cf17
	D. Černá and V. Finěk, Sparse wavelet representation of differential operators with piecewise polynomial coefficients. Axioms. 6 (2017), no. 4, 1–21.
	
	cf20
	D. Černá and V. Finěk, Galerkin method with new quadratic spline wavelet for integral and integro-differential equations, J. Comput. Appl. Math. 63 (2020), 426–443.
	
cq92
	C. K. Chui and E. Quak, Wavelets on a bounded interval. Numerical methods in approximation theory, Vol. 9, 53–75, Internat. Ser. Numer. Math., 105, Birkhäuser, Basel, 1992.

cw92	
C. K. Chui and J.-Z. Wang, On compactly supported spline wavelets and a duality principle, Trans. Amer. Math. Soc. 330 (1992) 903–915.

	d01
	W. Dahmen, Wavelet methods for PDEs–some recent developments. J. Comput. Appl. Math. 128 (2001), 133-185.
	
	dhjk00
	W. Dahmen, B. Han, R.-Q. Jia, and A. Kunoth, Biorthogonal multiwavelets on the interval: cubic Hermite splines. Constr. Approx. 16 (2000), no. 2, 221–259.
	
	dk92
	W. Dahmen and A. Kunoth, Multilevel preconditioning. Numer. Math. 63 (1992), 315–344.
	
	dku99
	W. Dahmen, A. Kunoth and K. Urban, Biorthogonal spline wavelets on the interval—stability and moment conditions. Appl. Comput. Harmon. Anal. 6 (1999), no. 2, 132–196.
	
	ds98
	W. Dahmen and R. Schneider, Wavelets with complementary boundary conditions–function spaces on the cube. Result. Math. 34 (1998), 255–293.
	
	ds99
	W. Dahmen and R. Schneider, Wavelets on manifolds I: construction and domain decomposition, SIAM J. Math. Anal. 31 (1999), no. 1, 184–230.
	
	d88
	I. Daubechies, Orthonormal bases of compactly supported wavelets. Comm. Pure Appl. Math. 41 (1988), no. 7, 909-996.
	
	ds10
	T. J. Dijkema and R. Stevenson. A sparse Laplacian in tensor product wavelet coordinates. Numer. Math. 115 (2010), no. 3, 433–449.
	
	dss09
	T. J. Dijkema, C. Schwab, and R. Stevenson. An adaptive wavelet method for solving high-dimensional elliptic PDEs. Constr. Approx. 30 (2009), 423–455.
	
	dms19
	G. C. Diwan, A. Moiola, E. A. Spence, Can coercive formulations lead to a fast and accurate solution of the Helmholtz equation? J. Comput. Appl. Math. 352 (2019), 110–131.
	
	gz19
	M. J. Gander and H. Zhang. A class of iterative solvers for the Helmholtz equation: factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz method. SIAM Rev. 61 (2019), no. 1, 3–76.
	
	ggm84
	P. Goupillaud, A. Grossman, and J. Morlet, Cycle-octave and related transforms in seismic signal analysis. Geoexploration 23 (1984), no. 1, 85–102.

han01
B. Han, Approximation properties and construction of Hermite interpolants and biorthogonal multiwavelets. J. Approx. Theory 110 (2001), no. 1, 18–53.

han03jat
B. Han, Vector cascade algorithms and refinable function vectors in Sobolev space, J. Approx. Theory. 124(2003) 44–88.

han10
B. Han, Pairs of frequency-based nonhomogeneous dual wavelet frames in the distribution space. Appl. Comput. Harmon. Anal. 29 (2010), 330–353.

han12
B. Han, Nonhomogeneous wavelet systems in high dimensions. Appl. Comput. Harmon. Anal. 32 (2012), 169–196.

hanbook
	B. Han, Framelets and wavelets: Algorithms, analysis, and applications. Applied and Numerical Harmonic Analysis. Birkhäuser/Springer, Cham, 2017. xxxiii + 724 pp.
	
	hm18
	B. Han and M. Michelle, Construction of wavelets and framelets on a bounded interval. Anal. Appl. 16 (2018), no. 6, 807-849.
	
	hm19
	B. Han and M. Michelle, Derivative-orthogonal Riesz wavelets in Sobolev spaces with applications to differential equations. Appl. Comp. Harmon. Anal. 47 (2019), no. 3, 759-794.
	
	hmw19
	B. Han, M. Michelle, and Y. S. Wong, Biorthogonal wavelets on intervals with vanishing moments and polynomial reproduction. Preprint. (2019).
	
	hs09
	B. Han and Z. Shen, Dual wavelet frames and Riesz bases in Sobolev spaces. Constr. Approx. 29 (2009), no. 3, 369–406.
	
hm99
D. P. Hardin and S. A. Marasovich, Biorthogonal multiwavelets on [-1,1].
Appl. Comput. Harmon. Anal. 7 (1999), no. 1, 34–53.

	hmp16
	R. Hiptmair, A. Moiola, and I. Perugia, A survey of Trefftz methods for the Helmholtz equation. Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations. Lecture Notes in Computational Science and Engineering. Springer, Cham, 2016.
	
	h99
	M. Holmström, Solving hyperbolic PDEs using interpolating wavelets. SIAM J. Sci. Comput. 21 (1999), no. 2, 405–420.
	
	j09
	R.-Q. Jia, Spline wavelets on the interval with homogeneous boundary conditions. Adv. Comput. Math. 30 (2009), no. 2, 177–200.
	
	jl06
	R.-Q. Jia and S.-T. Liu, Wavelet bases of Hermite cubic splines on the interval. Adv. Comput. Math. 25 (2006), no. 1-3, 23–39.
	
	jz11
	R.-Q. Jia and W. Zhao, Riesz bases of wavelets and applications to numerical solutions of elliptic equations. Math. Comp. 80 (2011), no. 275, 1525–1556.
	
	ks06
	A. Kunoth and J. Sahner, Wavelets on manifolds: an optimized construction. Math. Comp 75 (2006), no. 255, 1319–1349.
	
	m89
	S. G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans. Patter Anal. Mach. Intell. (1989), no. 7, 674–693.
	
	m90
	Y. Meyer, Ondelettes et opérateurs. I, II, and III. Actualités Mathématiques. Hermann, Paris, (1980).
	
	m91
	Y. Meyer, Ondelettes sur l'intervalle. Rev. Mat. Iberoamericana 7 (1991), no. 2, 115–133.
	
	ms14
	A. Moiola and E. A. Spence, Is the Helmholtz equation really sign-indefinite? SIAM Rev. 56 (2014), no. 2, 274–312.
	
	pst95
	G. Plonka, K. Selig, and M. Tasche, On the construction of wavelets on a bounded interval. Adv. Comput. Math. 4 (1995), no. 4, 357–388.

rs92	
S.D. Riemenschneider and Z. Shen, Wavelets and pre-wavelets in low dimensions, J. Approx. Theory 71 (1992), 18–38.

	ss09
	C. Schwab and R. Stevenson, Space-time adaptive wavelet methods for parabolic evolution problems. Math. Comp. 78 (2009), no. 267, 1293–1318.
	
	ww14
	K. Wang and Y. S. Wong, Pollution-free finite difference schemes for non-homogeneous Helmholtz equation. Int. J. Numer. Anal. Mod. 11 (2014), no. 4, 787–815.


	



dgh96
G. Donovan, J. Geronimo, and D. Hardin, Intertwining multiresolution analyses and the construction of piecewise-polynomial wavelets. SIAM J. Math. Anal. 27 (1996), no. 6, 1791–1815.

dau88
I. Daubechies, Orthonormal bases of compactly supported wavelets. Comm. Pure Appl. Math. 41 (1988), no. 7, 909–996.

ds10 T. J. Dijkema and R. Stevenson. A sparse Laplacian in tensor product wavelet corrdinates. Numer. Math. 115 (2010), no. 3, 433–449.

ghm94
J. S. Geronimo, D. P. Hardin, and P. R. Massopust, Fractal functions and wavelet expansions based on several scaling functions. J. Approx. Theory 78 (1994), no. 3, 373–401.


hm99
D. P. Hardin and S. A. Marasovich, Biorthogonal multiwavelets on [-1,1].
Appl. Comput. Harmon. Anal. 7 (1999), no. 1, 34–53.

jj03
R.-Q. Jia and Q. T. Jiang, Spectral analysis of the transition operator and its applications to smoothness analysis of wavelets. SIAM J. Matrix Anal. Appl. 24 (2003), no. 4, 1071–1109.

jrz99
R. Q. Jia, S. D. Riemenschneider, D.-X. Zhou, Smoothness of multiple refinable functions and multiple wavelets. SIAM J. Matrix Anal. Appl. 21 (1999), 1–28.

mad97
W. R. Madych, Finite orthogonal transforms and multiresolution analyses on intervals. J. Fourier Anal. Appl. 3 (1997), no. 3, 257–294.

pst95
G. Plonka, K. Selig, and M. Tasche, On the construction of wavelets on a bounded interval. Adv. Comput. Math. 4 (1995), no. 4, 357–388.

pz03
G. Plonka and D.-X. Zhou, Properties of locally linearly independent refinable function vectors, J. Approx. Theory, 122 (2003), 24–41.



