Abstract and Applied Analysis

Volume 2013 (2013), Article ID 938167, 20 pages

http://dx.doi.org/10.1155/2013/938167

Research Article

Solving Integral Equations on Piecewise Smooth Boundaries Using the RCIP Method: A Tutorial

Johan Helsing

Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden

Received 31 October 2012; Accepted 21 December 2012

Academic Editor: Juan J. Nieto

Copyright © 2013 Johan Helsing. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Recursively compressed inverse preconditioning (RCIP) is a numerical method for obtaining highly accurate solutions to integral equations on piecewise smooth surfaces. The method originated in 2008 as a technique within a scheme for solving Laplace’s equation in two-dimensional domains with corners. In a series of subsequent papers, the technique was then refined and extended as to apply to integral equation formulations of a broad range of boundary value problems in physics and engineering. The purpose of the present paper is threefold: first, to review the RCIP method in a simple setting; second, to show how easily the method can be implemented in MATLAB; third, to present new applications of RCIP to integral equations of scattering theory on planar curves with corners.

1. Introduction

This paper is about an efficient numerical solver for elliptic boundary value problems in piecewise smooth domains. Such a solver is useful for applications in physics and engineering, where computational domains of interest often have boundaries that contain singular points such as corners and triple junctions. Furthermore, elliptic boundary value problems in nonsmooth domains are difficult to solve irrespective of what numerical method is used. The reason is that the solution, or the quantity representing the solution, often exhibits a nonsmooth behavior close to boundary singularities. That behavior is hard to resolve by polynomials, which underly most approximation schemes. Mesh refinement is needed. This is costly and may lead to artificial ill-conditioning and the loss of accuracy.

The numerical solver we propose takes its starting point in an integral equation reformulation of the boundary value problem at hand. We assume that the problem can be modeled as a Fredholm second kind integral equation with integral operators that are compact away from singular boundary points and whose solution is a layer density representing the solution to the original problem. We then seek a discrete approximation to the layer density using a modification of the Nyström method [1, Chapter 4]. At the heart of the solver lie certain integral transforms whose inverses modify the kernels of the integral operators in such a way that the layer density becomes piecewise smooth and simple to resolve by polynomials. The discretizations of these inverses are constructed recursively on locally refined temporary meshes. Conceptually, this corresponds to applying a fast direct solver [2] locally to regions with troublesome geometry. Finally, a global iterative method is applied. This gives us many of the advantages of fast direct methods, for example, the ability to deal with certain classes of operators whose spectra make them unsuitable for iterative methods. In addition, the approach is typically much faster than using only a fast direct solver.

Our method, or scheme, has previously been referred to as recursive compressed inverse preconditioning [3–6], and there is a good reason for that name. The scheme relies on applying a relieving right inverse to the integral equation, on compressing this inverse to a low-dimensional subspace, and on carrying out the compression in a recursive manner. Still, the name recursive(ly) compressed inverse preconditioning is a bit awkward, and we will here simply use the acronym RCIP.

A strong motivation for writing the present paper is that the original references [3–6] are hard to read for nonexperts and that efficient codes, therefore, could be hard to implement. Certain derivations in [3–6] use complicated intermediary constructions, application specific issues obscure the general picture, and the notation has evolved from paper to paper. In the present paper, we focus on the method itself, on how it works and how it can be implemented, and refer to the original research papers for details. Demo codes, in Matlab, are part of the exposition and can be downloaded from http://www.maths.lth.se/na/staff/helsing/Tutor/.

Section 2 provides some historical background. The basics of the RCIP method are then explained by solving a simple model problem in Sections 3–6. Sections 7–15 review various extensions and improvements. Some of this material is new, for example, the simplified treatment of composed operators in Section 14. The last part of the paper, Sections 16–18, contains new applications to scattering problems.

2. Background

The line of research on fast solvers for elliptic boundary value problems in piecewise smooth domains, leading up to the RCIP method, grew out of work in computational fracture mechanics. Early efforts concerned finding efficient integral equation formulations. Corner singularities were either resolved with brute force or by using special basis functions; see [7–9] for examples. Such strategies, in combination with fast multipole [10] accelerated iterative solvers, work well for simple small-scale problems.

Real world physics is more complicated and, for example, the study [11] on a high-order time-stepping scheme for crack propagation (a series of biharmonic problems for an evolving piecewise smooth surface) shows that radically better methods are needed. Special basis functions are too complicated to construct, and brute force is not economical—merely storing the discretized solution becomes too costly in a large-scale simulation.

A breakthrough came in 2007, when a scheme was created that resolves virtually any problem for Laplace's equation in piecewise smooth two-dimensional domains in a way that is fully automatic, fast, stable, memory efficient, and whose computational cost scales linearly with the number of corners in the computational domain. The resulting paper [5] constitutes the origin of the RCIP method. Unfortunately, however, there are some flaws in [5]. The expressions in [5, Section 9] are not generally valid, and the paper fails to apply RCIP in its entirety to the biharmonic problem of [5, Section 3], which was the ultimate goal.

The second paper on RCIP [6] deals with elastic grains. The part [6, Appendix ], on speedup, is particularly useful.

The third paper on RCIP [3] contains improvement relative to the earlier papers, both in the notation and in the discretization of singular operators. The overall theme is mixed boundary conditions, which pose similar difficulties as piecewise smooth boundaries do.

The paper [4], finally, solves the problem of [5, Section 3] in a broad setting, and one can view the formalism of [4] as the RCIP method in its most general form, at least in two dimensions. Further work on RCIP has been devoted to more general boundary conditions [12], to problems in three dimension [13], and to solving elliptic boundary value problems in special situations, such as for aggregates of millions of grains, where special techniques are introduced to deal with problem specific issues [14, 15].

