An Adaptive Selection Method for Shape Parameters in MQ-RBF Interpolation for Two-Dimensional Scattered Data and Its Application to Integral Equation Solving

: The paper proposes an adaptive selection method for the shape parameter in the multiquadratic radial basis function (MQ-RBF) interpolation of two-dimensional (2D) scattered data and achieves good performance in solving integral equations (O-MQRBF). The effectiveness of MQ-RBF interpolation for 2D scattered data largely depends on the choice of the shape parameter. However, currently, the most appropriate parameter is chosen by empirical techniques or trial and error, and there is no widely accepted method. Fourier transform can linearly represent 2D scattering data as a combination of sine and cosine functions. Therefore, the paper employs an improved stochastic walk optimization algorithm to determine the optimal shape parameters for sine functions and their linear combinations, generating a dataset. Based on this dataset, the paper trains a particle swarm optimization backpropagation neural network (PSO-BP) to construct an optimal shape parameter selection model. The adaptive model accurately predicts the ideal shape parameters of the Fourier expansion of 2D scattering data, signiﬁcantly reducing computational cost and improving interpolation accuracy. The adaptive method forms the basis of the O-MQRBF algorithm for solving one-dimensional integral equations. Compared with traditional methods, this algorithm signiﬁcantly improves the precision of the solution. Overall, this study greatly facilitates the development of MQ-RBF interpolation technology and its widespread use in solving integral equations.


Introduction
The prevalence of scattered data problems is increasing in industries such as engineering design and financial analysis due to the consistent evolution of science and technology. One well-known mesh-free method that is typically utilized to handle this data is radial basis function interpolation. Frank [1] conducted numerous scattered data experiments comparing the accuracy of 29 interpolation methods and concluded that the MQ-RBF interpolation method is the most accurate. It has been pointed out in numerous studies that the accuracy of interpolation is heavily influenced by the shape parameter of the MQ-RBF [2,3].
The MQ-RBF interpolation method was initially proposed by Hardy [4] with the selection of a shape parameter, c = 0.815d. Here, d = 1 N ∑ N i=1 d i where d i represents the distance between the ith data point and its closest neighbor. This method provides an adequate fitting effect on terrain problems, leading to numerous researchers exploring the selection of shape parameters in MQ-RBF interpolation. The leave-one-out cross-validation (LOOCV) method, introduced by Rippa [5], has proven to be most influential and effective among these methods. Rippa proposed a cost function to represent the root mean square error (RMSE) between the interpolation function and the original function. The mnbrak and brent in [6] were then used to determine more suitable shape parameters that minimize the cost function. The effectiveness of the method was verified from multiple perspectives, method to identify suitable shape parameters for MQ-RBF and applies it successfully to one-dimensional integral equations. The performance of this method is evaluated through numerical simulation of various one-dimensional integral equations.

MQ-RBF Interpolation
According to E.M. Stein and G. Weiss [16], a radial basis function (ϕ(x)) is a real-valued function whose value is solely determined by the distance from the origin. If |x 1 | = |x 2 |, then ϕ(x 1 ) would be equal to ϕ(x 2 ). Table 1 illustrates the commonly employed radial basis functions.

Name ϕ(r)
Gaussian The function is defined as follows [26]: Here, λ j is the jth weight of the sample point, and N refers to the number of sample points. ϕ j (r) represents the basis function, as given by c is a shape parameter. When employing MQ-RBF interpolation, it determines the efficacy of the interpolation. x j is the jth sample point. As the basis functionf (x) passes through the sample points, we obtain the equation: The basis function matrix Ψ N×N can be defined as shown below: where ϕ ij is the basis function about the sample points and the jth sample point: Furthermore, let W denote the weight vector λ j and F denote the vector of f (x j ); we end up obtaining

