# Mikhlin’s Integral Equation and the Integral Equation with the Generalized Neumann Kernel on Simply Connected Domains

Authors: Samir Naqos; Ali H. M. Murid; Mohamed M. S. Nasser
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

Mikhlin’s integral equation is a classical integral equation for solving boundary value problems for Laplace’s equation. The kernel of the integral equation is known as the Neumann kernel. Recently, an integral equation for solving the Riemann–Hilbert problem was derived. The kernel of the new integral equation is a generalization of the Neumann kernel, and hence, it is called the generalized Neumann kernel. The objective of this paper is to present a detailed comparison between these two integral equations with emphasis on their similarities and differences. This comparison is done through applying both equations to solve Laplace’s equation with Dirichlet boundary conditions in simply connected domains with smooth and piecewise smooth boundaries.


---

## Body

## 1. Introduction

Fredholm integral equations of the second kind provide an important and useful tool for solving linear boundary value problems of the elliptic type [1]. Numerical methods based on these equations have a variety of applications ranging from mathematical physics to electrostatics, materials science, and electromagnetism (see, e.g., [1–4]). In the context of solving the Dirichlet problem for 2D Laplace’s equation, a standard boundary integral equation method is to represent the solution as a double-layer potential [1–5]. The integral equation that arises in this method is equivalent to the integral equation of Mikhlin whose kernel is referred to as the Neumann kernel [6, p. 282]. In the latter method, the solution is represented as the real part of a Cauchy-type integral with a real density [2, §29]. Indeed, Mikhlin’s integral equation can be seen as the complex counterpart of the integral equation where the solution is sought as a double-layer potential, since the real part of a Cauchy-type integral is the equivalent complex representation of the logarithmic double-layer potential [3, p. 74]. Integral equations of Mikhlin type have been used to solve the interior and exterior Dirichlet problem for Laplace’s equation both in simply and multiply connected domains (see, e.g., [2, 5, 7]). Other applications of Mikhlin’s integral equation include solving Cauchy problems [8] and computing conformal mappings [6].Recently, an integral equation for solving the Riemann–Hilbert boundary value problem was derived in [9, 10]. The kernel of the derived integral equation is a generalization of the classical Neumann kernel. Solvability has been studied for simply connected domains in [10] and for multiply connected domains in [11]. The integral equation with the generalized Neumann kernel (briefly, gNk) has been used to solve Laplace’s equation in simply and multiply connected domains [12, 13]. Other applications include computing conformal mappings [14], the logarithmic capacity of compact sets [15], and the capacity of generalized condensers [16].It may seem that the integral equation with the gNk approach might be equivalent to Mikhlin’s integral equation approach with no additional advantage. The purpose of this paper is to compare both integral equations through solving Laplace’s equation with Dirichlet boundary condition. The main motive behind that is twofold. First, to show the close connection between the two integral equations. Second, to give numerical evidence that the integral equation with the gNk is not only as accurate as Mikhlin’s integral equation for smooth and piecewise smooth boundaries but also converges faster in many cases.Depending on the domain, our results show that there is a difference with respect to their efficiency as measured by the number of discretization points required to attain the same level of accuracy. This is very important from a numerical point of view and reflects directly on the convergence rate. In this sense and for a given level of accuracy, Mikhlin’s integral equation method is more efficient particularly for elongated ellipses. The integral equation with the gNk method is more efficient for boundaries with highly varying curvature and domains with several corners.The paper is organized as follows. Section2 contains some preliminary material and known solvability conditions for both integral equations. Section 3 is devoted to the Mikhlin’s integral equation method. A detailed numerical treatment from derivation to discretization is presented. The details for numerically computing the solution at interior points are also highlighted. Section 4 covers the method based on the integral equation with the gNk. A simpler setup for deriving the integral equation is presented. Both methods are applied to solve the Dirichlet problem in domains bounded by ellipses in Section 5, and comparison of the numerical results of two test cases is presented. Section 6 treats domains with rapidly changing curvature. Further numerical test cases are compared. Domains with corners are treated in Section 7. The paper ends with some concluding remarks and a discussion in Section 8.


## 2. Notation and Preliminaries

In this section, we recall some preliminary material and establish solvability conditions.


### 2.1. The Domain and the Boundary

We consider in our study both smooth and piecewise smooth planar domains in the complex plane. LetG be a bounded simply connected boundary domain with smooth bounding Jordan curve Γ≔∂G. The boundary Γ is parameterized by a 2π-periodic twice continuously differentiable complex function ηt, 0≤t≤2π, with η′t=dηt/dt≠0 which traverses Γ in the positive orientation, i.e., η:0,2π⟶ℂ and Γ=ηt:t∈0,2π. The positive orientation in this context is that for which G is on the left. We let G− denote the domain exterior to Γ, i.e., G−=ℂ¯\G¯, where G¯≔G∪Γ and ℂ¯=ℂ∪∞. We assume for convenience that the origin of the coordinates is interior to Γ. See Figure 1 for an illustration.Figure 1
A simply connected domainG with smooth boundary Γ.We denote byH the space of Hölder continuous 2π-periodic functions on the boundary Γ. A Hölder continuous function defined on Γ can be viewed both as a function of position and a function of parameter depending on the argument. In this respect, no distinction shall be made between a function of position ϕηt defined on the boundary Γ and a function of parameterization ϕt of the same boundary defined on 0,2π. The case where Γ is piecewise smooth will be discussed in Section 7.


### 2.2. Cauchy-Type Integrals and the Sokhotski-Plemelj Formulae

An analytic function in a domain can be uniquely expressed with the help of an integral over the boundary of the domain. Iffz is an analytic function in G and continuous in the closure G¯, then according to the Cauchy integral formula [3, 17]
(1)12πi∫Γfηη−zdη=fz,z∈G,0,z∈G−.If, however,fz is analytic in G− and continuous in the closure G−¯≔G−∪Γ, then
(2)12πi∫Γfηη−zdη=f∞,z∈G,−fz+f∞,z∈G−.The integral on the left-hand sides of (1) and (2) is known as Cauchy’s integral. For a Hölder continuous function h on Γ, the Cauchy-type integral
(3)Ψz≔12πi∫Γhηη−zdη,z∈Γ, defines a function Ψ that is analytic in G and in G−, i.e., sectionally analytic [3, 17, 18]. It is well known from the theory of analytic functions [3] that the boundary values Ψ+ from inside and Ψ− from outside can be determined by the Sokhotski-Plemelj formulae
(4)Ψ±ζ=±12hζ+12πi∫Γhηη−ζdη,ζ∈Γ,where the Cauchy integral in (4) exists as a Cauchy principal value. The boundary functions Ψ± are both Hölder continuous on Γ [3, 18].


### 2.3. The Generalized Neumann Kernel

LetAζ be a given continuously differentiable complex function on the boundary Γ such that Aζ≠0 everywhere. The gNk Ns,t is defined for s,t∈0,2π×0,2π by [10]
(5)Ns,t≔1πImAsAtη′tηt−ηs,s≠t.The kernelNs,t is a generalization of the well-known Neumann kernel which is a classical kernel in potential theory obtained by setting A=1 [6, p. 282]. Remarkably, the kernel Ns,t is in fact continuous [10], although the denominator looks problematic. We have
(6)Nt,t=1πIm12η′′tη′t−A′tAt.Consequently, the integral operatorN defined on H by
(7)Nρs≔∫ΓNs,tρtdt,s∈0,2π,is compact. We consider also the following kernel:
(8)Ms,t≔1πReAsAtη′tηt−ηs,s,t∈0,2π×0,2π,s≠t,which is singular and its singular part involves the cotangent function [10]. The kernel Ms,t can be split as
(9)Ms,t=−Ks,t+M~s,t,where Ks,t is Hilbert’s singular kernel ([3, p. 46]; [2, p. 118]) known also as the conjugation kernel [19]
(10)Ks,t≔12πcots−t2.The kernelM~s,t is continuous and takes on the diagonal the values
(11)M~t,t=1π12Reη′′tη′t−ReA′tAt.The integral operatorM defined on H by
(12)Mρs≔∫ΓMs,tρtdt,s∈0,2π,is bounded on H [10].The possibility ofλ=±1 being eigenvalues of the gNk depends on the index of the function A, which is defined as the winding number of A with respect to 0,
(13)κ≔indA≔12πargA02π,i.e., the change of the argument of A over one period divided by 2π [3, 18].In this paper, we consider two special cases of the gNk, namely, the gNkN1s,t formed with
(14)A1t≔1,and the gNk N2s,t formed with
(15)A2t≔ηt.The kernelN1s,t is the classical Neumann kernel (see [6, p. 282], and [19, p. 371]). Ultimately, we define the kernels M1s,t and M2s,t from (8) with At replaced by A1t=1 and A2t=ηt, respectively.We consider two integral equations. The first integral equation is Mikhlin’s integral equation whose kernel isN1s,t [2]. The kernel of the second integral equation is N2s,t, and this integral equation is known as the boundary integral equation with the gNk [10, 20].Note thatκ1=indA1=0 and κ2=indA2=1. Thus, we have the following theorem (see, e.g., [1, p. 255], [12]).Theorem 1.
(a)λ=−1 is not an eigenvalue of N1 and λ=1 is a simple eigenvalue of N1 with the constant function as the corresponding eigenfunction.
(b)λ=1 is not an eigenvalue of N2 and λ=−1 is a simple eigenvalue of N2 with the constant function as the corresponding eigenfunction.We have also the following theorem from [10, 12].Theorem 2.
(a) Ifλ is an eigenvalue of N1, then λ∈−1,1.
(b) Ifλ is an eigenvalue of N2, then λ∈−1,1.
(c) Ifλ≠1(λ≠−1) is an eigenvalue of N1(N2), then −λ(−λ) is an eigenvalue of N1(N2).The following theorem follows from [12] (Theorem 8).Theorem 3.
Ifλ≠±1, then λ is an eigenvalue of the kernel N2 if and only if λ is an eigenvalue of N1.


### 2.4. The Dirichlet Problem