We end this retrospection by pointing out that several research groups in recent years have proposed numerical schemes for integral equations stemming from elliptic partial differential equations in piecewise smooth domains. See, for example, [16–21]. There is also a widespread notion that a slight rounding of corners is a good idea for numerics. While rounding may work in particular situations, we do not believe it is a generally viable method. For one thing, how does one round a triple junction?

3. A Model Problem

Let be the closed contour of Figure 1 with the parameterization

Let be the fundamental solution to Laplace's equation which in the plane can be written as

We will solve the integral equation

numerically for the unknown layer density . Here, is the exterior unit normal of the contour, is an element of arc length, is a parameter, and is a unit vector. Equation (3) models an electrostatic inclusion problem [5].

Figure 1: The contour of (1) with a corner at the origin. The solid curve corresponds to opening angle . The dashed curve has .

Using complex notation, where vectors , , and in the real plane correspond to points , , and in the complex plane , one can express (3) as

where the “bar” symbol denotes complex conjugation. Equation (4) is a simplification over (3) from a programming point of view.

From an algorithmic point of view, it is advantageous to write (3) in the abbreviated form

where is the identity. If is smooth, then (5) is a Fredholm second kind integral equation with a compact, non-self-adjoint, integral operator whose spectrum is discrete, bounded by one in modulus, and accumulates at zero.

We also need a way to monitor the convergence of solutions to (5). For this purpose, we introduce a quantity , which corresponds to dipole moment or polarizability [13]

Remark 1. Existence issues are important. Loosely speaking, the boundary value problem modeled by (3) has a unique finite energy solution for a large class of nonsmooth when is either off the real axis or when is real and . See [13] for sharper statements. The precise meaning of a numerical solution to an integral equation such as (3) also deserves comment. In this paper, a numerical solution refers to approximate values of at a discrete set of points . The values should, in a postprocessor, enable the extraction of quantities of interest including values of at arbitrary points , functionals of such as of (6), and the solution to the underlying boundary value problem at points in the domain where that problem was set.

4. Discretization on Two Meshes

We discretize (5) using the original Nyström method based on composite 16-point Gauss-Legendre quadrature on two different meshes: a coarse mesh with quadrature panels and a fine mesh which is constructed from the coarse mesh by times subdividing the panels closest to the corner in a direction toward the corner. The discretization is in parameter. The four panels on the coarse mesh that are closest to the corner should be equi-sized in parameter. These panels form a subset of called . See Figure 2.

(a)

(b)

Figure 2: (a) A coarse mesh with ten panels on the contour of (1) with opening angle . A subset of , called , covers four coarse panels as indicated by the dashed curve. (b) A fine mesh created from the coarse mesh by subdividing the panels closest to the corner times.

The linear systems resulting from discretization on the coarse mesh and on the fine mesh can be written formally as

where and are square matrices and and are column vectors. The subscripts fin and coa indicate what type of mesh is used. Discretization points on a mesh are said to constitute a grid. The coarse grid has points. The fine grid has points.

The discretization of (5) is carried out by first rewriting (4) as

Then Nyström discretization with points and weights on gives

The program demo1.m sets up the system (8), solves it using the GMRES iterative solver [22] incorporating a low-threshold stagnation avoiding technique [23, Section 8], and computes of (6). The user has to specify the opening angle , the parameter , the number of coarse panels on , the unit vector and the number of subdivisions . The opening angle should be in the interval . We choose , , , and . The quantity converges initially as is increased, but for the results start to get worse. See Figure 3. This is related to the fact, pointed out by Bremer [17], that the original Nyström method captures the behavior of the solution , while our is unbounded. See, further, Appendix D.

(a)

(b)

Figure 3: Convergence for of (6) using (8) and the program demo1b.m (a loop of demo1.m) with lambda = 0.999, theta = pi/2, npan = 10, and evec = 1. The reference value is . There are unknowns in the main linear system. (a) Convergence with . (b) The number of iterations needed to meet an estimated relative residual of .

5. Compressed Inverse Preconditioning

Let us split the matrices and of (7) and (8) into two parts each

Here, the superscript indicates that only entries of a matrix whose indices and correspond to points and that both belong to the boundary segment are retained. The remaining entries are zero.

Now, we introduce two diagonal matrices and which have the quadrature weights on the diagonal. Furthermore, we need a prolongation matrix which interpolates functions known at points on the coarse grid to points on the fine grid. The construction of relies on panelwise 15-degree polynomial interpolation in parameter using Vandermonde matrices. We also construct a weighted prolongation matrix via

The matrices and share the same sparsity pattern. They are rectangular matrices, similar to the the identity matrix, but with one full block. Let superscript denote the transpose. Then,

holds exactly. See Appendix A and [4, Section 4.3].

Equipped with and , we are ready to compress (8) on the fine grid to an equation essentially on the coarse grid. This compression is done without the loss of accuracy—the discretization error in the solution is unaffected and no information is lost. The compression relies on the variable substitution

Here, is the discretization of a piecewise smooth transformed density. The compression also uses the low-rank decomposition

which should hold to about machine precision.

The compressed version of (8) reads

where the compressed weighted inverse is given by

See Appendix B for details on the derivation. The compressed weighted inverse , for of (1), is a block diagonal matrix with one full block and the remaining entries coinciding with those of the identity matrix.

After having solved (17) for , the density can easily be reconstructed from in a postprocessor; see Section 9. It is important to observe, however, that is not always needed. For example, the quantity of (6) can be computed directly from . Let be a column vector which contains the absolute values of the boundary derivatives multiplied with weights , positions , and the complex scalar . Then,

6. The Recursion for R

