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

Authors: Johan Helsing
Journal: Abstract and Applied Analysis
 (2013)

---

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


---

## Body

## 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 asrecursive 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/.Section2 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 B], 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
(1)r(s)=sin(πs)(cos((s-0.5)θ),sin((s-0.5)θ)),(πs)(cos((s-0.5)θ),sin((s-0.5)θ))s∈[0,1].
Let G(r,r′) be the fundamental solution to Laplace's equation which in the plane can be written as
(2)G(r,r′)=-12πlog|r-r′|.
We will solve the integral equation
(3)ρ(r)+2λ∫Γ∂∂νrG(r,r′)ρ(r′)dσr′=2λ(e·νr),r∈Γ,
numerically for the unknown layer density ρ(r). Here, ν is the exterior unit normal of the contour, dσ is an element of arc length, λ is a parameter, and e 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  θ=π/3. The dashed curve has  θ=4π/3.Using complex notation, where vectorsr, r′, and ν in the real plane ℝ2 correspond to points z, τ, and n in the complex plane ℂ, one can express (3) as
(4)ρ(z)+λπ∫Γρ(τ)ℑ{nzn-τdττ-z}=2λℜ{e-nz},z∈Γ,
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
(5)(I+λK)ρ(z)=λg(z),z∈Γ,
where I  is the identity. If  Γ is smooth, then (5) is a Fredholm second kind integral equation with a compact, non-self-adjoint, integral operator K whose spectrum is discrete, bounded by one in modulus, and accumulates at zero.We also need a way to monitor the convergence of solutionsρ(z) to (5). For this purpose, we introduce a quantity q, which corresponds to dipole moment or polarizability [13]
(6)q=∫Γρ(z)ℜ{e-z}d|z|.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 λ∈[-1,1). 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 ρ(r) at a discrete set of points ri∈Γ. The values ρ(ri) should, in a postprocessor, enable the extraction of quantities of interest including values of ρ(r) at arbitrary points r∈Γ, functionals of ρ(r)  such as q 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 npan quadrature panels and a fine mesh which is constructed from the coarse mesh by nsub 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) A coarse mesh with ten panels on the contourΓ of (1) with opening angle θ=π/2. 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  nsub=3 times.
(a)
(b)The linear systems resulting from discretization on the coarse mesh and on the fine mesh can be written formally as(7)(Icoa+λKcoa)ρcoa=λgcoa,(8)(Ifin+λKfin)ρfin=λgfin,
where I and K are square matrices and ρ and g 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 np=16npan points. The fine grid has np=16(npan+2nsub) points.The discretization of (5) is carried out by first rewriting (4) as
(9)ρ(z(s))+λπ∫01ρ(τ(t))ℜ{nz(s)|τ′(t)|dtτ(t)-z(s)}=2λℜ{e-nz(s)},s∈[0,1].
Then Nyström discretization with np points zi and weights wi on Γ gives
(10)ρi+λπ∑j=1npρjℜ{ni|zj′|wjzj-zi}=2λℜ{e-ni},i=1,2,…,np.The programdemo1.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 q of (6). The user has to specify the opening angle θ, the parameter λ, the number npan of coarse panels on Γ, the unit vector e and the number of subdivisions nsub. The opening angle should be in the interval π/3≤θ≤5π/3. We choose λ=0.999, θ=π/2, npan=10, and  e=1. The quantity q converges initially as nsub is increased, but for  nsub>44 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 L∞ behavior of the solution ρ, while our ρ is unbounded. See, further, Appendix D.Convergence forq 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 q=1.1300163213105365. There are np=160+32nsub unknowns in the main linear system. (a) Convergence with nsub. (b) The number of iterations needed to meet an estimated relative residual of  ϵmach.
(a)
(b)


## 5. Compressed Inverse Preconditioning