The interior Dirichlet problem for Laplace’s equation is to find a harmonic functionuz which satisfies
(16a)Δu=0,inG,(16b)u=γ,onΓ,where Δ is the Laplacian operator Δu=∂2u/∂x2+∂2u/∂y2. The boundary data γ is a prescribed Hölder continuous real-valued function on the boundary Γ. The Dirichlet problem has a unique solution u ([17], p. 93) which can be regarded as the real part of a single-valued function F, analytic in G and continuous up to the boundary Γ ([2], p. 137), i.e.,
(17)uz=ReFz,inG∪Γ.To ensure uniqueness of the functionFz, we assume that [3, p. 208]
(18)ImF0=0.It follows from (16b) that Fz satisfies the boundary condition
(19)ReF+ζ=γζforζ∈Γ,where F+ is the restriction of the function F on Γ. The problem (19) is known in [3] as Schwarz problem. In the following two sections, we discuss two integral equation methods for solving the Dirichlet problem.


## 2.1. The Domain and the Boundary

We consider in our study both smooth and piecewise smooth planar domains in the complex plane. LetG be a bounded simply connected boundary domain with smooth bounding Jordan curve Γ≔∂G. The boundary Γ is parameterized by a 2π-periodic twice continuously differentiable complex function ηt, 0≤t≤2π, with η′t=dηt/dt≠0 which traverses Γ in the positive orientation, i.e., η:0,2π⟶ℂ and Γ=ηt:t∈0,2π. The positive orientation in this context is that for which G is on the left. We let G− denote the domain exterior to Γ, i.e., G−=ℂ¯\G¯, where G¯≔G∪Γ and ℂ¯=ℂ∪∞. We assume for convenience that the origin of the coordinates is interior to Γ. See Figure 1 for an illustration.Figure 1
A simply connected domainG with smooth boundary Γ.We denote byH the space of Hölder continuous 2π-periodic functions on the boundary Γ. A Hölder continuous function defined on Γ can be viewed both as a function of position and a function of parameter depending on the argument. In this respect, no distinction shall be made between a function of position ϕηt defined on the boundary Γ and a function of parameterization ϕt of the same boundary defined on 0,2π. The case where Γ is piecewise smooth will be discussed in Section 7.


## 2.2. Cauchy-Type Integrals and the Sokhotski-Plemelj Formulae

An analytic function in a domain can be uniquely expressed with the help of an integral over the boundary of the domain. Iffz is an analytic function in G and continuous in the closure G¯, then according to the Cauchy integral formula [3, 17]
(1)12πi∫Γfηη−zdη=fz,z∈G,0,z∈G−.If, however,fz is analytic in G− and continuous in the closure G−¯≔G−∪Γ, then
(2)12πi∫Γfηη−zdη=f∞,z∈G,−fz+f∞,z∈G−.The integral on the left-hand sides of (1) and (2) is known as Cauchy’s integral. For a Hölder continuous function h on Γ, the Cauchy-type integral
(3)Ψz≔12πi∫Γhηη−zdη,z∈Γ, defines a function Ψ that is analytic in G and in G−, i.e., sectionally analytic [3, 17, 18]. It is well known from the theory of analytic functions [3] that the boundary values Ψ+ from inside and Ψ− from outside can be determined by the Sokhotski-Plemelj formulae
(4)Ψ±ζ=±12hζ+12πi∫Γhηη−ζdη,ζ∈Γ,where the Cauchy integral in (4) exists as a Cauchy principal value. The boundary functions Ψ± are both Hölder continuous on Γ [3, 18].


## 2.3. The Generalized Neumann Kernel

LetAζ be a given continuously differentiable complex function on the boundary Γ such that Aζ≠0 everywhere. The gNk Ns,t is defined for s,t∈0,2π×0,2π by [10]
(5)Ns,t≔1πImAsAtη′tηt−ηs,s≠t.The kernelNs,t is a generalization of the well-known Neumann kernel which is a classical kernel in potential theory obtained by setting A=1 [6, p. 282]. Remarkably, the kernel Ns,t is in fact continuous [10], although the denominator looks problematic. We have
(6)Nt,t=1πIm12η′′tη′t−A′tAt.Consequently, the integral operatorN defined on H by
(7)Nρs≔∫ΓNs,tρtdt,s∈0,2π,is compact. We consider also the following kernel:
(8)Ms,t≔1πReAsAtη′tηt−ηs,s,t∈0,2π×0,2π,s≠t,which is singular and its singular part involves the cotangent function [10]. The kernel Ms,t can be split as
(9)Ms,t=−Ks,t+M~s,t,where Ks,t is Hilbert’s singular kernel ([3, p. 46]; [2, p. 118]) known also as the conjugation kernel [19]
(10)Ks,t≔12πcots−t2.The kernelM~s,t is continuous and takes on the diagonal the values
(11)M~t,t=1π12Reη′′tη′t−ReA′tAt.The integral operatorM defined on H by
(12)Mρs≔∫ΓMs,tρtdt,s∈0,2π,is bounded on H [10].The possibility ofλ=±1 being eigenvalues of the gNk depends on the index of the function A, which is defined as the winding number of A with respect to 0,
(13)κ≔indA≔12πargA02π,i.e., the change of the argument of A over one period divided by 2π [3, 18].In this paper, we consider two special cases of the gNk, namely, the gNkN1s,t formed with
(14)A1t≔1,and the gNk N2s,t formed with
(15)A2t≔ηt.The kernelN1s,t is the classical Neumann kernel (see [6, p. 282], and [19, p. 371]). Ultimately, we define the kernels M1s,t and M2s,t from (8) with At replaced by A1t=1 and A2t=ηt, respectively.We consider two integral equations. The first integral equation is Mikhlin’s integral equation whose kernel isN1s,t [2]. The kernel of the second integral equation is N2s,t, and this integral equation is known as the boundary integral equation with the gNk [10, 20].Note thatκ1=indA1=0 and κ2=indA2=1. Thus, we have the following theorem (see, e.g., [1, p. 255], [12]).Theorem 1.
(a)λ=−1 is not an eigenvalue of N1 and λ=1 is a simple eigenvalue of N1 with the constant function as the corresponding eigenfunction.
(b)λ=1 is not an eigenvalue of N2 and λ=−1 is a simple eigenvalue of N2 with the constant function as the corresponding eigenfunction.We have also the following theorem from [10, 12].Theorem 2.
(a) Ifλ is an eigenvalue of N1, then λ∈−1,1.
(b) Ifλ is an eigenvalue of N2, then λ∈−1,1.
(c) Ifλ≠1(λ≠−1) is an eigenvalue of N1(N2), then −λ(−λ) is an eigenvalue of N1(N2).The following theorem follows from [12] (Theorem 8).Theorem 3.
Ifλ≠±1, then λ is an eigenvalue of the kernel N2 if and only if λ is an eigenvalue of N1.


## 2.4. The Dirichlet Problem

The interior Dirichlet problem for Laplace’s equation is to find a harmonic functionuz which satisfies
(16a)Δu=0,inG,(16b)u=γ,onΓ,where Δ is the Laplacian operator Δu=∂2u/∂x2+∂2u/∂y2. The boundary data γ is a prescribed Hölder continuous real-valued function on the boundary Γ. The Dirichlet problem has a unique solution u ([17], p. 93) which can be regarded as the real part of a single-valued function F, analytic in G and continuous up to the boundary Γ ([2], p. 137), i.e.,
(17)uz=ReFz,inG∪Γ.To ensure uniqueness of the functionFz, we assume that [3, p. 208]
(18)ImF0=0.It follows from (16b) that Fz satisfies the boundary condition
(19)ReF+ζ=γζforζ∈Γ,where F+ is the restriction of the function F on Γ. The problem (19) is known in [3] as Schwarz problem. In the following two sections, we discuss two integral equation methods for solving the Dirichlet problem.


## 3. Mikhlin’s Integral Equation

### 3.1. The Integral Equation

In this section, we review the well-known Mikhlin’s integral equation which is used to compute the values of the analytic functionFz. We seek the function Fz in the form of a Cauchy-type integral [2, §29]
(20)Fz=12πi∫Γhηη−zdη,z∈G,where hη is an unknown Hölder continuous real-valued density function defined on the boundary Γ. The task is now reduced to find the density hη. Using the Sokhotski-Plemelj formulae (4), we obtain
(21)F+ζ=12hζ+12πi∫Γhηη−ζdη.Using the boundary condition (19), we get
(22)Re12hζ+12πi∫Γhηη−ζdη=γζ.Further simplification yields(23)12hζ+12π∫ΓImdηη−ζhη=γζ.The resulting integral equation (23) is a Fredholm integral equation of the second kind known as Mikhlin’s integral equation [2, §29].We use the parameterizationη=ηt,0≤t≤2π, of the boundary and set ζ=ηs,0≤s≤2π to obtain
(24)hηs+∫02πN1ηs,ηthηtdt=2γηs.Equation (24) can be written as
(25)hs+∫02πN1s,thtdt=2γs,where the kernel N1 is the classical Neumann kernel defined in Section 2.3. By Theorem 1, λ=−1 is not an eigenvalue of the kernel N1. Consequently, by the Fredholm alternative, the integral equation (25) is uniquely solvable. Since ∫02πN1s,tdt=1, using regularization, the integral equation (25) can be written as
(26)2hs+∫02πN1s,tht−hsdt=2γs.In particular, using regularization before discretization into a system of linear equations is advantageous and provides an alternative to using limits for the diagonal elements (formula (6)), since the entire integrand in equation (9) vanishes for s=t ([21, p. 101]).


### 3.2. The Nyström Method

