**1. Introduction**

Recently, the meshless approach has raised extensive attention due to its computational efficiency as well as simple collocation scheme. Many varieties of the radial basis functions (RBFs) have been developed for dealing with partial differential equations (PDEs) [1–3]. Most popular RBFs, such as the Gaussian [4–6], multiquadric (MQ) [7,8], and inverse multiquadric (IMQ) [9–11], require the shape parameter. Among them, the Kansa method [12] is recognized as one of the most popular domain-type meshfree approaches for solving PDEs. The MQ RBF adopted by the Kansa method becomes the well-known RBF, which has been successfully adopted for solving numerous engineering problems [13,14]. Despite the success of the Kansa method, limitations regarding to the accuracy affecting by the shape parameter still remain. The MQ RBF depends on the shape parameter that plays an important role for remaining the RBF as a smooth and non-singular function for solving PDEs. Attempts regarding for identifying proper value for the shape parameter of the MQ RBF have been widely studied, such as the LOOCV optimization technique [15–17]. The question of finding the optimal shape parameter in the MQ RBF, however, is still very challenging.

In this study, we propose radial polynomials (RPs) rooted in the MQ RBF for solving PDEs. Formulated from the binomial series using the Taylor series expansion of the MQ RBF, the new global RPs include only even order radial terms. The proposed RPs may be regarded as the equivalent expression of the MQ RBF in series form. Not only are the RPs infinitely smooth and differentiable in nature, but the proposed RPs do not require any extra shape parameters. Therefore, the challenging task for finding the optimal shape parameter in the Kansa method is avoided. Several numerical implementations, including two- and three-dimensional problems, are conducted to verify the accuracy and robustness of the proposed RPs. The structure of this article is organized as follows: In Section 2, formulation of the radial polynomial basis function is presented. To verify the proposed RPs, we conduct a convergence analysis in Section 3. Section 4 is devoted to present several numerical examples in two and three dimensions. The discussion of this paper is addressed in Section 5. Conclusions are given finally.

### **2. Formulation of the Radial Polynomials**

Considering a region, Ω, with the boundary, ∂Ω, the governing equation for the three-dimensional PDE can be expressed as follows.

$$
\Delta u(\mathbf{x}) + D \frac{\partial u(\mathbf{x})}{\partial \mathbf{x}} + E \frac{\partial u(\mathbf{x})}{\partial y} + F \frac{\partial u(\mathbf{x})}{\partial z} + Gu(\mathbf{x}) = H \text{ in } \Omega,\tag{1}
$$

$$u(\mathbf{x}) = \mathcal{g}(\mathbf{x}) \text{ on } \partial \Omega^D,\tag{2}$$

$$\frac{\partial u(\mathbf{x})}{\partial n} = f(\mathbf{x}) \text{ on } \partial \Omega^N,\tag{3}$$

in which Δ represents Laplace operator, **x** = (*x*, *y*, *<sup>z</sup>*), *u*(**x**) is the unknown, *D*, *E*, *F*, *G* and *H* are given functions. Ω is a bounded domain with boundary ∂Ω*<sup>D</sup>* and ∂Ω*N*. ∂Ω*<sup>D</sup>* denotes boundary subjected to Dirichlet data, ∂Ω*<sup>N</sup>* denotes boundary subjected to Neumann data, *g*(**x**) and *f*(**x**) represent given boundary data. The meshless method using the MQ RBF is often named the Kansa method, where the RBFs are directly implemented for the approximation of the solution of partial differential equations. We may express the unknown by the RBF as follows.

$$\mu(\mathbf{x}) = \sum\_{j=1}^{M\_{\xi}} \lambda\_{j} \wp(r\_{j}),\tag{4}$$

where *rj* is the radial distance, *rj* = **x** − **s***j*, ϕ(*rj*) represents the RBF which is the distance of **x** and **<sup>s</sup>***j*, **s***j* is the center, **x** denotes an arbitrary point inside the domain, λ*j* is the coefficient to be solved and *Mc* is the number of the center points. The MQ RBF may be expressed as follows.