The compressed weighted inverse is costly to compute from its definition (18). As we saw in Section 4, the inversion of large matrices on highly refined grids could also be unstable. Fortunately, the computation of can be greatly sped up and stabilized via a recursion. This recursion is derived in a roundabout way and using a refined grid that differs from that of the present tutorial in [5, Section 7.2]. A better derivation can be found in [4, Section 5], but there the setting is more general so that text could be hard to follow. Here, we focus on results.

6.1. Basic Prolongation Matrices

Let be a prolongation matrix, performing panelwise 15-degree polynomial interpolation in parameter from a 64-point grid on a four-panel mesh to a 96-point grid on a six-panel mesh as shown in Figure 4. Let be a weighted prolongation matrix in the style of (13). If T16 and W16 are the nodes and weights of 16-point Gauss-Legendre quadrature on the canonical interval , then and can be constructed as

                                                         end            %                           %                         

(a)

(b)

Figure 4: (a) The prolongation operator performs panelwise interpolation from a grid on a four-panel mesh to a grid on a six-panel mesh. (b) The sparsity pattern of .

See [23, Appendix ] for an explanation of why high-degree polynomial interpolation involving ill-conditioned Vandermonde systems gives accurate results for smooth functions.

6.2. Discretization on Nested Meshes

Let , , be a sequence of subsets of with and . Let there also be a six-panel mesh and a corresponding 96-point grid on each . The construction of the subsets and their meshes should be such that if , , is a local parameterization of , then the breakpoints (locations of panel endpoints) of its mesh are at , and the breakpoints of the mesh on are at . We denote this type of nested six-panel meshes type b. The index is the level. An example of a sequence of subsets and meshes on is shown in Figure 5 for . Compare [3, Figure 2] and [4, Figure 5.1].

Figure 5: The boundary subsets , , and along with their corresponding type b meshes for .

Let denote the discretization of on a type b mesh on . In the spirit of (11) and (12), we write

where the superscript indicates that only entries with both indices corresponding to points on the four inner panels are retained.

6.3. The Recursion Proper

Now, let denote the full diagonal block of . The recursion for is derived in Appendix C, and it reads

where the operator expands its matrix argument by zero-padding (adding a frame of zeros of width 16 around it). Note that the initializer (22) makes the recursion (21) take the first step

The program demo2.m sets up the linear system (17), runs the recursion of (21) and (22), and solves the linear system using the same techniques as demo1.m; see Section 4. In fact, the results produced by the two programs are very similar, at least up to . This supports the claim of Section 5 that the discretization error in the solution is unaffected by compression.

Figure 6 demonstrates the power of RCIP: fewer unknowns and faster execution, better conditioning (the number of GMRES iterations does not grow), and higher achievable accuracy. Compare Figure 3. We emphasize that the number of recursion steps (levels) used in (21) corresponds to the number of subdivisions used to construct the fine mesh.

(a)

(b)

Figure 6: Same as Figure 3, but using (17) and the program demo2b.m (a loop of demo2.m). There are only unknowns in the main linear system.

7. Schur-Banachiewicz Speedup of the Recursion

The recursion (21) can be sped up using the Schur-Banachiewicz inverse formula for partitioned matrices [24], which in this context can be written [6, Appendix ] as

where plays the role of , and are submatrices of and , and , , and refer to blocks of .

The program demo3.m is based on demo2.m, but has (24) incorporated. Besides, the integral equation (3) is replaced with

which has the same solution but is more stable for close to one. For the discretization of (25) to fit the form (17), the last term on the left hand side of (25) is added to the matrix of (17).

The execution of demo3.m is faster than that of demo2.m. Figure 7 shows that a couple of extra digits are gained by using (25) rather than (3) and that full machine accuracy is achieved for .

(a)

(b)

Figure 7: Same as Figure 6, but the program demo3b.m is used.

8. Various Useful Quantities

Let us introduce a new discrete density via

Rewriting (17) in terms of gives

which resembles the original equation (7). We see that , which is discretized using Gauss-Legendre quadrature, acts on . Therefore, one can interpret as pointwise values of the original density , multiplied with weight corrections suitable for integration against polynomials. We refer to as a weight-corrected density. Compare [6, Section 5.4].

Assume now that there is a square matrix which maps to discrete values of the original density on the coarse grid

The matrix allows us to rewrite (17) as a system for the original density

We can interpret the composition as a matrix of multiplicative weight corrections that compensate for the singular behavior of on when Gauss-Legendre quadrature is used.

Let denote the rectangular matrix

and let be a restriction operator which performs panelwise 15-degree polynomial interpolation in parameter from a grid on the fine mesh to a grid on a the coarse mesh. We see from (15) that is the mapping from to . Therefore the columns of can be interpreted as discrete basis functions for . It holds by definition that

The quantities and interpretations of this section come in handy in various situations, for example, in 3D extensions of the RCIP method [13]. An efficient scheme for constructing will be presented in Section 10.

9. Reconstruction of from

The action of on , which gives , can be obtained by, in a sense, running the recursion (21) backwards. The process is described in detail in [3, Section 7]. Here, we focus on results.

The backward recursion on reads

Here, is a column vector with elements. In particular, is the restriction of to , while are taken as elements of for . The elements and of are the reconstructed values of on the outermost panels of a type mesh on . Outside of , coincides with .

When the recursion is completed, the reconstructed values of on the four innermost panels are obtained from

Should one wish to interrupt the recursion (33) prematurely, say at step , then

gives values of a weight-corrected density on the four innermost panels of a type mesh on . That is, we have a part-way reconstructed weight-corrected density on a mesh that is times refined. This observation is useful in the context of evaluating layer potentials close to their sources.

If the memory permits, one can store the matrices and in the forward recursion (21) and reuse them in the backward recursion (33). Otherwise, they may be computed afresh.

