Galerkin Scheme Using Biorthogonal Wavelets on Intervals for 2D Elliptic Interface Problems

We introduce a new second-order Galerkin scheme using the tensor product of biorthogonal wavelets on intervals for the model problem in (1.1d). As mentioned earlier, to achieve optimal convergence rates (i.e., those consistent with the approximation order of the scheme), special treatments are required to handle the interface. Our method involves adding extra wavelet elements, which touch the interface and belong to higher scale levels, to our approximated solution. This simple approach conveniently handles the complicated geometry of the interface, effectively captures the singularity along the interface, and easily handles high-contract coefficients $a$, thereby enabling us to achieve near optimal convergence rates $\mathscr{O}(h^{2}|\log(h)|^{2})$ in the $L_{2}(\Omega)$-norm and $\mathscr{O}(h|\log(h)|)$ in the $H^{1}(\Omega)$-norm, as stated in Theorem 2.2. More specifically, the convergence rates of our method get arbitrarily close to second-order in the $L_{2}(\Omega)$-norm and first-order in the $H^{1}(\Omega)$-norm as the scale level increases. We establish these near optimal convergence rates by extensively using the dual part of the biorthogonal wavelet basis, relying on the weighted Bessel property and results of wavelets in fractional Sobolev spaces, and employing standard FEM arguments. Not only does our method easily handle the interface geometry, but it can also manage high-contrast coefficients effectively as demonstrated by our numerical experiments later.

It is also worth noting that the added/augmented wavelet elements have scale levels that are at most double the maximum scale level of the other regular basis functions positioned throughout the domain. Consequently, the number of terms used in the approximated solution with these extra functions is only a fixed constant multiple of the number of terms without them. This fixed constant depends on the shape of the interface curve and the support of the wavelet elements. The new linear system corresponding to the coefficients is also a fixed constant multiple of the previous one.

Our method, in a sense, can be interpreted as a meshfree method in that we do not need to generate a mesh that depends on the domain and the interface. To obtain a more accurate solution, there is no need for re-meshing of the entire domain with a smaller mesh size. Instead, to increase accuracy, we raise the scale level of the approximated solution, which entails adding more wavelet elements throughout the domain and additional wavelet elements near the interface. Furthermore, same as the XFEM, our method is considered to be conforming, since each wavelet element is continuous.

Coefficient matrices of many FEMs are known to have condition numbers that are growing proportionally to $h^{-2}$, where $h$ is the mesh size (e.g., [1, 6, 28]). Our wavelet Galerkin method produces coefficient matrices whose condition numbers are relatively small and uniformly bounded regardless of its size. More precisely, we prove in (2.19) of Theorem 2.2 that the condition numbers are bounded by $C_{w}\|a\|_{L_{\infty}(\Omega)}\|a^{-1}\|_{L_{\infty}(\Omega)}$, where the constant $C_{w}$ only depends on the wavelet basis and the domain $\Omega$, but $C_{w}$ is independent of the interface $\Gamma$. Additionally, the smallest singular values of the coefficient matrices are uniformly bounded away from zero. This is an advantage that we inherit directly from the fact that our 2D wavelet basis is a Riesz basis of $H^{1}_{0}(\Omega)$. Having such nice condition number properties is beneficial, when an iterative solver is employed to solve the linear system, as it often leads to a relatively small number of iterations required to reach a given tolerance level.

At present, we solely aim to lay the groundwork of our wavelet Galerkin method for solving the 2D elliptic interface problem in (1.1d). Other equally important problems like its high-order version and extensions to the 3D setting are left as a future work, since their implementations and effective calculation of quadratures are much more demanding.

In Section 2.1, we revisit some basic concepts and definitions of wavelets. In Section 2.2, we present the biorthogonal wavelet basis derived from the bilinear function, which we shall use throughout the paper. In Section 2, we formally present our wavelet Galerkin method for the model problem (1.1d) and state our main result in Theorem 2.2 on convergence rates and uniform boundedness of condition numbers. In Section 3, we discuss how we handle nonhomogeneous first jump conditions and/or Dirichlet boundary conditions, and present some numerical experiments to demonstrate the performance of our proposed method. Finally, we present the proofs of our theoretical findings in Theorem 2.2 on convergence rates in Section 4.

