A well-conditioned and e ﬃ cient implementation of dual reciprocity method for Poisson equation

: One of the attractive and practical techniques to transform the domain integrals to equivalent boundary integrals is the dual reciprocity method (DRM). The success of DRM relies on the proper treatment of the non-homogeneous term in the governing di ﬀ erential equation. For this purpose, radial basis functions (RBFs) interpolations are performed to approximate the non-homogeneous term accurately


Introduction
The boundary element method (BEM) is a well-known, attractive and easy tool to deals with partial differential equations (PDEs). The idea is to convert the PDEs to equivalent boundary integral equations (BIEs). The BEM is attractive to its competitor's methods like the finite volume method (FVM) and the finite element method (FEM). Firstly, the reduction of dimensions, secondly, to deal the problems with an infinite domain other than introducing an extra condition. The BEM faces difficulties when there lies a non-homogeneous term in the PDEs.
There exists a dual reciprocity method (DRM), which handles the PDEs with non-homogeneous terms. The idea is to transform the domain integrals into equivalent BIEs while using radial basis function (RBFs), the well-known interpolation technique. In the literature, efforts are carried out to better interpolate the non-homogeneous terms while using different RBFs, see, Agnantiaris et al., [1] for thin-plate spline (TPS), Karur and Ramachandran [16] for augmented thin-plate spline (ATPS), Golberg et al., [10] for multiquadrics (MQs), Yadama et al., [30] for Gaussians (GA) and, Jumarhon et al., [14] for higher order splines.
In the work of Wendland [27] and Schaback [24], compactly supported positive definite radial basis functions (CSRBFs) have been developed and used for multivariate surface reconstruction. It was Chen et al., [4,6], who applied CSRBFs to DRM for solving the Poisson equation in 2D and 3D while using the method of fundamental solution (MFS). A detailed survey on the CSRBFs for the Laplace operator and Helmholtz operator and particular solutions of CSRBFs can be found in the work of Chen et al., [5]. Further, Fasshauer [9] implemented CSRBFs along with the Hermite collocation method for solving different types of PDEs.
In the literature, meshless methods have got much attention for the evaluation of BIEs. The method of Zhang et al., known as the dual interpolation boundary face method (DiBFM) [31,32,34,35], is one of the attractive method in terms of accuracy and efficiency, in contrast to the BEM. The DiBFM has been successfully implemented by Zhang et al., [32] to potential problems, Zhang et al., [33] to elasticity problems and to Poisson equation by Khan et al., [17] in connection to the DRM using ATPS. For complex geometries with sharp edges, the DiBFM with Hermite-type moving least squares (HMLS) approximation is adopted and applied to Potential problems by Zhang et al., [35], to elasticity problems by Zhang et al., [36], and to Poisson equation by He et al., [13] with DRM using ATPS. It was Zhou et al., [37], who implemented the variable-shaped RBFs with connection to the dual reciprocity boundary faced method to solve boundary value problems governed by the Poisson equations, and have got stable and accurate algorithms. Along with, Cheng et al., [7], have developed an iterative DRM based on CSRBFs for the solution of the Poisson equation without the need of assembling a matrix. Additionally, Loeffle et al., [18,19], have applied a new BEM technique, the direct interpolation boundary element method (DIBEM), the method is simpler and closely resembles to an interpolation procedure despite the use of primitive radial functions in the transformation of the domain integrals into boundary integrals. The DIBEM has been successfully applied with both the global RBFs and CSRBFs and has attained better accuracy and efficiency.
The RBFs have many exciting and attractive features in general, but most of them are globally defined functions, which means that the interpolation matrix can be dense and ill-conditioned. These global RBFs behave severe stability issues and computational cost, especially for large-scale data interpolation.
Keeping in view the above limitations of the global RBFs, we were looking for basis functions with local support. Fortunately, Wendland's [27] CSRBFs exist, which have the local support and can produce a sparse interpolation matrix, especially for large-scale problems. The CSRBFs is a powerful interpolation technique and having a solid mathematical basis for treating the large-scale data interpolation, see [24,27]. The implementation of CSRBFs for large-scale problems leads to a sparse matrix due to the compact support.
Therefore, the interest of this paper is on the implementation of DRM with CSRBFs for large-scale interpolation points. Moreover, the method is also compared to the DRM with the well-known MQ RBFs. Also, the convergence analysis of DRM-MQ and DRM-CSRBF is obtained. Furthermore, the CSRBFs produce a sparse interpolation matrix compared to MQ RBFs for the same interpolation points.
The remaining paper is outlined as follows. In Section 2 a straightforward implementation of DRM for the Poisson equation is introduced, and different types of interpolation functions are described. Furthermore, different types of global RBFs and CSRBFs are introduced, the error estimates and stability analysis are presented. In Section 3, the detailed convergence analysis of DRM-MQ and DRM-CSRBF are drawn. In Section 4, few numerical experiments are conducted to compare MQ and CSRBFs results for small and large interpolation points. Finally, Section 5 deals with the conclusions of the paper.