The program demo4.m builds on the program demo3.m, using (17) for (25). After the main linear system is solved for , a postprocessor reconstructs via (33). Then, a comparison is made with a solution obtained by solving the uncompressed system (8). Figure 8 shows that for , the results are virtually identical. This verifies the correctness of (33). For , the results start to deviate. That illustrates the instabilities associated with solving (8) on a highly refined mesh. Compare Figure 3.

(a)

(b)

Figure 8: Output from demo4.m and demo5.m. (a) A comparison of from the unstable equation (8) and reconstructed from of (17) via (33). (b) Relative accuracy in of (6) from part-way reconstructed solutions .

The program demo5.m investigates the effects of premature interruption of (33). The number of recursion steps is set to , and the recursion is interrupted at different levels. The density is reconstructed on outer panels up to the level of interruption. Then, a weight-corrected density is produced at the innermost four panels according to (35). Finally, of (6) is computed from this part-way reconstructed solution. Figure 8(b) shows that the quality of is unaffected by the level of interruption.

10. The Construction of

This section discusses the construction of and other auxiliary matrices. Note that in many applications, these matrices are not needed.

The entries of the matrices , , , , , and can only differ from those of the identity matrix when both indices correspond to discretization points on . For example, the entries of only differ from the identity matrix for the block denoted in (21). In accordance with this notation, we introduce , , , , and for the restriction of , , , , and to . In the codes of this section, we often use this restricted type of matrices, leaving the identity part out.

We observe that is a square matrix, , and are rectangular matrices, and is a rectangular matrix. Furthermore, is very sparse for large . All columns of with column indices corresponding to points on panels that result from more than eight subdivisions are identically zero.

The program demo6.m sets up , , and , shows their sparsity patterns, and verifies the identities (14) and (31). The implementations for and rely on repeated interpolation from coarser to finer intermediate grids. The implementation of relies on keeping track of the relation between points on the original coarse and fine grids. Output from demo6.m is depicted in Figure 9(a). Note that the matrices and are never needed in applications.

(a)

(b)

Figure 9: (a) The identities (14) and (31) hold to high accuracy in our implementation, irrespective of the degree of mesh refinement. (b) The solution to (25) on (1) with parameters as specified in Section 4. The solution with RCIP (17) and (28), shown as blue stars, agrees with the solution from (8), shown as a red solid line. The solution diverges in the corner.

We are now ready to construct . Section 9 presented a scheme for evaluating the action of on discrete functions on the coarse grid on . The matrix , itself, can be constructed by applying this scheme to a identity matrix. The matrix was set up in demo6.m. Composing these two matrices gives ; see (32). This is done in the program demo7.m, where the identity part is added as to get the entire matrix . In previous work on RCIP, we have found use for in complex situations where (29) is preferable over (17); see [13, Section 9]. If one merely needs from in a postprocessor, setting up and using (28) are not worthwhile. It is cheaper to let act on and then let act on the resulting vector. Anyhow, demo7.m builds on demo4.m and gives as output computed via (28); see Figure 9(b). For comparison, , computed from the uncompressed system (8), is also shown.

11. Initiating R Using Fixed-Point Iteration

It often happens that is wedge-like. A corner of a polygon, for example, has wedge-like at all levels. If is merely piecewise smooth, then the are wedge-like to double precision accuracy for .

Wedge-like sequences of open up for simplifications and speedup in the recursion of (21) and (22). Particularly so if the kernel of the integral operator of (5) is scale invariant on wedges. Then the matrix becomes independent of . It can be denoted and needs only to be constructed once. Furthermore, the recursion of (21) and (22) assumes the form of a fixed-point iteration

The iteration (36) can be run until converges. Let be the converged result. One needs not worry about predicting the number of levels needed. Choosing the number of levels needed, in order to meet a beforehand given tolerance in , is otherwise a big problem in connection with (21) and (22) and non-wedge-like . This number has no general upper bound.

Assume now that the kernel of is scale invariant on wedges. If all are wedge-like, then (36) and (37) replace (21) and (22) entirely. If is merely piecewise smooth, then (36) and (37) can be run on a wedge with the same opening angle as , to produce an initializer to (21). That initializer could often be more appropriate than (22), which is plagued with a very large discretization error whenever (10) is used. This fixed-point recursion initializer is implemented in the program demo8b.m, which is a simple upgrading of demo3b.m, and produces Figure 10. A comparison of Figure 10 with Figure 7 shows that the number of levels needed for full convergence is halved compared to when using the initializer (22).

(a)

(b)

Figure 10: Same as Figure 7, but the program demo8b.m is used.

There are, generally speaking, several advantages with using the fixed-point recursion initializer, rather than (22), in (21) on a non-wedge-like . First, the number of different matrices and needed in (21) and in (33) is reduced as the recursions are shortened. This means savings in storage. Second, the number of levels needed for full convergence in (21) seems to always be bounded. The hard work is done in (36). Third, Newton's method can be used to accelerate (36). That is the topic of Section 12.

12. Newton Acceleration

When solving integral equations stemming from particularly challenging elliptic boundary value problems with solutions that are barely absolutely integrable, the fixed-point iteration of (36) and (37) on wedge-like may need a very large number of steps to reach full convergence. See [15, Section 6.3] for an example where steps are needed.

Fortunately, (36) can be cast as a nonlinear matrix equation

where , as in Section 11, is the fixed-point solution and

The nonlinear equation (38), in turn, can be solved for with a variant of Newton’s method. Let be a matrix-valued perturbation of , and expand to first order in . This gives a Sylvester-type matrix equation

for the Newton update . One can use the Matlab built-in function dlyap for (40), but GMRES seems to be more efficient and we use that method. Compare [15, Section 6.2].