Algorithm Selection
The present study employs optimization algorithms to determine the shape parameters of MQ-RBF interpolation for sine functions. This approach offers several advantages, including cost reductions as well as improved accuracy of interpolation. The optimization problem that represents the selection of shape parameters in MQ-RBF is as follows: The parameters of the interpolation include the interpolation basis function s(x, c), the interpolation primitive function f (x), the maximum error Emax(c), and the required optimal shape parameter copt for the interpolation.
To determine the ideal optimization algorithm for determining MQ-RBF's interpolation shape parameters in sine functions, we utilized different optimization approaches [27][28][29][30][31], such as Gradient Descent (GD), Newton-Raphson method (NR), Genetic Algorithm (GA), Tabu Search (TS), and Random Walk (RW), for MQ-RBF interpolation of numerous sine functions. The leave-one-out cross-validation method [5] was utilized in our study to select the initial shape parameter of Function (9), which yielded a value of 0.4133. After a series of experiments, we identified the optimal settings for the initial shape parameter optimization using different algorithms. These settings consist of a learning rate of 0.1 for GD, a population size of 10 for GA, a taboo length of 10 for TS, and 10 walks for RW. We compared the performance of these algorithms in terms of interpolation accuracy, computation time, and the number of iterations required to reach the optimal shape parameters, with all algorithms being set to a maximum iteration of 20. The results provided us with important insights regarding the optimization algorithms' capability to identify the optimal shape parameters. Table 2 highlights the results of Function (9). According to our experimental results, GA and TS were found to produce relatively small shape parameter errors and to require fewer iterations to identify the optimal shape parameters compared to other algorithms when initially configured with the same number of iterations and shape parameters. However, these algorithms require more computation time. In contrast, GD and NR exhibit faster training but produce larger shape parameter errors when the maximum iteration is reached. On the other hand, while ensuring high interpolation accuracy, RW requires a relatively short computation time. Therefore, we recommend using the RW to determine the shape parameter of the MQ-RBF interpolation function for the sine function.
The selection of the optimal shape parameter (c opt ) using the Random Walk (RW) algorithm involves the following steps: Step 1: Define i(i = 1, 2 . . . , M) as the number of walks, k(k = 1, 2 . . . , N) as the number of current iterations, the accuracy θ for step control, and the accuracy ε for error control. Set k equal to 1 and establish the initial parameter c 0 .
Step 3: Compute the value of Emax(c i ): (1) If Emax(c i ) < Emax(c i−1 ), the ith step is completed. Take c i as the new initial parameter, reset k to 1, and begin the next walk. The walk process is repeated until Emax(c) < ε or i = M, at which point the algorithm ends.
(2) If Emax(c i ) > Emax(c i−1 ), it indicates that no better parameter than the present one has been found. If k < N, return to step 2 to regenerate the random vectors u ik+1 , · · · , u iN−1 and continue the search. When k = N and no better parameter is found, the optimal parameter c opt is regarded as in the sphere with a center c i−1 and a radius λ. If λ < θ, end the algorithm; otherwise, set λ = λ 0 /2, go back to step 1, and initiate a new round of walking.

Improved Random Walk Algorithm
The Random Walk algorithm, however, exhibits some issues in finding parameters. If a superior parameter is discovered in the initial parameter's neighborhood, the algorithm will advance to the next walk, regardless of whether the iteration meets the specified N times or not. As a result, the outcome may regress into local optimization.
We have made enhancements to the Random Walk algorithm and labeled it the Improved Random Walk Algorithm (IRW). The enhancements are as follows: every walk is iterated for N times, and the parameter registered with the corresponding minimum error in this walk is taken as the starting parameter for the next walk. By incorporating this improvement, the algorithm covers a wider parameter range and offers more directions. Figure 1 displays a flowchart of IRW for determining the best shape parameters. Table 3 presents the interpolation data for the optimal shape parameter selected by the IRW algorithm based on Equation (9) for the purpose of comparison with the results in Table 2. Our validation process has repeatedly shown that the IRW algorithm can identify the optimal shape parameter with minimal iterations, thereby enhancing the interpolation accuracy without significantly increasing the time cost. Furthermore, Figures 2 and 3 demonstrate the impact and absolute error of the MQ-RBF interpolation based on the ideal shape parameters established by the IRW algorithm in Equation (9). These figures indicate that the chosen shape parameters have excellent interpolation effects. The low cost and high precision of the IRW algorithm make it the most suitable choice for our problem as we need to accumulate a large number of data.

