Inﬁnitely Smooth Polyharmonic RBF Collocation Method for Numerical Solution of Elliptic PDEs

: In this article, a novel inﬁnitely smooth polyharmonic radial basis function (PRBF) collocation method for solving elliptic partial differential equations (PDEs) is presented. The PRBF with natural logarithm is a piecewise smooth function in the conventional radial basis function collocation method for solving governing equations. We converted the piecewise smooth PRBF into an inﬁnitely smooth PRBF using source points collocated outside the domain to ensure that the radial distance was always greater than zero to avoid the singularity of the conventional PRBF. Accordingly, the PRBF and its derivatives in the governing PDEs were always continuous. The seismic wave propagation problem, groundwater ﬂow problem, unsaturated ﬂow problem, and groundwater contamination problem were investigated to reveal the robustness of the proposed PRBF. Comparisons of the conventional PRBF with the proposed method were carried out as well. The results illustrate that the proposed approach could provide more accurate solutions for solving PDEs than the conventional PRBF, even with the optimal order. Furthermore, we also demonstrated that techniques designed to deal with the singularity in the original piecewise smooth PRBF are no longer required.


Introduction
A large number of multidisciplinary problems may require solutions involving complicated mathematical models and numerical techniques [1,2]. Among these approaches, the meshfree method has gained the attention of researchers from different scientific fields due to its capability of dealing with partial differential equations (PDEs) with complicated and irregular geometry [3][4][5][6][7]. The radial basis function collocation method (RBFCM) is one meshfree approach for analyzing governing equations, where the unknowns are represented by a function approximation [8][9][10]. The RBFCM is also capable of solving linear and nonlinear PDEs requiring solutions involveing sophisticated mathematical models and numerical algorithms [8,9]. Function approximation using the interpolation method of the radial basis function (RBF) was proposed in 1971 by Hardy [11], who used the multiquadric (MQ) RBF for scattered data interpolation to solve problems. In addition to the MQ RBF, several RBFs may be found, such as the inverse multiquadric (IMQ), Gaussian, and polyharmonic spline (PS) functions [12][13][14]. Among them, the PS and the MQ RBFs may provide much more accurate solutions than other RBFs [15,16]. The RBFs are usually categorized into piecewise and infinitely smooth functions for solving PDEs in the RBFCM. For example, the PS is piecewise smooth. On the other hand, the MQ is infinitely smooth, with the shape parameter [17,18]. For the PS, a minimum order of two must be adopted to keep it smooth for solving Laplace-type equations [19]. Despite the great success of the above RBFs as effective numerical techniques for dealing with several kinds of PDEs, there is still growing interest in the application and development of new and advanced RBFs [20]. A significant number of modifications to RBFs have been proposed, such as the pseudo-spectral RBF [21,22], Gaussian RBF [23], RBF QR alternative basis method [24], finite difference RBF [25,26], partition of unity RBF [27,28], stabilized expansion of the Gaussian RBF [29], rational RBF [30,31], and RBF based on partition of unity of Taylor series expansion [32].
An adaptive greedy technique for the Kansa method is adopted for general PDE problems to be solved on complicated domains [5]. The adaptive greedy algorithm prefers exterior centers for large scaling parameters. Similar to the idea of exterior centers, Ku et al. presented a fictitious source collocation scheme where the centers are replaced by fictitious sources collocated outside the domain as well [33,34]. However, different from the greedy algorithm, the role of fictitious sources is emphasized on the prevention of singularity without the shape parameter for the MQ. Since the distance between the source and the interior point is always greater than zero, the MQ and its derivatives are always smooth and globally infinitely differentiable without using the shape parameter.
Motivated by this concept, a novel infinitely smooth polyharmonic radial basis function (PRBF) collocation method for solving elliptic partial differential equations (PDEs) is presented. The seismic wave propagation problem, groundwater flow problem, unsaturated flow problem, and groundwater contamination problem were investigated to reveal the robustness of the proposed PRBF. Comparisons of the conventional PRBF with the proposed method were carried out as well. The formulation of the PRBF is introduced in Section 2. The convergence analysis is presented in Section 3. Section 4 details the results of the investigation into the geomechanical problems listed above. The findings of this research are finally summarized in Section 5.