Figure 11 shows a comparison between the fixed-point iteration of (36) and (37) and Newton’s method for computing the fixed-point solution to (38) on a wedge-like . The program demo9.m is used, and it incorporates Schur-Banachiewicz speedup in the style of Section 7. The wedge opening angle is , the integral operator is the same as in (5), and . The relative difference between the two converged solutions is . Figure 11 clearly shows that (36) and (37) converge linearly (in 68 iterations), while Newton's method has quadratic convergence. Only four iterations are needed. The computational cost per iteration is, of course, higher for Newton's method than for the fixed-point iteration, but it is the same at each step. Recall that the underlying size of the matrix , that is inverted according to (18), grows linearly with the number of steps needed in the fixed-point iteration. This example therefore demonstrates that one can invert and compress a linear system of the type (18) in sublinear time.

Figure 11: Output from the program demo9.m. The fixed-point iteration (36) and (37) is compared to Newton’s method for (38).

13. On the Accuracy of the “Solution”

The integral equation (3) comes from a boundary value problem for Laplace's equation where the potential field at a point in the plane is related to the layer density via

see [5, Section 2.1].

Figure 12 shows how the field solution of (41) converges with mesh refinement at a point inside the contour of (1). We see that the accuracy in , typically, is one digit better than the accuracy of the dipole moment of (6). One can therefore say that measuring the field error at a point some distance away from the corner is more forgiving than measuring the dipole moment error. It is possible to construct examples where the differences in accuracy between field solutions and moments of layer densities are more pronounced, and this raises the question of how accuracy should be best measured in the context of solving integral equations. We do not have a definite answer, but we think that higher moments of layer densities give more fair estimates of overall accuracy than field solutions do at points some distance away from the boundary. See, further, Appendix E.

Figure 12: Similar as Figure 7(a) but with the program demo3c.m. The potential at is evaluated via (41), and npan = 14 is used.

14. Composed Integral Operators

Assume that we have a modification of (5) which reads

Here, and are as in (5), is a new, bounded integral operator, and is an unknown layer density to be solved for. This section shows how to apply RCIP to (42) using a simplified version of the scheme in [4].

Let us, temporarily, expand (42) into a system of equations by introducing a new layer density . Then,

and after discretization on the fine mesh,

Standard RCIP gives

where the compressed inverse is partitioned into four equi-sized blocks.

Now, we replace and with a single unknown via

The change of variables (46) is chosen so that the second block-row of (45) is automatically satisfied. The first block-row of (45) becomes

When (47) has been solved for , the weight-corrected version of the original density can be recovered as

Figure 13 shows results for (42) with being the double layer potential

Figure 13(a) shows the convergence of of (6) with using the inner product preserving discretization scheme of Appendix D for (42) as implemented in demo10.m. Figure 13(b) shows produced with RCIP according to (47) and (48) as implemented in demo10b.m. The reference value for is computed with the program demo10c.m, which uses inner product preserving discretization together with compensated summation [25, 26] in order to enhance the achievable accuracy. One can see that, in addition to being faster, RCIP gives and extra digit of accuracy. Actually, it seems as if the scheme in demo10.m converges to a that is slightly wrong.

(a)

(b)

Figure 13: Convergence for of (6) with from (42). The curve is as in (1), and theta = pi/2, npan = 11, and evec = 1. The reference value is taken as . (a) Results from the inner product preserving scheme of Appendix D produced with demo10.m. (b) Results with RCIP according to (47) and (48) produced with demo10b.m.

In conclusion, in this example and in terms of stability, the RCIP method is better than standard inner product preserving discretization and on par with inner product preserving discretization enhanced with compensated summation. In terms of computational economy and speed, RCIP greatly outperforms the two other schemes.

15. Nyström Discretization of Singular Kernels

The Nyström scheme of Section 4 discretizes (5) using composite 16-point Gauss-Legendre quadrature. This works well as long as the kernel of the integral operator and the layer density are smooth on smooth . When the kernel is not smooth on smooth , the quadrature fails and something better is needed. See [27] for a comparison of the performance of various modified high-order accurate Nyström discretizations for weakly singular kernels and [28] for a recent high-order general approach to the evaluation of layer potentials.

We are not sure which modified discretization is optimal in every situation. When logarithmic- and Cauchy-singular operators need to be discretized in Sections 16–18, we use two modifications to composite Gauss-Legendre quadrature called local panelwise evaluation and local regularization. See [3, Section 2] for a description of these techniques.

16. The Exterior Dirichlet Helmholtz Problem

Let be the domain enclosed by the curve , and let be the exterior to the closure of . The exterior Dirichlet problem for the Helmholtz equation

has a unique solution under mild assumptions on and [29] and can be modeled using a combined integral representation [30, Chapter 3]

where is the fundamental solution to the Helmholtz equation in two dimensions

Here, is a Hankel function of the first kind. Insertion of (53) into (51) gives the combined field integral equation

where

Figure 14 shows the performance of RCIP applied to (55) for 1000 different values of . The program demo11.m is used. The boundary is as in (1) with , and the boundary conditions are chosen as with inside . The error in of (53) is evaluated at outside . Since the magnitude of varies with , peaking at about unity, the absolute error is shown rather than the relative error. The number of panels on the coarse mesh is chosen as npan = 0.6 * omega + 18 rounded to the nearest integer.

(a)

(b)

Figure 14: The exterior Dirichlet problem for Helmholtz equation with RCIP applied to (55). The program demo11.m is used with as in (1), and . The boundary condition of (51) is generated by a point source at . (a) The absolute error in at . (b) The number of GMRES iterations needed to meet an estimated relative residual of .

17. The Exterior Neumann Helmholtz Problem

The exterior Neumann problem for the Helmholtz equation

has a unique solution under mild assumptions on and [29] and can be modeled as an integral equation in several ways. We will consider two options: an “analogy with the standard approach for Laplace’s equation,” which is not necessarily uniquely solvable for all , and a “regularized combined field integral equation,” which is always uniquely solvable. See, further, [16, 31].