Let us split the matricesKcoa and Kfin of (7) and (8) into two parts each
(11)Kcoa=Kcoa⋆+Kcoa∘,(12)Kfin=Kfin⋆+Kfin∘.
Here, the superscript ⋆ indicates that only entries of a matrix Kij whose indices i and j correspond to points zi and zj that both belong to the boundary segment Γ⋆ are retained. The remaining entries are zero.Now, we introduce two diagonal matricesWcoa and Wfin which have the quadrature weights wi on the diagonal. Furthermore, we need a prolongation matrix P which interpolates functions known at points on the coarse grid to points on the fine grid. The construction of P relies on panelwise 15-degree polynomial interpolation in parameter using Vandermonde matrices. We also construct a weighted prolongation matrix  PW via
(13)PW=WfinPWcoa-1.
The matrices P and PW share the same sparsity pattern. They are rectangular matrices, similar to the the identity matrix, but with one full (4+2nsub)16×64 block. Let superscript  T denote the transpose. Then,
(14)PWTP=Icoa
holds exactly. See Appendix A and [4, Section 4.3].Equipped withP and PW, 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
(15)(Ifin+λKfin⋆)ρfin=Pρ~coa.
Here, ρ~coa is the discretization of a piecewise smooth transformed density. The compression also uses the low-rank decomposition
(16)Kfin∘=PKcoa∘PWT,
which should hold to about machine precision.The compressed version of (8) reads
(17)(Icoa+λKcoa∘R)ρ~coa=λgcoa,
where the compressed weighted inverse  R is given by
(18)R=PWT(Ifin+λKfin⋆)-1P.
See Appendix B for details on the derivation. The compressed weighted inverse R, for Γ of (1), is a block diagonal matrix with one full  64×64  block and the remaining entries coinciding with those of the identity matrix.After having solved (17) for ρ~coa, the density ρfin can easily be reconstructed from ρ~coa in a postprocessor; see Section 9. It is important to observe, however, that ρfin is not always needed. For example, the quantity q of (6) can be computed directly from ρ~coa. Let ζcoa be a column vector which contains the absolute values of the boundary derivatives zi′ multiplied with weights wi, positions zi, and the complex scalar e-. Then,
(19)q=ℜ{ζcoa}TRρ~coa.


## 6. The Recursion for R

The compressed weighted inverseR is costly to compute from its definition (18). As we saw in Section 4, the inversion of large matrices (I+K) on highly refined grids could also be unstable. Fortunately, the computation of R 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

LetPbc 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 PWbc 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 [-1,1], then Pbc and  PWbc can be constructed as T32=[T16-1;T16+1]/2W32=[W16;W16]/2A=ones(16)AA=ones(32,16)fork=2:16A(:,k)=A(:,k-1).*T16AA(:,k)=AA(:,k-1).*T32endIP=AA/AIPW=IP.*(W32*(1./W16)′)%Pbc=zeros(96,64)Pbc(1:16,1:16)=eye(16)Pbc(17:48,17:32)=IPPbc(49:96,33:64)=flipud(fliplr(Pbc(1:48,1:32)))%PWbc=zeros(96,64)PWbc(1:16,1:16)=eye(16)PWbc(17:48,17:32)=IPWPWbc(49:96,33:64)=flipud(fliplr(PWbc(1:48,1:32))).(a) The prolongation operatorPbc performs panelwise interpolation from a grid on a four-panel mesh to a grid on a six-panel mesh. (b) The sparsity pattern of  Pbc.
(a)
(b)See [23, Appendix A] 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Γi⋆, i=1,2,…,nsub, be a sequence of subsets of Γ⋆ with Γi-1⋆⊂Γi⋆ and Γnsub⋆=Γ⋆. Let there also be a six-panel mesh and a corresponding 96-point grid on each Γi⋆. The construction of the subsets and their meshes should be such that if z(s), s∈[-2,2], is a local parameterization of Γi⋆, then the breakpoints (locations of panel endpoints) of its mesh are at s∈{-2,-1,-0.5,0,0.5,1,2}, and the breakpoints of the mesh on Γi-1⋆ are at s={-1,-0.5,-0.25,0,-0.25,0.5,1}. We denote this type of nested six-panel meshes type b. The index i is the level. An example of a sequence of subsets and meshes on Γ⋆ is shown in Figure 5 for nsub=3. Compare [3, Figure 2] and [4, Figure 5.1].Figure 5
The boundary subsetsΓ3⋆, Γ2⋆, and Γ1⋆ along with their corresponding type b meshes for nsub=3.LetKib denote the discretization of K on a type b mesh on Γi⋆. In the spirit of (11) and (12), we write
(20)Kib=Kib⋆+Kib∘,
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, letRnsub denote the full 64×64 diagonal block of R. The recursion for Rnsub is derived in Appendix C, and it reads
(21)Ri=PWbcT(𝔽{Ri-1-1}+Ib∘+λKib∘)-1Pbc,i=1,…,nsub,(22)𝔽{R0-1}=Ib⋆+λK1b⋆,
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
(23)R1=PWbcT(Ib+λK1b)-1Pbc.The programdemo2.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 nsub=40. This supports the claim of Section 5 that the discretization error in the solution is unaffected by compression.Figure6 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 nsub of recursion steps (levels) used in (21) corresponds to the number of subdivisions nsub used to construct the fine mesh.Same as Figure3, but using (17) and the program  demo2b.m (a loop of  demo2.m). There are only  np=160  unknowns in the main linear system.
(a)
(b)


