An algorithm for choosing best shape parameter for numerical solution of partial differential equation via inverse multiquadric radial basis function

Radial Basis Function (RBF) is a real valued function whose value rests only on the distance from some other points called a center, so that a linear combination of radial basis functions are typically used to approximate given functions or differential equations. Radial Basis Function (RBF) approximation has the ability to give an accurate approximation for large data sites which gives smooth solution for a given number of knots points; particularly, when the RBFs are scaled to the nearly flat and the shape parameter is chosen wisely. In this research work, an algorithm for solving partial differential equations is written and implemented on some selected problems, inverse multiquadric (IMQ) function was considered among other RBFs. Preference is given to the choice of shape parameter, which need to be wisely chosen. The strategy is written as an algorithm to perform number of interpolation experiments by changing the interval of the shape parameters and consequently select the best shape parameter that give small root means square error (RMSE). All the computational work has been done using Matlab. The interpolant for the selected problems and its corresponding root means square errors (RMSEs) are tabulated and plotted.


Introduction
R BF is an important method for approximating and interpolating functions and solving differential equations (both ordinary and partial) for data with scattered node locations and in computation domain in higher dimensions as discussed in [1]. RBF methods were first introduced by [2], who was a geodesist at Iowa State University. In his proposed method, he simplified the computation of scattered data problem relative to the previously used polynomial interpolation. In many cases, RBF methods have been proved to be effective, where polynomial interpolation fails. RBF method are often used in topographical representation and other intricate 3D shapes.
The use of RBF to interpolate scattered data problem emanated from the fact that the interpolation problem becomes insensitive to the dimension (d) of the space in which the data set lies. Recently, different RBF based methods have gained attention in scientific, computing and engineering applications such as functions interpolation [3,4], numerical solutions to partial differential equation (PDEs) [5][6][7], Multivariate scattered data processing [8] and so on. The main advantages of this method are spectral convergence rates that can be achieved using infinitely smooth basis functions, geometrical flexibility and ease of implementation [3,9].

Idea of radial basis functions
The origin of RBFs can be traced back to [2], when Hardy introduced RBF multi-quadratics (MQ) to solve surface fitting on topography and irregular surfaces. As a field engineer from 1947 to 1951, Hardy was first interested in stream and ridge lines. Hardy believed that an interpolation method containing an exact fit of data to a topographical region should exist. After thorough investigation, he discovered what later become multiquadric (MQ). Prior to the knowledge of MQ, trigonometric and algebraic polynomials were employed.
Due to the quadric surface, Hardy termed his method multiquadric. The MQ transformed scattered data to a very accurate fit model of graph or surface. Subsequently he considered the MQ as a consistent solution of biharmonic solution [10] to physically demonstrate its bravery. Other researchers like [11] extensively tested twenty-nine different algorithms on the typical benchmark function on interpolation problems and ranked the MQ-RBF and thin-plate splines TPS-RBF (introduced in [12] via employing the minimum bending energy theory of the surface of a thin plate) as top two candidates based on the following criteria: timing, storage, accuracy, visual pleasantness of the surface, and ease of implementation. Following the success recorded on researches, RBFs have become popular in the scientific computing world, such as computer graphics, data processing, and economics.
In 1986, Charles Micchelli developed the theory behind the MQ method, and proved that the MQ method matrix system is invertible [13]. Few years later, [14,15] a collocating method was derived for solving partial differential equations using MQ. The breakthrough by Kansa triggered a research boom in RBF. The Kansa method is meshless and superior to the classical method due to some advantages such as: integration-free, ease of implementation and superior convergence. In addition, it may be of interest that even before Kansa's pioneer work, Nardini and Brebbia (1983) [16] applied the function 1 + r without prior knowledge of RBFs as the basis function in the Dual Reciprocity Method (DRM) in the context of the boundary element methods (BEM) to get rid of domain integral effectively. Over the years, several researchers proposed different methods which include:

Radial basis interpolation
RBF interpolation method use linear combinations of N (where N is the number of data sites), radial functions ϕ(r), where r = ||x − x i || and x i are the centers in R d , where d is the dimension of the problem. In this section, we are considering radial basis expansion to solve a scattered data interpolation problem in R d . It is further assumed that The coefficients a k can be obtained by activating the interpolation conditions and when N = Υ, (2) becomes N × N system of linear equations, that is: where A is the following interpolation matrix: At evaluation point Φ, on the set Ξ := {λ j } Φ j=1 , the Equation (1) becomes:

Roles of shape parameter ( )
Shape parameter play a very important role in the determination of the shape of basis functions [3,[22][23][24]. It is obvious from its name (shape parameter) that if the value of is moved over the interval [0, ∞], the shape of the radial function will change effectively, (see Figure 1).

Figure 1. Shapes of 1D Gaussian and IMQ RBF of different shape parameters
To have a better knowledge of how the parameter changes the shape of a function, we present the following theorem of strictly positive definite (SPD) functions [25,26]: Observe that we choose ϕ so that it can be defined over the entire region R.
In summary, it can be deduced from the above theorem that −ϕ(0) Since we are considering basis functions that are strictly positive definite (SPD): IMQ, Wedland and Gaussian, the Theorem 1 gives new bounds as seen in Equation (6).
Note that φ( ; r) = φ(1; ; r). If we take these radial basis as φ( ; r), as a result, the bounds explain how would behave on the shape of the function. i.e, increment in the value of leads to a continuous alter of the shape of the function to a spiky shape as shown in Figure 1, where the spike is denoted by φ(0). Conversely decrements in the value of an leads to a continuous change in the shape of the function to a flat shape, having limit at φ(0) as well, see Figure 1 for explanation.
Earlier in this section, we mentioned the role of the shape parameter ( ). Several researchers worked in this direction [3,[22][23][24] but only few have discussed the role of shape parameter on the accuracy and shape of the interpolant. Numerous articles and books have been dedicated to explain how vital is to determine the shape of a radial function therefore it is obvious that the shape parameter play a vital role in the accuracy and the shape of an interpolant, therefore, the shape parameter need to be chosen wisely.
The question is: how can we choose the optimal shape parameter * with least root mean square error (RMSE)? The answer of this question will be discussed later in this section.
In an attempt to determine the interval in which the optimal shape parameter exists, different algorithms have been provided by different researchers but the purpose of this research work is to provide a Trial and Error algorithm. This algorithm explain explicitly the effect of shape parameter on the the solution of PDE problems selected for interpolation and narrows the interval in which the optimal shape parameter * is obtain. This algorithm will be discussed in Subsection 3.2.

Algorithm for selecting interval of the shape parameter
The determination of best shape parameter * received attention from various researchers like [23] used trigonometric shape parameter to generalize Multiquadric RBF. Similarly, [22] showed that random variable shape parameter yield the most accurate results provided that the centers are uniformly spaced. An algorithm that determine the interval for Multiquadric shape parameter was developed in [3]. This research work, will concentrate on writing Trial and Error algorithm for choosing optimal shape parameter and its implementation on PDE problems. In [11,27], authors described the Leave-One-Out Cross Validation and while in [28,29], the authors introduced the Contour-Pade algorithm.

Algorithm 1: The Algorithm
Input: f : problem to be interpolated, Ω : Data sites, Φ: the set of evaluation points, d: dimension data sites and the problem Output: * : best (optimal) shape parameter 1 Define: the data sites Γ from Ω; 2 Define: the range for the choice of ∈ [a, b]; 3 Define: number shape parameters to be tested, namely n; 4 Define: number of data sites to be use, namely n; 5 for j = 1, . . . , n do 6 ν j = 2 j + 1 It can be observed that the lower values of result in a ( f latter) basis function and the higher values of result in a (spikier) basis functions and = 0 result to constant basis functions. Accuracy and stability are major concerns in numerical computations. It is well known that instability increases as decreases. The significance of this statement lies in the fact that highest accuracy is often found at some small shape parameters, which may be in an unstable region. "Trade-off principle" is the term used to denote the conflict between accuracy and numerical stability [25]. The choice of the basis function is linked with this "trade-off principle". Recently, different methods have been introduced to stably compute interpolants with small parameters by using an alternate basis, such as Contour-Pade and RBF-QR, as in [30] and [31].
Trial and Error Algorithm: To use Trial and Error algorithm, the PDE problems must be transform to matrix form i.e., Ab = f , the generating data must be known [29]. If the function to be interpolate is unknown, it becomes very difficult to choose the best shape parameter. One of the easiest way to choose ( * ) is to perform a series of interpolation experiments on a large range of shape parameters then narrow it down to some selected shape parameters relative to least RMSE. Trial and Error algorithm enhances the clarification effectively on how shape parameter acts on solution of the PDE problems and gives an insight on how to narrow the interval of .
This scheme works only when the problem to be interpolated is known. We have the Algorithm 1:

Root mean square error (RMSE)
Define where α j ∈ M, and M is the numbers of evaluation points. From Equation (7), the RMSE is said to be the square root of mean squared error, where α j , j = 1, · · · , M are evaluation points and should not be mistaken for data sites; data sites errors usually disappears since we are considering an interpolant.
Due to the fact that the RMSE allows mapping of every interpolant to a single value which consequently allows comparison of results in a single way. The RMSE is in some way superior and more practically than the point-wise error.

Illustrative examples
In this section, the Trial and Error algorithm has been employed on some selected PDE problems in R 2 using inverse multiquadric (IMQ) RBF with evaluation points M = 1600 and with different knots N by letting Ω = [0, 1] 2 . The Root Mean Square Error is calculated by the formula given in (7).
with boundary conditions u(x, 0) = u(x, 1) = u(0, y) = u(1, y). The exact solution is u(x, y) = 16x(1 − x)y(1 − y). Table 5 presents the RMSE for different data sites which correspond to Figure 8. The shape parameter that best fits the problem is tabulated in Table 6 and its corresponding figure is shown in Figure 9.

Discussion of Results
Tables 2-5 are the product of the errors obtained from interpolating Examples 1-4. Table 6 summarizes all the RMSE for the examples considered with their respective best shape parameters. The data sites used in this research work are N = [9, 25, 81, 289, 1089, 4225] (gotten from 2 k + 1 d , k = 1, 2, . . .). It is visually obvious from the tables that the errors follow a pattern irrespective of the examples (that's the errors decline as the data sites increases) as shown in Table 6. In Example 3, it can be noticed that with small data sites; [9,25,81] (that is 2 k + 1 d , k = 1, 2, 3 and d = 2) the RMSE are very bad, irrespective of the choice of shape parameter, but as the data sites increases (that is, 2 k + 1 d , k = 4, 5, . . . and d = 2) the RMSE improves. In fact, the best shape parameter exists at ( * = 0.72).

Conclusion
In this research work, Trial and Error scheme is proposed to determine the range in which the best shape parameter * lies via IMQ radial function by performing interpolation experiments and selecting which experiment produces the least RMSE for different data sites. We have illustrated the formulation and implementation of this process on some selected PDE problems. Obviously, numerical results demonstrated that the scheme alongside the Matlab codes has immensely reduced time and ease the complexity involved in the computation, thus, making this scheme an effective and efficient technique for solving partial differential equations.