17.1. An Analogy with the Standard Laplace Approach

Let be the adjoint to the double-layer integral operator of (56)

Insertion of the integral representation

into (59) gives the integral equation

Figure 15 shows the performance of RCIP applied to (63). The program demo12.m is used, and the setup is the same as that for the Dirichlet problem in Section 16. A comparison between Figures 15 and 14 shows that the number of GMRES iterations needed for full convergence now grows much faster with . Furthermore, the relative error in the solution to the Neumann problem is larger, particularly when happens to be close to values for which the operator in (63) has a nontrivial nullspace. Recall that (55) is always uniquely solvable, while (63) is not.

(a)

(b)

Figure 15: The exterior Neumann problem for Helmholtz equation with RCIP applied to (63). The program demo12.m is used with as in (1), and . The boundary condition of (59) is generated by a point source at . (a) The absolute error in at . (b) The number of GMRES iterations needed to meet an estimated relative residual of .

17.2. A Regularized Combined Field Integral Equation

The literature on regularized combined field integral equations for the exterior Neumann problem is rich, and several formulations have been suggested. We will use the representation [31]

which after insertion into (59) gives the integral equation

where

The hypersingular operator of (66) can be expressed as a sum of a simple operator and an operator that requires differentiation with respect to arc length only [32]

This makes it possible to write (65) in the form

where , , and the action of the operators , , and is given by

All integral operators in (68) are such that their discretizations admit the low-rank decomposition (16). We use the temporary expansion technique of Section 14 for (68), with two new layer densities that are later eliminated, to arrive at a single compressed equation analogous to (47). That equation involves nine equi-sized blocks of the compressed inverse .

Solving the problem in the example of Section 17.1 again, we now take the number of panels on the coarse mesh as npan = 0.6 * omega + 48 rounded to the nearest integer. Figure 16 shows results from the program demo13b.m The resonances, visible in Figure 15, are now gone. It is interesting to observe in Figure 16 that, despite the presence of several singular operators and compositions in (68), the results produced with RCIP are essentially fully accurate, and the number of GMRES iterations needed for convergence grows very slowly with .

(a)

(b)

Figure 16: The same exterior Neumann problem for Helmholtz equation as in Figure 15, but RCIP is now applied to (68). The program demo13b.m is used.

The program demo13c.m differs from demo13b.m in that it uses local regularization for the Cauchy-singular operators of (70) and (71) rather than local panelwise evaluation. The results produced by the two programs are virtually identical; so, we do not show yet another figure.

18. Field Evaluations

Strictly speaking, a boundary value problem is not properly solved until its solution can be accurately evaluated in the entire computational domain. The program demo11b.m is a continuation of demo11.m which, after solving (55) for with RCIP and forming via (26), computes the solution via (53) using three slightly different discretizations.

(i)When is away from , 16-point Gauss-Legendre quadrature is used in (53) on all quadrature panels. (ii)When is close to , but not close to a panel neighboring a corner, 16-point Gauss-Legendre quadrature is used in (53) on panels away from , and local panelwise evaluation is used for panels close to . (iii)When is close to a panel neighboring a corner, the density is first used to reconstruct according to Section 9. Then 16-point Gauss-Legendre quadrature is used in (53) on panels away from , and local panelwise evaluation is used for panels close to .

The first two discretizations only use the coarse grid on . The third discretization needs a grid on a partially refined mesh on .

The program demo13d.m is a continuation of demo13b.m which, after solving (65) with RCIP as described in Section 17.2, computes the solution via (64) using the three discretizations of the previous paragraph.

Figures 17 and 18 show that RCIP in conjunction with the modified quadrature rules of [3, Section 2] is capable of producing very accurate solutions to exterior Helmholtz problems in, essentially, the entire computational domain.

Figure 17: The error in the solution to the exterior Dirichlet Helmholtz problem of Section 16 evaluated at 62392 points on a cartesian grid using demo11b.m with . The source at is shown as a blue star.

Figure 18: Same as Figure 17, but the exterior Neumann Helmholtz problem is solved using demo13d.m. The accuracy is even higher than in Figure 17.

Appendices

A. Proof That

Let and be two column vectors, corresponding to the discretization of two panelwise polynomials with panelwise degree 15 on the coarse mesh of . Then,

because composite 16-point Gauss-Legendre quadrature has panelwise polynomial degree 31. The diagonal matrix has siz e .

Since there are linearly independent choices of and of , it follows from (A.1) that

which, using (13), can be rewritten as

B. Derivation of the Compressed Equation

The compression of (8), leading up to (17), was originally described in [5, Section 6.4]. Here, we give a summary.

The starting point is (5) which, using the operator split analogous to (11) and (12)

and the variable substitution

gives the right preconditioned equation

Now, let us take a close look at (B.3). We observe that is an operator whose action on any function gives a function that is smooth on the innermost two panels of the coarse grid on . This is so since is constructed so that its action on any function gives a function that is smooth on the innermost two panels of the coarse grid on . Furthermore, the right hand side of (B.3) is assumed to be panelwise smooth. Using an argument of contradiction, we see that has to be a panelwise smooth function on the innermost two panels of the coarse grid on .

Having concluded that is panelwise smooth close to the corner, we can write

We also have

the discrete version of (B.2) on the fine grid

and the relations (12) and (16) which we now repeat as

Substitution of (B.4), (B.5), (B.6), and (B.7) into (8), which we now repeat as

gives

Applying (or ) to the left in (B.9) and using the identities (14) (or (31)) give the compressed equation (17).

C. Derivation of the Recursion