The integral equation (26) is discretized into a linear system Sx=y, where S is a square matrix, by means of the Nyström method and the trapezoidal quadrature rule. Note that since the integrand is smooth and periodic, the trapezoidal rule is an optimal choice [22]. The Nyström method has the property of preserving both the stability of the original integral equation ([17, p. 282]) ([1, p. 383]) and the convergence order of the underlying quadrature rule ([17, p. 282]).Given a positive integern, the integral in (26) is discretized using the trapezoidal rule with equal weights wj=2π/n and equally spaced nodes tj=j−12π/n, j=1,⋯,n. We obtain the semidiscrete equation
(27)2hs+2πn∑j=1nN1s,tjhtj−hs=2γs,s∈0,2π.To get a fully discrete equation, we require that (27) should hold at the quadrature points. We set s=tk, k=1,⋯,n. This results in the system of equations
(28)2htk+2πn∑j=1j≠knN1tk,tjhtj−htk=2γtk,k=1,⋯,n,where the term under the summation sign is zero when k=j since the kernel N1 is continuous. Hence, (28) can be written as
(29)2−2πn∑j=1j≠knN1tk,tjhtk+2πn∑j=1j≠knN1tk,tjhtj=2γtk,k=1,⋯,n.We can write (29) in matrix form as
(30)Sx=y,where S≔2I−diagB1+B, tk=tk, k=1,⋯,n, x=ht, y=2γt, I is the n×n identity matrix, 1 is an n×1 vector of ones, diagB1 is a diagonal matrix whose diagonal elements are the elements of the vector B1, and
(31)Bkj=0,k=j,2πnN1tk,tj,k≠j.Note that matrixB has zeros in the main diagonal. Matrix S of the linear system (30) is nonsymmetric, invertible, and dense.


### 3.3. Computing the Functionuz

Once we solve the linear system (30) and obtain an approximation of the density function h, we use it to compute Fz, according to
(32)Fz=12πi∫02πhtηt−zη′tdt≈1ni∑k=1nhtkη′tkηtk−z.The solutionu in the domain G can be evaluated as uz=ReFz.It is worth mentioning that as long asz does not lie on Γ, the integrand in the middle of equation (32) is a smooth periodic function on 0,2π. However, when z gets closer to the boundary Γ, the integrand in (32) is nearly singular and the accuracy of the quadrature in (32) is lost [23]. There are several techniques to overcome such loss (see, e.g., [5, 23–25]). An accurate method has been presented in [5]. The idea in this method is to use the numerical solution of the integral equation to first approximate the boundary values F+ζ of the analytic function Fz and then use the Cauchy integral formula with singularity subtraction to compute the values of Fz for z∈G. This method gives accurate results even when z is very close to the boundary Γ [5, 26].For convenience, we present here the details of the method. Let us first split the second term on the right-hand side in (21) into real and imaginary parts as
(33)F+ζ=12hζ+Re12πi∫Γhηη−ζdη+iIm12πi∫Γhηη−ζdη,which in view of (22) can be written as
(34)F+ζ=γζ+iIm12πi∫Γhηη−ζdη.Settingζ=ηs in (34) and using the definition of M1s,t, we get
(35)F+ζs=γs−i2∫02πM1s,thtdt.Since∫02πM1s,tdt=0 [12], equation (35) can be written as
(36)F+ζs=γs−i2∫02πM1s,tht−hsdt.Note that the kernelM1 is singular. As in (9), it can be split as
(37)M1s,t=−Ks,t+M~1s,t,where M~1s,t is continuous and Ks,t is the Hilbert kernel given by (10).Equation (34) is the same as equation (24) in [5] (Eq. (24)). In [5], the values of F+ζs are approximated by discretization of the integral in (35) by the trapezoidal rule where the integrand is rewritten such that it is continuous even when t=s. However, for t=s, the integrand involves the derivative of ht which is computed numerically via n-point polynomial interpolation ([5, p. 2904]). In this paper, we follow the approach used in [20] which is based on using Wittich’s method and does not require the computation of the derivative of the function ht.We substitutes=tk, k=1,⋯,n, to obtain
(38)F+tk=γtk+i2∫02πKtk,t−M~1tk,tht−htkdt.We proceed by discretizing the continuous kernelM~1 using the trapezoidal rule and the Hilbert kernel K by Wittich’s method [19] to get the fully discrete scheme for assembling F+tk, i.e.,
(39)F+tk=γtk+i2∑j=1j≠knKkj−Ckjhtj−htk,k=1,⋯,n, where the term under the summation sign is zero when k=j since Kk,j=0 whenever k=j and the kernel M~1 is continuous. The matrices C and K are given by
(40)Ck,j=0,k=j,2nReη′tkηtk−ηtj+1ncotk−jπn,k≠j.Kk,j=0,k−jeven,2ncotk−jπn,k−jodd. Or
(41)−Kkj+Ckj=Lkj+Dkj, where
(42)Lkj=0,k=j,−1k−j1ncotk−jπn,k≠j,(43)Dkj=0,k=j,2πnM1tk,tj,k≠j.In light of (41), (42), and (43), equation (39) can be written as
(44)F+tk=γtk−i2∑j=1nLkjhtj−htk−i2∑j=1nDkjhtj−htk or equivalently in matrix form as
(45)F+=y−i2L+D−diagL1+diagD1x, where F+=F+t, x=ht, y=γt, and 1 is a n×1 vector of ones.Now that we have accurately computed an approximation of the boundary valuesF+, the function Fz can be evaluated at any point z in the domain G via the Cauchy integral formula
(46)Fz=12πi∫ΓF+ζζ−zdζ,z∈G.The integrand in (46) has a pole at z=ζ. This singularity can be removed by dividing the usual Cauchy integral formula for F by the same formula for 1, i.e., 1/2πi∫Γ1/ζ−zdζ=1, and rearranging to obtain [5, 26, 27]
(47)∫ΓF+ζ−Fzζ−zdζ=0,z∈G.Contrary to the integrand in (46), the integrand in (47) has no pole at z=ζ and is by consequence an analytic function of z∈G whose Cauchy integral must be equal to zero [26]. Using the trapezoidal rule to discretize this integral yields
(48)∑j=1nF+ηtj−Fzηtj−zη′tj=0,z∈G.Solving forFz in (48) results in the barycentric formula [5, 26, 27]
(49)Fz=∑j=1nF+ηtjη′tj/ηtj−z∑j=1nη′tj/ηtj−z.Formula (49) is numerically stable even when the evaluation point z approaches the boundary Γ arbitrarily closely (see [5, 26]). It is referred to in [26] as discrete Cauchy integral of the second kind.The above method for the numerical solution of the Dirichlet problem can be summarized in the following algorithm:Algorithm 1: Solving the Dirichlet problem with Mikhlin’s integral equation method
Step 1. Solve the integral equation (25) for h.Step 2. Find the boundary values ofF given in equation (35).Step 3. Compute the values ofFz for z∈G using (49).Step 4. Finally, evaluate the solutionuz through uz=ReFz.


## 3.1. The Integral Equation

In this section, we review the well-known Mikhlin’s integral equation which is used to compute the values of the analytic functionFz. We seek the function Fz in the form of a Cauchy-type integral [2, §29]
(20)Fz=12πi∫Γhηη−zdη,z∈G,where hη is an unknown Hölder continuous real-valued density function defined on the boundary Γ. The task is now reduced to find the density hη. Using the Sokhotski-Plemelj formulae (4), we obtain
(21)F+ζ=12hζ+12πi∫Γhηη−ζdη.Using the boundary condition (19), we get
(22)Re12hζ+12πi∫Γhηη−ζdη=γζ.Further simplification yields(23)12hζ+12π∫ΓImdηη−ζhη=γζ.The resulting integral equation (23) is a Fredholm integral equation of the second kind known as Mikhlin’s integral equation [2, §29].We use the parameterizationη=ηt,0≤t≤2π, of the boundary and set ζ=ηs,0≤s≤2π to obtain
(24)hηs+∫02πN1ηs,ηthηtdt=2γηs.Equation (24) can be written as
(25)hs+∫02πN1s,thtdt=2γs,where the kernel N1 is the classical Neumann kernel defined in Section 2.3. By Theorem 1, λ=−1 is not an eigenvalue of the kernel N1. Consequently, by the Fredholm alternative, the integral equation (25) is uniquely solvable. Since ∫02πN1s,tdt=1, using regularization, the integral equation (25) can be written as
(26)2hs+∫02πN1s,tht−hsdt=2γs.In particular, using regularization before discretization into a system of linear equations is advantageous and provides an alternative to using limits for the diagonal elements (formula (6)), since the entire integrand in equation (9) vanishes for s=t ([21, p. 101]).


## 3.2. The Nyström Method

The integral equation (26) is discretized into a linear system Sx=y, where S is a square matrix, by means of the Nyström method and the trapezoidal quadrature rule. Note that since the integrand is smooth and periodic, the trapezoidal rule is an optimal choice [22]. The Nyström method has the property of preserving both the stability of the original integral equation ([17, p. 282]) ([1, p. 383]) and the convergence order of the underlying quadrature rule ([17, p. 282]).Given a positive integern, the integral in (26) is discretized using the trapezoidal rule with equal weights wj=2π/n and equally spaced nodes tj=j−12π/n, j=1,⋯,n. We obtain the semidiscrete equation
(27)2hs+2πn∑j=1nN1s,tjhtj−hs=2γs,s∈0,2π.To get a fully discrete equation, we require that (27) should hold at the quadrature points. We set s=tk, k=1,⋯,n. This results in the system of equations
(28)2htk+2πn∑j=1j≠knN1tk,tjhtj−htk=2γtk,k=1,⋯,n,where the term under the summation sign is zero when k=j since the kernel N1 is continuous. Hence, (28) can be written as
(29)2−2πn∑j=1j≠knN1tk,tjhtk+2πn∑j=1j≠knN1tk,tjhtj=2γtk,k=1,⋯,n.We can write (29) in matrix form as
(30)Sx=y,where S≔2I−diagB1+B, tk=tk, k=1,⋯,n, x=ht, y=2γt, I is the n×n identity matrix, 1 is an n×1 vector of ones, diagB1 is a diagonal matrix whose diagonal elements are the elements of the vector B1, and
(31)Bkj=0,k=j,2πnN1tk,tj,k≠j.Note that matrixB has zeros in the main diagonal. Matrix S of the linear system (30) is nonsymmetric, invertible, and dense.


## 3.3. Computing the Functionuz