Basic formulation of DRM
Consider the PDEs: with a source function f , and domain Ξ is enclosed by Υ = Υ u ∪ Υ q , along with boundary conditions defined by whereū andq are the given values of potential and normal flux along Υ u and Υ q , respectively, and n represent normal to the boundary Υ.
In integral form the boundary value problems govern by the equations (2.1)-(2.3) can be written as follows: [2,22] cu where u * represents the fundamental solution of Laplace equation for 2D problems and can be written as: Here r is the distance among the collocation point and field point. The constant c has value from 0 to 1, it's value is 1 2 if the boundaries are smooth and 1 when the source point lies in the domain [2]. Now the target is to transform the domain integral in Eq (2.4) into an equivalent boundary integrals using the well known procedure DRM [22]. The idea is to expand f (x, y) by the following approximate functions: where ψ(r) represents the radial basis function, r i = |x − ξ i |, β i are the unknown coefficients and N + L are the collocation points (N boundary and L interior). While using Green's reciprocity identity [22] domain integral can be written as a series of surface integrals, which in boundary integral form can be represented as: (2.7) The discretized form of Eq (2.7) is The values in the matrices H jk and G jk can be get from the integration of q * and u * respectively, at each boundary element. In matrix form Eq (2.8) can be written as The coefficients β i can be determined from equation (2.5) by the collocation where B is a matrix with coefficients B jk = ψ jk and f = f i . Using β from Eq (2.5) into Eq (2.9) yields By utilizing the boundary conditions, the above equation can be reduced to where x contains N unknowns values of u or q.

