The wavelet-capacitance matrix method

To solve general nonseparable boundary value problems using the Wavelet-Galerkin method we have developed a Capacitance Matrix method. We will describe the method, as developed by Qian and Weiss [47,48], for the Harmonic Helmholtz equation

$$\left(-\Delta +\alpha\right)U=F$$

in a domain $D$ with boundary conditions $U=g$ on the boundary of $D$. One version of the direct method is equivalent to a numerical implementation the single layer potential [41]. A method based on the double layer potential is also a possibility [46]. The algorithm is based on the calculation of a numerical partial Green's Function [41].

The outline of our method is as follows. Regard the domain $D$ as contained (embedded) in a periodic cell, $S$. We extend $F$ from $D$ to $S$ in a smooth way. The extension $\hat F$ is periodic on $S$. We also define a periodic function $\hat \rho$ where $\hat \rho$ is zero except on the support of $\partial D\in S$. We determine $\hat \rho$ so that the periodic solution in $S$

$$\left(-\Delta + \alpha\right)U=\hat F + \hat \rho$$

will verify the boundary conditions $U=g$ on $\partial D$. By construction the equation $\left(-\Delta +\alpha\right)U=F$ is satisfied in $D$.

We have extended the method by allowing the support of $\hat \rho$ to be separate from the boundary of $D$, $\partial D$. When the equations are discretized by the Wavelet-Galerkin method, this extension eliminates the boundary residuals and defines a spectrally accurate method for non-separable domains. To our knowledge this algorithm is the first implementation of its' type. We have presented an extensive series of numerical calculations that support our conclusions about accuracy and convergence [47,48].

The numerical implementation is straight forward. In effect, we expand the solution in periodic, wavelet-Galerkin basis

$$U=\sum\sum U_{i,j}\varphi(x-i)\varphi(y-j)$$

where $\varphi$ is a scaling function. To calculate the Green's Function we resolve the delta function in the space of translates of scaling function

$$\lambda_{x_0,y_0}(x,y)=\sum\sum\varphi(x_0-i)\varphi(x-i) \varphi(y_0-j)\varphi(y-j).$$

Since the translates of the scaling function are orthogonal and complete in $L^2$, the above expression implies that for a square integrable function $f$

$$f(x_0,y_0)=\int\int dxdy \lambda_{x_0,y_0}(x,y) f(x,y),$$

which is the definition of the delta function.

Therefore, we solve, by the wavelet-Galerkin method [47,48], the equation

$$\left(-\Delta + \alpha\right)G(x,x_0;y,y_0)=\lambda_{x_0,y_0}(x,y)$$

for the partial Green's Function, $G$. To find the usual Capacitance Matrix, $C$, we discretize the boundary into a series of points $\hat x_j$ and form the matrix whose $(i,j)$ component is $G(\hat x_i,\hat x_j)$. The evaluation of $G$ requires only one solution of the periodic, fast, wavelet-Galerkin solver [46].

In our formulation of the algorithm, we discretize the boundary by the points $\hat x_j$ and the support of $\hat \rho$ in $S$ by the points $\hat y_j$. The definition of the capacitance matrix is then

$$C_{i,j}=G(\hat x_i,\hat y_j).$$

Depending on the cardinality of the sets $\hat x$ and $\hat y$, the system of equations for the discrete potential $\hat \rho$ are determined, overdetermined or underdetermined. We have examined these possibilities and present the results in ref. [47,48]. In general, if $\hat y$ is exterior to $\hat x$, we obtain excellent numerical results that depend stably on the choice of $\hat y$.

In terms of the (extended) Capacitance Matrix, the discrete potential of a single layer is a solution of the system

$$\hat g=C\hat\rho.$$

For non-determined systems we use a singular value decomposition of $C$ to find the least square or minimal norm solution [25].

The Capacitance Matrix is a fast and general method for solving boundary value problems in nonseparable domains. It uses fast periodic solvers based on the FFT to drive direct or iterative (Conjugate Gradient) algorithms. The geometry at the boundary is enforced by potentials with singular support on the boundary. The use of functions with singular support effectively restricts the Capacitance Matrix method to low order solvers, requiring a high level of discretization to produce accurate results. Due to boundary residuals, the introduction of higher order solvers can cause the rate of convergence to become worse. For problems with complicated geometries this fact limits the applicability of the method.

By combining a reformulation of the Capacitance Matrix method with a wavelet discretization, we have defined a Wavelet-Capacitance Matrix method. This allows the use of higher order approximations with rapid (even spectral) convergence and produces highly accurate solutions for low to moderate levels of discretization. In effect, we cure the Capacitance Matrix method of its' most serious limitation, while retaining the method's advantages.

The method applies equally to equations with three space dimensions and problems with a time dependence. For instance, we have already applied the method to the long time integration of Euler flow, with excellent results [60].

A preliminary comparison of wavelet methods to finite difference and spectral methods is reported in [31,47,58]. In general, the wavelet methods were found to be more stable and more accurate. For instance, Euler flow could not be solved as described above using Spectral methods (became unstable) [58] and the Wavelet-Capacitance method is more accurate than Capacitance methods with comparable, higher-order finite difference methods [47].

Figure 2 shows two frames from the time evolution of Navier-Stokes flow in a L-Shape region. The noslip boundary conditions are imposed by the Capacitance Matrix method. The Reynolds number is $21,000$.

Figure 2: Navier-Stokes flow in a L-Shape region at Reynolds number $21,000$