Let us first review some basic concepts of wavelets, which follow a similar presentation as in [28]. Let $\phi:=\{\phi^{1},\dots,\phi^{r}\}^{\mathsf{T}}$ and $\psi:=\{\psi^{1},\dots,\psi^{s}\}^{\mathsf{T}}$ be square integrable functions in $L_{2}(\mathbb{R})$. Define a wavelet affine system by

where $J_{0}\in\mathbb{Z}$, $\phi^{\ell}_{J_{0};k}:=2^{J_{0}/2}\phi^{\ell}(2^{J_{0}}\cdot-k)$, and $\psi^{\ell}_{j;k}:=2^{j/2}\psi^{\ell}(2^{j}\cdot-k)$. We say that $\operatorname{\mathsf{AS}}_{J_{0}}(\phi;\psi)$ is a Riesz basis for $L_{2}(\mathbb{R})$ if (1) the linear span of $\operatorname{\mathsf{AS}}_{J_{0}}(\phi;\psi)$ is dense in $L_{2}(\mathbb{R})$, and (2) there exist $C_{1},C_{2}>0$ such that

for all finitely supported sequences $\{c_{\eta}\}_{\eta\in\operatorname{\mathsf{AS}}_{J_{0}}(\phi;\psi)}$. The relation in (2.2) holds for some $J_{0}\in\mathbb{Z}$ if and only if it holds for all $J_{0}\in\mathbb{Z}$ with identical positive constants $C_{1}$ and $C_{2}$ (see for example [23, Theorem 6]). As a result, we simply refer to $\{\phi;\psi\}$ as a Riesz multiwavelet in $L_{2}(\mathbb{R})$ if $\operatorname{\mathsf{AS}}_{0}(\phi;\psi)$ is a Riesz basis for $L_{2}(\mathbb{R})$. Further, let $\tilde{\phi}:=\{\tilde{\phi}^{1},\dots,\tilde{\phi}^{r}\}^{\mathsf{T}}$ and $\tilde{\psi}:=\{\tilde{\psi}^{1},\dots,\tilde{\psi}^{s}\}^{\mathsf{T}}$ be vectors of square integrable functions in $L_{2}(\mathbb{R})$. We say that $(\{\tilde{\phi};\tilde{\psi}\},\{\phi;\psi\})$ is a biorthogonal multiwavelet in $L_{2}(\mathbb{R})$ if $(\operatorname{\mathsf{AS}}_{0}(\tilde{\phi};\tilde{\psi}),\operatorname{\mathsf{AS}}_{0}(\phi;\psi))$ is a biorthogonal basis in $L_{2}(\mathbb{R})$, i.e., (1) $\operatorname{\mathsf{AS}}_{0}(\tilde{\phi};\tilde{\psi})$ and $\operatorname{\mathsf{AS}}_{0}(\phi;\psi)$ are Riesz bases in $L_{2}(\mathbb{R})$, and (2) $\operatorname{\mathsf{AS}}_{0}(\tilde{\phi};\tilde{\psi})$ and $\operatorname{\mathsf{AS}}_{0}(\phi;\psi)$ are biorthogonal to each other in $L_{2}(\mathbb{R})$.

The wavelet function $\psi$ has $m$ vanishing moments if $\int_{\mathbb{R}}x^{j}\psi(x)dx=0$ for all $j=0,\ldots,m-1$. By convention, we define $\operatorname{vm}(\psi):=m$ with $m$ being the largest of such an integer.