Once we solve the linear system (30) and obtain an approximation of the density function h, we use it to compute Fz, according to
(32)Fz=12πi∫02πhtηt−zη′tdt≈1ni∑k=1nhtkη′tkηtk−z.The solutionu in the domain G can be evaluated as uz=ReFz.It is worth mentioning that as long asz does not lie on Γ, the integrand in the middle of equation (32) is a smooth periodic function on 0,2π. However, when z gets closer to the boundary Γ, the integrand in (32) is nearly singular and the accuracy of the quadrature in (32) is lost [23]. There are several techniques to overcome such loss (see, e.g., [5, 23–25]). An accurate method has been presented in [5]. The idea in this method is to use the numerical solution of the integral equation to first approximate the boundary values F+ζ of the analytic function Fz and then use the Cauchy integral formula with singularity subtraction to compute the values of Fz for z∈G. This method gives accurate results even when z is very close to the boundary Γ [5, 26].For convenience, we present here the details of the method. Let us first split the second term on the right-hand side in (21) into real and imaginary parts as
(33)F+ζ=12hζ+Re12πi∫Γhηη−ζdη+iIm12πi∫Γhηη−ζdη,which in view of (22) can be written as
(34)F+ζ=γζ+iIm12πi∫Γhηη−ζdη.Settingζ=ηs in (34) and using the definition of M1s,t, we get
(35)F+ζs=γs−i2∫02πM1s,thtdt.Since∫02πM1s,tdt=0 [12], equation (35) can be written as
(36)F+ζs=γs−i2∫02πM1s,tht−hsdt.Note that the kernelM1 is singular. As in (9), it can be split as
(37)M1s,t=−Ks,t+M~1s,t,where M~1s,t is continuous and Ks,t is the Hilbert kernel given by (10).Equation (34) is the same as equation (24) in [5] (Eq. (24)). In [5], the values of F+ζs are approximated by discretization of the integral in (35) by the trapezoidal rule where the integrand is rewritten such that it is continuous even when t=s. However, for t=s, the integrand involves the derivative of ht which is computed numerically via n-point polynomial interpolation ([5, p. 2904]). In this paper, we follow the approach used in [20] which is based on using Wittich’s method and does not require the computation of the derivative of the function ht.We substitutes=tk, k=1,⋯,n, to obtain
(38)F+tk=γtk+i2∫02πKtk,t−M~1tk,tht−htkdt.We proceed by discretizing the continuous kernelM~1 using the trapezoidal rule and the Hilbert kernel K by Wittich’s method [19] to get the fully discrete scheme for assembling F+tk, i.e.,
(39)F+tk=γtk+i2∑j=1j≠knKkj−Ckjhtj−htk,k=1,⋯,n, where the term under the summation sign is zero when k=j since Kk,j=0 whenever k=j and the kernel M~1 is continuous. The matrices C and K are given by
(40)Ck,j=0,k=j,2nReη′tkηtk−ηtj+1ncotk−jπn,k≠j.Kk,j=0,k−jeven,2ncotk−jπn,k−jodd. Or
(41)−Kkj+Ckj=Lkj+Dkj, where
(42)Lkj=0,k=j,−1k−j1ncotk−jπn,k≠j,(43)Dkj=0,k=j,2πnM1tk,tj,k≠j.In light of (41), (42), and (43), equation (39) can be written as
(44)F+tk=γtk−i2∑j=1nLkjhtj−htk−i2∑j=1nDkjhtj−htk or equivalently in matrix form as
(45)F+=y−i2L+D−diagL1+diagD1x, where F+=F+t, x=ht, y=γt, and 1 is a n×1 vector of ones.Now that we have accurately computed an approximation of the boundary valuesF+, the function Fz can be evaluated at any point z in the domain G via the Cauchy integral formula
(46)Fz=12πi∫ΓF+ζζ−zdζ,z∈G.The integrand in (46) has a pole at z=ζ. This singularity can be removed by dividing the usual Cauchy integral formula for F by the same formula for 1, i.e., 1/2πi∫Γ1/ζ−zdζ=1, and rearranging to obtain [5, 26, 27]
(47)∫ΓF+ζ−Fzζ−zdζ=0,z∈G.Contrary to the integrand in (46), the integrand in (47) has no pole at z=ζ and is by consequence an analytic function of z∈G whose Cauchy integral must be equal to zero [26]. Using the trapezoidal rule to discretize this integral yields
(48)∑j=1nF+ηtj−Fzηtj−zη′tj=0,z∈G.Solving forFz in (48) results in the barycentric formula [5, 26, 27]
(49)Fz=∑j=1nF+ηtjη′tj/ηtj−z∑j=1nη′tj/ηtj−z.Formula (49) is numerically stable even when the evaluation point z approaches the boundary Γ arbitrarily closely (see [5, 26]). It is referred to in [26] as discrete Cauchy integral of the second kind.The above method for the numerical solution of the Dirichlet problem can be summarized in the following algorithm:Algorithm 1: Solving the Dirichlet problem with Mikhlin’s integral equation method
Step 1. Solve the integral equation (25) for h.Step 2. Find the boundary values ofF given in equation (35).Step 3. Compute the values ofFz for z∈G using (49).Step 4. Finally, evaluate the solutionuz through uz=ReFz.


## 4. The Integral Equation with the gNk

The integral equation with the gNk has been derived for solving the Riemann–Hilbert problem in simply and multiply connected domains (see, e.g., [9–11]). As stated in [10] (Theorem 11), the Dirichlet problem in simply connected domains can be solved using the integral equation with the gNk. However, as the objective of this paper is to compare both integral equation methods, we proceed, for completeness, with deriving the integral equation first. The approach used here is slightly different, but simpler than the one used in [10].


### 4.1. The Riemann–Hilbert Problem

Since we are interested in computinguz=ReFz, we can assume, without loss of generality, that F0=c for some real constant c. We define a function f in G as
(50)fz=Fz−cz,so that
(51)Fz=zfz+c.The functionfz is analytic in G and its boundary values satisfy the boundary condition
(52)ReA2tf+ηt=γt−cforηt∈Γ,with the function A2t=ηt. The problem (52) is a Riemann–Hilbert problem with the coefficient A2t=ηt which has the index κ2=1. This problem is solvable only if the right-hand side satisfies one solvability condition [10]. If this condition is satisfied, then the problem has a unique solution. The undetermined real constant c will be chosen so that the solvability condition is satisfied.


### 4.2. Derivation of the Integral Equation

Letfz be the unique solution of the Riemann–Hilbert problem (52) and
(53)μη≔ImA2ηf+ηforη∈Γ.Thus, the boundary values of the functionfz satisfy
(54)A2ηf+η=γη−c+iμηforη∈Γ.Letz∈G−, since fz is analytic in G, then by the Cauchy integral formula (1)
(55)12πi∫Γf+ηη−zdη=0.Solving (54) for f+ and substituting in (55) yields
(56)12πi∫Γγη−c+iμηA2ηdηη−z=0,which implies
(57)12πi∫Γγη+iμηA2ηdηη−z=c12πi∫Γ1/ηη−zdη.Since the functiongz=1/z in the numerator of the right-hand side in (57) is analytic for z∈G− and Γ is counterclockwise oriented, then using the Cauchy integral formula (2) yields
(58)12πi∫Γγη+iμηA2ηdηη−z=−cz,z∈G−.We proceed by taking the limitG−∋z⟶ζ∈Γ on both sides of equation (58) and applying the Sokhotski-Plemelj formula (4) to the left-hand side to obtain
(59)−12γζ+iμζA2ζ+12πi∫Γγη+iμηA2ηdηη−ζ=−cζ.SinceA2ζ=ζ, multiplying (59) by −2A2ζ=−2ζ gives
(60)γζ+iμζ−1πi∫Γγη+iμηA2ζA2ηdηη−ζ=2c.Using the parameterizationη=ηt,0≤t≤2π, and ζ=ηs,0≤s≤2π, we get
(61)γs+iμs−1i∫02πγt+iμt1πA2sA2tη′tηt−ηsdt=2c.Then, using the definitions of the kernelsN2 and M2 (see Section 2.3), we obtain
(62)γs+iμs+∫02π−M2s,tμt−N2s,tγt+iM2s,tγt−N2s,tμtdt=2c.Taking the imaginary parts of both sides in (62) yields
(63)μs−∫02πN2s,tμtdt=−∫02πM2s,tγtdt,which is the integral equation in [10] (Eq. (98)). This equation is known as the integral equation with the gNk. The integral equation (63) is uniquely solvable (Theorem 1). Note that the integral equation (63) does not involve the undetermined real constant c. Manifestly, computing the boundary values of F requires only computing the function μ, the solution of the integral equation (63). Indeed, from (51) and (54), we have
(64)F+ζ=γζ+iμζ,ζ∈Γ.The values ofFz for z∈G can be computed using (49). Note that the boundary values in formulas (34) and (64) can differ from each other only by an imaginary constant [2].


### 4.3. Discretization of the Integral Equation

The regularized form of the integral equation with the gNk (63) is discretized by the Nyström method and the trapezoidal quadrature rule which gives the linear system [20]
(65)S^x^=−y^,where S^=2I+diagB^1−B^, x^=μt, and matrix B^ is defined by
(66)B^kj=0,k=j,2πnN2tk,tj,k≠j.The right-hand side in (65) is given by [20]
(67)y^=D^−diagD^1+Ly, where y=γt, the matrix L is defined in equation (42), 1 is a n×1 vector of ones, and matrix D^ is defined by
(68)D^kj=0,k=j,2πnM2tk,tj,k≠j.We omit the details here and refer the reader to [20] for a thorough description of the method.In [20], the linear system (65) was solved by the GMRES iterative method accelerated by the Fast Multipole Method (FMM). Since our objective here is the comparison between the two integral equations, we shall solve both linear systems using the MATLAB \ operator.Algorithm 2: Solving the Dirichlet problem using the integral equation with the gNk method
Step 1. Solve the integral equation (63) for μ.Step 2. Compute the values ofFz for z∈G using formula (49), where the boundary values of F are given by (64).Step 3. Finally, evaluate the solutionuz as uz=ReFz.


## 4.1. The Riemann–Hilbert Problem