The recursion (21) for the rapid construction of the diagonal blocks of the compressed weighted inverse was originally derived in [5, Section 7] using different notations and different meshes than in the present tutorial. The recursion was derived a second time in [6, Section 7] using new meshes. Better notation was introduced in [3, Section 6]. A third derivation, in a general setting, takes place in [4, Section 5], and it uses the same notation and meshes as in the present tutorial.

A problem when explaining the derivation of (21) is that one needs to introduce intermediate meshes and matrices whose appearance may cause enervation at a first glance. Particularly, since these meshes and matrices are not needed in the final expression (21), we emphasize that the underlying matrix property that permits the recursion is the low rank of certain off-diagonal blocks in discretizations of of (B.1) on nested meshes.

The recursion (21) only uses one type of mesh explicitly—the type b mesh of Figure 5. On each , there is a type b mesh and a corresponding discretization of denoted . Here, we need two new types of meshes denoted, type a and type c, along with corresponding discrete operators. For example, is the discretization of on a type a mesh on . The three types of meshes are depicted in Figure 19. Actually, a straight type c mesh was already introduced in Figure 4.

Figure 19: Meshes of type a, type b, and type c on the boundary subset for and . The type a mesh has panels. The type b mesh has six panels. The type c mesh has four panels. The type a mesh is the restriction of the fine mesh to . For , the type c mesh is the restriction of the coarse mesh to . The type a mesh and the type b mesh coincide for .

Now, we define as

where and are prolongation operators (in parameter) from a grid on a type c mesh on to a grid on a type a mesh on . Note that for , according to the definition (C.1), is identical to the full diagonal block of of (18). Note also that comes cheaply. The rest of this appendix is about finding an expression for in terms of that is cheap to compute.

Let us split into two parts

where , and is such that

holds to about machine precision; compare (16). The prolongation operators and act from a grid on a type b mesh to a grid on a type a mesh. It holds that

Summing up, we can rewrite (C.1) as

The subsequent steps in the derivation of (21) are to expand the term within parenthesis in (C.5) in a Neumann series, multiply the terms in this series with from the left and with from the right, and bring the series back in closed form. The result is

which, in fact, is (21) in disguise. To see this, recall from (C.1) that

Then,

where the second equality uses ; see Appendix A. Substitution of (C.8) in (C.6) gives the recursion in the familiar form

D. An Inner Product Preserving Scheme

In [17], Bremer describes a scheme that stabilizes the solution to the discretized system (8) on the fine mesh. The scheme can be interpreted as an inner product preserving discretization. In practice, it corresponds to making a similarity transformation of the system matrix. While inner product preserving Nyström discretization elegantly solves problems related to stability (the condition number of the system matrix is improved), it does not reduce the number of discretization points (unknowns) needed to achieve a given precision in the solution. Neither does it affect the spectrum of the system matrix (similarity transformations preserve eigenvalues), and, hence, it does not in any substantial way improve the convergence rate of the GMRES iterative method [33, Lecture 35].

For completeness, we have implemented inner product preserving Nyström discretization in the program demo1d.m. The program is a continuation of demo1b.m where we also have replaced (3) with the more stable integral equation (25). This should facilitate comparison with the program demo3b.m and the results shown in Figure 7.

Figure 20 shows results produced by demo1d.m. Beyond , one now achieves essentially full machine precision in of (6). Despite this success, inner product preserving Nyström discretization can perhaps not quite compete with the RCIP method in this example. The differences in performance relate to issues of memory and speed. The RCIP method uses a much smaller linear system ( unknowns) than does inner product preserving Nyström discretization ( unknowns). Besides, the RCIP method converges in only eight GMRES iterations, irrespective of . See Figure 7.

(a)

(b)

Figure 20: Same as Figure 7, but the program demo1d.m is used.

E. The Nature of Layer Densities in Corners

It is considered difficult to solve integral equations on piecewise smooth boundaries. One reason being that solutions (layer densities) often have complicated asymptotics close to corner vertices. Also, stable numerical schemes, such as that of Appendix D, are burdened with the task of resolution. It takes a lot of mesh refinement to resolve nonpolynomial-like layer densities with polynomial basis functions.

Sometimes, layer densities diverge, and sometimes they stay bounded, but they are seldom polynomial-like close to corner vertices. While we refrain from giving a characterization of what properties of Fredholm second kind integral equations give rise to unbounded solutions, we observe the following.

(i)The solution to (3) diverges in corners whenever and is such that exists. Figure 9 illustrates this for . (ii)The size of the zone where the asymptotic shape of a layer density is visible depends on the smoothness of the solution to the PDE that underlies the integral equation. See [3, Figure 3]. (iii)The divergence of layer densities close to corner vertices can be very strong. In classical materials science applications, it is not unusual with layer densities that barely lie in [5, Section 2]. In metamaterial applications, layer densities may barely lie in [15, Section 3]. Furthermore, it is likely that spaces are not the most natural function spaces for the characterization of layer densities [13, Section 5].

Variable separation is often used to predict the asymptotic behavior of layer densities close to corner vertices. See, for example, [5, Section 2] and [7, Section 2]. But perhaps asymptotic studies of layer densities are not that important? Assume that there are two different Fredholm second kind integral equations for the same underlying PDE: one where the layer density is bounded and one where it is unbounded. If the RCIP method is to be used, it really does not matter which equation is chosen. The recursion of (21) and (22) might need fewer steps for the equation with a bounded solution. The achievable accuracy for functionals of the layer density, on the other hand, is often slightly higher if computed from the unbounded solution. The reason being, loosely speaking, that more information can be contained in a rapidly varying function than in a slowly varying function.

Acknowledgments

The idea to make this tutorial came up during a discussion with Alex Barnett and Adrianna Gillman on 5/20/12 at the FACM’12 conference at NJIT. Input and feedback from these colleagues and from Rikard Ojala have been of great value. The work was supported by the Swedish Research Council under Contract 621-2011-5516.