The Relationship between ω and c opt
Sine and cosine functions are collectively referred to as sine functions in practical applications. Their general function expression is as follows: The expression for trigonometric functions includes four parameters, amplitude, offset, initial phase, and angular frequency, denoted as A, B, ϕ, and ω respectively. These parameters determine the basic shape of the trigonometric curve. As stated by [7], the basic shape of the MQ-RBF is determined by its parameter c. Our numerous experimental results indicate that A, B, and ϕ have negligible influence on the copt, while ω exerts a profound impact on the copt [32]. Consequently, the IRW algorithm is employed to explore the relationship between ω and the c opt .
Let ω = kπ(k = 2, · · · , 10) in Equation (11); i.e., expand the angular frequency of Equation (9) by a factor of k. Selected experimental results are presented in Table 4. According to the experimental results, changes in ω had a minor impact on the interpolation accuracy. When the angular frequency is multiplied by a factor of k, the corresponding c opt decreases by approximately k. Further verification was performed by dividing the ω by a factor of k (k = 2, · · · , 10), and some of the experimental results are presented in Table 5. Table 5. ω = π k (k = 2, · · · , 10) experimental results. Numerous numerical experiments have demonstrated an approximate inverse proportionality between variables ω and c opt for trigonometric functions. The IRW algorithm is employed to select parameters for every individual trigonometric function. Subsequently, the MQ-RBF interpolation shape parameter selection formula of a trigonometric function is fitted using the least-square method [33] based on a large number of data points that correspond to a one-to-one relationship with respect to the results. Figure 4 presents the fitting image of some data, and the resulting formula is Equation (12).

Establishment of the data set and selection of regression model
A linear combination of trigonometric functions can be expressed mathematically as follows [34]: The number of terms in the linear combination of trigonometric functions is denoted by N, and the amplitude, angular frequency, and initial phase of the kth term's trigonometric function are represented by A k , ω k , and θ k , respectively. The c opt in the MQ-RBF interpolation for linear combinations of trigonometric functions is determined by the IRW algorithm. Experimental research has shown that the angular frequency is the primary determinant of the corresponding c opt . Specifically, the highest angular frequency term significantly affects c opt for MQ-RBF interpolation in linear combinations of trigonometric functions. Based on these results, a dataset is proposed that includes the angular frequency of the linear combination of trigonometric functions along with their corresponding c opt .
Using the Pandas library in Python 3.10, the angular frequencies and corresponding optimal shape parameters for 1 million linear combinations of sine functions were divided into segments. To generate a training and testing dataset, we conducted three train-test splits with ratios of 7:3, 6:4, and 9:1, respectively. A 6:4 ratio is more suitable for smaller datasets since it can help prevent overfitting. A 7:3 ratio ensures model accuracy while avoiding overfitting and underfitting, making it best suited for moderate-sized datasets. A 9:1 ratio allocates more data for model training, improving model accuracy by allowing for a better understanding of the dataset's characteristics and patterns. Given the large size of our dataset, we validated and compared the different ratios, ultimately selecting the 9:1 ratio as the most appropriate for our needs.
We trained five models [35][36][37][38][39], namely, Back Propagation Neural Network (BP), Multiple Linear Regression (MLR), Gated Recurrent Unit (GRU) networks, Support Vector Machine (SVM), and Long Short-Term Memory (LSTM), using 900,000 data points as the training set. We compared and evaluated the models using 100,000 data points as the test set. We used three evaluation indices, namely training time, mean square error (MSE), and prediction accuracy. Refer to Table 6 for results. Our experimental results show that, despite its shorter training time, MLR exhibits the poorest predictive accuracy, suggesting that there is no clear linear relationship between the data. The performance of SVR in handling large-scale samples results in average training and predictive accuracy. Compared to SVR and LSTM, the BP predicts the shape parameters of linear combinations of trigonometric functions with the highest accuracy, with the same amount of training time. After a comprehensive comparison, we selected the BP to perform shape parameter prediction.