$$
\varphi(r\_j) = \sqrt{r\_j^2 + c^2}.\tag{5}
$$

With the introduction of the shape parameter, the MQ RBF becomes a smooth and non-singular function. Because the Kansa method is a domain-type method, it has to discretize the governing equation inside the domain using the MQ RBF. We may insert the above equation into Equation (1). After obtaining the MQ RBF derivatives, we may obtain the following equation in two-dimensions.

$$\sum\_{j=1}^{M\_c} \lambda\_j \frac{r\_j^2 + 2c^2}{\left(r\_j^2 + c^2\right)^{1.5}} + \sum\_{j=1}^{M\_c} \lambda\_j \frac{D\left(\mathbf{x} - \mathbf{x}\_j\right) + E\left(\mathbf{y} - \mathbf{y}\_j\right)}{\left(r\_j^2 + c^2\right)^{0.5}} + G \sum\_{j=1}^{M\_c} \lambda\_j \left(r\_j^2 + c^2\right)^{0.5} = H \text{ in } \Omega. \tag{6}$$

The above equation demonstrates that the derivatives of the MQ basis function may become singular at the center point (*rj* = 0) if the shape parameter is zero. It is obvious that the MQ RBF is infinitely differentiable depending on the shape parameter. To avoid the singularity, the shape parameter must not be equal to zero. In this study, we propose RPs based on the MQ RBF without the shape parameter. For the mathematical formulation of the RPs, we may start from the MQ RBF. Equation (5) can be rewritten as follows.

$$
\varphi(r\_{\dot{\jmath}}) = c \sqrt{\left(r\_{\dot{\jmath}}/c\right)^2 + 1}.\tag{7}
$$

Using the binomial series from the Taylor series of Equation (7), we have

$$c\sqrt{(r\_{\dot{\jmath}}/c)^2 + 1} = c\sum\_{k=0}^{\infty} \binom{\alpha}{k} \left( (r\_{\dot{\jmath}}/c)^2 \right)^k,\tag{8}$$

where α*k* := <sup>α</sup>(<sup>α</sup>−<sup>1</sup>)(<sup>α</sup>−<sup>2</sup>)···(<sup>α</sup>−*k*+<sup>1</sup>) *k*! and α = 1/2.

Using the finite terms, *Mn*, to approximate the solution, we may express the MQ RBF in series form as follows.

$$\varphi(r\_j) = c \sum\_{k=0}^{M\_\pi} \binom{0.5}{k} \left(\frac{1}{c}\right)^{2k} r\_j^{2k} \,. \tag{9}$$

where *Mn* is the order of the radial polynomials. In this study, we propose a novel meshless method to approximate the solution in terms of the RPs as follows.

$$u(\mathbf{x}) = \sum\_{j=1}^{M\_c} a\_j \varphi(r\_j),\tag{10}$$

where *Mc* represents the center point number. The above equation proves that the MQ RBF can be expressed as a radial polynomial with only even order terms. Equation (8) can be regarded as the equivalent series form of the MQ RBF. Inserting Equation (8) into Equation (10), we have

$$u(\mathbf{x}) \approx \sum\_{j=1}^{M\_{\mathcal{L}}} a\_j c \sum\_{k=0}^{M\_{\mathcal{R}}} \binom{0.5}{k} \binom{1}{c}^2 r\_j^{2k}.\tag{11}$$

Combining the constants in the above equation, we obtain

$$u(\mathbf{x}) \approx \sum\_{j=1}^{M\_c} \sum\_{k=0}^{M\_n} b\_{j,k} r\_j^{2k} \,. \tag{12}$$

in which *bj*,*<sup>k</sup>* are the coefficients to be solved. Using Equation (12) for the discretization of Equation (1), we may obtain the following equation:

$$\sum\_{j=1}^{M\_{\mathcal{L}}} \sum\_{k=0}^{M\_{\mathcal{R}}} b\_{j,k} L\_{1c} r\_j^{2k-2} + \sum\_{j=1}^{M\_{\mathcal{L}}} \sum\_{k=0}^{M\_{\mathcal{R}}} b\_{j,k} 2k L\_{2c} r\_j^{2k-2} + G \sum\_{j=1}^{M\_{\mathcal{L}}} \sum\_{k=0}^{M\_{\mathcal{R}}} b\_{j,k} r\_j^{2k} = H \text{ in } \Omega,\tag{13}$$

where *L*1*c* = 4*k*2, *L*2*c* = (*D*(*x* − *xj*) + *<sup>E</sup>*(*y* − *yj*)) and *L*1*c* = 4*k*<sup>2</sup> + 2*k*, *L*2*c* = (*D*(*x* − *xj*) + *<sup>E</sup>*(*y* − *yj*) + *<sup>F</sup>*(*z* − *zj*)) are in two and three dimensions, respectively. To determine the unknown coefficients, we apply the approximate solution with the boundary data at collocation points to satisfy the governing equation. We may ge<sup>t</sup> the system of simultaneous equations.

$$\mathbf{A}\mathbf{b} = \mathbf{R},\tag{14}$$

where **b** is the unknown coefficient with the size of *N* × 1 to be evaluated, **R** is the known function with the size of *M* × 1, **A** is an *M* × *N* matrix where *M* = *Mi* + *Mb* and *N* = *Mc* × *Mn*. The above equation can be written as follows:

$$
\begin{bmatrix}
\mathbf{A}\_{\mathrm{I}} \\
\mathbf{A}\_{\mathrm{B}}
\end{bmatrix}
\begin{bmatrix}
\mathbf{b}
\end{bmatrix} = 
\begin{bmatrix}
\mathbf{R}\_{\mathrm{I}} \\
\mathbf{R}\_{\mathrm{B}}
\end{bmatrix}.
\tag{15}
$$

In the preceding equations, **A**I represents the *Mi* × *N* submatrix from the inner collocation points, **A**B represents the *Mb* × *N* submatrix from the boundary collocation points, **R**I is the vector of function values at the inner points which is a *Mi* × 1 vector, **R**B is the data at the boundary points which is an *Mb* × 1 vector, *Mb* is the boundary point number, *Mi* is the inner point number. The root mean square error (RMSE) is adopted to evaluate the accuracy which is defined by

$$\text{Root mean square error} = \sqrt{\frac{1}{M\_m} \sum\_{i=1}^{M\_m} \left( \vec{u}(\mathbf{x}\_i) - \boldsymbol{u}(\mathbf{x}\_i) \right)^2} \tag{16}$$

in which *Mm* represents the number of the measuring points with uniform distribution; *u*(*xi*) and *u*<sup>ˆ</sup>(*xi*) are the exact and approximate solutions at the *ith* collocation point, respectively.

### **3. Accuracy and Convergence Analysis**

We first investigate a Laplacian problem in two dimensions enclosed by an irregular domain. The governing equation is

$$
\Delta u(\mathbf{x}) = 0, (\mathbf{x}) \in \Omega. \tag{17}
$$

The star–like object boundary in two dimensions can be expressed in the following form:

$$\partial \Omega = \left\{ (\mathbf{x}, y) \Big| \mathbf{x} = \rho(\theta) \cos \theta, y = \rho(\theta) \sin \theta, \rho(\theta) = \sec \left( 3\theta \right)^{\sin(6\theta)}, 0 \le \theta \le 2\pi \right\}.\tag{18}$$

The exact solution of Equation (17) is designated as

$$u(\mathbf{x}) = \varepsilon^x \cos(y) + \varepsilon^y \sin(\mathbf{x}).\tag{19}$$