## 6.1. Basic Prolongation Matrices

LetPbc 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 PWbc 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 [-1,1], then Pbc and  PWbc can be constructed as T32=[T16-1;T16+1]/2W32=[W16;W16]/2A=ones(16)AA=ones(32,16)fork=2:16A(:,k)=A(:,k-1).*T16AA(:,k)=AA(:,k-1).*T32endIP=AA/AIPW=IP.*(W32*(1./W16)′)%Pbc=zeros(96,64)Pbc(1:16,1:16)=eye(16)Pbc(17:48,17:32)=IPPbc(49:96,33:64)=flipud(fliplr(Pbc(1:48,1:32)))%PWbc=zeros(96,64)PWbc(1:16,1:16)=eye(16)PWbc(17:48,17:32)=IPWPWbc(49:96,33:64)=flipud(fliplr(PWbc(1:48,1:32))).(a) The prolongation operatorPbc performs panelwise interpolation from a grid on a four-panel mesh to a grid on a six-panel mesh. (b) The sparsity pattern of  Pbc.
(a)
(b)See [23, Appendix A] 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Γi⋆, i=1,2,…,nsub, be a sequence of subsets of Γ⋆ with Γi-1⋆⊂Γi⋆ and Γnsub⋆=Γ⋆. Let there also be a six-panel mesh and a corresponding 96-point grid on each Γi⋆. The construction of the subsets and their meshes should be such that if z(s), s∈[-2,2], is a local parameterization of Γi⋆, then the breakpoints (locations of panel endpoints) of its mesh are at s∈{-2,-1,-0.5,0,0.5,1,2}, and the breakpoints of the mesh on Γi-1⋆ are at s={-1,-0.5,-0.25,0,-0.25,0.5,1}. We denote this type of nested six-panel meshes type b. The index i is the level. An example of a sequence of subsets and meshes on Γ⋆ is shown in Figure 5 for nsub=3. Compare [3, Figure 2] and [4, Figure 5.1].Figure 5
The boundary subsetsΓ3⋆, Γ2⋆, and Γ1⋆ along with their corresponding type b meshes for nsub=3.LetKib denote the discretization of K on a type b mesh on Γi⋆. In the spirit of (11) and (12), we write
(20)Kib=Kib⋆+Kib∘,
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, letRnsub denote the full 64×64 diagonal block of R. The recursion for Rnsub is derived in Appendix C, and it reads
(21)Ri=PWbcT(𝔽{Ri-1-1}+Ib∘+λKib∘)-1Pbc,i=1,…,nsub,(22)𝔽{R0-1}=Ib⋆+λK1b⋆,
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
(23)R1=PWbcT(Ib+λK1b)-1Pbc.The programdemo2.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 nsub=40. This supports the claim of Section 5 that the discretization error in the solution is unaffected by compression.Figure6 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 nsub of recursion steps (levels) used in (21) corresponds to the number of subdivisions nsub used to construct the fine mesh.Same as Figure3, but using (17) and the program  demo2b.m (a loop of  demo2.m). There are only  np=160  unknowns in the main linear system.
(a)
(b)


