On Solving Modified Helmholtz Equation in Layered Materials Using the Multiple Source Meshfree Approach

: This paper presents a study for solving the modified Helmholtz equation in layered materials using the multiple source meshfree approach (MSMA). The key idea of the MSMA starts with the method of fundamental solutions (MFS) as well as the collocation Trefftz method (CTM). The multiple source collocation scheme in the MSMA stems from the MFS and the basis functions are formulated using the CTM. The solution of the modified Helmholtz equation is therefore approximated by the superposition theorem using particular nonsingular functions by means of multiple sources located within the domain. To deal with the two-dimensional modified Helmholtz equation in layered materials, the domain decomposition method was adopted. Numerical examples were carried out to validate the method. The results illustrate that the MSMA is relatively simple because it avoids a complicated procedure for finding the appropriate position of the sources. Additionally, the MSMA for solving the modified Helmholtz equation is advantageous because the source points can be collocated on or within the domain boundary and the results are not sensitive to the location of source points. Finally, compared with other methods, highly accurate solutions can be obtained using the proposed method. the two-dimensional modified Helmholtz equation over arbitrary geometries with doubly and multiply connected domains.


Introduction
The solution of the modified Helmholtz equation plays a crucial role in the science and engineering fields, such as for boundary detection problems [1], water wave problems [2], Cauchy problems [3,4], diffusion equations [5], topological sensitivity analyses [6], and boundary value problems [7]. Over the past 10 years, conventional approaches have been widely adopted for solving the modified Helmholtz equation [8][9][10]. In contrast to the conventional approaches, several boundary collocation meshfree methods that have the basis functions satisfying the partial differential equation have also been proposed, for instance, the boundary knot method (BKM) [11,12], the boundary node method [13,14], the singular boundary method (SBM) [15][16][17], and the regularized meshless method [18,19]. In 2003, Chen and Hon [12] adopted the BKM to solve Helmholtz, modified Helmholtz, and the combination of diffusion and convection problems. The convergence of the BKM for the above equations was also investigated. Later, Chen et al. [15] applied the SBM to solve the modified Helmholtz equation. To examine the accuracy of the SBM, several wave numbers were considered in the numerical model [15]. For the modeling of Helmholtz-type equations in composite layered materials, Bin-Mohsin and Lesnic [20] adopted the method of fundamental solutions (MFS).
The MFS [21,22] is categorized into the indirect boundary element method or the boundary collocation meshfree approach, since only the boundary collocation points are used in the solution process. To avoid the singularity of the fundamental solutions, the MFS must collocate the source points outside the physical domain. Since it is a challenge to identify the location of source points, several attempts regarding this issue have been studied [23,24]. In 2007, Young et al. [25] proposed the modified MFS for finding the solution to the Laplace equation with an arbitrary domain. Later, Liu [26] adopted the multiple-length MFS (MLMFS) to investigate the appropriate placement of the sources for the Laplace equation. By using the MLMFS, the condition number can be reduced and the accuracy of the MLMFS may be increased. In recent years, Shigeta et al. [27] presented the multilayer method of fundamental solutions in conjunction with the weighted greedy QR decomposition (WGQRD) to solve the Laplace equation. As for the WGQRD, it can be used to find the necessary source points and unnecessary ones are removed in the numerical model.
Recently, Ku et al. [28] proposed a novel boundary-type meshfree method which is a hybrid of the MFS as well as the collocation Trefftz method (CTM). It adopts the collocation scheme from the MFS and the matrix form can be constructed in terms of superpositioning the nonsingular basis function from the CTM. This approach was then used for solving the Laplace-type subsurface flow problem. Xiao et al. [29] further applied the method to the problems of heat transfer in heterogeneous multilayer materials in two dimensions. Since the hybrid boundary-type meshless methods are only found to solve Laplace-type problems, it is especially interesting to explore this newly developed method for solving other non-Laplace problems. To the authors' knowledge, pioneering work using the new boundary-type meshfree method to solve the modified Helmholtz equation in layered materials has not yet been proposed. We therefore propose a new boundary-type meshfree method using the basis functions based on the first-kind modified Bessel function capable of solving the twodimensional modified Helmholtz equation over arbitrary geometries with doubly and multiply connected domains.
Starting from previous studies [28,29], we conducted pioneering work to solve the modified Helmholtz equation in layered materials. Since it is necessary to derive particular nonsingular basis functions for the two-dimensional modified Helmholtz equation, the specific method in this study is named the multiple source meshfree approach (MSMA). The key idea of the MSMA starts with the MFS and the CTM. The multiple source collocation scheme in the MSMA stems from the MFS, and the basis functions are formulated from the CTM. The solution of the modified Helmholtz equation is therefore approximated by the superposition theorem using particular nonsingular functions by means of multiple sources located within the domain. To deal with the two-dimensional modified Helmholtz equation in layered materials, the domain decomposition method (DDM) [30,31] was adopted. Numerical examples were carried out to validate the method.
The article is organized as follows: A brief introduction of the two-dimensional modified Helmholtz equation and the formulation of the MSMA considering the domain decomposition method are given in Section 2. In Section 3, several numerical examples are included to evaluate the accuracy of the proposed approach by comparing it with exact solutions. In Section 4, we discuss the research. Findings are summarized in Section 5.