To verify the accuracy and convergence, we conduct a series of testing cases for the radial polynomial terms in which all cases adopt the same configurations of the boundary, center and inner points as shown in Figure 1. In the analysis, *Mb*, *Mi* and *Mc* are 1208, 151, and 151, respectively. The number of the RPs terms, *Mn*, needs to be given for the proposed method. As shown in Figure 2, for the RPs, it is found that the RMSE decreases with the increase in the number of RPs terms in which solutions with high accuracy may be found with the radial polynomial terms from 6 to 12.

**Figure 1.** Layout of the collocation points.

**Figure 2.** The convergence analysis for the Kansa method and the RPs.

On the other hand, other testing cases using the Kansa method for considering different shape parameters are conducted. Figure 2 shows different shape parameters versus the RMSE. The optimal shape parameter is found within a narrow range of 0.5 to 1. We may also observe that the shape parameter is very sensitive to the accuracy in the Kansa method, such that attempts regarding for identifying proper values for the shape parameter of the Kansa method may be important. It is apparent that the minimums of the RMSE for the Kansa method and the RPs are in the order of 10−<sup>9</sup> and <sup>10</sup>−12, respectively.

In addition, to investigate the accuracy, another convergence analysis for investigating the boundary and inner point number is carried out. Table 1 shows the comparison of this study with the Kansa method. We find that very high accurate results may be obtained using the proposed RPs.


**Table 1.** Results comparison between this study and the Kansa method with the optimal shape parameter.

### **4. Numerical Examples**

To investigate the applicability of the proposed RPs, four numerical examples are conducted, in which Sections 4.1 and 4.2 are steady-state linear two-dimensional PDEs, Section 4.3 is the three-dimensional modified Helmholtz equation, and Section 4.4 is the three-dimensional Poisson equation.

### *4.1. A Two-Dimensional Ameoba-Shaped Problem*

We first consider the following two-dimensional PDEs.

$$
\Delta u(\mathbf{x}) + D \frac{\partial u(\mathbf{x})}{\partial \mathbf{x}} + E \frac{\partial u(\mathbf{x})}{\partial y} + F \frac{\partial u(\mathbf{x})}{\partial z} + Gu(\mathbf{x}) = H, \ \mathbf{x} \in \Omega,\tag{20}
$$

where *D* = *y* cos(*x*), *E* = −sin *h*(*x*), *F* = 0, *G* = *x*2 + *y*2. The function, *H*, can be directly derived from the exact solution as follows:

$$\begin{array}{l} H = \left(-\pi^2 + 1\right) \left[\sin(\pi x)\cosh(y)\right] + \left(\pi^2 - 1\right) \left[\cos(\pi x)\sinh(y)\right] \\ + \left(\pi y\cos(x) - \sinh(x)\right) \left[\sin(\pi x)\sinh(y)\right] + \left(\pi y\cos(x) + \sinh(x)\right) \left[\cos(\pi x)\cosh(y)\right] \\ + \left(\pi^2 + y^2\right) \left[\sin(\pi x)\cosh(y) - \cos(\pi x)\sinh(y)\right] \end{array} \tag{21}$$

The amoeba-like object boundary in two dimensions is defined as

$$\partial \Omega = \left\langle (\mathbf{x}, y) \middle| \mathbf{x} = \rho(\theta) \cos \theta, y = \rho(\theta) \sin \theta, \rho(\theta) = e^{\left(\sin \theta \sin \theta \mathbf{s}\right)^2} + e^{\left(\cos \theta \cos \theta \mathbf{c}\right)^2}, 0 \le \theta \le 2\pi \right\}. \tag{22}$$

Both Dirichlet and Neumann boundary conditions are considered as follows:

$$\ln(\mathbf{x}) = \sin(\pi \mathbf{x})\cosh(y) - \cos(\pi \mathbf{x})\sinh(y), (\mathbf{x}, y) \in \partial \Omega^{D},\tag{23}$$