Construction of the c opt Selection Model Based on PSO-BP
The BP is composed of three layers: the input layer, the hidden layer, and the output layer. The signal transmission in the BP progresses forward sequentially through the input layer, hidden layer, and output layer, while its error is propagated backward, starting from the output layer, then the hidden layer, and finally the input layer. The learning ability of the neural network is directly impacted by the number of nodes in the hidden layer, and the formula used to calculate it is as follows: where h represents the number of nodes in the hidden layer, l represents the number of nodes in the input layer, and e is a constant in the range of [1,10] and takes an integer value. The architecture of our BP neural network includes two hidden layers, with 80 and 30 neurons in the first and second hidden layers, respectively. We selected Rectified Linear Unit (ReLU) as our activation function while setting the learning rate to 0.05 and the momentum to 0.9. Particle Swarm Optimization (PSO) [40] is an optimization algorithm that imitates the predation behavior of birds. Unlike the random gradient descent method used in BP training, PSO is a global optimization algorithm that can find the optimal solution to a problem in the entire region. The algorithm generates a set of random solutions and then updates the particle velocity and position in each iteration to find the optimal solution. The rules for updating the particle velocity and position are as follows: The variables V k i and X k i denote the velocity and position of particle i during the kth iteration. The variable ω corresponds to the inertia factor. The variable c 1 corresponds to the individual learning factor. The variable c 2 corresponds to the social learning factor. The variables r 1 and r 2 correspond to random numbers in the range [0,1]. The variables p best and g best indicate the best positions found so far by the current single particle and by all particles, respectively.
The BP can perform the nonlinear mapping from input to output, but it is prone to reaching a local minimum after a specific number of iterations. PSO can fully utilize the nonlinear application of BP and overcome the issue of the slow convergence of weights in the BP training neural network, which can easily cause it to fall into local optima. Table 7 presents the evaluation results of the Particle Swarm Optimization Backpropagation model (PSO-BP). Without imposing significant time costs, the model demonstrates significant improvements in MSE and prediction accuracy. In our study, we propose using Particle Swarm Optimization (PSO) as the optimizer to enhance the effectiveness of our model. The PSO algorithm is configured with the following parameters: swarm size, maximum number of iterations, inertia weight, cognitive learning factor, and social learning factor. We set the swarm size to 64, the maximum number of iterations to 100, the inertia weight to 0.8, the cognitive learning factor to 1.5, and the social learning factor to 2.0; these parameters were selected based on prior research and our own experimentation to ensure optimal model performance. Utilizing the configured PSO optimizer aims to maximize the accuracy and efficiency of the model and achieve improved results. As shown in Figure 5, a comparison between predicted and actual values of the partial data reveals that the model's predicted values closely align with the actual values, thus demonstrating strong predictive performance. Figure 6 illustrates the loss function of the PSO-BP that reveals the gradual approach of both the training and testing loss functions to zero with an increase in training iterations. This indicates enhanced predictive accuracy of the model as the number of iterations is increased.

Verification Experiment
This section focuses on determining the optimal shape parameters for the MQ-RBF interpolation applied to various sine functions (given as test functions in Table 8). We accomplish this by utilizing the Formula (12) presented in Section 3.1 and the model in Section 3.2.2. We also provide a direct comparison of the obtained results with those obtained through IRW to ensure the accuracy of our method. Detailed comparison results are shown in Table 9.
According to the experimental findings, the predicted c opt and the corresponding MaxError are relatively consistent with the algorithm's direct outcomes. Figure 7 displays the MQ-RBF interpolation effect for functions in Table 8 employing optimal parameters selected by the model. The diagram indicates that the interpolation effect of the predicted optimal shape parameters is satisfactory, no matter what kind of linear combination of sine functions is chosen.

Fourier Expansion of 2-D Scattered Data
Assuming f (x) is a periodic function with a period of 2L and satisfies the Dirichlet convergence condition, its Fourier expansion can be expressed as follows: In this equation, a q and b q represent the Fourier coefficients: Moreover, assuming a non-periodic function f (x) is defined within the interval [−L, L] and satisfies the Dirichlet convergence condition [41], it can be expanded into a Fourier series by employing periodic continuation. Based on the aforementioned theory, it can be concluded that, under certain conditions, any given 2D scattered data or continuous function over a specified interval can be represented as a linear combination of sine and cosine functions through the Fourier transform.