The Fourier transform is defined by $\widehat{f}(\xi):=\int_{\mathbb{R}}f(x)e^{-\textsf{i}x\xi}dx,\xi\in\mathbb{R}$ for $f\in L_{1}(\mathbb{R})$ and is naturally extended to square integrable functions in $L_{2}(\mathbb{R})$. Meanwhile, the Fourier series of $u=\{u(k)\}_{k\in\mathbb{Z}}\in(l_{0}(\mathbb{Z}))^{r\times s}$ is defined by $\widehat{u}(\xi):=\sum_{k\in\mathbb{Z}}u(k)e^{-\textsf{i}k\xi}$ for $\xi\in\mathbb{R}$, which is an $r\times s$ matrix of $2\pi$-periodic trigonometric polynomials. By $\bm{\delta}$ we denote the sequence such that $\bm{\delta}(0)=1$ and $\bm{\delta}(k)=0$ if $k\neq 0$.

Now, we are ready to recall a result of biorthogonal wavelets in $L_{2}(\mathbb{R})$.

To solve the elliptic interface problem (1.1d), we take the tensor product of wavelets on $\mathcal{I}:=(0,1)$. Without explicitly involving the dual, the direct approach presented in [27] allows us to construct all possible locally compactly supported biorthogonal wavelets in $L_{2}(\mathcal{I}$) satisfying prescribed boundary conditions and vanishing moments from any compactly supported biorthogonal (multi)wavelets in $L_{2}(\mathbb{R})$. That is, our direct approach produces a biorthogonal wavelet $(\tilde{\mathcal{B}}^{1D}_{J_{0}},\mathcal{B}^{1D}_{J_{0}})$ in $L_{2}(\mathcal{I}$), where

the integer $J_{0}\in\mathbb{N}$ denotes the coarsest scale level, and

with $n_{l,\phi},n_{h,\phi},n_{l,\psi},n_{h,\psi}$ being known integers, $\phi^{L},\phi^{R}$ being boundary refinable functions, and $\psi^{L},\psi^{R}$ being boundary wavelets that are finite subsets of functions in $L_{2}(\mathcal{I})$. Recall that $\psi_{j;k}:=2^{j/2}\psi(2^{j}\cdot-k)$. We define $\tilde{\mathcal{B}}^{1D}_{J_{0}}$ the same way, except we add $\sim$ to each element in $\mathcal{B}^{1D}_{J_{0}}$ for a natural bijection.

Throughout the paper, for simplicity of presentation, we consider the domain $\Omega=(0,1)^{2}$. Though many biorthogonal wavelet bases in $H^{1}_{0}(\Omega)$ can be used for numerically solving the elliptic interface problems (e.g., [27, Section 7], [28, Section 3.2] and [8, 15]), we shall restrict our attention to one specific biorthogonal wavelet basis on the bounded interval $\mathcal{I}:=(0,1)$.

Interpolating functions play a critical role in numerical PDEs, wavelet analysis, and computer aided geometric design (e.g., see [25] and references therein). The simplest example of compactly supported interpolating functions is probably the hat function $\phi(x):=\max(1-|x|,0)$ for $x\in\mathbb{R}$, which is extensively used in numerical PDEs and approximation theory. The hat function $\phi$ satisfies the refinement equation $\phi=\frac{1}{2}\phi(2\cdot-1)+\phi(2\cdot)+\frac{1}{2}\phi(2\cdot+1)$ and $\widehat{\phi}(0)=1$. In what follows, we recall a biorthogonal wavelet basis in $L_{2}(\mathcal{I})$ derived from the hat function $\phi$ and discussed in [27, Example 7.5], which will be the only biorthogonal wavelet basis used in this paper.

Let $\phi$ be the hat function. Consider the scalar biorthogonal wavelet $(\{\tilde{\phi};\tilde{\psi}\},\{\phi;\psi\})$ in $L_{2}(\mathbb{R})$ with $\widehat{\phi}(0)=\widehat{\tilde{\phi}}(0)=1$ and a biorthogonal wavelet filter bank $(\{\tilde{a};\tilde{b}\},\{a;b\})$ given by

In other words, the refinable functions $\phi,\tilde{\phi}$ and the wavelet functions $\psi,\tilde{\psi}$ are determined through the equations in (2.3) and (2.4). Note that the analytic expression of the hat function is $\phi:=(x+1)\chi_{[-1,0)}+(1-x)\chi_{[0,1]}$. As discussed in [27, Example 7.5], the boundary refinable functions and boundary wavelet functions are defined to be