## 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 B] as (24)[PW⋆T00I][A-1UVD]-1[P⋆00I]=[PW⋆TAP⋆+PW⋆TAU(D-VAU)-1VAP⋆-PW⋆TAU(D-VAU)-1-(D-VAU)-1VAP⋆(D-VAU)-1],where A plays the role of  Ri-1, P⋆ and  PW⋆ are submatrices of  Pbc and  PWbc, and  U, V, and D refer to blocks of  λKib∘.The programdemo3.m is based on  demo2.m, but has (24) incorporated. Besides, the integral equation (3) is replaced with
(25)ρ(r)+2λ∫Γ∂∂νrG(r,r′)ρ(r′)dσr′+∫Γρ(r′)dσr′=2λ(e·νr),r∈Γ,
which has the same solution ρ(r) 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 λKcoa∘ of (17).The execution ofdemo3.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  nsub>60.Same as Figure6, but the program  demo3b.m is used.
(a)
(b)


## 8. Various Useful Quantities

Let us introduce a new discrete densityρ^coa   via
(26)ρ^coa=Rρ~coa.
Rewriting (17) in terms of  ρ^coa gives
(27)(R-1+λKcoa∘)ρ^coa=λgcoa,
which resembles the original equation (7). We see that Kcoa∘, which is discretized using Gauss-Legendre quadrature, acts on ρ^coa. Therefore, one can interpret ρ^coa as pointwise values of the original density ρ(r), multiplied with weight corrections suitable for integration against polynomials. We refer to ρ^coa as a weight-corrected density. Compare [6, Section 5.4].Assume now that there is a square matrixS  which maps  ρ~coa to discrete values  ρcoa  of the original density on the coarse grid
(28)ρcoa=Sρ~coa.
The matrix S allows us to rewrite (17) as a system for the original density
(29)(S-1+λKcoa∘RS-1)ρcoa=λgcoa.
We can interpret the composition RS-1 as a matrix of multiplicative weight corrections that compensate for the singular behavior of ρ(r) on Γ⋆ when Gauss-Legendre quadrature is used.LetY denote the rectangular matrix
(30)Y=(Ifin+λKfin⋆)-1P,
and let Q 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 Y is the mapping from ρ~coa to ρfin. Therefore the columns of Y can be interpreted as discrete basis functions for ρ(r). It holds by definition that
(31)QP=Icoa,(32)QY=S.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  S will be presented in Section 10.


## 9. Reconstruction ofρfin from ρ~coa

The action ofY on  ρ~coa,  which gives ρfin, 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
(33)ρ→coa,i=[Ib-λKib∘(𝔽{Ri-1-1}+Ib∘+λKib∘)-1]Pbcρ~coa,i,Ib-λKib∘(𝔽{Ri-1-1}+Ib∘+λKib∘)-1i=nsub,…,1.
Here, ρ~coa,i is a column vector with 64 elements. In particular, ρ~coa,nsub is the restriction of ρ~coa to Γ⋆, while ρ~coa,i are taken as elements {17:80} of ρ→coa,i+1 for i<nsub. The elements {1:16} and {81:96} of ρ→coa,i are the reconstructed values of ρfin on the outermost panels of a type b mesh on Γi⋆. Outside of Γ⋆, ρfin coincides with ρ~coa.When the recursion is completed, the reconstructed values ofρfin on the four innermost panels are obtained from
(34)R0ρ~coa,0.
Should one wish to interrupt the recursion (33) prematurely, say at step  i=j, then
(35)Rj-1ρ~coa,(j-1)
gives values of a weight-corrected density on the four innermost panels of a type b mesh on Γj⋆. That is, we have a part-way reconstructed weight-corrected density ρ^part  on a mesh that is nsub-j+1 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 matricesKib∘ and Ri in the forward recursion (21) and reuse them in the backward recursion (33). Otherwise, they may be computed afresh.The programdemo4.m builds on the program  demo3.m, using (17) for (25). After the main linear system is solved for ρ~coa, a postprocessor reconstructs ρfin via (33). Then, a comparison is made with a solution ρfin obtained by solving the uncompressed system (8). Figure 8 shows that for nsub<10, the results are virtually identical. This verifies the correctness of (33). For  nsub>10, the results start to deviate. That illustrates the instabilities associated with solving (8) on a highly refined mesh. Compare Figure 3.Output fromdemo4.m and  demo5.m. (a) A comparison of ρfin  from the unstable equation (8) and  ρfin reconstructed from  ρ~coa  of  (17) via (33). (b) Relative accuracy in  q  of (6) from part-way reconstructed solutions  ρ^part.
(a)
(b)The programdemo5.m investigates the effects of premature interruption of (33). The number of recursion steps is set to nsub=100, and the recursion is interrupted at different levels. The density ρfin 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, q of (6) is computed from this part-way reconstructed solution. Figure 8(b) shows that the quality of q is unaffected by the level of interruption.