$$\frac{\partial \mu(\mathbf{x})}{\partial \mathbf{n}} = \left[\nabla(\sin(\pi \mathbf{x})\cosh(y) - \cos(\pi \mathbf{x})\sinh(y))\right] \cdot \overleftarrow{n}\_{\prime}(\mathbf{x}, y) \in \partial \Omega^N. \tag{24}$$

In this example, the Kansa method and the proposed RPs are examined. Figure 3 depicts the configuration of the collocation points. The over-specified Dirichlet as well as Neumann boundary data are imposed on the whole boundary. In the analysis, *Mb*, *Mi* and *Mc* are 1750, 500 and 500, respectively. The analysis of convergence for the RPs terms is conducted, as shown in Figure 4. According to Figure 4, it is found that highly accurate solutions may be solved with the radial polynomial terms from 7 to 12. Consequently, the terms of the RPs are set to 9. The result comparison for the Kansa method and the proposed RPs is shown in Table 2. Table 2 demonstrates that highly accurate results are obtained in which the RMSE of the proposed method is within the order of 10−8. On the other hand, the minimum RMSE for the Kansa method with the optimal shape parameter can only reach to the order of 10−4. Figure 4 demonstrates results of the convergence analysis in which it is found that solutions with high accuracy may be obtained with the radial polynomial terms from 6 to 11. Moreover, it is clear that the number of terms is not very sensitive to the result. Figure 5 depicts the numerical solution is identical to the exact solution.

**Figure 3.** Layout of the collocation points.

**Figure 4.** The terms of the RPs versus the RMSE.



**Figure 5.** Result comparison between numerical and the analytical solutions.

### *4.2. A Two-Dimensional Star-Shaped Problem*

The second example is a problem enclosed by a star-shaped boundary in two dimensions.

$$
\Delta u(\mathbf{x}) + D \frac{\partial u(\mathbf{x})}{\partial \mathbf{x}} + E \frac{\partial u(\mathbf{x})}{\partial y} + F \frac{\partial u(\mathbf{x})}{\partial z} + Gu(\mathbf{x}) = H, \ \mathbf{x} \in \Omega,\tag{25}
$$

where *D* = *y*2 sin(*x*), *E* = *xey*, *F* = 0, *G* = sin(*x*) + cos(*y*). The function, *H*, can be directly derived from the exact solution as follows:

$$\begin{array}{l} H = \left( -\pi^2 y \sin(\pi x) - \pi^2 x \cos(\pi y) \right) + \left( y^2 \sin(\pi) \right) [\pi y \cos(\pi x) + \cos(\pi y)]\\ + (\pi e^y) [\sin(\pi x) - \pi x \sin(\pi y)] + (\sin(x) + \cos(y)) [y \sin(\pi x) + x \cos(\pi y)] \end{array} \tag{26}$$

The star-like object boundary in two dimensions is defined as

$$\partial \Omega = \left| (\mathbf{x}, y) \right| \mathbf{x} = \rho(\theta) \cos \theta,\\ y = \rho(\theta) \sin \theta,\\ \rho(\theta) = 1 + (\cos 4\theta)^2,\\ 0 \le \theta \le 2\pi \text{ } \right|. \tag{27}$$

The Dirichlet boundary conditions are considered as follows:

$$u(\mathbf{x}) = y \sin(\pi \mathbf{x}) + \mathbf{x} \cos(\pi y), (\mathbf{x}, y) \in \partial \Omega^{\mathcal{D}}.\tag{28}$$

In this example, the Dirichlet data are applied on the whole boundary using Equation (28). In the analysis, *Mb*, *Mi* and *Mc* are 1800, 200 and 200, respectively. We conduct the convergence analysis for the RPs terms. Figure 6 displays the configuration of the boundary, inner and center collocation points. Figure 7 displays the terms of the RPs versus the RMSE in which we may find that solutions with high accuracy in the order of 10−<sup>8</sup> may be found with the radial polynomial terms from 7 to 12. Consequently, the terms of the RPs are set to 9.