Adaptive Selection Method of the Shape Parameter in the MQ-RBF Interpolation for 2D Scattered Data
We propose an adaptive selection method for shape parameters in MQ-RBF interpolation of 2D scattered data by combining the theory of Fourier series and the optimal parameter selection model for the sine function and its linear combination constructed in Section 3.2.2. The steps are as follows: Step 1: Utilizing the Fourier series for fitting two-dimensional scattered data points, we acquire the corresponding Fourier expansion.
Step 2: We use the MQ-RBF interpolation shape parameter selection model we provide for the sine function and its linear combination, based on the Fourier expansion, to predict the corresponding optimal shape parameters.
Step 3: We use the MQ-RBF interpolation shape parameters predicted by the model to interpolate the original two-dimensional scattered data.
Our adaptive method eliminates the need for selecting initial shape parameters, resulting in reduced accuracy loss during iteration and greatly reduced operating costs. Instead, we only need to perform Fourier expansion on the sampling point data and use the established model to predict the appropriate shape parameters directly.
In this study, we generate 2D scattered data points from seven theoretical functions selected from [8,[10][11][12]16] (Table 10). We use the adaptive method mentioned above to calculate these data points and compare the results with those obtained using Rippa's algorithm (Table 11). f Interval Data Point 1+25x 2 x ∈ [5, 10] 67 x ∈ [0, 1] 72  Table 11 demonstrates an improvement in interpolation accuracy and a significant reduction in operation costs compared to Rippa's algorithm. These results provide strong evidence for the effectiveness of the adaptive method proposed in this paper.
The assessment of optimal shape parameters for MQ-RBF interpolation is critical to the accuracy of the approximation results. Given that the optimal shape parameters can vary depending on the nature of the functions, it is essential to propose suitable methods for each function. Thus, we utilize the proposed adaptive method, which can obtain optimal shape parameters regardless of any domain or range considerations of functions.
The formula (12) is effective in determining the optimal shape parameters from a geometric perspective. By analyzing the errors between the original and the interpolation functions, we can acquire information about the ideal shape parameters. The methodology proposed in Section 3.2.2 provides a more robust approach that considers the situations where the optimal parameters may not be evident. The proposed methodology expands on an arbitrary sine function using the Fourier series and determines the optimal parameters while minimizing the mean squared error between the original function and the interpolation function.
In summary, our approach combines the geometry of MQ-RBF, the adaptive strategy, and the Fourier series expansion method to obtain optimal shape parameters for sine functions. These methods ensure higher accuracy in function approximation and enable successful applications of the MQ-RBF to a wide range of fields.

Application of the Adaptive Method in Solving One-Dimensional Integral Equations
The procedure for solving an integral equation using the MQ-RBF method is similar to that of solving a differential equation using the same method. Firstly, we approximate the unknown function by a linear combination of MQ-RBF functions, which we then substitute into the integral equation. Next, we determine the weight coefficients using the collocation point method and obtain an approximate solution for the unknown function. Unlike differential equations, the collocation point equation is represented by integral formulas containing MQ-RBFs rather than the differential equation differentiated at collocation points. The solution process for linear and nonlinear integral equations with one variable using MQ-RBFs is presented in detail below.

MQ-RBF Collocation Approximation for One-Dimensional Linear Integral Equation
The one-dimensional linear integral equation in its general form can be expressed as where Selecting J collocation points results in a collocation equation represented in matrix and K ij represents a definite integral, which can be evaluated using the Gaussian quadrature formula: The integration interval is transformed to t = p(ξ) = (b+a) 2 + (b−a)ξ 2 to calculate the coefficient K ij , which can be approximated as The collocation points x j and RBF center points xi are assumed to be identical in this paper, resulting in a square matrix for Ψ.

MQ-RBF Collocation Approximation for One-Dimensional Nonlinear Integral Equation
The one-dimensional nonlinear integral equation can be expressed in general as shown below: To represent the unknown function f (x) in terms of MQ-RBF, we use MQ-RBF approximations, which results in the following equation: Next, consider a specific collocation point x j ∈ [a, b]; Equation (24) becomes By choosing J collocation points, we arrive at a nonlinear system of equations in matrix form: a h(t)dt is a non-linear function with respect to the weight coefficient [W]. Therefore, iterative techniques are required to solve it. In each iteration step, all terms in Equation (27) can be computed based on the current value of [W] = [λ 1 , λ 2 , · · · , λ N ] T . For instance, the coefficient K j can be expressed by using the Gaussian quadrature formula:

Solving One-Dimensional Integral Equations Using the MQ-RBF Method with an Optimal Shape Parameter
The MQ-RBF method with optimal shape parameters (O-MQRBF) utilizes the following steps to solve one-dimensional integral equations.
Step 1: Set the MQ-RBF and configuration point by selecting the center point x i of MQ-RBF and the configuration point x j on the domain of the function. The optimal shape parameter α can be determined using the method mentioned in Section 4.2.
Step 2: Construct collocation equations by substituting MQ-RBF interpolation format into integral equations and applying collocation conditions to form collocation equations that contain definite integrals. The matrix form of the linear integral equations is given by Step 3: Calculate the matrix elements. For one-dimensional rectangular areas, the Gaussian integral method can be applied to calculate definite integrals.
Step Step 5: Calculate the RBF approximation of the function

Numerical Example
The following section presents four examples of linear and nonlinear integral equations with a single variable that include both Fredholm and Volterra forms, making the examples more convincing. The setting of MQ-RBF parameters for each example is also provided. It is worth mentioning that we determined all shape parameters using the adaptive selection method developed in Section 5.3.

1.
One-dimensional Fredholm linear integral equation: The exact solution to this equation is f (x) = x. The integration of the interpolation coefficients is computed using a 20-point Gauss integration. Our adaptive method was utilized to select a shape parameter of c = 21.35, with a center distance of h = 0.1, and 11 points share the same center, h t = 0.001.

2.
One-dimensional Volterra linear integral equation: The exact solution of the function is f (x) = cos(x). The interpolation coefficients are integrated using a 60-point Gauss integration. The center points of MQ-RBF are set to h = 0.1. Following the application of our adaptive method, a suitable shape parameter is determined to be c = 1.88, with both the center and collation point taking position 11. We select a total of 501 measuring points at an interval of h t = 0.002. 3.
One-dimensional Fredholm nonlinear integral equation: The exact solution of this equation is f (x) = 2 − x 2 . The integration of the interpolation coefficients is computed using a 10-point Gauss integration. The center points of MQ-RBF are set to h = 0.1. Utilizing our adaptive method, we settled on a shape parameter of c = 2.75, with the center point and the collocation point being 10, and h t = 0.001. 4.
One-dimensional Volterra nonlinear integral equation: The exact solution for this equation is f (x) = e −x . The integration of the interpolation coefficients is computed using a 10-point Gauss integration. The shape parameter c = 1.25 was chosen using our adaptive strategy, and h t = 0.005 between the measuring points.
The high accuracy of the Haar wavelet method in [18] and the Maleknejad method in [20,21] commonly used in solving one-dimensional integral equations cannot be overlooked. Nonetheless, to illustrate the undeniable superiority of O-MQRBF, an accuracy comparison was conducted with those two methods for four different examples, as presented in Table 12.
Our experimental results demonstrate that the O-MQRBF method enhances the accuracy in solving one-dimensional integral equations, affirming its effectiveness. Simultaneously, the results also indicate that our proposed adaptive shape parameter selection method is both convenient and effective, further adding to its potential for widespread application.

Conclusions
A method utilizing MQ-RBF interpolation with adaptive shape parameters is proposed for solving one-dimensional integral equations. When 2D scattered data can be expanded into a linear combination of sine and cosine functions via a Fourier series, the angular frequency of the sine function and its linear combination is observed to significantly impact the optimal shape parameters. Thus, we developed an optimal shape parameter selection model for both the sine function and its linear combination. By comparing and verifying the results, our adaptive model accurately selects shape parameters for the Fourier expansion of 2D scattered data while also reducing running costs. We applied the adaptive method to solve one-dimensional integral equations and conducted comparative experiments, which demonstrate notable enhancements to interpolation accuracy. This paper offers a practical and effective technique for solving one-dimensional integral equations while providing a convenient approach for picking shape parameters in MQ-RBF interpolation of 2D scattered data. Future research can extend our approach to higher dimensional data.