The Methodology
We consider a two-dimensional domain for the problem of interest in 2 ℜ and assume that Ω is enclosed by a boundary Γ Ω = ∂ . The two-dimensional modified Helmholtz equation and the corresponding boundary conditions are stated as The Dirichlet and Neumann conditions are in which 2 ℜ denotes the two-dimensional domain, 2 ∇ denotes the Laplacian, u may represent temperature (or acoustic pressure), Ω represents the domain of the problem, λ is the wave number, g and f are the Dirichlet and Neumann boundary data, D Γ and f Γ denote the boundaries where the Dirichlet and Neumann conditions are given, and n represents the outward normal direction. The idea of this study was originally motivated by the MFS and the CTM, in which the collocation scheme is from the MFS and the basis functions are from the Trefftz formulation. To elaborate the proposed MSMA, the collocation scheme is introduced first. Figure 1 demonstrates the collocation points for the MFS, CTM, and MSMA bounded by simply and multiply connected domains. As illustrated in Figure 1a, a source point is located within the physical boundary enclosed by a simply connected domain in the CTM. Instead of using just a source point, the concept of the MFS is to approximate the solutions via a fundamental solution by means of locating the sources away from the physical boundary, as illustrated in Figure 1b. As for the MSMA, the solution is approximated through the linear combination of the particular nonsingular functions in terms of multiple sources collocated within the physical boundary, as illustrated in Figure 1c. For a multiply connected domain, which includes more than one cavity, multiple sources are placed within the physical boundary. The collocation scheme of the sources is demonstrated in Figure 1d-f for the CTM, MFS, and MSMA. From Figure 1, it is obvious that the proposed MSMA is advantageous because we may place the source and boundary points at the same position. Also, unlike the MFS, the source points can be placed at any position and are not sensitive to the results.  In the MSMA, the solution of the modified Helmholtz equation is approximated by the superposition theorem using nonsingular general solutions by means of multiple sources. Numerical solutions of the modified Helmholtz equation in two dimensions may be expressed as the combination of nonsingular general solutions. We may rewrite the two-dimensional modified Helmholtz equation in the polar coordinate as follows: where r represents the radial coordinate and θ is the angular coordinate. For Equation (4), particular nonsingular functions can be derived using the separation of variables. Since the solutions are expressed as the linear superposition of the nonsingular functions, we may obtain the numerical solution of the two-dimensional modified Helmholtz equation represented by the linear superposition of the nonsingular basis functions: where k I is the modified Bessel function of order k for the first kind; j ρ denotes the distance between the boundary collocation point and the source point, in which Q l 1,..., = , Q represents the boundary point number, and ) ( l g x represents the Dirichlet boundary data given on the boundary points. From Equation (3), the boundary of the Neumann condition is rewritten as where x n and y n at l x denote the outward normal direction vector. Using Equation (5) Substituting Equations (9) and (10) into Equation (8), the following equation is listed below: Neumann boundary data are also imposed on the domain boundary using Equation (8): where Q l 1,..., = , and ) ( l f x represents the Neumann boundary data given on the boundary points. Using Equations (7) and (14), the following equation can be obtained as II  ,  I  ,  II  1  ,  I  1  ,   II  ,  I  ,  II  1  ,  I  1  ,   II  ,  I  ,  II  1  ,  I 1 , and  = = + To deal with the two-dimensional modified Helmholtz equation in layered materials, the DDM [30,31] was adopted. It should be noted that the DDM was also used to deal with the problems bounded by doubly and multiply connected domains in this study because the proposed MSMA approximates the solution of the two-dimensional modified Helmholtz equation as the combination of only particular nonsingular basis functions.
Unlike the CTM or the MFS, the MSMA does not consider the negative basis functions or the socalled singular basis functions in the proposed method. As a result, we are not able to tackle the singular solution in the cavity by using only the nonsingular basis functions. It is therefore necessary to split the domain into several simply connected subregions, such that we may guarantee subregions existing only in a nonsingular solution.
In the DDM, the physical boundary is split into several simply connected subregions which are connected at the interface between two consecutive subregions, as depicted in Figure 2. At the interface, the overspecified boundary data need to be satisfied. A doubly connected domain, Ω , is split into two subregions, 1 Ω Ω , we may use both the Dirichlet and Neumann boundary conditions to satisfy the continuity condition on the interface. The form of the continuity equations is described as follows: Using the DDM, we can obtain the matrix form, which is expressed as follows:

Numerical Examples
To validate the MSMA for solving the modified Helmholtz equation, six examples were conducted to verify the proposed meshfree approach. The absolute error and relative root-meansquare error were used to evaluate the accuracy of the numerical solution as follows: where nt is the number of inner points in the domain, and

Modeling of Modified Helmholtz Equation Bounded by a Simply Connected Region
This example is the modeling of the modified Helmholtz equation bounded by a simply connected region with a gear-shaped domain, as depicted in Figure 3. The shape parameter of the physical boundary is given by where ( ) In this case, the Dirichlet and Neumann boundary data are imposed using the exact solution.
The wave number λ is set as 2 . Figure 3 illustrates the maximum absolute error (MAE) versus the terms of the particular solutions, k . The MSMA can obtain highly accurate results in the order of -11 10 while 15 = k . Figure 4 demonstrates the relationship of the MAE with the source point number. When the source points number is more than 100, the proposed MSMA can obtain accurate results with the order of -10 10 . As shown in Figure 5, we found that the accuracy of the results decreased with the increase of the wave number λ . To obtain the MAE, 1211 inner points were uniformly located within the domain. Figure 6 shows the comparison of the numerical and analytical solutions. The results demonstrate that the solutions of the MSMA agree very well with the exact solution.

Accuracy Comparison of the Proposed Method
The second case is the modeling of the modified Helmholtz equation bounded by an amoebashaped boundary, as illustrated in Figure 7. In this example, the results computed by MSMA were compared with those from the SBM [15]. The shape parameter of the irregular physical boundary is described by In the numerical implementation, the terms of the particular solutions, k , was 15. The Dirichlet boundary data were imposed on the domain boundary using three different analytical solutions, as depicted in Equations (27)- (29). The wave number λ was set as π . The computed results using the MSMA were compared with the SBM [15]. We obtained highly accurate numerical solutions while the boundary and source point number was greater than 100, as shown in Figure 8a. Compared with the results using the SBM [15] as shown in Figure 8b, we found that the accuracy of the MSMA was much more accurate than the SBM. Figure 9a shows the accuracy for the relative root-mean-square error versus the wave number. It was found that the accuracy of the increased wave number may decrease the accuracy of the computed results. Figure 9b shows the relative root-mean-square error versus the wave number computed by the SBM [15]. Figure 9 shows that the accuracy of the MSMA is better than that from the SBM.

Investigation of the Wave Number
The third case is the modeling of the modified Helmholtz equation bounded by a peanut-shaped boundary, as depicted in Figure 10. To examine the effect of the wave number for the numerical results, different values of the wave number were considered in the example. The shape parameter of the irregular physical boundary is described by To solve the modified Helmholtz equation, the exact solution can be found as (31) The Dirichlet boundary values were imposed on the domain boundary using the exact solution. We investigated the accuracy of the proposed method by using different values of the wave number ranging from 1 to 50. Figure 10 shows the accuracy for the MAE versus the wave number using the proposed meshfree approach. As illustrated in Figure 10, the accuracy can be found in the order of

Solution of the Modified Helmholtz Equation Bounded by a Doubly Connected Region
The fourth case is the modeling of the modified Helmholtz equation bounded by a doubly connected region, as depicted in Figure 11. The shape parameter of the concentric circles is described by The exact solution of this example can be found as follows: where 3 I is the modified Bessel function of the third order for the first kind, and 3 K is the modified Bessel function of the third order for the second kind. In this example, the terms of the particular solutions were set to be 15. The wave number was set as 2 . For this example, the physical boundary, enclosed by a doubly connected region, was split into two simply connected subregions using the DDM [30,31], as depicted in Figure 11. Four hundred boundary and source points were interactively arranged on the physical boundary. Figure 12 shows that the sub-boundaries at the interface contained four sub-boundaries: 1 Γ , 3 Γ , 5 Γ , and 7 Γ . The overspecified boundary conditions enforced on the boundary and source points at the interface are described by To examine the accuracy, 2415 inner points were placed within the domain to obtain the numerical results. We also compared the numerical solutions with the above exact solution, as illustrated in Figure 13. It seems that the numerical solutions perfectly agree with the exact solution. Furthermore, the MAE of the MSMA can be obtained in the order of -6 10 .  y (m)

Solution of Modified Helmholtz Equation Bounded by a Multiply Connected Region
The fifth case is the modified Helmholtz equation bounded by a multiply connected region, as depicted in Figure 14. For the multiply connected domain including more than one cavity, we must place multiple source points within the physical boundary. Because the proposed MSMA does not consider the negative basis functions or the so-called singular basis functions, we cannot tackle the singular solution in the cavity by using only the nonsingular basis functions. Accordingly, to deal with problems bounded by multiply connected domains, the DDM was adopted, which split the domain into several simply connected subregions, such that we may guarantee subregions existing only in a nonsingular solution. The shape parameter of the physical boundary is given by The exact solution can be found as follows: where 1 I is the modified Bessel function of the first order for the first kind, and 1 K is the modified Bessel function of the first order for the second kind. In this example, the term of the particular solutions was 15. The wave number was set as 2 . For this example, the physical boundary, enclosed by a doubly connected region, was split into two simply connected subregions using the DDM [30,31], as depicted in Figure 14. Four hundred boundary and source points were interactively arranged on the physical boundary. Figure 15 shows that the sub-boundaries at the interface contained six sub-boundaries: 2 Γ , 4 Γ , 6 Γ , 8 Γ , 10 Γ , and 12 Γ . The overspecified boundary conditions enforced on the boundary and source points at the interface are described by To examine the accuracy, 2737 inner points were placed within the domain to obtain the numerical results. We also compared the numerical solutions with the exact solution, as illustrated in Figure 16. The results demonstrate that the numerical solutions perfectly agree with the exact solution. Furthermore, the MAE of the MSMA can be obtained in the order of

Solution of Modified Helmholtz Equation in Two Layered Materials
The last example is the modeling of the modified Helmholtz equation in two layered materials, as depicted in Figure 17. The wave number in two subdomains were set to be 1 λ and 2 λ , in which 2 1 = λ and 1 2 = λ , respectively. The shape parameter of the boundary is described by The exact solutions are as follows: In this example, the terms of the particular solutions were 10. For modeling the modified Helmholtz equation in two layered materials, the physical boundary was split into two simply connected subregions using the DDM [30,31], as illustrated in Figure 17. For each subregion, 600 boundary and source points were interactively arranged on the physical boundary. Figure 17 also shows that the sub-boundaries at the interface contained two sub-boundaries: 1 Γ and 5 Γ . The overspecified boundary conditions enforced on the boundary and source points at the interface are described by To evaluate the accuracy of the problem, 900 inner points were arranged within the domain to obtain the numerical results. We also compared the numerical solutions with the exact solution, as depicted in Figure 18. The results demonstrate that our results agree well with the above exact solution. Furthermore, we found that the MAE of the MSMA can be obtained in the order of -7 10 , as depicted in Figure 19.    Figure 19. The absolute error.

Discussion
In this paper, we proposed a multiple source meshfree approach for solving a two-dimensional modified Helmholtz equation in layered materials. The novel conception of the proposed MSMA started with the MFS and the CTM. The advantages of the proposed method are as follows.
The proposed method for solving the modified Helmholtz equation avoids a complicated procedure for finding the location of the sources. By adopting the nonsingular functions, the source points can be collocated on or within the domain boundary. Compared with other methods, highly accurate solutions can be obtained regardless of the location of source points. We also demonstrated that the proposed method incorporated with the DDM can be used to solve the modified Helmholtz equation in layered materials as well as the domain bounded by simply and multiply connected regions.
We investigated the accuracy of the proposed method by using different values of the wave number ranging from 1 to 50 for the modified Helmholtz equation. The accuracy of the MSMA was found in the order of -9 10 , while the wave number was 50. The computed results revealed that the MSMA may yield highly accurate results even when a high wave number is considered. The proposed method is, however, limited for solving linear governing equations due to the approximation of the solution using the addition theorem.

Conclusions
In this study, we presented a meshfree approach for solving the modified Helmholtz equation bounded by simply, doubly, and multiply connected regions using the MSMA in two dimensions. The modified Helmholtz equation in two layered materials was also considered. The method was verified and numerical examples were also performed. We may point out the two major achievements of the paper: 1. The key idea of the MSMA stems from the indirect boundary element method, which adopts multiple source points. Because of the adoption of nonsingular functions, the sources can be collocated on or within the domain boundary without using a complicated searching algorithm to find the appropriate location of the source points. 2. To the best of our knowledge, the MSMA using nonsingular basis functions is newly developed.
A pioneering work for solving the modified Helmholtz equation bounded by a multiply connected region was conducted using the MSMA in this study. Furthermore, the MAE of the proposed approach for the modified Helmholtz equation in two layered materials can reach up to the order of -7 10 . From the computed results, we conclude that the MSMA is relatively simple because it avoids a complicated procedure for finding the appropriate location of the source points. Moreover, the MSMA has advantages of highly accurate and boundary collocation only for solving problems with complex geometry.