**Figure 6.** Layout of the collocation points.

**Figure 7.** The terms of the RPs versus the RMSE.

Figures 8 and 9 demonstrate the RMSE versus the boundary and inner point numbers, respectively. We may find that promising solutions may be found while the boundary and inner point numbers are greater than 500 and 100, respectively. Figure 10 demonstrates the comparison of the analytical and the numerical solutions. It can be found that the results agree with the analytical solutions.

**Figure 8.** The convergence analysis for the boundary point number.

**Figure 9.** The convergence analysis for the inner point number.

**Figure 10.** Result comparison between numerical and the analytical solutions.

### *4.3. A Three-Dimensional Modified Helmholtz Problem*

Consider a three-dimensional modified Helmholtz equation. The equation is written as follows.

$$
\Delta u(\mathbf{x}) + Gu(\mathbf{x}) = H, \; \mathbf{x} \in \Omega,\tag{29}
$$

where *D* = *E* = *F* = 0, *G* = −<sup>λ</sup>2, *H* = (1 − <sup>λ</sup><sup>2</sup>)(*e<sup>x</sup>* + *ey* + *e<sup>z</sup>*) + *xyz*, λ represents wave number and λ = 100. The domain in three dimensions can be expressed in the following form.

$$\delta\Omega = \left| (\mathbf{x}, y, z) \right| \mathbf{x} = \rho(\theta) \sin \theta \cos \varphi,\\ \mathbf{y} = \rho(\theta) \sin \theta \sin \varphi,\\ z = \rho(\theta) \cos \varphi \Big|\_{\prime} \tag{30}$$

where ρ(<sup>θ</sup>,ϕ) = 1 + 1/8 sin(10 θ) sin(9 ϕ), 0 ≤ θ ≤ 2, 0 ≤ ϕ ≤ π. The Dirichlet boundary data are applied on ∂Ω using the following exact solution.

$$
\mu(\mathbf{x}) = \mathbf{e}^{\mathbf{x}} + \mathbf{e}^{y} + \mathbf{e}^{z} - \mathbf{x}yz/\lambda^{2}, (\mathbf{x}, y) \in \partial \Omega^{D}.\tag{31}
$$

In this example, the layout of the domain is depicted in Figure 11. The Dirichlet data are applied on the whole boundary using Equation (31). In the analysis, *Mb*, *Mi* and *Mc* are 7569, 800 and 800, respectively. Figure 12 demonstrates solutions with high accuracy in the order of 10−<sup>8</sup> may be found with the radial polynomial terms from 8 to 11. Consequently, the terms of the RPs are set to 9. The result comparison for the Kansa method and the proposed RPs is shown in Table 3. Table 3 demonstrates that highly accurate results are obtained in which the RMSE is within the order of 10−11. On the other hand, the best RMSE for the Kansa method with the optimal shape parameter can only reach to the order of 10−7.

**Figure 11.** Layout of the three-dimensional modified Helmholtz equation.



### *4.4. A Three-Dimensional Poisson Problem*

The last example under consideration is a three-dimensional Poisson equation enclosed by an irregular domain. The governing equation is expressed as follows.

$$
\Delta u(\mathbf{x}) = H, \; \mathbf{x} \in \Omega,\tag{32}
$$

where *D* = *E* = *F* = *G* = 0 and *H* = −sin *x* − cos *y* − sin *z*. The domain in three dimensions can be expressed in the following parametric equation.

$$\delta\Omega = \left\{ (\mathbf{x}, y, z) \middle| \mathbf{x} = \rho(\theta) \cos \theta, y = \rho(\theta) \sin \theta \sin \phi, z = \rho(\theta) \sin \theta \cos \phi \right\}, \tag{33}$$