Since we are interested in computinguz=ReFz, we can assume, without loss of generality, that F0=c for some real constant c. We define a function f in G as
(50)fz=Fz−cz,so that
(51)Fz=zfz+c.The functionfz is analytic in G and its boundary values satisfy the boundary condition
(52)ReA2tf+ηt=γt−cforηt∈Γ,with the function A2t=ηt. The problem (52) is a Riemann–Hilbert problem with the coefficient A2t=ηt which has the index κ2=1. This problem is solvable only if the right-hand side satisfies one solvability condition [10]. If this condition is satisfied, then the problem has a unique solution. The undetermined real constant c will be chosen so that the solvability condition is satisfied.


## 4.2. Derivation of the Integral Equation

Letfz be the unique solution of the Riemann–Hilbert problem (52) and
(53)μη≔ImA2ηf+ηforη∈Γ.Thus, the boundary values of the functionfz satisfy
(54)A2ηf+η=γη−c+iμηforη∈Γ.Letz∈G−, since fz is analytic in G, then by the Cauchy integral formula (1)
(55)12πi∫Γf+ηη−zdη=0.Solving (54) for f+ and substituting in (55) yields
(56)12πi∫Γγη−c+iμηA2ηdηη−z=0,which implies
(57)12πi∫Γγη+iμηA2ηdηη−z=c12πi∫Γ1/ηη−zdη.Since the functiongz=1/z in the numerator of the right-hand side in (57) is analytic for z∈G− and Γ is counterclockwise oriented, then using the Cauchy integral formula (2) yields
(58)12πi∫Γγη+iμηA2ηdηη−z=−cz,z∈G−.We proceed by taking the limitG−∋z⟶ζ∈Γ on both sides of equation (58) and applying the Sokhotski-Plemelj formula (4) to the left-hand side to obtain
(59)−12γζ+iμζA2ζ+12πi∫Γγη+iμηA2ηdηη−ζ=−cζ.SinceA2ζ=ζ, multiplying (59) by −2A2ζ=−2ζ gives
(60)γζ+iμζ−1πi∫Γγη+iμηA2ζA2ηdηη−ζ=2c.Using the parameterizationη=ηt,0≤t≤2π, and ζ=ηs,0≤s≤2π, we get
(61)γs+iμs−1i∫02πγt+iμt1πA2sA2tη′tηt−ηsdt=2c.Then, using the definitions of the kernelsN2 and M2 (see Section 2.3), we obtain
(62)γs+iμs+∫02π−M2s,tμt−N2s,tγt+iM2s,tγt−N2s,tμtdt=2c.Taking the imaginary parts of both sides in (62) yields
(63)μs−∫02πN2s,tμtdt=−∫02πM2s,tγtdt,which is the integral equation in [10] (Eq. (98)). This equation is known as the integral equation with the gNk. The integral equation (63) is uniquely solvable (Theorem 1). Note that the integral equation (63) does not involve the undetermined real constant c. Manifestly, computing the boundary values of F requires only computing the function μ, the solution of the integral equation (63). Indeed, from (51) and (54), we have
(64)F+ζ=γζ+iμζ,ζ∈Γ.The values ofFz for z∈G can be computed using (49). Note that the boundary values in formulas (34) and (64) can differ from each other only by an imaginary constant [2].


## 4.3. Discretization of the Integral Equation

The regularized form of the integral equation with the gNk (63) is discretized by the Nyström method and the trapezoidal quadrature rule which gives the linear system [20]
(65)S^x^=−y^,where S^=2I+diagB^1−B^, x^=μt, and matrix B^ is defined by
(66)B^kj=0,k=j,2πnN2tk,tj,k≠j.The right-hand side in (65) is given by [20]
(67)y^=D^−diagD^1+Ly, where y=γt, the matrix L is defined in equation (42), 1 is a n×1 vector of ones, and matrix D^ is defined by
(68)D^kj=0,k=j,2πnM2tk,tj,k≠j.We omit the details here and refer the reader to [20] for a thorough description of the method.In [20], the linear system (65) was solved by the GMRES iterative method accelerated by the Fast Multipole Method (FMM). Since our objective here is the comparison between the two integral equations, we shall solve both linear systems using the MATLAB \ operator.Algorithm 2: Solving the Dirichlet problem using the integral equation with the gNk method
Step 1. Solve the integral equation (63) for μ.Step 2. Compute the values ofFz for z∈G using formula (49), where the boundary values of F are given by (64).Step 3. Finally, evaluate the solutionuz as uz=ReFz.


## 5. Domains Bounded by Ellipses

In this section, we consider the simply connected domainG interior to the ellipse Γ with the parameterization
(69)ηt=acost+isint,0≤t≤2π,for a>0, i.e., the minor radius is 1 and the major radius is a (see Figure 2 for a=5). We present in the following subsections a comparison between Mikhlin’s integral equation method and the integral equation with the gNk method in domains bounded by the ellipse in (69). We start by comparing the eigenvalues of the coefficient matrices.Figure 2
The ellipse with parameterization in equation (69) with a=5.


### 5.1. Eigenvalues

For the ellipse (69), the explicit form of the eigenvalues of the Neumann kernel are known and are given by (see [28], Lemma 1.1, and [1, p. 255])
(70)λ^0=1,λ^2k−1=−θk,λ^2k=θk,k=1,2,3,⋯,where
(71)θ=a−1a+1,0≤θ<1.This is in agreement with Theorem2(c). Further, since θk⟶0 as k⟶∞, the point 0 is an accumulated point of the eigenvalues of the Neumann kernel. In view of Theorems 1 and 3, the eigenvalues of the gNk are
(72)λ~0=−1,λ^2k−1=θk,λ^2k=−θk,k=1,2,3,⋯.Consequently, in light of Theorem2, the coefficient matrix 2I−diagB1+B of Mikhlin’s integral equation (25) and the coefficient matrix 2I+diagB^1−B^ of the integral equation with the gNk (63) have the same eigenvalues and these eigenvalues lie in the interval 0,2 for sufficiently large values of n. The n approximate eigenvalues of the coefficient matrices, sorted decreasingly, are given by
(73)λ1=2,λ2=1+θ,λ3=1+θ2,⋯,λn−1=1−θ2,λn=1−θ.Sinceθk≈0 for small θ and large k, these eigenvalues are clustered in a remarkably symmetric way around 1. In fact, for small θ and for sufficiently large values of n, most of these eigenvalues are equal to 1. For elongated domains, for which a>>1 or a<<1 and hence θ near to 1, we can notice that more eigenvalues are different from 1. In all cases, the largest eigenvalue of the coefficient matrices is λmax≈2 and the smallest eigenvalue is λmin≈1−θ. For a=1.5, we have θ=0.2 and the eigenvalues are too accumulated around 1 (see Figure 3(a)), where the computed smallest eigenvalue is λmin=0.799999999999999. For the elongated ellipse with a=19, the computed eigenvalues are shown in Figure 3(b). In this case, θ=0.9, and hence, we have several eigenvalues slightly away from 1 where the computed smallest eigenvalue is λmin=0.1.Figure 3
The eigenvalues of the coefficient matrices for (a)a=1.5 and (b) a=19.In Figure3, we present only the eigenvalues of the matrix 2I+diagB^1−B^. However, we chose the values of n such that the maximum norm between the approximate eigenvalues of the two matrices 2I−diagB1+B and 2I+diagB^1−B^ is less than 10−13.


### 5.2. Closed-Form Expressions for the KernelsN1, M1, N2, and M2