## 10. The Construction ofS

This section discusses the construction ofS and other auxiliary matrices. Note that in many applications, these matrices are not needed.The entries of the matricesP, PW, Q, R, S, and Y can only differ from those of the identity matrix when both indices correspond to discretization points on  Γ⋆. For example, the entries of  R only differ from the identity matrix for the 64×64 block denoted Rnsub in (21). In accordance with this notation, we introduce Pnsub, PWnsub, Qnsub, Snsub, and Ynsub for the restriction of P, PW, Q, S, and Y to Γ⋆. In the codes of this section, we often use this restricted type of matrices, leaving the identity part out.We observe thatSnsub is a square 64×64 matrix, Pnsub, PWnsub and Ynsub are rectangular 16(4+2nsub)×64 matrices, and Qnsub is a rectangular 64×16(4+2nsub) matrix. Furthermore, Qnsub is very sparse for large nsub. All columns of Qnsub with column indices corresponding to points on panels that result from more than eight subdivisions are identically zero.The programdemo6.m sets up Pnsub, PWnsub, and Qnsub, shows their sparsity patterns, and verifies the identities (14) and (31). The implementations for Pnsub and PWnsub rely on repeated interpolation from coarser to finer intermediate grids. The implementation of Qnsub 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 Pnsub and PWnsub are never needed in applications.(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.
(a)
(b)We are now ready to constructS. Section 9 presented a scheme for evaluating the action of Ynsub on discrete functions on the coarse grid on Γ⋆. The matrix  Ynsub, itself, can be constructed by applying this scheme to a 64×64 identity matrix. The matrix Qnsub was set up in  demo6.m. Composing these two matrices gives Snsub; see (32). This is done in the program  demo7.m, where the identity part is added as to get the entire matrix S. In previous work on RCIP, we have found use for S in complex situations where (29) is preferable over (17); see [13, Section 9]. If one merely needs  ρcoa from ρ~coa in a postprocessor, setting up S and using (28) are not worthwhile. It is cheaper to let Y act on ρ~coa and then let Q act on the resulting vector. Anyhow,  demo7.m builds on  demo4.m and gives as output ρcoa computed via (28); see Figure 9(b). For comparison,  ρfin, computed from the uncompressed system (8), is also shown.


## 11. Initiating R Using Fixed-Point Iteration

It often happens thatΓi⋆ is wedge-like. A corner of a polygon, for example, has wedge-like Γi⋆ at all levels. If Γ⋆ is merely piecewise smooth, then the Γi⋆ are wedge-like to double precision accuracy for i≪nsub.Wedge-like sequences ofΓi⋆ open up for simplifications and speedup in the recursion of (21) and (22). Particularly so if the kernel of the integral operator K of (5) is scale invariant on wedges. Then the matrix  Kib∘ becomes independent of i. It can be denoted Kb∘ and needs only to be constructed once. Furthermore, the recursion of (21) and (22) assumes the form of a fixed-point iteration
(36)Ri=PWbcT(𝔽{Ri-1-1}+Ib∘+λKb∘)-1Pbc,i=1,…(37)𝔽{R0-1}=Ib⋆+λKb⋆.
The iteration (36) can be run until  Ri converges. Let R* be the converged result. One needs not worry about predicting the number of levels needed. Choosing the number nsub of levels needed, in order to meet a beforehand given tolerance in R*, 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 ofK is scale invariant on wedges. If all Γi⋆ 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 nsub of levels needed for full convergence is halved compared to when using the initializer (22).Same as Figure7, but the program  demo8b.m is used.
(a)
(b)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  Ri and Kib∘ needed in (21) and in (33) is reduced as the recursions are shortened. This means savings in storage. Second, the number nsub 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ρ(r) 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  2·105 steps are needed.Fortunately, (36) can be cast as a nonlinear matrix equation
(38)G(R*)≡PWbcTA(R*)Pbc-R*=0,
where R*, as in Section 11, is the fixed-point solution and
(39)A(R*)=(𝔽{R*-1}+Ib∘+λKb∘)-1.
The nonlinear equation (38), in turn, can be solved for R* with a variant of Newton’s method. Let X be a matrix-valued perturbation of  R*, and expand  G(R*+X)=0 to first order in X. This gives a Sylvester-type matrix equation
(40)X-PWbcTA(R*)𝔽{R*-1XR*-1}A(R*)Pbc=G(R*)
for the Newton update X. 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].Figure11 shows a comparison between the fixed-point iteration of (36) and (37) and Newton’s method for computing the fixed-point solution  R*  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  θ=π/2, the integral operator K is the same as in (5), and λ=0.999. The relative difference between the two converged solutions R* is 5.6·10-16. 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 (Ifin+λKfin⋆), 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 programdemo9.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 U(r) at a point r in the plane is related to the layer density ρ(r) via
(41)U(r)=(e·r)-∫ΓG(r,r′)ρ(r′)dσr′;
see [5, Section 2.1].Figure12 shows how the field solution U(r) of (41) converges with mesh refinement at a point r=(0.4,0.1) inside the contour Γ of (1). We see that the accuracy in U(r), typically, is one digit better than the accuracy of the dipole moment q 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 Figure7(a) but with the program  demo3c.m. The potential U(r) at r=(0.4,0.1) is evaluated via (41), and  npan = 14 is used.


