Numerical Methods

4.1 Quadrature

We begin this chapter with a brief review of quadrature. Recall a quadrature formula is one such that

where \(w\) is a weight function and \(g\) is a continuous function defined on \(G\). We assume that \(G\) is compact and Jordan measurable. A quadrature rule is then a discrete approximation to this continuous formula.

Here, \(x_j\) are the quadrature nodes where we evaluate the function and \(\alpha_j = w(x_j)\) are the weights. For convenience, the order of a quadrature rule \(n\) is the number of quadrature nodes we use in the approximation.

That is a sequence of quadrature rules is convergent if the linear functionals \(Q_n\) converge pointwise to the formula \(Q\). This review will be useful as we consider the numerical solutions to integrals and the Nyström method below.

4.2 Degenerate Kernel Approximations

If we recall our treatment of integral operators with weakly singular kernels in Theorem 1.1.21, we had a kernel we wanted to work with but instead approximated it be a series of more well-behaved kernels. This is the same idea for degenerate kernel approximations.

Let \(\langle X, X \rangle\) be a dual system generated by a non-degenerate bilinear form and \(A_n: X \to X\) be a bounded linear operator with

where \(\{ a_i \}_{i=1}^{n}\), \(\{ b_i \}_{i=1}^{n} \subset X\) and the \(a_i\)'s are linearly independent.

We call \(A_n\) a degenerate operator and clearly its action can be expressed a finite linear combination of basis elements for a space \(X\). In this way, solutions to the second kind equation

can be found by solving a finite dimensional linear system.

The justification for this approach is given in the following two theorems.

This result tells us that when the integral equation with the convergent operator is uniquely solvable, then the integral equations with the degenerate operators are uniquely solvable.

Analogously to the previous theorem, this theorem tells us that if the operators in the sequence of degenerate operators yield unique solutions to \((I - A_n)\varphi_n = f\) then the convergent operator \(A\) has a unique solution to \((I-A)\varphi = f\).

We can therefore use degenerate approximations of the kernels we have seen for second kind integral equations resulting from the boundary integral formulations of the Laplace and Helmholtz equations. To do so with these integral operators, we can approximate the continuous and weakly singular kernels by a sequence of continuous kernels converging to to the kernel in question.

Concretely, if the dual system we are looking at is \(\langle C(G), C(G) \rangle\) then the finite dimensional operator \(A_n: C(G) \to C(G)\) has the form

where this integral operator has a degenerate kernel given by

In this scenario, as we see from Theorem 4.2.3, the integral equation of the second time becomes

and the solution follows from the same theorem. Thus, we can frame the solution \((I-A)\varphi = f\) as the solution of a linear system. We just need to find these degenerate kernels.

4.2.1 Finding Degenerate Kernels

One way we can find these degenerate kernels is through interpolation. Much like the interpolative quadrature methods of classic numerical analysis, we replace the kernel by a piecewise-linear or global periodic interpolation. Then we solve the linear systems with these known interpolants. Explicitly, we interpolate the kernel \(K\) with respect to the non-integrating variable

where \(L_j\) is the interpolating basis. We have from the analysis of the interpolation method used that \(K_n \to K\) as \(n \to \infty\). The linear system we need to solve then becomes

In practice the evaluation of these integrals would require the use of quadrature scheme. With this addition the numerical computation of solutions to integral equations is complete. The error analysis of these degenerate approximations found via interpolation can then be done with the standard tools of approximation theory for interpolation and the quadrature methods used.

Alternatively, we can find the degenerate kernels through expanding the kernel in an orthonormal basis. Let \((\cdot, \cdot)\) be scalar product on \(C(G)\) and \(\{u_1, u_2, \dots \}\) be a complete orthonormal basis. Then for each fixed \(y \in G\) the kernel is expanded with respect to \(x\). Giving us the degenerate kernel

From approximation theory we have that \(K_n \to K\) as \(n \to \infty\). The resulting linear system of this degenerate kernel is then

Note that the scalar product given is often an integral. So for each matrix coefficient and element of the solution vector in the linear system we need to compute two integrals. This is not very efficient, and is the tradeoff for the relatively easy error analysis we get by using interpolation and basis expansion techniques.

4.3 Nyström Method

The Nyström method is another way to approximate integral equations with continuous or weakly singular kernels. In this scheme, we force the value of the unknown density layer \(\varphi\) to equal the known boundary data \(f\) through the boundary integral operator precisely at the quadrature nodes we choose. Whereas before, we were interpolating or expanding with respect to the non-integrating variable, now we discretize this variable so that the integral operator

is approximated by

The validity of approximating solutions to \((I - A) \varphi = f\) with \((I-A_n)\varphi_n = f_n\) is given by the following theorem.

The remarkable property of the Nyström method is that the error convergence for the solution to the integral equation is of the same order as the underlying quadrature used to build the discretization [8].

We see one immediate issue with the Nyström method: we will evaluate the kernel at \(K(x_j, x_n)\) when \(j = n\). We have seen for the Laplace and Helmholtz kernels that they are weakly singular. Thus we have a problem not only when \(j=n\) but also for the quadrature nodes near \(j=n\). Since here the kernel will be blowing up at a rate that the smooth quadrature will not be able to keep up with. In this neighborhood we see a degradation of the convergence rate. This shows up on the diagonal of the discretized matrix. To combat this, we introduce the idea of correction schemes in the style of Kapur and Rokhlin [5].

The idea of these correction schemes is that we “fix” a small number of quadrature weights near the singularity to get a better convergence rate. For the purposes of this discussion, assume that our scatterer is closed and smooth, so it is parameterized by a global periodic function \(f(t)\) and we can use the trapezoidal rule for the quadrature. This means all weights will be the same and we denote them \(\alpha = \alpha_1 = \alpha_2 = \dots = \alpha_n\) As stated above, our strategy will be to modify a selection of the quadrature weights. Say we allow ourselves to modify four total weights. Owing to the symmetry of the problem, we really only have two total free variables. Borrowing the notation of [9] we denote the order of the Kapur-Rokhlin correction scheme as the number of modified weights on each side of the singularity \(Q_n^{KR(2)}\). We can write the scheme as

Our job is then to choose these \(\eta_1, \eta_2\) so that the singularity blow up can resolved and we get an improved order of convergence. Kapur and Rokhlin show how to do this in [5] so that error converges like

where \(h\) is the distance between the equispaced point quadrature nodes of the trapezoid rule. In the case of the log singular kernels of Helmholtz and Laplace, this comes from subtracting the singularity and Taylor expanding. The coefficients of the Taylor expansion give the correction weights.

This is not the only way to deal with the singularity in the Nyström method and indeed this is an active area of research.