Interpolation functions
Let the value of function g(x) be given, on nodes we can say that g(x) interpolates the given data {(x j , g(x j ))} N j=1 . When the interpolation conditions are imposed on g(x), we get (2.11) This produce a system of N linear equations in N unknowns Aβ = b, i.e.
The matrix A of order N × N is called the interpolation matrix. To ensure the solvability of the system for g(x i ), it is necessary and sufficient that the interpolation matrix be nonsingular. We will need the following definitions.
where r = ||x|| and ||.|| is some norm on R d , usually the Euclidean norm.
Definition 2. A continuous function ψ : R d → R is said to be positive definite (PD) iff for all N ∈ ℵ, all sets of pairwise distinct centers X = {x 1 , ..., is a real valued strictly positive definite function. Then the real native Hilbert space of ψ on R d is introduced as Definition 6. The condition number K 2 (A) of the interpolation matrix A is given by with s > d/2 and two positive constants c 1 ≤ c 2 . Then the native space ℵ ψ (R d ) corresponding to ψ coincides with the Sobolev space H s (R d ), and the native space norm and the Sobolev norm are equivalent.

Theorem 1. [29]
Suppose ψ is a symmetric positive definite kernel on a compact set Ξ ⊆ R d . Then its native space is given by , and the inner product has the representation Theorem 2. [29] Let ψ d,k = ψ d,k (||.|| 2 ) denote the compactly supported radial basis function of minimal degree that is positive definite on R d and in C 2k . Let d ≥ 3 if k = 0. Then there exist constants c 1 , c 2 > 0 depending only on d and k such that for all ξ ∈ R d . This means in particular that i.e. the native space for these basis functions is a classical Sobolev space.
The proper choice of RBFs in DRM is much concerned because the accuracy of DRM depends on the best approximation of a particular solution. For more details of the choice of RBFs, one may refer to the review articles [11,12,23]. We have adopted DRM-MQ with multiquadric RBFs (ψ(r) = √ r 2 + 2 , where is a positive parameter) in our computation. Due to the dense interpolation matrix of MQ, we were looking for basis functions that can result in a sparse matrix. Fortunately, Wendland's CSRBFs defined in Table 1 which produces a sparse matrix, even for large-scale data. We adopted Wendland's CSRBFs defined in Table 1 as a basis function in DRM-CSRBF for evaluation of the Poisson equation, especially for large-scale problems. Table 1. Wendland's CS-RBFs [27].
After the invention of the CSRBFs, some researchers hoped that using these basis functions was the solution to the trade-off problem. For the practical usefulness of CSRBFs, it is observed that for fixed support radius CSRBFs is not recommended. To have a sparse system matrix, the support needs to be kept so small that the resulting fit is unacceptable (i.e., the errors are large), or the support is chosen to be large enough to make the error acceptable, but which produce a dense matrix, and so computationally expensive. As the trade-off principle [26], a large bandwidth of A results in small errors, but computational inefficiency, whereas a small bandwidth means O(N) complexity, but large errors.
In numerical experiments we have adopted the following Wendland's CSRBFs: and α is the scaling parameter of CSRBFs. As accuracy of MQ-RBFs rely on the optimal value of shape parameter, just like that the accuracy and efficiency of CSRBFs are affected with the scaling parameter.
For more details about CSRBFs and scaling parameter the readers are referred to the articles [8,24,27].
We have adopted specific values of shape parameter and scaling parameter according to the problem and have shown its effect on accuracy and efficiency.

Error estimates and stability analysis of the interpolation functions
Here, we have described the most famous error estimates and stability analysis of the interpolation functions for global RBFs and CSRBFs. In his work, Schaback [25] pointed out the error estimate in the native space, represented by: where P(x) is the power function, which is just norm of the error functional on ℵ computed at x, where The generalized Fourier transform of ψ is denoted withψ. The power function P(x) is bounded above by F(h(x)), which has different values for different RBFs. In Schaback work [25], the stability analysis was given by a lower bound of the eigenvalues λ of the interpolation matrix. The lower bound were represented by G(h(x)), which are evaluated for different global RBFs. The exact value of the upper bound F(h(x)) and lower bound G(h(x)) for global RBFs can be found in [25]. Here, we will mention the lower bound of MQ, Gaussian (GA) and inverse multiquadrics (IMQ) along with shape parameter .
Here, h = sup y∈Ω min x∈X ||x − y||, d is the dimension, is the shape parameter and X is the interpolation set. Moreover, the error estimation of CSRBFs with minimal degree was evaluated by Wendland [28], which is given as: where g is a function from H s (R d ) and h defined as above and C 1 is a constant. H s is the usual Sobolev space of functions with s derivatives bounded in L 2 . The regularity of the data is given by s = d 2 + k + 1 2 (k ≥ 1 for d = 0, 1), and k defines the smoothness of the functions ψ, i.e., ψ ∈ C 2k . Furthermore, the upper bound F(h(x)) were computed by Wendland and are given in his book [29] for the CSRBFs. These bounds for CSRBFs are: In addition, the upper bounds on the condition number of the interpolation matrix are given: where C 2 is a constant independent of q x = 1 2 min i j ||x i − x j ||.

Convergence of DRM
This section describes the convergence analysis of DRM in the context of MQ-RBFs and CSRBFs.
[28] suppose ψ l (r) = (1 − r) l + and let ψ d,k (r) = I k ψ l (r) ∈PD∩C 2k be the function of minimal degree with l = d 2 + k + 1 > 1. Let X = {x 1 , ..., x N } ⊆ R d be a set of centers with separation distance 2q x = min i j ||x i − x j || > 0 and denote the interpolation matrix with A X,ψ d, The norm of the inverse of the interpolation matrix then satisfies: Here, I is an operator defined by: Now for the convergence of DRM consider: Let v b represents the boundary solution of v. The DRM solution of boundary value problems in Eqs (2.1) and (2.2) is u h , which can be written as The error e can be represented as: It follows from the results of Miranda [20], that Poisson equation satisfies: From the BEM theory (see Ref. [3]), we have: It remains to bound ||û||, it can be followed from [11] that: here N represents the number of interpolation points and A is the interpolation matrix: The bounds on ||A −1 || are available in [25] for different RBFs, for MQ it can be bounded as ||A −1 || ≤ C 8 e σ δ , where δ = sup y∈Ω min x∈X ||x − y||. Asf is MQ interpolant, it follows from the results of Madych [21] that || f −f || ≤ C 9 δ n , n ≥ 0. Following the above results, equation (3.4) becomes Moreover, the error estimation of CSRBFs with minimal degree was evaluated by Wendland [28], which is given by: where s = d 2 + k + 1 2 and g is a function from H s (R d ) and δ defined as in Eq (3.6) and C 11 is a constant. Now the convergence of DRM-CSRBF can be drawn while substituting Eqs (3.5) and (3.7) in Eq (3.4) and using Theorem 3, we obtain the following: (3.8)

Numerical examples
This section describes the implementation of the DRM for solving some benchmark problems. Three problems with a known exact solution are evaluated, and numerical results of DRM-MQ and DRM-CSRBF are compared in each problem. To judge the error and convergence analysis of the method, a Root Mean Square Error (RMSE) is adopted to measure the accuracy, which can be defined as Here For the numerical evaluation we have adopted constant element in the DRM. The spatial distribution of the non-homogeneous term in the right hand side of Eq (4.1) is presented in Figure 1 (Left), while the exact solution of Eq (4.1) is shown in Figure 1 (Right). For varying boundary nodes and constant number of inner nodes the convergence analysis of the normal flux q and potential u I are presented in Figure 2 (Left) and Figure 2 (Right), respectively. For the same normal flux q and potential u I condition number and computational time is presented in Figure 3      with the domain Ξ (0 ≤ x ≤ 1 and 0 ≤ y ≤ 1) and BC, u = 0 on δΞ. The analytical solution of the problem is: u(x, y) = sin(2πx) cos(2πy).
The spatial distribution of the non-homogeneous term in the Eq (4.2) is presented in Figure 4 (Left), while the exact solution of Eq (4.2) is shown in Figure 4 (Right). For varying boundary nodes and fixed inner nodes the convergence analysis of the normal flux q and potential u I are presented in Figure 5 (Left) and Figure 5 (Right), respectively. For the same normal flux q and potential u I condition number and computational time is presented in Figure 6        Here, we adopted l = 1. The analytical solution of u is described in [18] u = − l π 2 sinh(π) sinh πx l − x π 2 cos πy l .
The solution of normal flux q is defined with where n x and n y are components of the normal on the boundary. The spatial distribution of the nonhomogeneous term in the Eq (4.3) is presented in Figure 9 (Left), while the exact solution of Eq (4.3) is shown in Figure 9 (Right). For varying boundary nodes and fixed inner nodes the convergence analysis of the potential on the boundary u b and normal flux q are presented in Figure 10 (Left) and Figure 10 (Right), respectively, while the convergence of potential u I is given in Figure 11 (Left). For the potential u b , normal flux q, and potential u I , the condition number and computational time is presented in Figure 11 (Right) and Figure 12 (Left), respectively. For fixed boundary nodes and varying inner nodes the convergence analysis of potential on the boundary u b , normal flux q, and potential u I , are shown in Figure 12 (Right), Figure 13 (Left) and Figure 13 (Right), respectively. While the condition number and computational behavior is presented in Figure 14 (Left) and Figure 14 (Right), respectively. For the numerical evaluation of example 3, we have used ψ 2 r α = 1 − r α 4 + 4 r α + 1 with α = 0.8 in DRM-CSRBF and = 0.8 in DRM-MQ. The effects of the shape parameter in case of MQ-RBFs by the potential on the boundary u b , normal flux q, and inner potential u I , is given in Figure 15 (Left), while the effects of the scaling parameter in case of CSRBFs by the potential on the boundary u b , normal flux q, and inner potential u I , is presented in Figure 15 (Right). From the above figures, we can claim that the DRM-CSRBF is an accurate, efficient and well-conditioned method for evaluation of the Poisson equation with large scale data.

Conclusions
This work demonstrated the implementation of DRM with MQ RBFs and CSRBFs to evaluate the Poisson equation, especially for large-scale data, while using constant boundary elements in the BEM. As stated in the introduction, that the interpolation matrix of the MQ RBFs is dense and has a large condition number, especially for large-scale problems, so its computation is expensive and poses serious stability issues. For this purpose, we have applied CSRBFs, which produce a sparse matrix, even for large-scale problems, that can reduce the computational cost and risk of stability. Each example is evaluated with DRM-MQ and DRM-CSRBF, and the comparison of accuracy, efficiency, and condition numbers are presented. The numerical experiments have demonstrated that the DRM-CSRBF is an accurate, efficient, and well-conditioned method for evaluating the Poisson equation, especially for large-scale problems.