## 14. Composed Integral Operators

Assume that we have a modification of (5) which reads
(42)(I+MK)ρ1(z)=g(z),z∈Γ.
Here, K and g are as in (5), M is a new, bounded integral operator, and ρ1 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 ρ2(z)=Kρ1(z). Then,
(43)ρ1(z)+Mρ2(z)=g(z),-Kρ1(z)+ρ2(z)=0,
and after discretization on the fine mesh,
(44)([Ifin0fin0finIfin]+[0finMfin-Kfin0fin])[ρ1finρ2fin]=[gfin0].
Standard RCIP gives
(45)([Icoa0coa0coaIcoa]+[0coaMcoa∘-Kcoa∘0coa][R1R3R2R4])[ρ~1coaρ~2coa]=[gcoa0],
where the compressed inverse  R  is partitioned into four equi-sized blocks.Now, we replaceρ~1coa and  ρ~2coa  with a single unknown ρ~coa via
(46)ρ~1coa=ρ~coa-R1-1R3Kcoa∘R1ρ~coa,ρ~2coa=Kcoa∘R1ρ~coa.
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
(47)[Icoa+Mcoa∘(R4-R2R1-1R3)Kcoa∘R1+Mcoa∘R2-R1-1R3Kcoa∘R1]ρ~coa=gcoa.When (47) has been solved for ρ~coa, the weight-corrected version of the original density  ρ1 can be recovered as
(48)ρ^1coa=R1ρ~coa.Figure13 shows results for (42) with M being the double layer potential
(49)Mρ(z)=-2∫Γ∂∂νr′G(r,r′)ρ(r′)dσr′=1π∫Γρ(τ)ℑ{dττ-z}.
Figure 13(a) shows the convergence of q of (6) with nsub using the inner product preserving discretization scheme of Appendix D for (42) as implemented in  demo10.m. Figure 13(b) shows q produced with RCIP according to (47) and (48) as implemented in  demo10b.m. The reference value for q 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 q that is slightly wrong.Convergence forq of (6) with  ρ=ρ1 from (42). The curve Γ is as in (1), and  theta = pi/2,  npan = 11, and  evec = 1. The reference value is taken as q=1.95243329423584. (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.
(a)
(b)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 Section4 discretizes (5) using composite 16-point Gauss-Legendre quadrature. This works well as long as the kernel of the integral operator K 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 Sections16–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

LetD be the domain enclosed by the curve Γ, and let E be the exterior to the closure of  D. The exterior Dirichlet problem for the Helmholtz equation
(50)ΔU(r)+ω2U(r)=0,r∈E,(51)limr∋E→r∘U(r)=g(r∘),r∘∈Γ,(52)lim|r|→∞|r|(∂∂|r|-iω)U(r)=0,
has a unique solution U(r) under mild assumptions on Γ and  g(r)  [29] and can be modeled using a combined integral representation [30, Chapter 3]
(53)U(r)=∫Γ∂∂νr′Φω(r,r′)ρ(r′)dσr′-iω2∫ΓΦω(r,r′)ρ(r′)dσr′,r∈ℝ2∖Γ,
where Φω(r,r′) is the fundamental solution to the Helmholtz equation in two dimensions
(54)Φω(r,r′)=i4H0(1)(ω|r-r′|).
Here, H0(1)(·) is a Hankel function of the first kind. Insertion of (53) into (51) gives the combined field integral equation
(55)(I+Kω-iω2Sω)ρ(r)=2g(r),r∈Γ,
where
(56)Kωρ(r)=2∫Γ∂∂νr′Φω(r,r′)ρ(r′)dσr′,(57)Sωρ(r)=2∫ΓΦω(r,r′)ρ(r′)dσr′.Figure14 shows the performance of RCIP applied to (55) for 1000 different values of ω∈[1,103]. The program  demo11.m is used. The boundary Γ is as in (1) with θ=π/2, and the boundary conditions are chosen as g(r)=H0(1)(r-r′) with r′=(0.3,0.1) inside Γ. The error in U(r) of (53) is evaluated at r=(-0.1,0.2) outside Γ. Since the magnitude of U(r) 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.The exterior Dirichlet problem for Helmholtz equation with RCIP applied to (55). The program  demo11.m is used with Γ as in (1), and θ=π/2. The boundary condition g(r) of (51) is generated by a point source at (0.3,0.1). (a) The absolute error in U(r) at r=(-0.1,0.2). (b) The number of GMRES iterations needed to meet an estimated relative residual of  ϵmach.
(a)
(b)


## 17. The Exterior Neumann Helmholtz Problem

The exterior Neumann problem for the Helmholtz equation(58)ΔU(r)+ω2U(r)=0,r∈E,(59)limr∋E→r∘∂U(r)∂νr=g(r∘),r∘∈Γ,(60)lim|r|→∞|r|(∂∂|r|-iω)U(r)=0,
has a unique solution U(r)  under mild assumptions on  Γ and  g(r) [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

LetKω′ be the adjoint to the double-layer integral operator  Kω of (56)
(61)Kω′ρ(r)=2∫Γ∂∂νrΦω(r,r′)ρ(r′)dσr′.
Insertion of the integral representation
(62)U(r)=∫ΓΦω(r,r′)ρ(r′)dσr′,r∈ℝ2∖Γ
into (59) gives the integral equation
(63)(I-Kω′)ρ(r)=-2g(r),r∈Γ.Figure15 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  I-Kω′ in (63) has a nontrivial nullspace. Recall that (55) is always uniquely solvable, while (63) is not.The exterior Neumann problem for Helmholtz equation with RCIP applied to (63). The program  demo12.m is used with Γ as in (1), and θ=π/2. The boundary condition g(r) of (59) is generated by a point source at (0.3,0.1). (a) The absolute error in U(r) at r=(-0.1,0.2). (b) The number of GMRES iterations needed to meet an estimated relative residual of  ϵmach.
(a)
(b)


### 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]
(64)U(r)=∫ΓΦω(r,r′)ρ(r′)dσr′+i∫Γ∂∂νr′Φω(r,r′)(Siωρ)(r′)dσr′,r∈ℝ2∖Γ,
which after insertion into (59) gives the integral equation
(65)(I-Kω′-iTωSiω)ρ(r)=-2g(r),r∈Γ,
where
(66)Tωρ(r)=2∂∂νr∫Γ∂∂νr′Φω(r,r′)ρ(r′)dσr′.The hypersingular operatorTω 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]
(67)Tωρ(r)=2ω2∫ΓΦω(r,r′)(νr·νr′)ρ(r′)dσr′+2ddσr∫ΓΦω(r,r′)dρ(r′)dσr′dσr′.
This makes it possible to write (65) in the form
(68)(I+A+B1B2+C1C2)ρ(r)=-2g(r),r∈Γ,
where A=-K′, B2=Siω, and the action of the operators B1, C1, and C2 is given by
(69)B1ρ(r)=-2iω2∫ΓΦω(r,r′)(νr·νr′)ρ(r′)dσr′,(70)C1ρ(r)=-2iddσr∫ΓΦω(r,r′)ρ(r′)dσr′,(71)C2ρ(r)=2ddσr∫ΓΦiω(r,r′)ρ(r′)dσr′.
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 R.Solving the problem in the example of Section17.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 ω.The same exterior Neumann problem for Helmholtz equation as in Figure15, but RCIP is now applied to (68). The program  demo13b.m is used.
(a)
(b)The programdemo13c.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.


## 17.1. An Analogy with the Standard Laplace Approach

LetKω′ be the adjoint to the double-layer integral operator  Kω of (56)
(61)Kω′ρ(r)=2∫Γ∂∂νrΦω(r,r′)ρ(r′)dσr′.
Insertion of the integral representation
(62)U(r)=∫ΓΦω(r,r′)ρ(r′)dσr′,r∈ℝ2∖Γ
into (59) gives the integral equation
(63)(I-Kω′)ρ(r)=-2g(r),r∈Γ.Figure15 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  I-Kω′ in (63) has a nontrivial nullspace. Recall that (55) is always uniquely solvable, while (63) is not.The exterior Neumann problem for Helmholtz equation with RCIP applied to (63). The program  demo12.m is used with Γ as in (1), and θ=π/2. The boundary condition g(r) of (59) is generated by a point source at (0.3,0.1). (a) The absolute error in U(r) at r=(-0.1,0.2). (b) The number of GMRES iterations needed to meet an estimated relative residual of  ϵmach.
(a)
(b)


## 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]
(64)U(r)=∫ΓΦω(r,r′)ρ(r′)dσr′+i∫Γ∂∂νr′Φω(r,r′)(Siωρ)(r′)dσr′,r∈ℝ2∖Γ,
which after insertion into (59) gives the integral equation
(65)(I-Kω′-iTωSiω)ρ(r)=-2g(r),r∈Γ,
where
(66)Tωρ(r)=2∂∂νr∫Γ∂∂νr′Φω(r,r′)ρ(r′)dσr′.The hypersingular operatorTω 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]
(67)Tωρ(r)=2ω2∫ΓΦω(r,r′)(νr·νr′)ρ(r′)dσr′+2ddσr∫ΓΦω(r,r′)dρ(r′)dσr′dσr′.
This makes it possible to write (65) in the form
(68)(I+A+B1B2+C1C2)ρ(r)=-2g(r),r∈Γ,
where A=-K′, B2=Siω, and the action of the operators B1, C1, and C2 is given by
(69)B1ρ(r)=-2iω2∫ΓΦω(r,r′)(νr·νr′)ρ(r′)dσr′,(70)C1ρ(r)=-2iddσr∫ΓΦω(r,r′)ρ(r′)dσr′,(71)C2ρ(r)=2ddσr∫ΓΦiω(r,r′)ρ(r′)dσr′.
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 R.Solving the problem in the example of Section17.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 ω.The same exterior Neumann problem for Helmholtz equation as in Figure15, but RCIP is now applied to (68). The program  demo13b.m is used.
(a)
(b)The programdemo13c.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 programdemo11b.m is a continuation of  demo11.m which, after solving (55) for ρ~coa with RCIP and forming  ρ^coa via (26), computes the solution U(r) via (53) using three slightly different discretizations.(i)
Whenr is away from  Γ, 16-point Gauss-Legendre quadrature is used in (53) on all quadrature panels.(ii)
Whenr 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 r, and local panelwise evaluation is used for panels close to r.(iii)
Whenr is close to a panel neighboring a corner, the density  ρ~coa is first used to reconstruct   ρ^part  according to Section 9. Then 16-point Gauss-Legendre quadrature is used in (53) on panels away from r, and local panelwise evaluation is used for panels close to  r.The first two discretizations only use the coarse grid onΓ. The third discretization needs a grid on a partially refined mesh on Γ.The programdemo13d.m is a continuation of  demo13b.m which, after solving (65) with RCIP as described in Section 17.2, computes the solution U(r) via (64) using the three discretizations of the previous paragraph.Figures17 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 solutionU(r) to the exterior Dirichlet Helmholtz problem of Section 16 evaluated at 62392 points on a cartesian grid using  demo11b.m with ω=10. The source at r′=(0.3,0.1) is shown as a blue star.Figure 18
Same as Figure17, but the exterior Neumann Helmholtz problem is solved using  demo13d.m. The accuracy is even higher than in Figure 17.



---

*Source: 938167-2013-04-16.xml*