In the case of the ellipse (69), it is well known that the Neumann kernel can be expressed in closed form as (see [21, p. 135–136])
(74)N1s,t≔1πImη′tηt−ηs=1πa1+a2+1−a2coss+t.The explicit closed-form expressions for the kernelsM1, N2, and M2 defined in Section 2.3 are given in the following theorem.Theorem 4.
LetΓ be the ellipse parameterized by ηt≔acost+isint, where 0≤t≤2π and a>0. The closed-form expressions for the kernels M1, N2, and M2 are, respectively, given by
(a)(75)M1s,t=−12πcots−t2−12π1−a2sint+s1+a2+1−a2cost+s,s≠t,
(b)(76)N2s,t=−1π2a1+a2−1−a2cos2t+1πa1+a2+1−a2cost+s,
(c)(77)M2s,t=−12πcots−t2−1π1−a2sin2t1+a2−1−a2cos2t−12π1−a2sint+s1+a2+1−a2cost+s,s≠t.Proof.
(a) Fors≠t, we have
(78)η′tηt−ηs=−asint+icostacost−coss+isint−sins.
Applying the sum-to-product trigonometric identitiescost−coss=−2sint+s/2sint−s/2 and sint−sins=2cost+s/2sint−s/2, (78) becomes
(79)η′tηt−ηs=asint−icost2sint−s/2asint+s/2−icost+s/2,=a2sintsint+s/2+costcost+s/2+iasintcost+s/2−costsint+s/22sint−s/2a2sin2t+s/2+cos2t+s/2.
Applying the half-angle identitiessin2θ/2=1−cosθ/2 and cos2θ/2=1+cosθ/2 to the bracket term in the denominator and taking the real parts yields
(80)Reη′tηt−ηs=a2sintsint+s/2+costcost+s/2sint−s/21+a2+1−a2cost+s.
Applying the product-to-sum trigonometric identities2sinAsinB=−cosA+B+cosA−B and 2cosAcosB=cosA+B+cosA−B to the numerator part of (80) and rearranging, we get
(81)Reη′tηt−ηs=1/21+a2cost−s/2+1−a2cos3t+s/2sint−s/21+a2+1−a2cost+s.
However,(82)cos3t+s2=cost−s2+t+s=cost−s2cost+s−sint−s2sint+s.
Substituting this result into (71) and rearranging gives
(83)Reη′tηt−ηs=1/2cost−s/21+a2+1−a2cost+s−1/21−a2sint−s/2sint+ssint−s/21+a2+1−a2cost+s.
Finally, multiplying both sides of (83) by 1/π and simplifying gives (75).
(b) Fors≠t, we have
(84)ηsηtη′tηt−ηs=ηs−ηt+ηtηtη′tηt−ηs=−η′tηt+η′tηt−ηs.
Observe that(85)η′tηt=1−a2sin2t+2ia1+a2−1−a2cos2t.
So(86)N2s,t=1πImηsηtη′tηt−ηs=−1πImη′tηt+1πImη′tηt−ηs.
In light of (74) and (85), we get
(87)N2s,t=−1π2a1+a2−1−a2cos2t+1πa1+a2+1−a2cost+s.
(c) Fors≠t, we have
(88)M2s,t=1πReηsηtη′tηt−ηs=−1πReη′tηt+1πReη′tηt−ηs.
Applying (85) and the result in (a) yields
(89)M2s,t=−12πcots−t2−1π1−a2sin2t1+a2−1−a2cos2t−12π1−a2sint+s1+a2+1−a2cost+s.Example 5.
To illustrate the accuracy of the two methods, consider the Dirichlet problem inside a domain with a boundary curve parameterized by the ellipse (69) and boundary condition
(90)uηt=γt≔Rea+1+ηtforηt∈Γ.
Note that−a is outside the domain and we can always choose any analytic branch of the function a+1+z.For the sake of comparison, we first solve Mikhlin’s integral equation (25) using the method described in Section 3.2, and then, for a=1.5 and α+iβ∈G with real numbers α and β, we use formula (32) to compute the values at α+iβ (i.e., without computing F+ and without using the barycentric formula (49)). Table 1 shows the estimated relative error un−u∞/u∞ of the approximate solution un at two given interior points. We observe that accurate results are obtained for α=β=0 even for small values of n (see Table 1). For α=0 and β=0.999, the point α+iβ is close to the boundary and the obtained results are very inaccurate even for large values of n. To get an accurate approximation, we compute the boundary values F+ according to equation (35) and then compute the solution at α=0 and β=0,0.999 as described in Section 3.3. The values of the relative error of the computed solution un are presented also in Table 1. It is observed that accurate results are obtained for both cases of β even for small values of n.Table 1
The relative errorun−u∞/u∞ of the approximate solution un of the Dirichlet problem for a=1.5 at the points z=α+iβ with α=0 and two values of β for different values of n using Mikhlin’s integral equation and the integral equation with the gNk (IE with the gNk).
nMikhlin (CIF withh)Mikhlin (CIF withF+)IE with the gNkβ=0β=0.999β=0β=0.999β=0β=0.99983.58×10−031.92×10+028.65×10−058.20×10−071.25×10−041.10×10−06165.71×10−069.59×10+011.00×10−072.64×10−081.62×10−072.76×10−08321.45×10−114.77×10+012.25×10−136.14×10−113.84×10−136.14×10−11642.22×10−162.36×10+011.40×10−165.51×10−162.22×10−165.51×10−161281.40×10−161.15×10+011.40×10−161.38×10−161.40×10−161.38×10−16We solve Example5 using the integral equation with the gNk method, as described in Section 4.3 and summarized in Algorithm 2. The estimated relative error of the approximate solutions at α+iβ for α=0 and β=0,0.999 is presented in Table 1. These results are also highly accurate at both points.In Figure4(a), the relative error of the approximate solution at 1000 random interior points is plotted against the number of discretization points for a=2. Both methods achieve high accuracy, exhibit good efficiency, and converge equally fast. However, in the case of the elongated ellipse with a=19 (see Figure 4(b)), although both methods achieve the same accuracy, there is a difference in terms of the cost. Mikhlin’s integral equation method converges faster and achieves high accuracy at half the number of discretization points needed using the integral equation with the gNk method. This computational advantage in favour of Mikhlin’s integral equation method has been noticed for other values of a for which the ellipse (69) is elongated.Figure 4
The relative errorun−u∞/u∞ of the approximate solution un of the Dirichlet problem at 1000 interior points as a function of the number of discretization points for (a) a=2 and (b) a=19.


## 5.1. Eigenvalues

For the ellipse (69), the explicit form of the eigenvalues of the Neumann kernel are known and are given by (see [28], Lemma 1.1, and [1, p. 255])
(70)λ^0=1,λ^2k−1=−θk,λ^2k=θk,k=1,2,3,⋯,where
(71)θ=a−1a+1,0≤θ<1.This is in agreement with Theorem2(c). Further, since θk⟶0 as k⟶∞, the point 0 is an accumulated point of the eigenvalues of the Neumann kernel. In view of Theorems 1 and 3, the eigenvalues of the gNk are
(72)λ~0=−1,λ^2k−1=θk,λ^2k=−θk,k=1,2,3,⋯.Consequently, in light of Theorem2, the coefficient matrix 2I−diagB1+B of Mikhlin’s integral equation (25) and the coefficient matrix 2I+diagB^1−B^ of the integral equation with the gNk (63) have the same eigenvalues and these eigenvalues lie in the interval 0,2 for sufficiently large values of n. The n approximate eigenvalues of the coefficient matrices, sorted decreasingly, are given by
(73)λ1=2,λ2=1+θ,λ3=1+θ2,⋯,λn−1=1−θ2,λn=1−θ.Sinceθk≈0 for small θ and large k, these eigenvalues are clustered in a remarkably symmetric way around 1. In fact, for small θ and for sufficiently large values of n, most of these eigenvalues are equal to 1. For elongated domains, for which a>>1 or a<<1 and hence θ near to 1, we can notice that more eigenvalues are different from 1. In all cases, the largest eigenvalue of the coefficient matrices is λmax≈2 and the smallest eigenvalue is λmin≈1−θ. For a=1.5, we have θ=0.2 and the eigenvalues are too accumulated around 1 (see Figure 3(a)), where the computed smallest eigenvalue is λmin=0.799999999999999. For the elongated ellipse with a=19, the computed eigenvalues are shown in Figure 3(b). In this case, θ=0.9, and hence, we have several eigenvalues slightly away from 1 where the computed smallest eigenvalue is λmin=0.1.Figure 3
The eigenvalues of the coefficient matrices for (a)a=1.5 and (b) a=19.In Figure3, we present only the eigenvalues of the matrix 2I+diagB^1−B^. However, we chose the values of n such that the maximum norm between the approximate eigenvalues of the two matrices 2I−diagB1+B and 2I+diagB^1−B^ is less than 10−13.


## 5.2. Closed-Form Expressions for the KernelsN1, M1, N2, and M2

In the case of the ellipse (69), it is well known that the Neumann kernel can be expressed in closed form as (see [21, p. 135–136])
(74)N1s,t≔1πImη′tηt−ηs=1πa1+a2+1−a2coss+t.The explicit closed-form expressions for the kernelsM1, N2, and M2 defined in Section 2.3 are given in the following theorem.Theorem 4.
LetΓ be the ellipse parameterized by ηt≔acost+isint, where 0≤t≤2π and a>0. The closed-form expressions for the kernels M1, N2, and M2 are, respectively, given by
(a)(75)M1s,t=−12πcots−t2−12π1−a2sint+s1+a2+1−a2cost+s,s≠t,
(b)(76)N2s,t=−1π2a1+a2−1−a2cos2t+1πa1+a2+1−a2cost+s,
(c)(77)M2s,t=−12πcots−t2−1π1−a2sin2t1+a2−1−a2cos2t−12π1−a2sint+s1+a2+1−a2cost+s,s≠t.Proof.
(a) Fors≠t, we have
(78)η′tηt−ηs=−asint+icostacost−coss+isint−sins.
Applying the sum-to-product trigonometric identitiescost−coss=−2sint+s/2sint−s/2 and sint−sins=2cost+s/2sint−s/2, (78) becomes
(79)η′tηt−ηs=asint−icost2sint−s/2asint+s/2−icost+s/2,=a2sintsint+s/2+costcost+s/2+iasintcost+s/2−costsint+s/22sint−s/2a2sin2t+s/2+cos2t+s/2.
Applying the half-angle identitiessin2θ/2=1−cosθ/2 and cos2θ/2=1+cosθ/2 to the bracket term in the denominator and taking the real parts yields
(80)Reη′tηt−ηs=a2sintsint+s/2+costcost+s/2sint−s/21+a2+1−a2cost+s.
Applying the product-to-sum trigonometric identities2sinAsinB=−cosA+B+cosA−B and 2cosAcosB=cosA+B+cosA−B to the numerator part of (80) and rearranging, we get
(81)Reη′tηt−ηs=1/21+a2cost−s/2+1−a2cos3t+s/2sint−s/21+a2+1−a2cost+s.
However,(82)cos3t+s2=cost−s2+t+s=cost−s2cost+s−sint−s2sint+s.
Substituting this result into (71) and rearranging gives
(83)Reη′tηt−ηs=1/2cost−s/21+a2+1−a2cost+s−1/21−a2sint−s/2sint+ssint−s/21+a2+1−a2cost+s.
Finally, multiplying both sides of (83) by 1/π and simplifying gives (75).
(b) Fors≠t, we have
(84)ηsηtη′tηt−ηs=ηs−ηt+ηtηtη′tηt−ηs=−η′tηt+η′tηt−ηs.
Observe that(85)η′tηt=1−a2sin2t+2ia1+a2−1−a2cos2t.
So(86)N2s,t=1πImηsηtη′tηt−ηs=−1πImη′tηt+1πImη′tηt−ηs.
In light of (74) and (85), we get
(87)N2s,t=−1π2a1+a2−1−a2cos2t+1πa1+a2+1−a2cost+s.
(c) Fors≠t, we have
(88)M2s,t=1πReηsηtη′tηt−ηs=−1πReη′tηt+1πReη′tηt−ηs.
Applying (85) and the result in (a) yields
(89)M2s,t=−12πcots−t2−1π1−a2sin2t1+a2−1−a2cos2t−12π1−a2sint+s1+a2+1−a2cost+s.Example 5.
To illustrate the accuracy of the two methods, consider the Dirichlet problem inside a domain with a boundary curve parameterized by the ellipse (69) and boundary condition
(90)uηt=γt≔Rea+1+ηtforηt∈Γ.
Note that−a is outside the domain and we can always choose any analytic branch of the function a+1+z.For the sake of comparison, we first solve Mikhlin’s integral equation (25) using the method described in Section 3.2, and then, for a=1.5 and α+iβ∈G with real numbers α and β, we use formula (32) to compute the values at α+iβ (i.e., without computing F+ and without using the barycentric formula (49)). Table 1 shows the estimated relative error un−u∞/u∞ of the approximate solution un at two given interior points. We observe that accurate results are obtained for α=β=0 even for small values of n (see Table 1). For α=0 and β=0.999, the point α+iβ is close to the boundary and the obtained results are very inaccurate even for large values of n. To get an accurate approximation, we compute the boundary values F+ according to equation (35) and then compute the solution at α=0 and β=0,0.999 as described in Section 3.3. The values of the relative error of the computed solution un are presented also in Table 1. It is observed that accurate results are obtained for both cases of β even for small values of n.Table 1
The relative errorun−u∞/u∞ of the approximate solution un of the Dirichlet problem for a=1.5 at the points z=α+iβ with α=0 and two values of β for different values of n using Mikhlin’s integral equation and the integral equation with the gNk (IE with the gNk).
nMikhlin (CIF withh)Mikhlin (CIF withF+)IE with the gNkβ=0β=0.999β=0β=0.999β=0β=0.99983.58×10−031.92×10+028.65×10−058.20×10−071.25×10−041.10×10−06165.71×10−069.59×10+011.00×10−072.64×10−081.62×10−072.76×10−08321.45×10−114.77×10+012.25×10−136.14×10−113.84×10−136.14×10−11642.22×10−162.36×10+011.40×10−165.51×10−162.22×10−165.51×10−161281.40×10−161.15×10+011.40×10−161.38×10−161.40×10−161.38×10−16We solve Example5 using the integral equation with the gNk method, as described in Section 4.3 and summarized in Algorithm 2. The estimated relative error of the approximate solutions at α+iβ for α=0 and β=0,0.999 is presented in Table 1. These results are also highly accurate at both points.In Figure4(a), the relative error of the approximate solution at 1000 random interior points is plotted against the number of discretization points for a=2. Both methods achieve high accuracy, exhibit good efficiency, and converge equally fast. However, in the case of the elongated ellipse with a=19 (see Figure 4(b)), although both methods achieve the same accuracy, there is a difference in terms of the cost. Mikhlin’s integral equation method converges faster and achieves high accuracy at half the number of discretization points needed using the integral equation with the gNk method. This computational advantage in favour of Mikhlin’s integral equation method has been noticed for other values of a for which the ellipse (69) is elongated.Figure 4
The relative errorun−u∞/u∞ of the approximate solution un of the Dirichlet problem at 1000 interior points as a function of the number of discretization points for (a) a=2 and (b) a=19.