Furthermore, we define

$\mathcal{B}^{1D}_{J_{0}}:=\Phi_{J_{0}}\cup\cup_{j=J_{0}}^{\infty}\Psi_{j}$, and $\tilde{\mathcal{B}}^{1D}_{J_{0}}:=\tilde{\Phi}_{J_{0}}\cup\cup_{j=J_{0}}^{\infty}\tilde{\Psi}_{j}$. Then, $(\tilde{\mathcal{B}}^{1D}_{J_{0}},\mathcal{B}^{1D}_{J_{0}})$, where $J_{0}\geqslant 3$, is a biorthogonal wavelet in $L_{2}(\mathcal{I})$. We shall use the tensor product of this one-dimensional biorthogonal wavelet in $L_{2}(\mathcal{I})$ throughout this paper. Due to item (3) of Theorem 2.1 and the relations stated in (2.7), there exist well-defined (refinability) matrices $A_{j,j^{\prime}}$ and $B_{j,j^{\prime}}$ such that

which may be convenient to use in the numerical implementation.

We now discuss how to obtain two-dimensional biorthogonal wavelets in $L_{2}(\Omega)$ with $\Omega=(0,1)^{2}$ using the tensor product of the one-dimensional biorthogonal wavelet in $L_{2}(\mathcal{I})$. Given 1D functions $f_{1}$ and $f_{2}$, define $(f_{1}\otimes f_{2})(x,y):=f_{1}(x)f_{2}(y)$ for $x,y\in\mathbb{R}$. If $F_{1},F_{2}$ are sets containing 1D functions, we define $F_{1}\otimes F_{2}:=\{f_{1}\otimes f_{2}:f_{1}\in F_{1},f_{2}\in F_{2}\}$. Also, define

where

where $\Phi_{j},\Psi_{j},\tilde{\Phi}_{j},\tilde{\Psi}_{j}$ are defined as in (2.8), and $J_{0}\geqslant 3$. By using an argument identical to [28, Theorem 1.2], we conclude that $(\tilde{\mathcal{B}}^{2D}_{J_{0}},\mathcal{B}^{2D}_{J_{0}})$ is a biorthogonal wavelet in $L_{2}(\Omega)$ and its properly scaled version defined below

is a Riesz basis of the Sobolev space $H^{1}_{0}(\Omega)$, see [28, Theorem 1.2]. That is, (1) the linear span of $\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0}}$ is dense in $H^{1}_{0}(\Omega)$, and (2) there exist positive constants $C_{\mathcal{B},1},C_{\mathcal{B},2}>0$ such that

for all finitely supported sequences $\{c_{\eta}\}_{\eta\in\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0}}}$.

We now describe our proposed method. As mentioned earlier, we shall always use the biorthogonal wavelet basis presented in Section 2.2.

For $J_{0}=3$ and $J\geqslant J_{0}$, we define the traditional finite-dimensional wavelet element space truncated at the scale level $J$ as follows:

Obviously, $\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0},J}$ is a finite subset of $\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0}}$, where $\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0}}$ is defined as in (2.11). Using a uniform grid, the standard FEM only uses the basis $\Phi^{2D}_{J}$ and its finite element space $V_{J}:=\mbox{span}(\Phi^{2D}_{J})$. It is very important to notice that $\mbox{span}(\mathcal{B}^{2D}_{J_{0},J})=V_{J}$ and $\mbox{span}(\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0},J})=V_{J}$. In other words, both $\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0},J}$ (or $\mathcal{B}^{2D}_{J_{0},J}$) and $\Phi^{2D}_{J}$ span the same (finite element) space $V_{J}$. The numerical solution to (1.1d) obtained by the traditional wavelet Galerkin method using only $\mathcal{B}^{2D}_{J_{0},J}$ is the same as the solution obtained by using $\Phi^{2D}_{J}$ (the standard bilinear FEM). In the context of the model problem (1.1d), using only the traditional wavelet basis $\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0},J}$ inevitably suffers from the same convergence issue faced in the standard FEM. More specifically, the observed convergence rate will typically be well below two because the exact solution $u\not\in H^{2}(\Omega)$ and has discontinuous gradients along the interface $\Gamma$.

