Solving Fredholm second-kind integral equations with singular right-hand sides on non-smooth boundaries
SSolving Fredholm second-kind integral equations with singularright-hand sides on non-smooth boundaries
Johan Helsing a , Shidong Jiang b a Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden b Department of Mathematics Sciences, New Jersey Institute of Technology,Newark, New Jersey 07102, USA
Abstract
A numerical scheme is presented for the solution of Fredholm second-kind boundaryintegral equations with right-hand sides that are singular at a finite set of boundarypoints. The boundaries themselves may be non-smooth. The scheme, which builds onrecursively compressed inverse preconditioning (RCIP), is universal as it is independentof the nature of the singularities. Strong right-hand side singularities, such as 1 / | r | α with α close to 1, can be treated in full machine precision. Adaptive refinement is usedonly in the recursive construction of the preconditioner, leading to an optimal numberof discretization points and superior stability in the solve phase. The performance ofthe scheme is illustrated via several numerical examples, including an application to anintegral equation derived from the linearized BGKW kinetic equation for the steadyCouette flow. Keywords: integral equation method, singular right-hand side, non-smooth domain,RCIP method, linearized BGKW equation
1. Introduction
Fredholm second-kind integral equations (SKIEs) have become standard tools forsolving boundary value problems of elliptic partial differential equations [9, 18, 21, 25,27]. Advantages include dimensionality reduction in the solve phase, elimination ofthe need to impose artificial boundary conditions for exterior problems, easily achievedhigh-order discretization, and optimal complexity when coupled with fast algorithmssuch as the fast multipole method [12].The present work is about the numerical solution to SKIEs of the form( I + K ) ρ ( r ) = f ( r ) , r ∈ Γ . (1)Here r ∈ R is a point in the plane; I is the identity operator; ρ is an unknown layerdensity to be solved for; Γ is a piecewise smooth closed contour (boundary) with a Email addresses: [email protected] (Johan Helsing), [email protected] (Shidong Jiang) a r X i v : . [ m a t h . NA ] F e b arameterization r ( s ) and with a finite number of corners; K is an integral operator onΓ which is compact away from the corners; and f is a right-hand side which is singularat a finite number of boundary points, which may or may not coincide with the cornervertices, but otherwise smooth. The union of corner vertices and boundary points where f is singular is referred to as singular points . These points are denoted γ j , j = 1 , , . . . .The kernel of K is denoted K ( r, r (cid:48) ).We shall construct an efficient scheme for the numerical solution of (1). The difficultyin this undertaking is that the singularities in f and the non-compactness of K at the γ j may require a very large number of unknowns for the resolution of ρ . This, in turn,may lead to high computing costs and also to artificial ill-conditioning and reducedachievable precision in quantities computed from ρ .The nature of the singularity of f can be rather arbitrary in the applications weconsider. There is no analysis of f involved in our work and neither is there anyfurther analysis of K . We only assume that f can be evaluated everywhere at Γ exceptfor at the γ j . If, however, it is known whether the leading singular behavior of f ishomogeneous on Γ in the shape of a wedge and whether K is scale invariant on such Γ,that information can be used to further improve the performance of our scheme.We shall solve (1) numerically using Nystr¨om discretization based on underlyingcomposite 16-point Gauss-Legendre quadrature and the parameterization r ( s ) and thenaccelerate and stabilize the solution process using an extended version of the recursivelycompressed inverse preconditioning (RCIP) method [14]. The discretization is chieflydone on a coarse mesh with quadrature panels of approximately the same size. Thecoarse quadrature panels are chosen so that the following holds: all singular points γ j coincide with panel endpoints; for r on a panel close to γ j and r (cid:48) away from γ j , K ( r ( s ) , r (cid:48) ) is smooth; for r away from γ j and r (cid:48) on a panel close to γ j , K ( r, r (cid:48) ( s )) issmooth.The RCIP method assumes that the SKIE has a panelwise smooth right-hand sideand consists of the following steps: • Transform the SKIE at hand into a form where the layer density to be solved foris panelwise smooth. • Use Nystr¨om discretization to discretize the transformed SKIE on a grid on a finemesh, obtained from the coarse mesh by repeated subdivision of the panels closestto each γ j . • Compress the transformed and discretized SKIE so that it can be solved on a gridon the coarse mesh without the loss of information. This involves the use of aforward recursion. • Solve the compressed equation. • Reconstruct the solution to the original SKIE from the solution to the compressedequation. This involves the use of a backward recursion.2n this paper, we extend the RCIP method to treat singular right-hand sides. Themethod is universal since algorithmic steps are completely independent of the natureof the singularities in the right-hand side function. The only information needed isthe locations of singularities. Very strong singularities can be treated in full machineprecision. Indeed, we have studied singularities of the form 1 / | r | α for a wide range of α and our method works very well even for α = 1 + 0 . nearly singular right-hand sides close to acorner, which can be treated easily by the method developed in this paper.The paper is organized as follows. In Section 2, a new transformation is introduced totreat the singular right-hand side f . Sections 3, 4, 5, and 6 present detailed modificationsto the other steps in the list that follow from the introduction of the new transform.Numerical examples are presented in Section 7. Finally, we discuss the application ofour new method to the integral equation derived from the linearized Bhatnagar-Gross-Krook-Welander (BGKW) kinetic equation [5, 30] for the steady Couette flow.In order to keep the presentation concise, we concentrate on new development andrefer the reader to the recently updated compendium [14] for details on the RCIPmethod. This compendium, in turn, contains references to original journal papers.
2. Transformed equation with smooth density
This section reviews the original RCIP transformation, discusses its shortcomingsfor singular right hand sides f , and presents a new and better transformation. Wefrequently use the concept of panelwise smooth functions . By this we mean functionswhich can be well approximated by polynomials of degree 15 in s on individual quadra-ture panels. We also introduce the boundary subsets Γ j(cid:63) , which refer to the four panelsthat are closest to a point γ j (two on each side), and the boundary subsets Γ j(cid:63)(cid:63) , whichrefer to the two panels that are closest to a point γ j (one on each side). The original RCIP transformation for SKIEs of the form (1) assumes that f is apanelwise smooth function and relies on a kernel split K ( r, r (cid:48) ) = K (cid:63) ( r, r (cid:48) ) + K ◦ ( r, r (cid:48) ) , r, r (cid:48) ∈ Γ , (2)and a corresponding operator split K = K (cid:63) + K ◦ . (3)In (3), the operator K (cid:63) denotes the part of K that accounts for self-interaction close tothe corner vertices γ j , and K ◦ is the compact remainder. The split (2) is determined3rom a geometric criterion: if r and r (cid:48) both are in Γ j(cid:63) for some j , then K (cid:63) ( r, r (cid:48) ) = K ( r, r (cid:48) ). Otherwise K (cid:63) ( r, r (cid:48) ) is zero.The change of variables ρ ( r ) = ( I + K (cid:63) ) − ˜ ρ ( r ) , r ∈ Γ , (4)makes (1) with (3) assume the form( I + K ◦ ( I + K (cid:63) ) − ) ˜ ρ ( r ) = f ( r ) , r ∈ Γ . (5)When f is panelwise smooth, the transformed layer density ˜ ρ in (5) will, loosely speak-ing, also be panelwise smooth thanks to that the action of K ◦ on any function resultsin a panelwise smooth function. In particular, ˜ ρ will be smooth on Γ j(cid:63)(cid:63) . The panelwisesmoothness of ˜ ρ on Γ j(cid:63)(cid:63) , inherited from f , is the key property which makes the trans-formed equation (5) efficient for the original problem (1). The efficiency comes fromthat a panelwise smooth unknown is easy to resolve by panelwise polynomials. When the right-hand side f in (1) is not panelwise smooth, but has singularities atthe corner vertices, neither ˜ ρ in the transformed equation (5) is panelwise smooth. Inorder to fix this problem we now propose a new transformation of (1) which, in additionto the split of K , also splits the right-hand side f and the unknown ρ as f ( r ) = f (cid:63) ( r ) + f ◦ ( r ) , r ∈ Γ , (6) ρ ( r ) = v ( r ) + g ( r ) , r ∈ Γ . (7)Here f (cid:63) ( r ) = f ( r ) if r ∈ Γ j(cid:63) for some j . Otherwise f (cid:63) ( r ) is zero. The functions v and g are given by v ( r ) = ( I + K (cid:63) ) − ˜ v ( r ) , (8) g ( r ) = ( I + K (cid:63) ) − f (cid:63) ( r ) , (9)where ˜ v is a new unknown transformed layer density.Use of (3), (6), (7), (8), and (9) makes (1) assume the form( I + K ◦ ( I + K (cid:63) ) − )˜ v ( r ) = f ◦ ( r ) − K ◦ ( I + K (cid:63) ) − f (cid:63) ( r ) , r ∈ Γ . (10)One can see, in (10), that ˜ v is panelwise smooth on Γ j(cid:63)(cid:63) . This is so since the right-handside of (10) is smooth on Γ j(cid:63) .We remark that if Γ is smooth so that there are no corners, only singularities in f ,then the local transformation (8) is not needed and (10) reduces to( I + K ) v ( r ) = f ◦ ( r ) − K ◦ ( I + K (cid:63) ) − f (cid:63) ( r ) , r ∈ Γ . (11)One can also imagine mixed situations where (8) is used only at those γ j which corre-spond to corner vertices. 4 igure 1: A contour Γ with a corner at γ of opening angle θ = π/ . Left: A coarse mesh with ten quadraturepanels on Γ . A subset of Γ , called Γ (cid:63) , covers the four coarse panels closest to γ . Right: A fine mesh createdfrom the coarse mesh by subdividing the panels closest to γ a number n sub = 3 of times.
3. Discretization of (5) and (10)This section summarizes the modifications needed in the RCIP method in orderfor it to apply to (10) rather than to (5). We first show how RCIP is applied to (5),with detailed references to [14], and then state the modifications needed for (10). Forsimplicity of presentation it is from now on assumed that there is only one singular point,denoted γ , with an associated four-panel neighboring zone denoted Γ (cid:63) as illustrated inFigure 1. We let Γ (cid:63)(cid:63) refer to the two panels closest to γ (one on each side).The discretization of (5) takes place on two different meshes: on the coarse meshand on a fine mesh. The fine mesh is constructed from the coarse mesh by n sub timessubdividing the panels closest to γ . See the right image of Figure 1 for an example.Discretization points on the coarse mesh constitute the coarse grid . Points on the finemesh constitute the fine grid . The number of refinement levels, n sub , is chosen so thatthe operator ( I + K (cid:63) ) − in (5) is resolved to a desired precision.As mentioned in Section 1, our Nystr¨om discretization relies on composite 16-pointGauss–Legendre quadrature. This means that an integral (cid:90) Γ h ( r ) d (cid:96) , where h is a smooth function and d (cid:96) is an element of arc length, can be approximatedby a sum on the coarse grid (cid:90) Γ h ( r ) d (cid:96) ≈ (cid:88) j h ( r ( s coa j )) | ˙ r ( s coa j ) | w coa j . (12)Here ˙ r ( s ) = d r ( s ) / d s and w coa j are appropriately scaled Gauss–Legendre weights. Aformula, analogous to (12) but with subscripts coa replaced with subscripts fin , holdswhen h is a singular function that needs the fine grid for resolution. It is shown in [14, Appendix B] that the discretization of (5) on the fine grid followedby RCIP-style compression leads to the system( I coa + K ◦ coa R ) ˜ ρ coa = f coa . (13)5ere R is a block-diagonal matrix defined as [14, Eq. (26)] R = P T W ( I fin + K (cid:63) fin ) − P , (14)and the subscripts coa and fin indicate what type of mesh is used for discretization. Theprolongation matrix P interpolates piecewise polynomial functions known at the coarsegrid to the fine grid and the matrix P W is a weighted prolongation matrix [14, Eq. (21)].The interpolation is done with respect to the parameter s . The superscript T denotesthe transpose. The discretization and compression of (10) can be constructed in a manner com-pletely analogous to that of (5) and leads to the system( I coa + K ◦ coa R )˜ v coa = f ◦ coa − K ◦ coa R f f (cid:63) coa , (15)where R is as in (14) and R f = P T W ( I fin + K (cid:63) fin ) − P f . (16)The prolongation matrix P f interpolates f from the coarse grid to the fine grid andcan easily be constructed from computed values of f on the two grids respectively. Thematrix P f is block-diagonal with blocks being either identity matrices or, for entriescorresponding to discretization points on Γ (cid:63)(cid:63) , the rectangular rank-one matrix P (cid:63)(cid:63)f = 1 f (cid:63)(cid:63) Hcoa f (cid:63)(cid:63) coa f (cid:63)(cid:63) fin f (cid:63)(cid:63) Hcoa . (17)Here f (cid:63)(cid:63) fin is a column vector with values of f on the fine grid on Γ (cid:63)(cid:63) , f (cid:63)(cid:63) coa is a columnvector with values of f on the coarse grid on Γ (cid:63)(cid:63) , and H denotes the conjugate transpose.Clearly, P (cid:63)(cid:63)f f (cid:63)(cid:63) coa = f (cid:63)(cid:63) fin .
4. Forward recursion for R and R f The matrices R and R f of (14) and (16) differ from the identity matrix only inthat they each contain a non-trivial 64 ×
64 diagonal block, associated with entriescorresponding to discretization points on Γ (cid:63) . These diagonal blocks can be efficientlyconstructed via recursions relying on the discretization of K on small meshes on ahierarchy of boundary subsets Γ (cid:63)i around γ . Central to the recursions are two types ofoverlapping meshes called type b and type c . The type b mesh contains six quadraturepanels. The type c mesh contains four quadrature panels. See Figure 2 for an illustrationand [14, Section 7] for more information.The recursions for R and R f , given below, use local prolongation matrices P bc , P W bc , and P fi bc similar to the global matrices P , P W and P f of Section 3. The matrix P bc performs polynomial interpolation from a grid on a type c mesh to a grid on a type b mesh on the same level. The matrix P fi bc interpolates f from a grid on a type c ype b type c Figure 2:
Top row: Meshes of type b and type c on the boundary subset Γ (cid:63) The type b mesh has six panels.The type c mesh has four panels. Bottom row: the boundary subsets Γ (cid:63) = Γ (cid:63) , Γ (cid:63) , and Γ (cid:63) along with theircorresponding type b meshes for n sub = 3 . mesh to a grid on a type b mesh on the same level. The subscript i indicates that P fi bc depends on the hierarchical level in which it appears.The matrices P bc and P W bc are level independent. Their construction is describedin [14, Section 7.1]. The matrix P fi bc is constructed analogously to the matrix P f ofSection 3. R It is shown in [14, Appendix D] that the non-trivial diagonal 64 ×
64 block of R canbe obtained from the recursion R i = P TW bc (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − P bc , i = 1 , . . . , n sub , (18)with the initializer F { R − } = I (cid:63) b + K (cid:63) . (19)Here K i b is the discretization of K on a type b mesh on level i in the hierarchy of localmeshes around γ . The operator F {·} expands its matrix argument by zero-padding(adding a frame of zeros of width 16 around it). The superscripts (cid:63) and ◦ denote matrixsplits analogous to the split (3).The forward recursion (18) starts at the finest refinement level, i = 1, and ascendsthrough the hierarchy of levels until it reaches the coarsest level i = n sub . The matrix R n sub is equal to the non-trivial 64 ×
64 block of R . R f A recursion for the non-trivial 64 ×
64 block of R f can be derived in complete analogywith the derivation of (18). See Appendix A for the key steps of the derivation. The7esult is R fi = P TW bc (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − (cid:0) F { R − i − R f ( i − } + I ◦ b (cid:1) P fi bc ,i = 1 , . . . , n sub . (20)The recursion (20) can be initialized with R f = R (21)and run in tandem with (18). From (15) it is evident that the matrix R f acts only on one particular known vector,namely on f (cid:63) coa . Therefore, rather than first finding R f via (20) and then computingthe vector r (cid:63)f = R f f (cid:63) coa , (22)one can instead modify (20) with (21) so that it produces r (cid:63)f directly r (cid:63)fi = P TW bc (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − (cid:16) F { R − i − r (cid:63)f ( i − } + f ◦ i b (cid:17) ,i = 1 , . . . , n sub , (23) r (cid:63)f = R f f (cid:63) . (24)The recursion (23) with (24) is a bit faster than (20) with (21).It is also worth noting that all three recursions (18), (20), and (23) can be imple-mented without the explicit inversion of R i − . The key to this is to use the Schur–Banachiewicz inverse formula for partitioned matrices [16, Eq. (8)]. Avoiding inversionis particularly important when R i − is ill conditioned. Details on an implementationfor (18), free from inversion of R i − , are given in [14, Section 8] (see also Appendix B). A minor problem with the recursions (18), (20), and (23) is that it could be difficultto determine a suitable recursion length n sub a priori. A too large n sub leads to unnec-essary work. A too small n sub fails to resolve the problem. Should it, however, happenthat K is scale-invariant on wedges, there may be an easy way around this problem.The key observation is that for large n sub and at levels i such that n sub − i (cid:29)
1, thematrices K ◦ i b have often converged to a matrix K ◦ b (in double precision arithmetic) thatis independent of i . This, typically, happens for n sub − i >
60 and means that (18)assumes the form of a fixed-point iteration R ∗ i = P TW bc (cid:16) F { R − ∗ ( i − } + I ◦ b + K ◦ b (cid:17) − P bc , i = 1 , , . . . . (25)It also means that all R i in (18) are the same for n sub − i (cid:29)
1. In view of the aboveone can replace the initializer R of (19) with the fixed-point matrix R ∗ obtained by8unning (25) until convergence. With the choice R = R ∗ , it is enough to take n sub = 60steps in the recursion (18). See, further, the discussion in [14, Sections 12–13].In a procedure similar to that just described, it is also possible to replace the initial-izers R f and r (cid:63)f of (21) and (24) with more efficient initializers. The requirements are,in addition to that K is scale invariant on wedges, that the leading singular behavior of f is homogeneous on wedges. If this holds, and if n sub − i (cid:29)
1, then (20) assumes theform of a linear fixed-point iteration R f ∗ i = P TW bc (cid:0) F { R − ∗ } + I ◦ b + K ◦ b (cid:1) − (cid:0) F { R − ∗ R f ∗ ( i − } + I ◦ b (cid:1) P f bc ,i = 1 , , . . . . (26)The fixed-point matrix R f ∗ can be found with direct methods solving a Sylvester equa-tion. With the choices R f = R f ∗ and r (cid:63)f = R f ∗ f (cid:63) it is enough to take n sub = 60steps in the recursions (20) and (23).We remark that efficient initializers can be found also under more general conditionson K than scale invariant on wedges. See [15, Section 5.3], for an example.
5. Computing integrals of ρ Often, in applications, one is not primarily interested in the solution ρ to (1) initself. Rather, one is interested in computing functionals of ρ of the type q = (cid:90) Γ h ( r ) ρ ( r ) d (cid:96) , (27)where h ( r ) is a smooth function. The RCIP method offers an elegant way to do thisthat only involves quantities appearing in the compressed equations (13) and (15).Assume that f is smooth and consider (13). The quantity q of (27) can then be wellapproximated by the sum on the coarse grid q ≈ (cid:88) j h ( r ( s coa j )) ˆ ρ coa j | ˙ r ( s coa j ) | w coa j , (28)where ˆ ρ coa j are elements of the weight-corrected density vectorˆ ρ coa = R ˜ ρ coa . (29)See [14, Appendix C] for a proof.The situation for a singular f and (15) is completely analogous. The expression (28)holds with (29) replaced byˆ ρ coa = R ˜ v coa + r (cid:63)f . (30)9 . Backward recursion for the reconstruction of ρ When the compressed equation (15) has been solved for ˜ v coa one might also beinterested in the reconstruction of the discretized solution ρ fin = v fin + g fin (31)to the original SKIE (1), compare (7). Such a reconstruction can be achieved by, looselyspeaking, running the recursions (18) and (20) backward on Γ (cid:63) . Outside of Γ (cid:63) , the coarsegrid and the fine grid coincide and it holds that ρ = v = ˜ v , see (7), (8), and (9). v fin We first review the reconstruction of ρ fin from the solution ˜ ρ coa to (13). The mech-anism for this reconstruction was originally derived in [13, Section 7] and is also sum-marized in [14, Section 10]. The backward recursion reads (cid:126) ρ coa ,i = (cid:104) I b − K ◦ i b (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − (cid:105) P bc ˜ ρ coa ,i , i = n sub , . . . , . (32)Here ˜ ρ coa ,i is a column vector with 64 elements. In particular, ˜ ρ coa ,n sub is the restrictionof ˜ ρ coa to Γ (cid:63) , while ˜ ρ coa ,i are taken as elements {
17 : 80 } of (cid:126) ρ coa ,i +1 for i < n sub . Theelements { } and {
81 : 96 } of (cid:126) ρ coa ,i are the reconstructed values of ρ fin on theoutermost panels of a type b mesh on Γ (cid:63)i .When the recursion is completed, there are no values assigned to ρ fin at points onthe four innermost panels (on Γ (cid:63) closest to γ ) on the fine grid. Reconstructed weight-corrected values of ρ fin on these panels can then be used, rather than true values, andare obtained from R ˜ ρ coa , . (33)We now observe that the reconstruction of v fin from the solution ˜ v coa to (15) isidentical to the reconstruction just described. This is so since both ˜ ρ of (5) and ˜ v of (10) are panelwise smooth functions. g fin The backward recursion for g fin from f (cid:63) coa is analogous to (32), but the vector (cid:126) g coa ,i ,corresponding to (cid:126) ρ coa ,i in (32), needs a split in a singular and a panelwise smooth parton a type b mesh on each Γ (cid:63)i (cid:126) g coa ,i = (cid:126) g smocoa ,i + f i b . (34)The backward recursion can then be written (cid:126) g smocoa ,i = (cid:104) I b − K ◦ i b (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − (cid:105) P bc ˜ g smocoa ,i − K ◦ i b (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − (cid:16) F { R − i − r (cid:63)f ( i − } + f ◦ i b (cid:17) , i = n sub , . . . , . (35)10ere ˜ g smocoa ,i is a column vector with 64 elements. In particular, ˜ g smocoa ,n sub = while ˜ g smocoa ,i are taken as elements {
17 : 80 } of (cid:126) g coa ,i +1 for i < n sub . The elements { } and {
81 : 96 } of (cid:126) g coa ,i in (34) are the reconstructed values of g fin on the outermost panelsof a type b mesh on Γ (cid:63)i . The reconstructed weight-corrected values of g fin on the fourinnermost panels on the fine grid are obtained from R ˜ g smocoa , + r (cid:63)f . (36)
7. Numerical examples
We now demonstrate the efficiency of our numerical scheme for (1). The schemeconsists of the compressed equation (15), the recursions (18), (23), (32), and (35), andinitializers obtained from (25) and (26). The code is implemented in
Matlab , release2020b, and executed on a 64 bit Linux laptop with a 2.10GHz Intel i7-4600U CPU.The implementations are standard and rely on built-in functions such as dlyap (SLI-COT subroutine SB04QD), for the Sylvester equation. Large linear systems are solvedusing GMRES, incorporating a low-threshold stagnation avoiding technique applicableto systems coming from discretizations of Fredholm integral equations of the secondkind. The GMRES stopping criterion is set to machine epsilon in the estimated relativeresidual.
We solve the SKIE (1) with an integral operator K defined by its action on ρ as Kρ ( r ) = 2 λ (cid:90) Γ ∂G∂ν ( r, r (cid:48) ) ρ ( r (cid:48) ) d (cid:96) (cid:48) , r ∈ Γ . (37)Here λ is a parameter set to λ = 0 . ν ( r ) is the exterior unit normal at position r on Γ, ∂/∂ν = ν ( r ) · ∇ , G ( r, r (cid:48) ) is the fundamental solution to Laplace’s equation in the plane G ( r, r (cid:48) ) = − π log | r − r (cid:48) | , (38)and Γ is the closed contour with a corner at γ = 0 parameterized as r ( s ) = sin( πs ) (cos(( s − . θ ) , sin(( s − . θ )) , s ∈ [0 , , (39)where θ corresponds to the opening angle of the corner. With the particular choice (37)for K , the SKIE (1) can model a transmission problem for Laplace’s equation [14,Section 4].In our experiments we shall use a coarse mesh that is sufficiently refined as toresolve (1) away from γ , vary the number of distinct recursion steps n sub , and monitorthe convergence of the scalar quantity q of (27) with h ( r ) = 1. For comparison we11ompute q in two ways: first on the coarse grid via (28) and then with ˆ ρ coa from (30),then on the fine grid via q ≈ (cid:88) j ρ fin j | ˙ r ( s fin j ) | w fin j (40)and with ρ fin from (31). We start with an example where the corner opening angle is set to θ = π . Then Γof (39) degenerates into a circle with a circumference of π . The right-hand side of (1)is set to f ( r ) = 1 (cid:96) α ( r ) + 1( π − (cid:96) ( r )) α , r ∈ Γ , (41)where α (cid:54) = 1 is a parameter controlling the strength of the singularity at γ and (cid:96) ( r )is the distance from γ to r measured in arc length along Γ as the circle is traversedin a counterclockwise fashion. The solution ρ to (1) can in this example be computedanalytically ρ ( r ) = f ( r ) + 2 λπ − α (1 − α )(1 − λ ) , r ∈ Γ , (42)as can the quantity q : q = 2 π (1 − α ) (1 − α )(1 − λ ) . (43)Note that ρ ( r ) of (42) diverges everywhere on Γ as α →
1, so there is no solution for α = 1. Neither is there a finite limit value for q . If α is complex, however, there is alimit solution as (cid:60) e { α } → α . For α < q computed on the coarse grid ( q coa ) and on the finegrid ( q fin ) agree completely – indicating that the reconstruction procedure of Section 6is stable.Note, in Figure 3(d), that for α = 1 + 0 . f that is not even in L , the scheme loses only about two digits of accuracy.Note also that, thanks to the use of initializers, a number n sub = 60 of distinct recursionsteps is more than enough to make q converge (to the achievable precision), irrespectiveof the value of α . 12
20 40 60 80 100 12010 -15 -10 -5 (a) -15 -10 -5 (b) -15 -10 -5 (c) -15 -10 -5 (d) Figure 3:
Convergence of q with the number of distinct recursion levels n sub in the example of Section 7.1.2.The singularity strengths α of (41) are taken as: (a) α = 0 . ; (b) α = 0 . ; (c) α = 0 . ; (d) α = 1 + 0 . . We now consider an example where the corner opening angle in (39) is set to θ = π/ f ( r ) = 1 | r | α + log | r | , r ∈ Γ . (44)Figure 4, analogous to Figure 3, shows results obtained with 160 discretizationpoints on the coarse grid on Γ and for various singularity strengths α . In the absenceof an analytical expression for q we use, as a reference, the value of q coa obtained with n sub = 500. For example, with α = 0 .
94 this value is q ≈ . q ≈ . ,
088 discretizationpoints on the fine grid on Γ using the L -norm-preserving Nystr¨om discretization ofAskham and Greengard [4, Section 4], which extends the L -inner-product-preservingdiscretization of Bremer [6, Section 2], in combination with compensated summation [17,20]. Since the relative difference between these two reference values is on the order ofour own error estimate in Figure 4(b), we believe that our error estimates are reliablealso in Figure 4(c,d), where the norm-preserving discretization cannot be used due tomemory constraints.We add that the execution of our code is very rapid. The computing time, per datapoint, in Figure 4 varies with n sub and also differs between q coa and q fin , but it is muchless than a second for all the data points shown.13
20 40 60 80 100 12010 -15 -10 -5 (a) -15 -10 -5 (b) -15 -10 -5 (c) -15 -10 -5 (d) Figure 4:
Convergence of q with the number of distinct recursion levels n sub in the example of Section 7.1.3.The singularity strengths α of (44) are taken as: (a) α = 0 . ; (b) α = 0 . ; (c) α = 0 . ; (d) α = 1 + 0 . . The exterior Dirichlet problem for the Helmholtz equation has been considered in[14, Section 18] in detail. Here we use this problem to illustrate that its solution U ( r )can be found accurately in the entire computational domain also when the right-handside f ( r ) is singular or nearly singular. The integral representation of U ( r ) and theresulting boundary integral equation (a combined field integral equation) are identicalto those in [14, Section 18]. The boundary Γ is given by (39) with θ = π/
2. We considertwo cases(a)
Singular right-hand side
The exact solution is U ( r ) = H (1)0 ( ω | r | ) (45)with f ( r ) being the restriction of U ( r ) to Γ. Here ω is the wavenumber and H (1)0 is the first-kind Hankel function of order zero, that is, the field is generated by anacoustic monopole right at the corner.(b) Nearly singular right-hand side
The exact solution is U ( r ) = H (1)1 ( ω | r − r (cid:48) | ) x − x (cid:48) | r − r (cid:48) | , (46)14 -15.6-15.4-15.2-15-14.8-14.6-14.4-14.2-14 -15.5-15-14.5-14-13.5-13 Figure 5: log of absolute error in the solution U ( r ) to the exterior Dirichlet Helmholtz problem. The bluecurve is the boundary Γ . Left: U ( r ) is an acoustic monopole field with the source at r (cid:48) = (0 , shown as ared star. Right: U ( r ) is an acoustic dipole field with the source at r (cid:48) = (10 − , shown as a green star. with f ( r ) being the restriction of U ( r ) to Γ, and r (cid:48) = (10 − , ω = 10 using the method in thispaper, the singular point is taken as γ = 0 both for (a) and for (b). For both cases, thecoarse grid on Γ has 56 panels, i.e., 896 discretization points; the number of refinementlevels is set to n sub = 112; and the number of GMRES iterations needed is 18. Thefield is evaluated using the scheme in [14, Section 20]. The solve phase takes about 4seconds, while the field evaluation takes about 100 seconds, which can be accelerated viathe fast multipole method [8] when necessary. A Cartesian grid of 200 ×
200 equispacedpoints is placed on the rectangle [ − . , . × [ − . , .
53] and evaluations are carriedout at those 27,760 grid points that are in the exterior domain. For case (a), 5,088target points activate local panelwise evaluation for close panels; 920 target points closeto the corner vertex require that the solution ρ fin on the fine grid is reconstructed usingthe backward recursion in Section 6. For case (b), the numbers are 6,364 and 1,326,respectively. Figure 5 shows the absolute error in the numerical solution, demonstratingthat our scheme achieves high accuracy in the entire computational domain for both singular and nearly singular right-hand sides.
8. Application to the linearized BGKW equation for the Couette flow
In this section, we revisit the integral equation that is derived from the linearizedBGKW equation for the steady Couette flow [7, 19, 22] u ( x ) − k √ π (cid:90) / − / J − (cid:18) | x − y | k (cid:19) u ( y ) dy = f ( x ) , (47a) f ( x ) = 12 √ π (cid:20) J (cid:18) / − xk (cid:19) − J (cid:18) / xk (cid:19)(cid:21) , (47b)15 able 1: The velocity u at x = 1 at sample values of the Knudsen number k . k u (1) in [19] u (1) (current) Error0 .
003 4 .
978 915 352 789 693 · − .
978 915 352 789 727 · − . · − .
01 4 .
930 697 807 742 208 · − .
930 697 807 742 217 · − . · − .
03 4 .
800 058 682 766 829 · − .
800 058 682 766 835 · − . · − . .
412 246 409 722 421 · − .
412 246 409 722 422 · − . · − . .
672 125 695 500 504 · − .
672 125 695 500 505 · − . · − . .
518 613 399 894 732 · − .
518 613 399 894 737 · − . · − . .
852 462 993 740 218 · − .
852 462 993 740 218 · − . . .
504 282 444 992 075 · − .
504 282 444 992 075 · − . . .
126 351 880 294 592 · − .
126 351 880 294 592 · − . · − . .
171 689 613 521 435 · − .
171 689 613 521 426 · − . · − . .
292 211 299 328 497 · − .
292 211 299 328 490 · − . · − . .
381 357 342 231 838 · − .
381 357 342 231 841 · − . · − . .
343 072 948 081 877 · − .
343 072 948 081 874 · − . · − where the parameter k is the Knudsen number and J n is the n th order Abramowitzfunction defined by J n ( x ) = (cid:90) ∞ t n e − t − x/t dt, n ≥ − . (48)See, for example, [11] and references therein for the properties of Abramowitz functionsand an accurate numerical scheme for their evaluation. The BGKW equation has beenstudied in [19], where it is shown that the solution of (47a) contains singular terms( x ln x ) n for n ∈ N := { , , . . . } at the end points. It is known that the kernel function J − ( x ) has both absolute value and logarithmic singularities at x = 0, and that theright hand side function (47b) has x ln x singularity at the end points. Some benchmarkcalculation has been carried out in [19], using dyadic refinement towards the end pointsto treat the singularities of the solution and the right hand side function, and generalizedGaussian quadrature [23] to treat the kernel singularities.We have implemented the method of this paper, based on (10), to solve (47a).We note that both end points are singular points and there is only one side to eachsingular point, since we are dealing with an open arc instead of a closed curve. Thisleads to straightforward modifications in the method and its implementation. Someimplementation details are as follows. First, the kernel-split quadrature is applied totreat the kernel singularity: the correction to the logarithmic singularity was done in[13]; the correction on diagonal blocks for the absolute value singularity can be derivedeasily in an identical manner; the splitting of the kernel into various parts is done usingeither the series expansion of the Abramowitz function [1] or the Chebyshev expansionfor each part as in [24]. When the Knudsen number is low, the kernel is sharply peakedat the origin. A modified version of the upsampling scheme in [2] is used to resolve thesharp peak of the kernel accurately and efficiently. Here the modification is that the16 able 2: The stress P xy at sample values of the Knudsen number k . k P xy in [19] P xy (current) Error0 . − .
490 909 702 131 201 · − − .
490 909 702 131 177 · − . · − . − .
900 405 009 657 547 · − − .
900 405 009 657 522 · − . · − . − .
413 798 601 526 842 · − − .
413 798 601 526 840 · − . · − . − .
155 607 782 558 620 · − − .
155 607 782 558 618 · − . · − . − .
344 983 511 356 682 · − − .
344 983 511 356 700 · − . · − . − .
694 625 753 368 226 · − − .
694 625 753 368 226 · − . . − .
083 322 536 749 375 · − − .
083 322 536 749 375 · − . . − .
266 437 497 658 084 · − − .
266 437 497 658 084 · − . . − .
446 632 678 455 994 · − − .
446 632 678 455 994 · − . . − .
536 943 539 674 479 · − − .
536 943 539 674 479 · − . . − .
611 624 603 488 405 · − − .
611 624 603 488 405 · − . . − .
743 853 873 277 227 · − − .
743 853 873 277 228 · − . · − . − .
796 682 147 138 912 · − − .
796 682 147 138 912 · − . · − exact centering in [2] is not enforced for local adaptive panels for each target. Instead,the so-called level-restricted property [10], i.e., sizes of adjacent panels can differ atmost by a factor of 2, is used to ensure the accuracy in the calculation of the integrals.Second, it is observed numerically that the condition number of the integral equation(47a) increases as k decreases, reaching about 2 · for k = 0 . u ( x ) approaches the asymptotic solution u asym ( x ) = x as k →
0. Thus, inorder to reduce the effect of the ill-conditioning of the integral equation for small valuesof k , we write u ( x ) = w ( x ) + x and solve the following equation for w ( x ) instead when k ≤ . w ( x ) − k √ π (cid:90) / − / J − (cid:18) | x − y | k (cid:19) w ( y ) dy = h ( x ) , (49a) h ( x ) = − k √ π (cid:20) J (cid:18) / − xk (cid:19) − J (cid:18) / xk (cid:19)(cid:21) . (49b)Since w ( x ) is a small perturbation when k is small, the inaccuracy in the calculationof w ( x ) has almost no effect in the overall accuracy of u ( x ). This allows us to achievenear machine precision for all physical quantities of interest for a much wider range ofthe Knudsen number k than reported in [19].We have repeated the calculations in [19] using the current method. As in [19],we have shifted the interval from [ − / , /
2] to [0 ,
1] in our calculation. For Knudsennumber k > .
1, eight panels are used for the whole interval [0 , k ≤ .
1, sixteenpanels are used. n sub is set to 40 for all values of the Knudsen number. Tables 1–3 listthe values of the velocity at x = 1, the stress P xy [19, Eq. (43)] and the half-channelmass flow rate Q [19, Eq. (42)] at sample values of the Knudsen number k , where thesecond column contains the values in [19], the third column contains the values using17 able 3: The half-channel mass flow rate Q at sample values of the Knudsen number k . k Q in [19] Q (current) Error0 .
003 1 .
242 445 655 299 172 · − .
242 445 655 299 170 · − . · − .
01 1 .
225 330 275 292 623 · − .
225 330 275 292 621 · − . · − .
03 1 .
180 147 037 188 893 · − .
180 147 037 188 892 · − . · − . .
057 028 408 172 292 · − .
057 028 408 172 292 · − . · − . .
560 111 699 820 618 · − .
560 111 699 820 629 · − . · − . .
804 708 735 555 459 · − .
804 708 735 555 461 · − . · − . .
281 659 776 113 917 · − .
281 659 776 113 916 · − . · − . .
489 298 506 190 833 · − .
489 298 506 190 833 · − . . .
627 042 060 967 383 · − .
627 042 060 967 383 · − . · − . .
147 460 412 330 841 · − .
147 460 412 330 841 · − . · − . .
714 449 048 590 649 · − .
714 449 048 590 649 · − . . .
043 009 085 700 258 · − .
043 009 085 700 261 · − . · − . .
226 757 181 742 400 · − .
226 757 181 742 400 · − . · − the current method, and the last column shows the relative difference between thesevalues. Here u (1) is obtained via the backward recursion in Section 6, P xy and Q arecalculated using the velocity on the coarse grid with special quadrature applied to treatthe logarithmic singularity at x = 1 / P xy . It is clear that the numerical resultsagree with those in [19] to near machine precision for all sample values of k .The current method is much more efficient as compared with the one used in [19].For k = 0 . . k = 10,the computation takes about 0 .
01 second of CPU time or less. In [19], the timing re-sults are 62 . . R also reduces the number of GMRES iterations. Indeed, for k = 0 . k = 10,only 7 iterations are needed to reach full double precision. All of these factors result insignificant reduction in the computational cost, i.e., a speedup of more than 300.
9. Conclusions
Finding efficient solvers for SKIEs on non-smooth boundaries is a field with substan-tial recent activity. However, almost all existing methods only work well for smoothright hand sides and may also involve a fair amount of operator-specific analysis andprecomputed quantities. In contrast, the method constructed in this paper for singu-lar and nearly singular right-hand sides computes all intermediary quantities neededon-the-fly and only involves a bare minimum of analysis. The method is both fast and18ccurate and is expected to have broad impacts and applications to the field of integralequation methods and beyond.
Acknowledgments
J. Helsing was supported by the Swedish Research Council under contract 2015-03780. S. Jiang was supported by the United States National Science Foundation undergrant DMS-1720405, and by the Flatiron Institute, a division of the Simons Foundation.The authors would like to thank Anders Karlsson at Lund University, Alex Barnett andLeslie Greengard at the Flatiron Institute for helpful discussions.
Appendix A. Derivation of the forward recursion formula (20) for R f First, we prove the following lemma.
Lemma 1.
Suppose that A is an invertible m × m matrix, B is an n × n matrix, U isan n × m matrix, and V is an m × n matrix with m ≥ n . Suppose further that UA − V and A + VBU are both invertible. Then U ( A + VBU ) − V = (cid:0) ( UA − V ) − + B (cid:1) − , (A.1) U ( A + VBU ) − = (cid:0) ( UA − V ) − + B (cid:1) − ( UA − V ) − UA − . (A.2) Proof.
Introduce a complex parameter λ and consider ( A + λ VBU ) − . When λ issufficiently small, the following Taylor expansion is valid( A + λ VBU ) − = A − − λ A − VBUA − + λ A − VBUA − VBUA − − . . . , (A.3)Mutiplying both sides of (A.3) with U from the left and with V from the right andregrouping, we obtain U ( A + λ VBU ) − V = UA − V − λ (cid:0) UA − V (cid:1) B (cid:0) UA − V (cid:1) + λ (cid:0) UA − V (cid:1) B (cid:0) UA − V (cid:1) B (cid:0) UA − V (cid:1) − . . . = (cid:16)(cid:0) UA − V (cid:1) − + λ B (cid:17) − , (A.4)where the second equality follows from the application of the Taylor expansion in thereverse order. By the cofactor formula of the matrix inverse and the so-called bigformula for the matrix determinant [29], both sides of (A.4) are rational functions of λ .Since these two rational functions are equal in a small neighborhood of the origin in thecomplex plane, they must be equal everywhere in the whole complex plane by analyticcontinuation [3]. Setting λ = 1 in (A.4), we establish (A.2) and the invertibility of( UA − V ) − + B simultaneously. (A.2) can be proved in an almost identical manner.We now recall the definition of R f in (16) R f = P T W ( I fin + K (cid:63) fin ) − P f , (A.5)19 ype a type b type clevel 1level 2level 3 Figure A.6:
Meshes of type a, type b, and type c at different levels. The type a mesh contains i dyadicfine panels at level i . The type b mesh always contains six panels, and the type c mesh always contains fourpanels. At level , the type a mesh is identical to the type b mesh. At level n sub , the type c mesh containsfour coarse panels on Γ (cid:63) in Figure 1. where K (cid:63) fin is the system matrix built on a fine mesh that is obtained via n sub level ofdyadic refinement of the four coarsest panels (with two panels on each side) near thesingular point. Let n gl be the number of Gauss-Legendre nodes on each panel. Thenthe total number of discretization points on the fine mesh is n gl (4 + 2 n sub ). Clearly,direct application of (A.5) is very expensive, inaccurate, and unrobust. Instead, theforward recursion is used to compute R f . The forward recursion starts from the finestsix panels around the singular point at level i = 1, adds one panel on each side as thelevel goes up, and reaches the full fine mesh at level i = n sub . See Figure A.6 for anillustration of three different types of meshes at different levels, where the type a meshis needed only in the derivation of the forward recursion. The actual recursion formulainvolves only type b and type c meshes. To be more precise, at any step in the forwardrecursion, one only needs to build a system matrix of size 6 n gl × n gl , i.e., on the sixpanels of a type b mesh, and the type c mesh is used only implicitly in the constructionof the prolongation matrix P and its weighted version P W . Similar to [14, Eq. (D.1)],we define R fi := P TW i ac ( I i a + K i a ) − P fi ac , i = 1 , · · · , n sub . (A.6)By the action of P fi ac and P TW i ac , R fi is always a matrix of size 4 n gl × n gl at any level i . A derivation similar to that of [14, Eq. (D.6)] leads to R fi = P TW bc P TW i ab ( I i a + F { K ( i − } + P i ab K ◦ i b P TW i ab ) − P fi ab P fi bc , (A.7)where we have assumed that the low-rank property K ◦ i a = P i ab K ◦ i b P TW i ab (A.8)holds to machine precision, as in [14, Eq. (D.3)]. Since the type a mesh differs fromthe type b mesh only on the two panels (i.e., Γ (cid:63)(cid:63)i b ) closest to the singular point in the20 igure A.7: Nonzero patterns of F { R − i − } and F { R f ( i − } (left), K ◦ i b (center), and I ◦ b (right). Note thatthe patterns depend on the fact that when constructing the system matrix in the forward recursion, thesources and targets are arranged in the order from the top panel to the bottom panel on a type b mesh, seeFigure A.6. type b mesh, the only diagonal blocks in P i ab and P W i ab that are not identity matricesare P i ab (I (cid:63)(cid:63)i a , I (cid:63)(cid:63)i b ) and P W i ab (I (cid:63)(cid:63)i a , I (cid:63)(cid:63)i b ), respectively. Here I (cid:63)(cid:63)i a and I (cid:63)(cid:63)i b contain indicescorresponding to discretization points in Γ (cid:63)(cid:63)i b on type a and type b meshes, respectively.Now it is straightforward to verify that (A.8) is equivalent to K i a (I (cid:63)(cid:63)i a , I ◦ i b ) = P i ab (I (cid:63)(cid:63)i a , I (cid:63)(cid:63)i b ) K i b (I (cid:63)(cid:63)i b , I ◦ i b ) , K i a (I ◦ i b , I (cid:63)(cid:63)i a ) = K i b (I ◦ i b , I (cid:63)(cid:63)i b ) P TW i ab (I (cid:63)(cid:63)i b , I (cid:63)(cid:63)i a ) , (A.9)where I ◦ i b contains indices corresponding to discretization points on Γ ◦ i b , that is, onthe two outermost panels from γ . Since Γ (cid:63)(cid:63)i b is always well-separated from Γ ◦ i b for anylevel i , (A.9) is equivalent to stating that the kernel function K ( r, r (cid:48) ) is smooth when r ∈ Γ (cid:63)(cid:63)i b and r (cid:48) ∈ Γ ◦ i b , or vice versa, and thus the interaction between Γ (cid:63)(cid:63)i b and Γ ◦ i b can bediscretized to machine precision with n gl points on each panel in Γ (cid:63)(cid:63)i b provided that n gl isnot too small. This holds for all kernels we have encountered in practice, including, say,highly oscillatory ones, if the coarse mesh is chosen in such a way that the oscillations ofthe kernel are well-resolved. We emphasize that this is a property of the kernel function K ( r, r (cid:48) ) and the type a and b meshes. It has nothing to do with the singularity of theright-hand side f .Applying (A.2) to part of the right-hand side of (A.7), we obtain P TW i ab ( I i a + F { K ( i − } + P i ab K ◦ i b P TW i ab ) − = (cid:104)(cid:0) P TW i ab ( I i a + F { K ( i − } ) − P i ab (cid:1) − + K ◦ i b (cid:105) − · (cid:0) P TW i ab ( I i a + F { K ( i − } ) − P i ab (cid:1) − P TW i ab ( I i a + F { K ( i − } ) − . (A.10)Recall that [14, Eq. (D.10)] states that F { R i − } + I ◦ b = P TW i ab ( I i a + F { K ( i − } ) − P i ab . (A.11)We observe that F { R i − } places R i − in the center 4 n gl × n gl block with zero paddingin a 6 n gl × n gl matrix, and the nonzero blocks of I ◦ b are the identity matrix of size n gl in the top and bottom diagonal blocks (Figure A.7). Thus, F { R − i − } + I ◦ b = ( F { R i − } + I ◦ b ) − = (cid:0) P TW i ab ( I i a + F { K ( i − } ) − P i ab (cid:1) − . (A.12)21imilarly, F { R f ( i − } = F { P TW ( i − ( I ( i − + K ( i − ) − P f ( i − } = P TW i ab ( I i a + F { K ( i − } ) − P fi ab − I ◦ b , (A.13)where the second equality uses (cid:0) P TW i ab P fi ab (cid:1) ◦ = I ◦ b . Combining (A.7), (A.10)–(A.13),we obtain R fi = P TW bc (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − (cid:0) F { R − i − } + I ◦ b (cid:1) (cid:0) F { R f ( i − } + I ◦ b (cid:1) P fi bc . (A.14)Finally, we arrive at (20) by observing that (cid:0) F { R − i − } + I ◦ b (cid:1) (cid:0) F { R f ( i − } + I ◦ b (cid:1) = F { R − i − R f ( i − } + I ◦ b , (A.15)due to the nonzero patterns of the involved matrices shown in Figure A.7. Appendix B. A synopsis of the RCIP demo codes
On the website , maintainedby the first author, a set of
Matlab demo codes are posted so that researchers who areinterested in the RCIP method can download and run them and gain first-hand experi-ence about the method’s efficiency, accuracy, and robustness for a few chosen problems.The codes are written so that they can easily be modified to solve other problems.The reader is expected to read the tutorial [14] in tandem with looking at the codes.Many demo codes use the one-corner contour Γ in Figure 1, which is parameterized by s ∈ [0 ,
1] with s = 0 corresponding to the corner at the origin, and going back to thecorner again in the counterclockwise direction as s increases to s = 1. Most demo codesconsist of the following building blocks or functions: • Boundary discretization functions. These include the discretization in parame-ter space using panel division and functions returning z ( s ), z (cid:48) ( s ), z (cid:48)(cid:48) ( s ) in com-plex form for a given parameter value s , e.g., zfunc , zpfunc , zppfunc , zinit , panelinit . • Integral operator discretization functions. These include splitting of the kernelinto various parts – smooth part, logarithmically singular part, Cauchy singularpart, hypersingular part, etc., corrections for singular and nearly singular in-tegrals (explicit kernel-split quadrature), e.g.,
LogCinit , WfrakLinit , wLCinit , StargClose , KtargClose . • System matrix construction functions such as
MAinit , MRinit , Soperinit , Koperinit . • Functions for the forward recursion formula for computing the nontrivial block of R . These include zlocinit , Rcomp , SchurBana . • Other utility functions such as myGMRES , Tinit16 , Winit16 that are self-explanatory.22e now provide some comments on the forward recursion functions. In the dis-cussions below, n gl = 16, as in most demo codes. The function zlocinit returns adiscretization on the type b mesh at level i , where the points are always arranged inthe order from the top panel to the bottom panel. Another important feature is thatin zlocinit , the singular points γ j is translated to the origin to reduce the effect ofround-off error. See Figure B.9 for detailed comments on zlocinit .Recall from (18) that the forward recursion formula for R is R i = P TW bc (cid:0) F { R − i − } + I ◦ b + K ◦ i b (cid:1) − P bc , i = 1 , . . . , n sub . (B.1)Let M i := I b + K i b and M ◦ i := I ◦ b + K ◦ i b . Then Figure A.7 shows that there is nooverlap in nonzero entries of M ◦ i and R i − . Direct implementation of (B.1) is shown inAlgorithm 1. Algorithm 1
Forward recursion via direct implementation of (B.1). for i = 1 , . . . , n sub do Obtain 96 discretization points on the type b mesh at level i . Construct 96 ×
96 system matrix M i at level i . if i == 1 then R = M (cid:63) − . That is, R is the inverse of the center 64 × block of M . end if Compute R − i − . Replace the center 64 ×
64 block of M i with R − i − to obtain (cid:102) M i . Compute (cid:102) M − i . Compress (cid:102) M − i to obtain R i . To be more precise, R i = P TW bc (cid:102) M − i P bc . end for Return the nontrivial 64 ×
64 block of R , i.e., R n sub .Algorithm 1 needs to invert two matrices, i.e., R i − of size 64 ×
64 and (cid:102) M i ofsize 96 ×
96. This may lead to loss of accuracy and even instability when R i − is ill-conditioned. As pointed out in Section 4.3, the algorithm can be improved via the blockmatrix inversion formula that avoids explicit inversion of R i − . Indeed, after columnand row permutations, we obtain the following block structure (Figure B.8) F { R − i − } + I ◦ b + K ◦ i b −→ (cid:20) A − UV D (cid:21) . (B.2)However, there is no need to carry out column/row permutation explicitly. Mat-lab ’s column/row indexing conveniently achieves this purpose without actual permu-tation. To be more precise, A = R i − is a 64 ×
64 block, D = M ◦ i ( circL , circL ) is a32 ×
32 block, U = M ◦ i ( starL , circL ) is a 64 ×
32 block, and V = M ◦ i ( circL , starL ) is23 igure B.8: Column/row permutation converting F { R − i − } + I ◦ b + K ◦ i b into a × block matrix. Blue blockis R − i − and red block is I ◦ b + K ◦ i b . Left: the original matrix F { R − i − } + I ◦ b + K ◦ i b . Right: the matrix aftercolumn/row permutation. a 32 ×
64 block, where circL =[1:16 81:96] and starL =17:80. And [14, Eq. (31)] followsfrom the following block matrix inversion formula (cid:20) A − UV D (cid:21) − = (cid:20) A + AU ( D − VAU ) − VA − AU ( D − VAU ) − − ( D − VAU ) − VA ( D − VAU ) − (cid:21) , (B.3)which can be verified via direct matrix multiplication. Thus, in order to achieve betterefficiency and stability, Algorithm 1 is replaced by function Rcomp (Figure B.10) thatcontains the main loop for the forward recursion, and function
SchurBana (Figure B.11)that replaces steps 8–11 in Algorithm 1 by [14, Eq. (31)]. We note that one only needsto invert a 32 ×
32 matrix in function
SchurBana . References [1] M. Abramowitz. Evaluation of the integral (cid:82) ∞ e − u − x/u du . J. Math. Phys. Camb. ,32:188–192, 1953.[2] L. af Klinteberg, F. Fryklund, and A.-K. Tornberg. An adaptive kernel-splitquadrature method for parameter-dependent layer potentials. arXiv preprintarXiv:1906.07713 , 2019.[3] L. V. Ahlfors.
Complex analysis: an introduction to the theory of analytic functionsof one complex variable . McGraw-Hill, New York, 1966.[4] T. Askham and L. Greengard. Norm-preserving discretization of integral equationsfor elliptic PDEs with internal layers I: the one-dimensional case.
SIAM Rev. ,56(4):625–641, 2014.[5] P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collision processes ingases. I. Small amplitude processes in charged and neutral one-component systems.
Phys. Rev. , 94(3):511–525, 1954.[6] J. Bremer. On the Nystr¨om discretization of integral equations on planar curveswith corners.
Appl. Comput. Harmon. Anal. , 32(1):45–64, 2012.24 unction [z,zp,zpp,nz,w,wzp]=zlocinit(theta,T,W,nsub,i,npan)% This function returns to the user the discretization of% the type b mesh at the ith level%% Input parameters:% theta - the parameter for the one-corner curve% T, W - Gauss-Legendre nodes and weights on [-1,1]% nsub - the total number of dyadic refinement in the% calculation of the preconditioner R% i - the level index% npan - the total number of coarse panels%% Output parameters:% z,zp,zpp,nz,w,wzp - complex column vectors of length 96% z - coordinates of discretization points (in complex% number format) from the top panel to the bottom panel% zp,zpp - z'(s), z''(s)% nz - unit normal vector% w - quadrature weights% wzp - w*z'(s)% 1/npan is the length of the coarse panel in the parameter space% h is the length of each panel on the type c mesh at level ih=1/npan/2ˆ(nsub-i);% s returns discretization points in the parameter space for% the bottom 3 panels in the type b mesh at level i.% That is, s contains scaled and shifted Gauss-Legendre nodes on% the three panels [0, 0.5h], [0.5h, h], [h, 2h].s=[T/4+0.25;T/4+0.75;T/2+1.5]*h;w=[W/4;W/4;W/2]*h; w=[flipud(w);w];z=zfunc(s,theta);% conj(flipud(z)) produces points on the top three panels of% the type b mesh by symmetry.z=[conj(flipud(z));z];zp=zpfunc(s,theta);zp=[-conj(flipud(zp));zp];zpp=zppfunc(s,theta);zpp=[conj(flipud(zpp));zpp];nz=-1i*zp./abs(zp);wzp=w.*zp;
Figure B.9:
Matlab function zlocinit . [7] C. Cercignani. Rarefied Gas Dynamics: From Basic Concepts to Actual Calcula-tions . Cambridge University Press, Cambridge, UK, 2000.[8] H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin,N. Yarvin, and J. Zhao. Remarks on the implementation of the wideband FMMfor the Helmholtz equation in two dimensions. In
Inverse problems, multi-scale unction R=Rcomp(theta,lambda,T,W,Pbc,PWbc,nsub,npan)% Forward recursion for computing the nontrivial 64x64 block% of the preconditioner R% indices for the center four panels on the type b meshstarL=17:80;% indices for the top and bottom panels on the type b meshcircL=[1:16 81:96];% forward recursion - level 1 is the finest level,% level nsub goes back to four coarse panels near the cornerfor level=1:nsub% discretization on the type b mesh at the current level[z,zp,zpp,nz,w,wzp]=zlocinit(theta,T,W,nsub,level,npan);% construct the interaction matrix on the type b meshK=MAinit(z,zp,zpp,nz,w,wzp,96);MAT=eye(96)+lambda*K; % the full system matrix% at the finest level, R is initialized by simply inverting% the "bad" part of the system matrixif level==1R=inv(MAT(starL,starL));end% use block inversion formula to carry out the forward recursionR=SchurBana(Pbc,PWbc,MAT,R,starL,circL);end Figure B.10:
Matlab function
Rcomp . analysis and effective medium theory , volume 408 of Contemp. Math. , pages 99–110. Amer. Math. Soc., Providence, RI, 2006.[9] D. Colton and R. Kress.
Inverse Acoustic and Electromagnetic Scattering Theory .Springer, New York, NY, 2012.[10] F. Ethridge and L. Greengard. A new fast-multipole accelerated Poisson solver intwo dimensions.
SIAM J. Sci. Comput. , 23(3):741–760, 2001.[11] Z. Gimbutas, S. Jiang, and L.-S. Luo. Evaluation of Abramowitz functions in theright half of the complex plane.
J. Comput. Phys. , 405:109169, 2020.[12] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations.
J. Comput.Phys. , 73(2):325–348, 1987.[13] J. Helsing. Integral equation methods for elliptic problems with boundary condi-tions of mixed type.
J. Comput. Phys. , 228(23):8892–8907, 2009.[14] J. Helsing. Solving integral equations on piecewise smooth boundaries using theRCIP method: a tutorial. arXiv preprint arXiv:1207.6737v9 , 2018.26 unction A=SchurBana(P,PW,K,A,starL,circL)% This function uses the block matrix inversion formula% to compute the forward recursion for the nontrivial% 64x64 block of the preconditioner R.%% Input parameters:% P - 96x64 prolongation matrix% PW - 96x64 weighted prolongation matrix% A - 64x64 preconditioner matrix R { i-1 } at level i-1% starL - 17:80 center four panels on the type b mesh% circL - [1:16 81:96] top and bottom panels on the type b mesh%% Output parameters:% A - 64x64 preconditioner matrix R { i } at level istarS=17:48; % center two panels on the type c meshcircS=[1:16 49:64]; % top and bottom panels on the type c mesh% Using the notation in Eq. (31) in the RCIP tutorial, we have% V = K(circL, starL);% U = K(starL, circL);% D = K(circL, circL);% A = R { i-1 } ;VA=K(circL,starL)*A; % V*APTA=PW(starL,starS)'*A; % P WˆT * APTAU=PTA*K(starL,circL); % P WˆT * A * UDVAUI=inv(K(circL,circL)-VA*K(starL,circL)); % (D-V A U)ˆ { -1 } DVAUIVAP=DVAUI*(VA*P(starL,starS)); % (D-V A U)ˆ { -1 } * (V A P)A(starS,starS)=PTA*P(starL,starS)+PTAU*DVAUIVAP; % (1,1) blockA(circS,circS)=DVAUI; % (2,2) blockA(circS,starS)=-DVAUIVAP; % (2,1) blockA(starS,circS)=-PTAU*DVAUI; % (1,2) block Figure B.11:
Matlab function
SchurBana . [15] J. Helsing and S. Jiang. On integral equation methods for the first dirichlet prob-lem of the biharmonic and modified biharmonic equations in nonsmooth domains. SIAM J. Sci. Comput. , 40(4):A2609–A2630, 2018.[16] H. V. Henderson and S. R. Searle. On deriving the inverse of a sum of matrices.
SIAM Rev. , 23(1):53–60, 1981.[17] N. J. Higham.
Accuracy and stability of numerical algorithms . SIAM, 2002.[18] G. C. Hsiao and W. L. Wendland.
Boundary integral equations , volume 164 of
Applied Mathematical Sciences . Springer-Verlag, Berlin, 2008.[19] S. Jiang and L.-S. Luo. Analysis and solutions of the integral equation derived27rom the linearized BGKW equation for the steady Couette flow.
J. Comput.Phys. , 316:416–434, 2016.[20] W. Kahan. Further remarks on reducing truncation errors.
Commun. Assoc. Com-put. Mach. , 8:40, 1965.[21] R. Kress.
Linear Integral Equations , volume 82 of
Applied Mathematical Sciences .Springer–Verlag, Berlin, third edition, 2014.[22] W. Li, L.-S. Luo, and J. Shen. Accurate solution and approximations of the lin-earized BGK equation for steady Couette flow.
Comput. Fluids , 111:18–32, 2015.[23] J. Ma, V. Rokhlin, and S. Wandzura. Generalized gaussian quadrature rules forsystems of arbitrary functions.
SIAM J. Numer. Anal. , 33(3):971–996, 1996.[24] A. J. MacLeod. Chebyshev expansions for Abramowitz functions.
Appl. Numer.Math. , 10:129–137, 1992.[25] P.-G. Martinsson.
Fast direct solvers for elliptic PDEs , volume 96 of
CBMS-NSFRegional Conference Series in Applied Mathematics . Society for Industrial andApplied Mathematics (SIAM), Philadelphia, PA, 2020.[26] Q. Nie and F.-R. Tian. Singularities in Hele–Shaw flows driven by a multipole.
SIAM J. Appl. Math. , 62(2):385–406, 2001.[27] C. Pozrikidis.
Boundary integral and singularity methods for linearized viscousflow . Cambridge Texts in Applied Mathematics. Cambridge University Press, Cam-bridge, 1992.[28] T. S. Rappaport, G. R. MacCartney, S. Sun, H. Yan, and S. Deng. Small-scale,local area, and transitional millimeter wave propagation for 5G communications.
IEEE Trans. Antennas Propag. , 65(12):6474–6490, 2017.[29] G. Strang.
Introduction to linear algebra . Wellesley-Cambridge Press, Wellesley,MA, 5th edition, 2016.[30] P. Welander. On the temperature jump in a rarefied gas.