The Governing Equation
The governing equation for elliptic partial differential equations in two dimensions can be expressed as where ∇ denotes the divergence operator; · denotes the dot product; u(x) denotes the unknown, x = (x, y); Ω denotes a domain with the boundary ∂Ω; g(x) denotes the boundary value; A is the velocity in two dimensions, A = A x , A y ; B is the given function; and C(x) is the given function value.
In the RBFCM, several RBFs can be implemented, such as the MQ, Gaussian, IMQ, and PS functions. However, the accuracy of the Gaussian, MQ, and IMQ RBFs may depend on the shape parameter. Thus, most applications of the Gaussian, MQ, and IMQ RBFs utilize optimization approaches or experimental techniques to determine the optimal shape parameter [17,18]. Since the PS may provide good agreement without providing an additional parameter to prevent singularity, we considered the PRBF in this study. The PRBF with natural logarithm is defined as ϕ(r) = r 2k ln(r), k = 1, 2, 3 . . .
where k is the order; r is the radial distance, r = |x − x s |; x denotes the interior point; and x s denotes the center point. The equation is a kth order PRBF.
Taking the derivative of (3) with respect to x gives where x s is the x-coordinate of the center point. We may then take the derivative of y and obtain ∂ϕ(r) ∂y = (y − y s )r 2k−2 [2kln(r) where y s is the y-coordinate of the center point. Similarly, we have Again, we take the derivative of y and obtain Utilizing the PRBF to build up the function approximation, we obtain where n c denotes the total number of center points; λ j denotes the unknown coefficient; ϕ(r j ) denotes the basis function; r j is the radial distance at the jth center point, r j = x − x s j ; and x s j denotes the jth center point, described as x s j = x s j , y s j . Inserting the aforementioned equations into (1), the function approximation of the governing equation is as follows: From (9), it is clear that ln(r j ) becomes singular when r j = 0, where the position of the center point exactly coincides with the position of the interior point. In this study, we propose a novel idea to convert the piecewise smooth PRBF into an infinitely smooth PRBF using source points collocated outside the domain to ensure that the radial distance is always greater than zero, thus avoiding the singularity of the conventional PRBF. Accordingly, the terms in (9) are always smooth, and no extra techniques are required to deal with the singularity. The concept of fictitious sources is elaborated in the following section.

The Fictitious Sources
In the conventional PRBF, the interior, center, and boundary points have to be placed where the positions of the interior and center points are usually the same, as displayed in Figure 1a. To prevent singularity, the PRBF with natural logarithm is usually implemented as follows [15]: Using (8), we may obtain the function approximation for the governing equation again where n s denotes the number of fictitious sources. To prevent singularity, a minimum order of two must be utilized to ensure that r 2k−3 j in the above equation is a smooth and completely monotone RBF. Accordingly, the PRBF with natural logarithm is a piecewise smooth function in the conventional RBFCM. In this study, fictitious sources are adopted. The center points in the PRBF are regarded as the fictitious sources, where the fictitious sources are placed outside the domain, as depicted in Figure 1b.
The following equation is utilized for placing the source points.
x s j = (ηρ s j cos θ s j , ηρ s j sin θ s j ), j = 1, . . . , n s (12) where η denotes the dilation parameter; ρ s j denotes the radius of the fictitious sources at the jth source point; and θ s j denotes the angle of the fictitious sources at the jth source point. Figure 1b depicts the proposed fictitious sources. As shown in Figure 1b, the positions of the source points can, therefore, be collocated with the consideration of different values of the dilation parameter by utilizing (12). Since the radial distance is always nonzero, the PRBF and its derivatives in the governing equation, as depicted in (3) and (9), are always continuous. Based on the proposed fictitious sources, the original piecewise smooth PRBF is then converted into an infinitely smooth PRBF.
To investigate the unknown coefficients, approximations are applied to satisfy (1) with the boundary conditions, as shown in (2), at each point. The following linear system is acquired: where A L denotes a n i × n s matrix for the interior points; A B denotes a n b × n s matrix for the boundary points; λ denotes a n s × 1 vector which includes unknown coefficients λ j to be solved; f denotes a n i × 1 vector of function values for the interior points, written as f = [C 1 , C 2 , . . . , C n i ]; g denotes an n b × 1 vector of boundary data, written as g = g 1 , g 2 , . . . , g n b ; n i denotes the number of interior points; and n b denotes the number of boundary points. The undetermined coefficients can be obtained by solving a system of simultaneous linear equations in matrix form as shown in (13). Accordingly, the solution of u(x) can be obtained using the function approximation of (8). To evaluate the error, the root-mean-square error (RMSE) and the maximum absolute error (MAE) are both adopted. The accuracy of the computed results is compared to the exact solutions using the following equations.
where I is the total number of interior points, RMSE is the root-mean-square error, MAE is the maximum absolute error, and u(x i ) and u E (x i ) represent the numerical and exact solutions at the ith interior point, respectively.

Convergence Analysis
To compare the proposed PRBF with the conventional PRBF, we firstly consider a steady-state groundwater flow problem [33]: An irregular boundary was defined as The boundary data, as depicted in (18), were assigned on the irregular boundaries.
We used 321 interior points, 396 fictitious sources, and 75 boundary points. In this example, we consider six cases with different shapes of fictitious sources using (19) to (24) as demonstrated in Figure 2a-f for Cases I, II, III, IV, V, and VI, respectively.  The fictitious sources were placed on the external domain by utilizing (25). The fictitious source boundaries were defined as [33,34] ∂Ω s = (x s j , y s j ) x s j = ηρ s j cos θ s j , y s j = ηρ s j sin θ s j (25) where ∂Ω s denotes the boundaries of fictitious sources. Figure 3 illustrates the MAE values of the conventional and proposed PRBFs with different orders (values of k). In Case I, the MAE tended to increase with an increase in the order k, as displayed in Figure 3a. In the other five cases, it was found that the MAE of the proposed PRBF remained lower than that of the conventional one with different orders, as displayed in Figure 3b-f. The accuracy is less sensitive to the order k. From Figure 3, accurate results with order k in the range of one to ten were acquired. Additionally, the results demonstrate that the proposed PRBF with an order of one provides more accurate results than the conventional one, as depicted in Table 1. The MAE of the conventional PRBF was in the order of 10 −1 to 10 −4 as the order k ranged from two to ten. It appears that the MAE for the conventional PRBF could only reach the order of 10 −4 when the order k was set to seven, as displayed in Figure 3.     The MAE of the proposed PRBF with different dilation parameters was investigated, as displayed in Figure 4. In Cases I, II, III, and IV, the MAE declined with an increase in the dilation parameters-thus, all cases could obtain accurate results. As depicted in Figure 4e,f these results show that accurate results in the order of 10 −6 to 10 −11 were yielded when an arbitrary dilation parameter in the range of zero to ten was utilized.    It is found that the error may grow with the increase of the dilation parameter. However, results from all cases reveal that we may obtain accurate results with the use of an arbitrary value in the range of 1.5 to 5 for the dilation parameter. Nevertheless, the error may increase while the larger dilation parameter is adopted. However, it reveals that we could still obtain better accuracy than original PS even in the worst case.

Numerical Examples
To demonstrate the applicability of the proposed method, several types of geomechanical problems-the seismic wave propagation problem, groundwater flow problem, unsaturated flow problem, and groundwater contamination problem-were investigated; the results are presented in this section.

The Seismic Wave Propagation Problem
The problem of seismic wave propagation [35] is considered as where λ 2 denotes the wave number depending on the property of the geomaterial. The boundary is Data are assigned on the boundaries by the exact solution: g(x) = e (λ/2)(x+ √ 3y) .
In this wave problem, we assigned λ 2 = 1 and λ 2 = 0.0001. The conventional PRBF and the proposed PRBF with fictitious sources were utilized to solve this problem. In the conventional PS RBF, the order k ≥ 2 was adopted to prevent the singularity. In the proposed PRBF with fictitious sources, the irregular boundaries of the exterior sources, ∂Ω s , defined by (24), were utilized in this example. Figure 5 shows the original collocation points for the conventional case. We used 480 center points and 75 boundary points. For the fictitious sources, demonstrated in Figure 6, the interior points, boundary points, and fictitious sources were placed within, on, and outside the region, respectively. The numbers of source, boundary, and interior points were 480, 75, and 405, respectively.   Figure 7a shows the order k versus the MAE when the wave number was 0.0001. From Figure 7a, the MAE of the conventional PRBF fluctuated between 10 −1 and 10 −4 as the order k ranged from two to ten. It appears that the proposed approach with order of one acquired the most accurate results, with MAE in the order of 10 −13 -much better than the conventional PS RBF. Figure 7b shows the order k versus the MAE when the wave number was 1. From Figure 7b, the MAE of the conventional PS RBF fluctuated between 10 −1 and 10 −4 as the order k ranged from two to ten. However, the MAE of the proposed PRBF utilizing fictitious sources drastically decreased from 10 2 to 10 −10 as the order k ranged from one to ten. It is clear that the proposed PRBF with fictitious sources can provide significantly more accurate results than the conventional PRBF. Figure 8 shows the MAE versus the dilation parameter. The results show that the locations of the fictitious sources are of low sensitivity with respect to the solutions.

The Groundwater Flow Problem
The Poisson equation [36] is referred to for groundwater flow with recharge, defined as where p denotes the parameter depending on the recharge; values of p of 1 and 0.0001 were considered in this problem. The boundary is The Dirichlet data are g(x) = y sin(px) + x cos(py).
The conventional PRBF and the proposed PRBF with fictitious sources were adopted to solve this problem. In the conventional PRBF, the order k ≥ 2 was adopted to prevent the singularity. In the proposed PRBF with fictitious sources, irregular boundaries of the exterior sources, ∂Ω s , defined by (24), were utilized in this example. Figure 9 depicts the collocation points for the conventional case. We considered 382 center points and 75 boundary points. For the proposed fictitious sources, as demonstrated in Figure 10, we assigned the interior, boundary, and fictitious sources within, on, and outside the region, respectively. The numbers of fictitious source, boundary, and interior points were 457, 75, and 382, respectively.   Figure 11a shows the order k versus the MAE when value of p was set to 0.0001. From Figure 11a, the MAE of the conventional PRBF fluctuated between 10 −1 and 10 −5 as the order k ranged from two to ten. However, it was found that the proposed method with order of one acquired significantly more accurate results than the conventional method. Figure 11b shows the order k versus the MAE when the value of p was set to 1. From Figure 11b, the MAE of the conventional PRBF fluctuated between 10 −1 and 10 −4 as the order k ranged from two to ten. However, the MAE of the proposed PRBF with fictitious sources drastically decreased from 10 1 to 10 −10 as the order k ranged from one to ten. It is clear that the proposed PRBF with fictitious sources can provide more accurate results than the conventional one. Figure 12 shows the MAE versus the dilation parameter. The results show that the locations of the fictitious sources may not be sensitive to the solutions.

The Unsaturated Flow Problem
The third example is the stationary convection-diffusion problem, usually referred to as unsaturated flow in two dimensions [1].
Here, K r (h) denotes the relative hydraulic conductivity, and h denotes the pressure head.
The following equation is obtained by utilizing the Gardner model [1]: where α g denotes the pore size distribution parameter;ĥ denotes the linearized pressure head, expressed asĥ = e α g h − e α g h d ; and h d denotes the dry pressure head. The irregular boundary is The boundary data are g(x) = xe −α g y .
In this example, soil pore size distribution values of α g = 0.0001 and α g = 1 were considered. The conventional PRBF and the proposed PRBF with fictitious sources were implemented to solve this problem. In the conventional PRBF, the order k ≥ 2 was adopted to prevent the singularity. In the proposed method with fictitious sources, irregular boundaries of the exterior sources, ∂Ω s , defined by (24), were utilized in this example. Figure 13 shows the collocation points for the conventional case. We used 502 center points and 75 boundary points. As demonstrated in Figure 14, we collocated the interior, boundary, and fictitious sources within, on, and outside the domain, respectively. The numbers of source, boundary, and interior points were 502, 75, and 427, respectively.   Figure 15a shows the order k versus the MAE when α g = 0.0001. From Figure 15a, the MAE of the conventional PRBF fluctuated between 10 −1 and 10 −4 as the order k ranged from two to ten. However, it was found that the proposed method with order of one acquired significantly more accurate results than the conventional method. Figure 15b shows the order k versus the MAE when the pore size distribution parameter was set to 1. From Figure 15b, the MAE of the conventional PS RBF fluctuated between 10 −1 and 10 −4 as the order k ranged from two to ten. The MAE of the proposed method with fictitious sources drastically decreased from 10 −1 to 10 −9 as the order k ranged from one to ten. It is clear that the proposed PRBF with fictitious sources can provide more accurate results than the conventional one. Figure 16 shows the MAE versus the dilation parameter. The results show that the positions of the fictitious sources may not be sensitive to the solutions.
The conventional PRBF and the proposed PRBF with fictitious sources were implemented to solve this problem. In the conventional PRBF, the order k ≥ 2 was utilized to prevent the singularity. In the proposed method, irregular boundaries of the exterior sources, ∂Ω s , defined by (24), were utilized in this example. Figure 17 shows the collocation points for the conventional case. We used 369 center points and 75 boundary points. For the fictitious sources, demonstrated in Figure 18, we assigned the interior, boundary, and source points within, on, and outside the domain, respectively. The numbers of source, boundary, and interior points were 444, 75, and 369, respectively.   Figure 19 shows the order k versus the MAE. From Figure 19, the MAE of the conventional PRBF fluctuated between 10 −1 and 10 −4 as the order k ranged from two to ten. However, the MAE of the proposed method utilizing fictitious sources drastically decreased from 10 −1 to 10 −8 as the order k ranged from one to ten. It is clear that the proposed PRBF with fictitious sources can provide more accurate results than the conventional PRBF. Figure 20 shows the MAE versus the dilation parameter. The results show that the positions of the fictitious sources may not be sensitive to the solutions.

Conclusions
An infinitely smooth polyharmonic radial basis function (PRBF) collocation method for solving elliptic partial differential equations (PDEs) was developed in this study. The significance of this study can be summarized as follows: (1) We convert the piecewise smooth PRBF into an infinitely smooth PRBF using source points collocated outside the domain. The distance between the interior point and fictitious source is, therefore, greater than zero. The PRBF and its derivatives in the governing equation always have the characteristics of global infinite differentiability and smoothness. Accordingly, techniques to deal with the singularity in the original PRBF are no longer required.
(2) To clarify the possible influences for the positions of the sources on the accuracy, an analysis was conducted to evaluate the accuracy versus the shapes and positions of sources. The results show that the error may grow with the increase of the dilation parameter. However, results from all cases reveal that we may obtain accurate results with the use of an arbitrary value in the range of 1.5 to 5 for the dilation parameter. Nevertheless, the error may increase while the larger dilation parameter is adopted. However, it reveals that we could still obtain better accuracy than original PS even in the worst case. Furthermore, it appears that the accuracy may be less sensitive to the order k from the numerical results. It is found that accurate results with order k in the range of 3 to 10 were acquired. (3) The seismic wave propagation problem, groundwater flow problem, unsaturated flow problem, and groundwater contamination problem were investigated to show the robustness and accuracy of the proposed approach. From the results, the proposed infinitely smooth PRBF could provide much more accurate solutions than the conventional PRBF even with the optimal order.