## 6. Domains with Complex Geometry

In this section, we perform numerical experiments, aiming to compare the two methods and highlight the differences in their accuracy in case the boundary of the simply connected domainG has a complex geometry. Using Mikhlin’s integral equation method, we follow the steps described in Algorithm 1, i.e., we solve the integral equation (25) first. Next, we use the computed approximate solution to find the boundary values F+ according to (35) and then we evaluate the solution according to (49). For the integral equation with the gNk method, we follow the steps described in Algorithm 2, i.e., we solve the integral equation (63). Then, we evaluate the solution according to (49).We consider two smooth Jordan curves, namely, the boundaryΓ1 parameterized by (see Figure 5(a))
(91)Γ1:ηt=ecostcos28t+esintsin28teit,0≤t≤2π,and the boundary Γ2 with the parameterization (see Figure 5(b))
(92)Γ2:ηt=1+0.5cos30teit,0≤t≤2π.Example 6.
LetG1 and G2 be the two simply connected domains bounded by the two smooth Jordan curves Γ1 and Γ2, respectively. The Dirichlet boundary conditions are constructed from a closed-form reference solution given by
(93)uz=logz−z0,where z is an interior point and z0 is an exterior point. We fix z0=2+1.5i for Γ=Γ1 and z0=1.5+i for Γ=Γ2.Figure 5
The boundaryΓ1 parameterized by (91) (a) and Γ2 parameterized by (92) (b).The approximate solutionun is computed for several values of n at one million randomly chosen interior points using both methods for Γ1 and Γ2. The estimated relative error un−u∞/u∞ is shown in Figure 6. For the curve Γ1, getting a relative error in the order of 10−14 requires around 1000 discretization points using the integral equation with the gNk method and around 1500 using Mikhlin’s integral equation method. For Γ2, reaching a relative error in the order of 10−14 requires around 1800 discretization points using the integral equation with the gNk method and around 4600 discretization points using Mikhlin’s integral equation method. The method based on the integral equation with the gNk is more efficient in both cases.Example 7.
LetG1 and G2 be the two domains described in Example 6. We assume that the exact solution is given by the harmonic function
(94)uz=Re∑j=131z−zj,where z is an interior point and the poles zj are located outside the domain, with z1=2.5−1.5i, z2=2.5+1.5i, and z3=−1.5+2.5i for Γ=Γ1 and z1=1.5+1.5i, z2=−1.5−1.5i, and z3=−1.5+1.5i for Γ=Γ2.Figure 6
The relative errorun−u∞/u∞ of the approximate solution un of the Dirichlet problem as a function of the number of discretization points n in Example 6 for (a) Γ1 and (b) Γ2 using both methods.Both methods are applied to find the approximate solutionun for different values of n at one million randomly chosen points in each interior domain for Γ1 and Γ2. The values of the relative error norm un−u∞/u∞ with the number of discretization points are plotted in Figure 7. Here again, we see that the integral equation with the gNk method requires less discretization points than Mikhlin’s integral equation method would require for the same level of accuracy.Figure 7
The relative errorun−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods in Example 7 for (a) Γ1 and (b) Γ2.The results of Examples6 and 7 show clearly that the method based on the integral equation with the gNk is more efficient for these types of curves with highly varying curvature. The convergence is faster, and the required number of discretization points for an accurate approximation of the solution is less.


## 7. Domains with Piecewise Smooth Boundaries

We consider in this section domains bounded by piecewise smooth curves. Now, in addition to the fact that the integral operatorsN1 and N2 are no longer compact, additional difficulties arise: the solutions of the associated integral equations exhibit a singular behaviour in the neighbourhood of the corner points, the trapezoidal quadrature rule with equidistant nodes looses its accuracy, and the Nyström method produces ill-conditioned linear systems [29, 30]. There are many successful approaches to overcome difficulties associated with corner points (see, e.g., [29–33]).


### 7.1. A Graded Mesh Quadrature

Since we are using the trapezoidal quadrature rule, we use the approach suggested by Kress in [29] which is based on replacing the equidistant nodes by a graded mesh with the same number of nodes constructed by a new variable substitution. This substitution has the particularity of making the derivatives of the new integrand vanish at the extremities of the integration interval. In particular, using this substitution in the parameterization of the curve renders the new transformed parameterization many times continuously differentiable along the whole curve. Kress [29] introduced a typical substitution based on the bijective and strictly monotonically increasing rational function ws:0,2π⟶0,2π defined as
(95)ws=vspvsp+v2π−sp,where p≥2 is the grading parameter and v is a cubic polynomial given by
(96)vs=1p−12π−sπ3+1ps−ππ+12.Notice thatv0=0 and v2π=1 and that w is infinitely differentiable.Assume that the boundaryΓ parameterized by ηt has n corner points. These corner points are at
(97)η0,η2π/n,η4π/n,⋯,η2n−1π/n.Define the functionω as [34]
(98)ωs=wns/n,s∈0,2π/n,wns−2π+2π/n,s∈2π/n,4π/n,wns−4π+4π/n,s∈4π/n,6π/n,⋮wns−2n−1π+2n−1π/n,s∈2n−1π/n,2π.Sincew has a zero of order p at s=0 and s=2π ([29], Thm.2.1), then ω∈Cp. We substitute t=ωs in the parameterization of the boundary and consequently obtain
(99)ηt=ηωs,η′tdt=ω′sη′ωsds.This substitution in the parameterization of the boundary eliminates the complexity arising from corner regions and the integral equations can be solved as in the case of smooth domains.


### 7.2. Numerical Experiments

We are now in a position to concretely use the described substitution technique in numerical examples for solving (2.12) in domains with corners. We consider four examples. The first example is a family of curves with one outward-pointing corner. The second example is a family of curves with one reentrant corner. The third example is a curve with 20 corners, half of them are inward-pointing and the other half are outward-pointing. The fourth example consists of two polygonal domains.Example 8.
LetGθ be the domain bounded by the curve Γθ parameterized by [33, §4.2]
(100)ηt=2sint/2−itanθ/2cost−1,t∈0,2π,where 0<θ<π. The boundary Γθ has a corner at t=0, with interior angle θ. The boundary data are constructed from the harmonic function (93), where z0 is a point outside Gθ. We consider two instances of Γθ, more precisely, Γπ/2 and Γπ/20 (see Figure 8). We set z0=−2 and use both methods to compute the solution of the Dirichlet problem at z=0 for different values of n. Figure 9 compares the convergence results for both methods. For Γπ/2, the two methods are almost equivalent. For Γπ/20, the difference is remarkable, the gNk method is more efficient, and the convergence rate decays linearly. It is similar to convergence rates seen in smooth boundaries. The convergence rate of Mikhlin’s integral equation on the other hand follows a nonlinear and slow pattern.Example 9.
LetΓϕ denote the family of curves parameterized by [33, §4.3]
(101)ηt=−2sin3t/2+itanϕ/2sint−1/2,where π<ϕ<2π is the interior angle. The boundary data are prescribed using the function (93) with z0=−3.5. We consider in this example the two curves Γ3π/2 and Γ2π−π/10 (see Figure 10). We use both methods to compute the solution of the Dirichlet problem at z=0 for different values of n.Figure 8
The curves (a)Γπ/2 and (b) Γπ/20 in Example 8.Figure 9
The relative errorun−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods in Example 8 with p=7 for (a) Γπ/2 and (b) Γπ/20.Figure 10
The curves (a)Γ3π/2 and (b) Γ2π−π/10 in Example 9.The relative error of the obtained solution is computed for different values ofn. The results appear in Figure 11. We notice that both methods are equally efficient for Γ3π/2 with comparable accuracy.Example 10.
We consider the boundary curve parameterized by ([33] §4.2)
(102)ηt=eit4+2cos10t−π20sin10t−π20.Figure 11
The relative errorun−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods in Example 9 with p=7 for (a) Γ3π/2 and p=6 for (b) Γ2π−π/10.The boundary data are constructed from the function (93) with z0=−6. Both methods are used to compute the solution of the Dirichlet problem at 1000 randomly chosen interior points (see Figure 12(a)).Figure 12
The quadrature nodes used to discretize the boundary in Example10 and the 1000 randomly chosen interior points (a). The relative error un−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods with different values of the grading parameter p in Example 10 (b).Figure12(a) compares the convergence results using both methods with different values of the grading parameter p and different numbers of discretization points n. As we clearly see, for all values of p, Mikhlin’s integral equation method stagnates at a lower accuracy. The gNk method is more efficient and accurate at, substantially, any given number of discretization points.Example 11.
In the last example, we consider two polygonal domains, namely, a hexagon and a five-armed star (see Figure13). The two polygonal domains were generated using PlgCirMap (see [14]). The Dirichlet boundary conditions are constructed from the harmonic function (93) with z0=−2 for the hexagon and z0=−1.5 for the five-armed star. The value of the grading parameter used is p=9. Figure 14 shows the estimated relative error of the computed solution at z=0 for different values of n.Figure 13
Two polygonal domains, a hexagon (a) and a five-armed star (b), in Example11.Figure 14
The relative errorun−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods in Example 11 for the hexagon (a) and the five-armed star (b).In both polygonal domains, the method based on the integral equation with the gNk is effectively more efficient, reaching highly accurate results in both domains. The difference between the two methods is more remarkable in the case of the five-armed star (see Figure14).