References

K. E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, UK, 1997.

W. Y. Kong, J. Bremer, and V. Rokhlin, “An adaptive fast direct solver for boundary integral equations in two dimensions,” Applied and Computational Harmonic Analysis, vol. 31, no. 3, pp. 346–369, 2011.

J. Helsing, “Integral equation methods for elliptic problems with boundary conditions of mixed type,” Journal of Computational Physics, vol. 228, no. 23, pp. 8892–8907, 2009.

J. Helsing, “A fast and stable solver for singular integral equations on piecewise smooth curves,” SIAM Journal on Scientific Computing, vol. 33, no. 1, pp. 153–174, 2011.

J. Helsing and R. Ojala, “Corner singularities for elliptic problems: integral equations, graded meshes, quadrature, and compressed inverse preconditioning,” Journal of Computational Physics, vol. 227, no. 20, pp. 8820–8840, 2008.

J. Helsing and R. Ojala, “Elastostatic computations on aggregates of grains with sharp interfaces, corners, and triple-junctions,” International Journal of Solids and Structures, vol. 46, pp. 4437–4450, 2009.

J. Helsing, “Corner singularities for elliptic problems: special basis functions versus ‘brute force’,” Communications in Numerical Methods in Engineering, vol. 16, pp. 37–46, 2000.

J. Helsing and A. Jonsson, “On the computation of stress fields on polygonal domains with V-notches,” International Journal for Numerical Methods in Engineering, vol. 53, pp. 433–454, 2002.

J. Helsing and G. Peters, “Integral equation methods and numerical solutions of crack and inclusion problems in planar elastostatics,” SIAM Journal on Applied Mathematics, vol. 59, no. 3, pp. 965–982, 1999.

L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” Journal of Computational Physics, vol. 73, no. 2, pp. 325–348, 1987.

J. Englund, “A higher order scheme for two-dimensional quasi-static crack growth simulations,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 21–24, pp. 2527–2538, 2007.

R. Ojala, “A robust and accurate solver of Laplace’s equation with general boundary conditions on general domains in the plane,” Journal of Computational Mathematics, vol. 30, pp. 433–448, 2012.

J. Helsing and K.-M. Perfekt, “On the polarizability and capacitance of the cube,” Applied and Computational Harmonic Analysis, vol. 34, no. 3, pp. 445–468, 2013.

J. Helsing, “The effective conductivity of random checkerboards,” Journal of Computational Physics, vol. 230, no. 4, pp. 1171–1181, 2011.

J. Helsing, “The effective conductivity of arrays of squares: large random unit cells and extreme contrast ratios,” Journal of Computational Physics, vol. 230, pp. 7533–7547, 2011.

J. Bremer, “A fast direct solver for the integral equations of scattering theory on planar curves with corners,” Journal of Computational Physics, vol. 231, no. 4, pp. 1879–1899, 2012.

J. Bremer, “On the Nyström discretization of integral equations on planar curves with corners,” Applied and Computational Harmonic Analysis, vol. 32, no. 1, pp. 45–64, 2012.

J. Bremer, Z. Gimbutas, and V. Rokhlin, “A nonlinear optimization procedure for generalized Gaussian quadratures,” SIAM Journal on Scientific Computing, vol. 32, no. 4, pp. 1761–1788, 2010.

J. Bremer and V. Rokhlin, “Efficient discretization of Laplace boundary integral equations on polygonal domains,” Journal of Computational Physics, vol. 229, no. 7, pp. 2507–2525, 2010.

J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” Journal of Computational Physics, vol. 229, no. 22, pp. 8259–8280, 2010.

O. P. Bruno, J. S. Ovall, and C. Turc, “A high-order integral algorithm for highly singular PDE solutions in Lipschitz domains,” Computing, vol. 84, no. 3-4, pp. 149–181, 2009.

Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 7, no. 3, pp. 856–869, 1986.

J. Helsing and R. Ojala, “On the evaluation of layer potentials close to their sources,” Journal of Computational Physics, vol. 227, no. 5, pp. 2899–2921, 2008.

H. V. Henderson and S. R. Searle, “On deriving the inverse of a sum of matrices,” SIAM Review, vol. 23, no. 1, pp. 53–60, 1981.

N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, Pa, USA, 1996.

W. Kahan, “Further remarks on reducing truncation errors,” Communications of the ACM, vol. 8, no. 1, p. 40, 1965.

S. Hao, A. H. Barnett, P. G. Martinsson, and P. Young, “High-order accurate Nyström discretization of integral equations with weakly singular kernels on smooth curves in the plane,” submitted, http://arxiv.org/abs/1112.6262v2

A. Klöckner, A. Barnett, L. Greengard, and M. O’Neil, “Quadrature by expansion: a new method for the evaluation of layer potentials,” submitted, http://arxiv.org/abs/1207.4461v1

M. Mitrea, “Boundary value problems and Hardy spaces associated to the Helmholtz equation in Lipschitz domains,” Journal of Mathematical Analysis and Applications, vol. 202, no. 3, pp. 819–842, 1996.

D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, vol. 93, Springer, Berlin, Germany, 2nd edition, 1998.

O. P. Bruno, T. Elling, and C. Turc, “Regularized integral equations and fast high-order solvers for sound-hard acoustic scattering problems,” International Journal for Numerical Methods in Engineering, vol. 91, pp. 1045–1072, 2012.

R. Kress, “On the numerical solution of a hypersingular integral equation in scattering theory,” Journal of Computational and Applied Mathematics, vol. 61, no. 3, pp. 345–360, 1995.

L. N. Trefethen and D. Bau, III, Numerical Linear Algebra, SIAM, Philadelphia, Pa, USA, 1997.