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.
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. Section 2 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. Let be a bounded simply connected boundary domain with smooth bounding Jordan curve . The boundary is parameterized by a -periodic twice continuously differentiable complex function , , with which traverses in the positive orientation, i.e., and . The positive orientation in this context is that for which is on the left. We let denote the domain exterior to , i.e., , where and . We assume for convenience that the origin of the coordinates is interior to . See Figure 1 for an illustration.
We denote by the space of Hölder continuous -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 defined on the boundary and a function of parameterization of the same boundary defined on . 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. If is an analytic function in and continuous in the closure , then according to the Cauchy integral formula [3, 17]
If, however, is analytic in and continuous in the closure , then
The integral on the left-hand sides of (1) and (2) is known as Cauchy’s integral. For a Hölder continuous function on , the Cauchy-type integral defines a function that is analytic in and in , 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 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
Let be a given continuously differentiable complex function on the boundary such that everywhere. The gNk is defined for by [10]
The kernel is a generalization of the well-known Neumann kernel which is a classical kernel in potential theory obtained by setting [6, p. 282]. Remarkably, the kernel is in fact continuous [10], although the denominator looks problematic. We have
Consequently, the integral operator defined on by is compact. We consider also the following kernel: which is singular and its singular part involves the cotangent function [10]. The kernel can be split as where is Hilbert’s singular kernel ([3, p. 46]; [2, p. 118]) known also as the conjugation kernel [19]
The kernel is continuous and takes on the diagonal the values
The integral operator defined on by is bounded on [10].
The possibility of being eigenvalues of the gNk depends on the index of the function , which is defined as the winding number of with respect to 0, i.e., the change of the argument of over one period divided by [3, 18].
In this paper, we consider two special cases of the gNk, namely, the gNk formed with and the gNk formed with
The kernel is the classical Neumann kernel (see [6, p. 282], and [19, p. 371]). Ultimately, we define the kernels and from (8) with replaced by and , respectively.
We consider two integral equations. The first integral equation is Mikhlin’s integral equation whose kernel is [2]. The kernel of the second integral equation is , and this integral equation is known as the boundary integral equation with the gNk [10, 20].
Note that and . Thus, we have the following theorem (see, e.g., [1, p. 255], [12]).
Theorem 1. (a) is not an eigenvalue of and is a simple eigenvalue of with the constant function as the corresponding eigenfunction.
(b) is not an eigenvalue of and is a simple eigenvalue of 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 , then .
(b) If is an eigenvalue of , then .
(c) If () is an eigenvalue of (), then () is an eigenvalue of ().
The following theorem follows from [12] (Theorem 8).
Theorem 3. If , then is an eigenvalue of the kernel if and only if is an eigenvalue of .
2.4. The Dirichlet Problem
The interior Dirichlet problem for Laplace’s equation is to find a harmonic function which satisfies where is the Laplacian operator . The boundary data is a prescribed Hölder continuous real-valued function on the boundary . The Dirichlet problem has a unique solution ([17], p. 93) which can be regarded as the real part of a single-valued function , analytic in and continuous up to the boundary ([2], p. 137), i.e.,
To ensure uniqueness of the function , we assume that [3, p. 208]
It follows from (16b) that satisfies the boundary condition where is the restriction of the function 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 function . We seek the function in the form of a Cauchy-type integral [2, §29] where is an unknown Hölder continuous real-valued density function defined on the boundary . The task is now reduced to find the density . Using the Sokhotski-Plemelj formulae (4), we obtain
Using the boundary condition (19), we get
Further simplification yields
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 , of the boundary and set to obtain
Equation (24) can be written as where the kernel is the classical Neumann kernel defined in Section 2.3. By Theorem 1, is not an eigenvalue of the kernel . Consequently, by the Fredholm alternative, the integral equation (25) is uniquely solvable. Since , using regularization, the integral equation (25) can be written as
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 ([21, p. 101]).
3.2. The Nyström Method
The integral equation (26) is discretized into a linear system , where 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 integer , the integral in (26) is discretized using the trapezoidal rule with equal weights and equally spaced nodes , . We obtain the semidiscrete equation
To get a fully discrete equation, we require that (27) should hold at the quadrature points. We set , . This results in the system of equations where the term under the summation sign is zero when since the kernel is continuous. Hence, (28) can be written as
We can write (29) in matrix form as where , , , , , is the identity matrix, is an vector of ones, is a diagonal matrix whose diagonal elements are the elements of the vector , and
Note that matrix B has zeros in the main diagonal. Matrix S of the linear system (30) is nonsymmetric, invertible, and dense.
3.3. Computing the Function
Once we solve the linear system (30) and obtain an approximation of the density function , we use it to compute , according to
The solution in the domain can be evaluated as .
It is worth mentioning that as long as does not lie on , the integrand in the middle of equation (32) is a smooth periodic function on . However, when 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 of the analytic function and then use the Cauchy integral formula with singularity subtraction to compute the values of for . This method gives accurate results even when 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 which in view of (22) can be written as
Setting in (34) and using the definition of , we get
Since [12], equation (35) can be written as
Note that the kernel is singular. As in (9), it can be split as where is continuous and is the Hilbert kernel given by (10).
Equation (34) is the same as equation (24) in [5] (Eq. (24)). In [5], the values of 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 . However, for , the integrand involves the derivative of which is computed numerically via -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 .
We substitute , , to obtain
We proceed by discretizing the continuous kernel using the trapezoidal rule and the Hilbert kernel by Wittich’s method [19] to get the fully discrete scheme for assembling , i.e., where the term under the summation sign is zero when since whenever and the kernel is continuous. The matrices and are given by Or where
In light of (41), (42), and (43), equation (39) can be written as or equivalently in matrix form as where , , , and is a vector of ones.
Now that we have accurately computed an approximation of the boundary values , the function can be evaluated at any point in the domain via the Cauchy integral formula
The integrand in (46) has a pole at . This singularity can be removed by dividing the usual Cauchy integral formula for by the same formula for 1, i.e., , and rearranging to obtain [5, 26, 27]
Contrary to the integrand in (46), the integrand in (47) has no pole at and is by consequence an analytic function of whose Cauchy integral must be equal to zero [26]. Using the trapezoidal rule to discretize this integral yields
Solving for in (48) results in the barycentric formula [5, 26, 27]
Formula (49) is numerically stable even when the evaluation point 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:
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 computing , we can assume, without loss of generality, that for some real constant . We define a function in as so that
The function is analytic in and its boundary values satisfy the boundary condition with the function . The problem (52) is a Riemann–Hilbert problem with the coefficient which has the index . 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 will be chosen so that the solvability condition is satisfied.
4.2. Derivation of the Integral Equation
Let be the unique solution of the Riemann–Hilbert problem (52) and
Thus, the boundary values of the function satisfy
Let , since is analytic in G, then by the Cauchy integral formula (1)
Solving (54) for and substituting in (55) yields which implies
Since the function in the numerator of the right-hand side in (57) is analytic for and is counterclockwise oriented, then using the Cauchy integral formula (2) yields
We proceed by taking the limit on both sides of equation (58) and applying the Sokhotski-Plemelj formula (4) to the left-hand side to obtain
Since , multiplying (59) by gives
Using the parameterization , and , we get
Then, using the definitions of the kernels and (see Section 2.3), we obtain
Taking the imaginary parts of both sides in (62) yields 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 . Manifestly, computing the boundary values of requires only computing the function , the solution of the integral equation (63). Indeed, from (51) and (54), we have
The values of for 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] where , , and matrix is defined by
The right-hand side in (65) is given by [20] where , the matrix is defined in equation (42), is a vector of ones, and matrix is defined by
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.
5. Domains Bounded by Ellipses
In this section, we consider the simply connected domain interior to the ellipse with the parameterization for , i.e., the minor radius is and the major radius is (see Figure 2 for ). 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.
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]) where
This is in agreement with Theorem 2(c). Further, since as , the point is an accumulated point of the eigenvalues of the Neumann kernel. In view of Theorems 1 and 3, the eigenvalues of the gNk are
Consequently, in light of Theorem 2, the coefficient matrix of Mikhlin’s integral equation (25) and the coefficient matrix of the integral equation with the gNk (63) have the same eigenvalues and these eigenvalues lie in the interval for sufficiently large values of . The approximate eigenvalues of the coefficient matrices, sorted decreasingly, are given by
Since for small and large , these eigenvalues are clustered in a remarkably symmetric way around . In fact, for small and for sufficiently large values of , most of these eigenvalues are equal to . For elongated domains, for which or and hence near to , we can notice that more eigenvalues are different from . In all cases, the largest eigenvalue of the coefficient matrices is and the smallest eigenvalue is . For , we have and the eigenvalues are too accumulated around (see Figure 3(a)), where the computed smallest eigenvalue is . For the elongated ellipse with , the computed eigenvalues are shown in Figure 3(b). In this case, , and hence, we have several eigenvalues slightly away from where the computed smallest eigenvalue is .
In Figure 3, we present only the eigenvalues of the matrix . However, we chose the values of such that the maximum norm between the approximate eigenvalues of the two matrices and is less than .
5.2. Closed-Form Expressions for the Kernels , , , and
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])
The explicit closed-form expressions for the kernels , , and defined in Section 2.3 are given in the following theorem.
Theorem 4. Let be the ellipse parameterized by , where and . The closed-form expressions for the kernels , , and are, respectively, given by
(a)
(b)
(c)
Proof. (a) For , we have
Applying the sum-to-product trigonometric identities and , (78) becomes
Applying the half-angle identities and to the bracket term in the denominator and taking the real parts yields
Applying the product-to-sum trigonometric identities and to the numerator part of (80) and rearranging, we get
However,
Substituting this result into (71) and rearranging gives
Finally, multiplying both sides of (83) by and simplifying gives (75).
(b) For , we have
Observe that
So
In light of (74) and (85), we get
(c) For , we have
Applying (85) and the result in (a) yields
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 Note that is outside the domain and we can always choose any analytic branch of the function .
For the sake of comparison, we first solve Mikhlin’s integral equation (25) using the method described in Section 3.2, and then, for and with real numbers and , we use formula (32) to compute the values at (i.e., without computing and without using the barycentric formula (49)). Table 1 shows the estimated relative error of the approximate solution at two given interior points. We observe that accurate results are obtained for even for small values of (see Table 1). For and , the point is close to the boundary and the obtained results are very inaccurate even for large values of . To get an accurate approximation, we compute the boundary values according to equation (35) and then compute the solution at and as described in Section 3.3. The values of the relative error of the computed solution are presented also in Table 1. It is observed that accurate results are obtained for both cases of even for small values of .
We solve Example 5 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 for and is presented in Table 1. These results are also highly accurate at both points.
In Figure 4(a), the relative error of the approximate solution at 1000 random interior points is plotted against the number of discretization points for . Both methods achieve high accuracy, exhibit good efficiency, and converge equally fast. However, in the case of the elongated ellipse with (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 for which the ellipse (69) is elongated.
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 domain 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 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 parameterized by (see Figure 5(a)) and the boundary with the parameterization (see Figure 5(b))
Example 6. Let and be the two simply connected domains bounded by the two smooth Jordan curves and , respectively. The Dirichlet boundary conditions are constructed from a closed-form reference solution given by where is an interior point and is an exterior point. We fix for and for .
The approximate solution is computed for several values of at one million randomly chosen interior points using both methods for and . The estimated relative error is shown in Figure 6. For the curve , getting a relative error in the order of requires around 1000 discretization points using the integral equation with the gNk method and around 1500 using Mikhlin’s integral equation method. For , reaching a relative error in the order of 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. Let and be the two domains described in Example 6. We assume that the exact solution is given by the harmonic function where is an interior point and the poles are located outside the domain, with , , and for and , , and for .
Both methods are applied to find the approximate solution for different values of at one million randomly chosen points in each interior domain for and . The values of the relative error norm 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.
The results of Examples 6 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 operators and 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 defined as where is the grading parameter and is a cubic polynomial given by
Notice that and and that is infinitely differentiable.
Assume that the boundary parameterized by has corner points. These corner points are at
Define the function as [34]
Since has a zero of order at and ([29], Thm.2.1), then . We substitute in the parameterization of the boundary and consequently obtain
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. Let be the domain bounded by the curve parameterized by [33, §4.2] where . The boundary has a corner at , with interior angle . The boundary data are constructed from the harmonic function (93), where is a point outside . We consider two instances of , more precisely, and (see Figure 8). We set and use both methods to compute the solution of the Dirichlet problem at for different values of . Figure 9 compares the convergence results for both methods. For , the two methods are almost equivalent. For , 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] where is the interior angle. The boundary data are prescribed using the function (93) with . We consider in this example the two curves and (see Figure 10). We use both methods to compute the solution of the Dirichlet problem at for different values of .
The relative error of the obtained solution is computed for different values of . The results appear in Figure 11. We notice that both methods are equally efficient for with comparable accuracy.
Example 10. We consider the boundary curve parameterized by ([33] §4.2)
The boundary data are constructed from the function (93) with . Both methods are used to compute the solution of the Dirichlet problem at 1000 randomly chosen interior points (see Figure 12(a)).
Figure 12(a) compares the convergence results using both methods with different values of the grading parameter and different numbers of discretization points . As we clearly see, for all values of , 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 Figure 13). The two polygonal domains were generated using PlgCirMap (see [14]). The Dirichlet boundary conditions are constructed from the harmonic function (93) with for the hexagon and for the five-armed star. The value of the grading parameter used is . Figure 14 shows the estimated relative error of the computed solution at for different values of .
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 Figure 14).
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 values of the analytic function where the unique solution of the Dirichlet problem in the domain is for . In both methods, once the boundary values 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 . For Mikhlin’s integral equation, note that computing by (35) requires computing , where 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 without any extra calculations. However, the right-hand side of the integral equation with the gNk is 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 computing . Both require solving a linear system and both require computing a singular integral ( or ). However, the function is a known function () for the integral equation with the gNk and is a computed approximate solution () 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 Figure 4). 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.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare no potential conflict of interests.
Authors’ Contributions
All the authors have contributed equally to this paper.
Acknowledgments
The authors wish to thank Universiti Teknologi Malaysia for supporting this work. This work was supported by the Ministry of Higher Education under Fundamental Research Grant Scheme (FRGS/1/2019/STG06/UTM/02/20). This support is gratefully acknowledged.