where ρ(θ) = cos(3<sup>θ</sup>) + 8 − sin<sup>2</sup>(3θ)1/3. The Dirichlet boundary data are assigned on ∂Ω*<sup>D</sup>* using the following exact solution.

$$
\mu(\mathbf{x}) = \sin(\mathbf{x}) + \cos(y) + \sin(z), (\mathbf{x}, y) \in \partial \Omega^D. \tag{34}
$$

The layout of the domain is depicted in Figure 13. The Dirichlet boundary conditions are given on the irregular domain in three dimensions using Equation (34). In the analysis, *Mb*, *Mi* and *Mc* are 7357, 756 and 756, respectively. Figure 14 shows the RMSE versus the terms of the RPs. It is apparent that the promising numerical solution in the order of 10−<sup>8</sup> may be obtained while the *Mn* is greater than 6. Consequently, the terms of the RPs is set to 9. Additionally, several cases for evaluating the number of the collocation points to the accuracy are conducted in Table 4. According to Table 4, it depicts that the accuracy can reach up to the order of 10−10.


**Table 4.** The RMSE of the RPs.

**Figure 13.** Layout of the three-dimensional Poisson equation.

**Figure 14.** The terms of the RPs versus the RMSE.

### **5. Discussion**

This study presents a collocation method using RPs which is regarded as the equivalent expression of the MQ RBF in series form for PDEs. The conception of the new global RPs includes only even order radial terms formulated from the binomial series using the Taylor series expansion of the MQ RBF. The discussions for this study are as follows.

The MQ RBF adopted by the Kansa method becomes one of the most successful RBFs for solving numerous problems. With the introduction of the shape parameter, the MQ RBF becomes a smooth and non-singular function. Even though the MQ RBF and its derivatives are smooth and global infinitely differentiable, the discretization of the governing equation may become singular while *c* = 0 at *rj* = 0. It is obvious the shape parameter plays a role for shifting from the singularity while the center point is coincided with the inner point. However, the near singular effects still remain. This may explain that the accuracy in the Kansa method is strongly affected by the shape parameter.

To deal with the issue, we adopt the RPs as the basis function in which the proposed RPs are the equivalent expression of the MQ RBF in series form. It is advantageous that the proposed RPs and their derivatives are infinitely smooth and differentiable in nature without using the shape parameter. Because RPs are a non-singular series function, there are no near singular or singular effects at all. Accordingly, the method may obtain more accurate solutions than the Kansa method in our numerical implementations. In addition, accurate results can be directly obtained without using the tedious procedure for finding the optimal shape parameter.

Even though the shape parameter is not required in the proposed method, the radial polynomial terms have to be decided in advance. From the numerical implementations, solutions with high accuracy in the order of 10−<sup>8</sup> may be found with radial polynomial terms from 6 to 12. The radial polynomial terms are selected to be 9 in our numerical examples. It demonstrates that the radial polynomial terms are considerably less significant to the accuracy than the shape parameter. In the numerical examples, it is found that satisfactory solutions could be obtained while the terms of the RPs are within the range of 6 to 12.

### **6. Conclusions**

A mathematical formulation of the RPs from the binomial series using the Taylor series expansion of the MQ RBF is presented. We prove that the proposed RPs are an equivalent expression of the MQ RBF in series form. Highly accurate results can be directly obtained without using the tedious procedure for finding the optimal shape parameter. Additionally, numerical comparisons reveal that the presented RPs could obtain better accurate solutions than those of the MQ RBF, even with the optimal shape parameter for solving multi-dimensional PDEs.

**Author Contributions:** Writing—review and editing, J.-E.X. and C.-Y.K.; conceptualization, C.-Y.K.; methodology, C.-Y.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially supported by the Ministry of Science and Technology of R.O.C.

**Acknowledgments:** The authors would like to thank the referees for invaluable comments.

**Conflicts of Interest:** The authors declare no conflict of interest.