To overcome such an issue, we propose incorporating of higher-resolution wavelets defined below

to capture the geometry of the interface $\Gamma$ and the singularity of the solution $u$ along the interface $\Gamma$ on top of the standard wavelet elements $\mathcal{B}^{2D}_{J_{0},J}$. More specifically, we shall use

where the superscript $S$ indicates that we add extra wavelet elements $S_{j},j=J,\ldots,2J-2$ to the traditional wavelet basis $\mathcal{B}^{2D}_{J_{0},J}$. Note that the cardinality of $\mathcal{B}^{2D}_{J_{0},J}$ is $\mathscr{O}(h^{-2})$, where $h:=2^{-J}$ is the mesh size. Because $\Gamma$ is a 1D closed curve inside $\Omega$, it is not difficult to observe that the cardinality of $\cup_{j=J}^{2J-2}\mathcal{S}_{j}$ is also $\mathscr{O}(h^{-2})$ (with the prefactor being independent of $J$). To put differently, the cardinality of $\mathcal{B}^{S}_{J_{0},J}$ is still $\mathscr{O}(h^{-2})$, which means that it is comparable to that of the original bases $\mathcal{B}^{2D}_{J_{0},J}$ and $\Phi^{2D}_{J}$.

By the weak formulation in (1.3) with $g=0$ and $g_{b}=0$ and considering the approximated function $u_{J}:=\sum_{\eta\in\mathcal{B}^{S,H^{1}_{0}(\Omega)}_{J_{0},J}}c_{\eta}\eta$ with to-be-determined unknown coefficients $\{c_{\eta}\}_{\eta\in\mathcal{B}^{S,H^{1}_{0}(\Omega)}_{J_{0},J}}$ (we shall also define $u_{h}:=u_{J}$ with $h:=2^{-J}$), our wavelet Galerkin method reduces to finding all the coefficients $c_{\eta}$ for $\eta\in\mathcal{B}^{H^{1}(\Omega)}_{J_{0},J}$ such that

Fig. 1 visualizes the basis functions used in our approximated solution $u_{J}$. Due to the symmetry in the biorthogonal wavelet basis in Section 2.2, each term of the approximated solution can be obtained by scaling, shifting, and rotating one of the functions in panels (c)-(h) of Fig. 1.

Meanwhile, Fig. 2 visualizes the overlapping supports of wavelet basis functions in $\mathcal{B}^{S,H^{1}_{0}(\Omega)}_{3,4}\cup\cup_{j=4}^{6}[2^{-j}\mathcal{S}_{j}]$, and gives us an insight as to how we add the wavelets along the interface in our approximated solution $u_{J}$.

Without loss of generality, we assume that $\Omega:=(0,1)^{2}$ is our domain. The next theorem is our main theoretical result on the convergence order of our proposed wavelet Galerkin method using $\mathcal{B}^{S,H^{1}_{0}(\Omega)}_{J_{0},J}$ in (2.14) as $J\to\infty$, and the uniform boundedness of the condition numbers of its coefficient matrices. We shall assume that $g=0$ on $\Gamma$ in the first jump condition (1.1b) for avoiding discontinuous $u$, and $g_{b}=0$ in $\partial\Omega$ for the homogeneous Dirichlet condition in accordance to the standard FEM argument. Since $\mathcal{B}^{S,H^{1}_{0}(\Omega)}_{J_{0},J}$ is a finite subset of the Riesz wavelet basis $\mathcal{B}^{H^{1}_{0}(\Omega)}_{J_{0}}$ in $H^{1}_{0}(\Omega)$, the condition numbers of coefficient matrices from (2.15) are uniformly bounded and independent of the mesh size $h$ and resolution level $J$. Due to the technicality, we defer the proof and its auxiliary results to Section 4.