## 7.1. A Graded Mesh Quadrature

Since we are using the trapezoidal quadrature rule, we use the approach suggested by Kress in [29] which is based on replacing the equidistant nodes by a graded mesh with the same number of nodes constructed by a new variable substitution. This substitution has the particularity of making the derivatives of the new integrand vanish at the extremities of the integration interval. In particular, using this substitution in the parameterization of the curve renders the new transformed parameterization many times continuously differentiable along the whole curve. Kress [29] introduced a typical substitution based on the bijective and strictly monotonically increasing rational function ws:0,2π⟶0,2π defined as
(95)ws=vspvsp+v2π−sp,where p≥2 is the grading parameter and v is a cubic polynomial given by
(96)vs=1p−12π−sπ3+1ps−ππ+12.Notice thatv0=0 and v2π=1 and that w is infinitely differentiable.Assume that the boundaryΓ parameterized by ηt has n corner points. These corner points are at
(97)η0,η2π/n,η4π/n,⋯,η2n−1π/n.Define the functionω as [34]
(98)ωs=wns/n,s∈0,2π/n,wns−2π+2π/n,s∈2π/n,4π/n,wns−4π+4π/n,s∈4π/n,6π/n,⋮wns−2n−1π+2n−1π/n,s∈2n−1π/n,2π.Sincew has a zero of order p at s=0 and s=2π ([29], Thm.2.1), then ω∈Cp. We substitute t=ωs in the parameterization of the boundary and consequently obtain
(99)ηt=ηωs,η′tdt=ω′sη′ωsds.This substitution in the parameterization of the boundary eliminates the complexity arising from corner regions and the integral equations can be solved as in the case of smooth domains.


## 7.2. Numerical Experiments

We are now in a position to concretely use the described substitution technique in numerical examples for solving (2.12) in domains with corners. We consider four examples. The first example is a family of curves with one outward-pointing corner. The second example is a family of curves with one reentrant corner. The third example is a curve with 20 corners, half of them are inward-pointing and the other half are outward-pointing. The fourth example consists of two polygonal domains.Example 8.
LetGθ be the domain bounded by the curve Γθ parameterized by [33, §4.2]
(100)ηt=2sint/2−itanθ/2cost−1,t∈0,2π,where 0<θ<π. The boundary Γθ has a corner at t=0, with interior angle θ. The boundary data are constructed from the harmonic function (93), where z0 is a point outside Gθ. We consider two instances of Γθ, more precisely, Γπ/2 and Γπ/20 (see Figure 8). We set z0=−2 and use both methods to compute the solution of the Dirichlet problem at z=0 for different values of n. Figure 9 compares the convergence results for both methods. For Γπ/2, the two methods are almost equivalent. For Γπ/20, the difference is remarkable, the gNk method is more efficient, and the convergence rate decays linearly. It is similar to convergence rates seen in smooth boundaries. The convergence rate of Mikhlin’s integral equation on the other hand follows a nonlinear and slow pattern.Example 9.
LetΓϕ denote the family of curves parameterized by [33, §4.3]
(101)ηt=−2sin3t/2+itanϕ/2sint−1/2,where π<ϕ<2π is the interior angle. The boundary data are prescribed using the function (93) with z0=−3.5. We consider in this example the two curves Γ3π/2 and Γ2π−π/10 (see Figure 10). We use both methods to compute the solution of the Dirichlet problem at z=0 for different values of n.Figure 8
The curves (a)Γπ/2 and (b) Γπ/20 in Example 8.Figure 9
The relative errorun−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods in Example 8 with p=7 for (a) Γπ/2 and (b) Γπ/20.Figure 10
The curves (a)Γ3π/2 and (b) Γ2π−π/10 in Example 9.The relative error of the obtained solution is computed for different values ofn. The results appear in Figure 11. We notice that both methods are equally efficient for Γ3π/2 with comparable accuracy.Example 10.
We consider the boundary curve parameterized by ([33] §4.2)
(102)ηt=eit4+2cos10t−π20sin10t−π20.Figure 11
The relative errorun−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods in Example 9 with p=7 for (a) Γ3π/2 and p=6 for (b) Γ2π−π/10.The boundary data are constructed from the function (93) with z0=−6. Both methods are used to compute the solution of the Dirichlet problem at 1000 randomly chosen interior points (see Figure 12(a)).Figure 12
The quadrature nodes used to discretize the boundary in Example10 and the 1000 randomly chosen interior points (a). The relative error un−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods with different values of the grading parameter p in Example 10 (b).Figure12(a) compares the convergence results using both methods with different values of the grading parameter p and different numbers of discretization points n. As we clearly see, for all values of p, Mikhlin’s integral equation method stagnates at a lower accuracy. The gNk method is more efficient and accurate at, substantially, any given number of discretization points.Example 11.
In the last example, we consider two polygonal domains, namely, a hexagon and a five-armed star (see Figure13). The two polygonal domains were generated using PlgCirMap (see [14]). The Dirichlet boundary conditions are constructed from the harmonic function (93) with z0=−2 for the hexagon and z0=−1.5 for the five-armed star. The value of the grading parameter used is p=9. Figure 14 shows the estimated relative error of the computed solution at z=0 for different values of n.Figure 13
Two polygonal domains, a hexagon (a) and a five-armed star (b), in Example11.Figure 14
The relative errorun−u∞/u∞ of the computed solution un of the Dirichlet problem versus the number of discretization points using both methods in Example 11 for the hexagon (a) and the five-armed star (b).In both polygonal domains, the method based on the integral equation with the gNk is effectively more efficient, reaching highly accurate results in both domains. The difference between the two methods is more remarkable in the case of the five-armed star (see Figure14).


## 8. Concluding Remarks and Discussion

In this work, we considered the comparison of two integral equation methods for solving the Dirichlet problem for Laplace’s equation in simply connected domains, namely, Mikhlin’s integral equation method and the integral equation with the gNk method. Both integral equations were discretized using the Nyström method and the trapezoidal rule. The resulting linear systems were solved using the MATLAB\ operator. The two integral equation methods are stable and highly accurate and have the same computational complexity.Solving the Dirichlet problem with both methods requires finding the boundary valuesF+ of the analytic function F where the unique solution of the Dirichlet problem in the domain G is uz=ReFz for z∈G. In both methods, once the boundary values F+ are determined, we use the barycentric formula (49) to evaluate the solution at given interior points. The main difference between the two methods though is how we calculate the boundary values F+. For Mikhlin’s integral equation, note that computing F+ by (35) requires computing M1h, where h is an approximate solution of the integral equation (25). On the other hand, solving the integral equation with the gNk provides us directly with the boundary values F+ without any extra calculations. However, the right-hand side of the integral equation with the gNk is −M2γ and needs to be computed first, i.e., we have only a computed approximation of the right-hand side of the integral equation, in contrast to Mikhlin’s integral equation where the right-hand side is given explicitly.To sum up, the two methods are computationally equivalent for computingF+. Both require solving a linear system and both require computing a singular integral Mkϕ (k=1 or 2). However, the function ϕ is a known function (=−γ) for the integral equation with the gNk and ϕ is a computed approximate solution (=h) for Mikhlin’s integral equation.In conclusion, for simply connected domains, the two methods based on these two integral equations are equivalent in terms of computational complexity and accuracy. The numerical examples show that for domains with simple geometry, both methods are highly accurate and exhibit good performance. The method based on Mikhlin’s integral equation is more efficient particularly for elongated ellipses (see Figure4). However, for boundary curves with rapidly varying curvature, the integral equation with the gNk method is more efficient (see Figures 6 and 7). In domains with corners, both methods have shown comparable accuracy (see Figures 9 and 11). However, for domains with several corners, the integral equation with the gNk method has shown better efficiency (see Figures 12 and 14).It is natural to devote further investigation on the comparison between the two integral equations for multiply connected domains. There are, actually, significant differences between the two integral equation methods in multiply connected domains. The integral equation with the gNk (63) is still uniquely solvable as it is [12]. In contrast, Mikhlin’s integral equation (25) is not uniquely solvable [27]. Other reformulations need to be considered in order to make it uniquely solvable (see [2, §31], [7]). A comparison between the two integral equation methods in both bounded and unbounded multiply connected domains will be the subject of a future work. This comparison will include the application of both integral equations in solving some other problems such as the computation of the Dirichlet-to-Neumann map and numerical conformal mappings.



---

*Source: 6488709-2022-04-20.xml*