The Wavelet-Capacitance Matrix Method

We present a numerical method for the solution of partial differential equations in nonseparable domains. This method was first developed by Qian and Weiss [20,21] for the Helmholtz equation in nonseparable geometry. The method uses a periodic wavelet-Galerkin solver with a specific adaptation of the standard capacitance matrix method.

We will describe the method 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 [18]. A method based on the double layer potential is also a possibility [19]. The algorithm is based on the calculation of a numerical partial Green's function [18]. The capacitance matrix method is developed in the references [1,2,11,15,18,19]. In particular, Martin develops the use of offset sources to fit airfoil geometry [15]. We found this procedure is to very useful for reducing the numerical residual error and accelerating convergence.

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

One advantage of this method is the use of fast and accurate periodic solvers to evaluate the solution. Another advantage is the efficient inclusion of general non-separable geometries. A possible disadvantage is the use of functions that are singularly supported in $S$. This can lower the accuracy of the numerical solution through Gibbs' phenomena and boundary residual errors.

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 will present an extensive series of numerical calculations that support our conclusions about accuracy and convergence.

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

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

where $\phi$ 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\phi (x_0-i)\phi (x-i) \phi (y_0-j)\phi (y-j).$$

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

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

which is the definition of the delta function in this subspace. As before, we assume, without loss of generality, that the integers are the finest scale of variation in these expressions. The scale factors are removed by a transform $x'=2^jx$, $y'=2^jy$.

Therefore, we solve, by the wavelet-Galerkin method [20,21,26], 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 [19].

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 will present the results in this report. 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.

For a specified geometry the capacitance matrix can be inverted once and for all and when $\hat \rho$ is extended from the discrete support to the periodic domain (equal to zero at non-support points), the solution of the Helmholtz equation with boundary data $g$ is found by one fast, periodic wavelet-Galerkin solution. As described in [19,18] the general inhomogeneous case reduces to this homogeneous problem.

The primary advantage of the direct method in comparison to the iterative methods is that the (discrete) boundary conditions are satisfied identically, and the method works near the resonance cases encountered in the solution of the Helmholtz equation. For domains that require many points for their discretization, the direct method could become impractical. We suggest that in these cases a fairly sparse selection of source points $\hat y$ in the direct method can effectively initialize an iterative (Conjugate-Gradient) method (even near resonance).