Solving Stochastic Nonlinear Poisson-Boltzmann Equations Using a Collocation Method Based on RBFs

Abstract

In this paper, we present a numerical scheme based on a collocation method to solve stochastic non-linear Poisson–Boltzmann equations (PBE). This equation is a generalized version of the non-linear Poisson–Boltzmann equations arising from a form of biomolecular modeling to the stochastic case. Applying the collocation method based on radial basis functions (RBFs) allows us to deal with the difficulties arising from the complexity of the domain. To indicate the accuracy of the RBF method, we present numerical results for two-dimensional models, we also study the stability of this method numerically. We examine our results with the RBF-reference value and the Chebyshev Spectral Collocation (CSC) method. Furthermore, we discuss finding the appropriate shape parameter to obtain an accurate numerical solution besides greatest stability. We have exerted the Newton–Raphson approach for solving the system of non-linear equations resulting from discretization by the RBF technique.

1. Introduction

There are numerous numerical methods for solving partial differential equations (PDEs). A popular category of numerical techniques is called mesh-free methods [1,2]. As a special case of the mesh-free methods, we recall radial basis functions (RBFs) collocation methods, which recently have been studied to solve a large number of PDEs related to science and engineering problems [3,4,5,6]. The implementation of RBF methods is very simple because the main property of the RBF method is that in order to discretize it does not require a mesh, indeed it only requires setting collocation points which allow us to deal with complicated domains. Furthermore, the implementation of this method is easy for higher dimensional problems [7,8,9]. So far, many authors have published papers on the improvement of the RBF method or on how to implement the method for various PDE problems [10,11,12]. It is worth noting that this method has been applied to the study of several biological systems [13,14]. We refer interested readers to [15] for new advances in RBF methods. Several radial basis functions have been defined so far, and some of them depend on a constant which is called a shape parameter $c>0$. Some popular RBFs are collected in

Table 1. Popular radial basis functions.

Name of Function	Abbreviation	Definition
Multiquadric	MQ	$\phi \left(r\right)=\sqrt{c^2+r^2}$
Inverse multiquadric	IMQ	$\phi \left(r\right)=\frac{1}{\sqrt{c^2+r^2}}$
Inverse quadratic	IQ	$\phi \left(r\right)=\frac{1}{c^2+r^2}$
Gaussian	GA	$\phi \left(r\right)=exp(-c^2{r}^2)$
Polyharmonic spline	PHS	$\phi \left(r\right)=\left\{\begin{array}{cc}{r}^{2n-1},\hfill & n\in \mathbb{N},\hfill \\ r^{2n}lnr,\hfill & n\in \mathbb{N}\hfill \end{array}\right.$
Thin plate spline	TPS	$\phi \left(r\right)=r^{2}logr$

.

The accuracy and stability of the numerical results obtained using the collocation method are highly dependent on the choice of the collocation points and the shape parameter. A smaller shape parameter generally improves the accuracy, but it may lead to instability. Thus, there exists a trade-off between the stability and accuracy of the approximation, and determining the best shape parameter remains a challenging problem. Consequently, the value of the shape parameter is included in the unknowns and needs to be determined. The shape parameter that leads to the lowest error is called the optimal shape parameter. Various techniques have been proposed by researchers to find the optimal shape parameter [16,17,18,19,20]. When describing some physical phenomena in the form of partial differential equations, adding stochastic factors enables the achieving of a more accurate model. In recent years, many researchers have become interested in studying stochastic models [21,22]. In this work, we apply a collocation method based on radial basis functions to find the discrete solutions of a stochastic non-linear Poisson–Boltzmann equation. This class of stochastic PDEs plays a role in describing various biological phenomena, making it imperative to solve these equations to obtain beneficial solutions. This strong motivation drives us to tackle these complex equations [23,24,25]. Here, we consider the model that is described in [26].

The main focus of this paper is to propose a formulation based on the RBF collocation method to solve stochastic non-linear Poisson–Boltzmann equations. Section 2 is devoted to describing a physical model for the non-linear stochastic non-linear Poisson–Boltzmann equation. In Section 3, we present the main results of this paper which is to implement the RBF collocation method to solve our model problem. In Section 4, we present numerical results, and in the final section we have summarized the conclusions.

2. The Stochastic Non-Linear Poisson–Boltzmann Equation

Let D be an open-bounded Lipschitz domain in ${\mathbb{R}}^d$ consisting of two non-overlapping sub-domains with different physical properties with the boundary $\partial D$. In this paper, we assume $D:=D_L\cup D_M$ where $D_L$ contains a solvent, such as an electrolyte, blood, or salty water, and $D_M$ is the impenetrable molecules region. The domain $D_L$ is filled with two free charge carriers, namely a positive and a negative carrier (see Figure 1). We also assume $\mathcal{P}:=(\Omega ,\Sigma ,\mathit{P})$ is a probability space.

The above process can be described by a deterministic Poisson–Boltzmann equation. This model could be generalized to a stochastic version when the molecules move randomly inside the solvent, which is controlled by Brownian motion. Here, the movements of the molecules are independent of each other. The governing equation for the stochastic Poisson–Boltzmann equation is considered as follows

$$\begin{array}{cccc}\hfill -\nabla ·\left(a\right(x,\omega )\nabla u(x,\omega \left)\right)+b(x,\omega ,u(x,\omega \left)\right)& =f(x,\omega )\hfill & & \mathrm{in}{D}_L,\hfill \end{array}$$

(1a)

$$\begin{array}{cccc}\hfill -\nabla ·\left(a\right(x,\omega )\nabla u(x,\omega \left)\right)& =g(x,\omega )\hfill & & \mathrm{in}{D}_M,\hfill \end{array}$$

(1b)

$$\begin{array}{cccc}\hfill u(x,\omega )& =u_D(x,\omega )\hfill & & \mathrm{on}\partial D,\hfill \end{array}$$

(1c)

where the elliptic operator depends on some coefficients $a(x,\omega )$ for $x\in D$ and the sample $\omega \in \Omega$. Additionally, the force terms $f(x,\omega )$ and $g(x,\omega )$ are assumed to be random. One of the interesting cases of (1a)–(1c) is the following stochastic non-linear Poisson–Boltzmann equation

$$\begin{array}{cccc}\hfill -\nabla ·\left(a\right(x,\omega )\nabla u(x,\omega \left)\right)+\kappa \left(x\right)sinh\beta u(x,\omega )& =0\hfill & & \mathrm{in}{D}_L,\hfill \end{array}$$

(2a)

$$\begin{array}{cccc}\hfill -\nabla ·\left(a\right(x,\omega )\nabla u(x,\omega \left)\right)& =\rho (x,\omega )\hfill & & \mathrm{in}{D}_M,\hfill \end{array}$$

(2b)

$$\begin{array}{cccc}\hfill u(x,\omega )& =u_D(x,\omega )\hfill & & \mathrm{on}\partial D,\hfill \end{array}$$

(2c)

where u describes the electrostatic potential which is due to free charges, $a(x,\omega )$ is the dielectric constant or permittivity function, and $\rho (x,\omega )$ represents the fixed charges of the molecules. Here, $\kappa \left(x\right)$ is the ionic concentration. $\kappa \left(x\right)$, $a(x,\omega )$ are the piecewise constant functions as below

$$a(x,\omega )=\left\{a_{M}x\in D_M,a_{L}x\in D_L,\right.\kappa \left(x\right)=\left\{\begin{array}{cc}0\hfill & x\in D_M,\hfill \\ \overline{\kappa }\hfill & x\in D_L.\hfill \end{array}\right.$$

$\beta$ is defined as $\beta :=\frac{q}{k_{B}T}$, where the constants ${\kappa }_{\mathcal{B}}$, T, and $q>0$ are the Boltzmann constant, the temperature, and the elementary charge, respectively.

This stochastic generalization describes physical systems where the fixed charges have probability distributions, and, hence, the electrostatic potential u is stochastic. The different configurations of the fixed charges are due to the movement of the molecules. In this situation, the input data vary randomly from one point of the physical domain D to another point. In reality, each point in the domain corresponds to a random variable, and the correlation between the values of any two distinct points in the domain is generally non-zero. In the model described here, the positions of subdomains are also stochastic.

3. Collocation Method Based on Radial Basis Functions for Solving

In this section, we present a collocation method based on radial basis functions to solve system (2a)–(2c). For the sake of simplicity, we rewrite the system (2a)–(2c) as

$$\begin{array}{cccc}\hfill \mathcal{L}u(·,\omega )+\mathcal{G}u(·,\omega )& =f(·,\omega )\hfill & \hfill & \mathrm{in}D,\hfill \end{array}$$

(3)

$$\begin{array}{cccc}\hfill u(x,\omega )& =u_D(x,\omega )\hfill & \hfill & \mathrm{on}\partial D,\hfill \end{array}$$

(4)

where the linear elliptic operator $\mathcal{L}$ is defined as

$$\mathcal{L}u(·,\omega ):=\left\{\begin{array}{cc}-\nabla ·(a_L(x,\omega )\nabla u(x,\omega ))\hfill & x\in D_L\hfill \\ -\nabla ·(a_M(x,\omega )\nabla u(x,\omega ))\hfill & x\in D_M\hfill \end{array}\right.,$$

the function $\mathcal{G}$ is a Lipschitz continuous function defined as

$$\mathcal{G}u(·,\omega ):=\left\{\begin{array}{cc}\kappa \left(x\right)sinh\beta u(x,\omega )\hfill & x\in D_L\hfill \\ 0\hfill & x\in D_M\hfill \end{array}\right.,$$

and the source term $f:=(x,\omega ):=\left\{\begin{array}{cc}0\hfill & x\in D_L\hfill \\ \rho (x,\omega )\hfill & x\in D_M\hfill \end{array}\right.$.

Let $\phi \left(r\right):[0,\infty )\to \mathbb{R}$ be a type of radial basis functions defined by ${\phi }_j\left(x\right):=\phi (\parallel x-x_j\parallel )$ at the collocation points ${\left\{x_j\right\}}_{j=1}^n$ belonging to $\overline{D}:=D\cup \partial D$. We approximate u by $\tilde{u}$ defined as a linear combination of the radial basis functions $\phi (\parallel ·\parallel )$, i.e.,

$$u\simeq \tilde{u}\left(x\right)=\sum _{j=1}^n{\lambda }_j{\phi }_j\left(x\right),x_j\in \overline{D},$$

(5)

where the first $n-1$ nodes belong to the interior of the domain and the last node in on the boundary of the domain. Coefficients ${\lambda }_j$ are unknown and should be determined. Thus, the collocation method to solve the system (3) and (4) can be written as find ${\lambda }_j\in \mathbb{R}$, such that

$$\begin{array}{cccc}\hfill \sum _{j=1}^{n-1}{\lambda }_j\mathcal{L}\phi (\parallel x_k-x_j\parallel ,\omega )+\mathcal{G}\left(\sum _{j=1}^n{\lambda }_j{\phi }_j\left(x\right),\omega \right)& =f(x_k,\omega )\mathrm{for}\hfill & \hfill & k=1,\dots ,n-1\hfill \\ \hfill \sum _{j=1}^{n-1}{\lambda }_j\mathcal{L}\phi (\parallel x_k-x_j\parallel ,\omega )+\mathcal{G}\left(\sum _{j=1}^n{\lambda }_j{\phi }_j\left(x\right),\omega \right)& =0\mathrm{for}\hfill & \hfill & k=n,\hfill \end{array}$$

which leads to the following system

$$\left[\begin{array}{cccc}\mathcal{L}\phi (\parallel x_1-x_1\parallel )& \mathcal{L}\phi (\parallel x_1-x_2\parallel )& \cdots & \mathcal{L}\phi (\parallel x_1-x_n\parallel )\\ \mathcal{L}\phi (\parallel x_2-x_1\parallel )& \mathcal{L}\phi (\parallel x_2-x_2\parallel )& \cdots & \mathcal{L}\phi (\parallel x_2-x_n\parallel )\\ ⋮& ⋮& \cdots & ⋮\\ \phi (\parallel x_n-x_1\parallel )& \phi (\parallel x_n-x_2\parallel )& \cdots & \phi (\parallel x_n-x_n\parallel )\end{array}\right]\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ ⋮\\ {\lambda }_n\end{array}\right]+G({\lambda }_1,\cdots ,{\lambda }_n)=\left[\begin{array}{c}f\left(x_1\right)\\ f\left(x_2\right)\\ ⋮\\ 0\end{array}\right],$$

(6)

where $\mathbf{G}({\lambda }_1,\cdots ,{\lambda }_n):=\left[\begin{array}{c}\mathcal{G}\left({\sum }_{j=1}^n{\lambda }_j{\phi }_j\left(x_1\right)\right)\\ \mathcal{G}\left({\sum }_{j=1}^n{\lambda }_j{\phi }_j\left(x_2\right)\right)\\ \cdots \\ 0\end{array}\right]$, and for the sake of simplicity we drop the stochastic variable $\omega$. We then solve the above system using Newton–Raphson’s method. Finally, we define the approximate solution of system (3) and (4) as

$$u\left(x\right)\simeq u_D\left(x\right)+\sum _{j=1}^n{\lambda }_j{\phi }_j\left(x\right).$$

4. Numerical Results

To indicate the effectiveness and efficiency of the RBF method, we present the achieved numerical results for the model equation in Section 3. Consider the 2D version of problem (2a)–(2c) that the forcing term on the right-hand side is considered a stochastic term. Then, by implementing the RBF method, the right-hand vector in system (6) would be random. For presenting numerical results consider the following assumptions and constants

$$u=0\mathrm{on}\partial D,$$

$$a_L=78.4,a_M=3.7,\kappa \left(x\right)=1,$$

$$q=1.6\times {10}^{-19},{\kappa }_B=1.38\times {10}^{-23}\mathrm{and}T=300.$$

Additionally, consider ${\left\{x_i\right\}}_{i=1}^n$ being the collocation points and $h={\displaystyle \underset{i=1\le i\le n}{max}}\parallel x-x_i\parallel \mathrm{for}\mathrm{all}x\in D\cup \partial D$. At first, we present our results for the problem consist of one molecule domain. We provide numerical results using some types of radial basis functions, such as the Multiquadric (MQ), Inverse Multiquadric (IMQ), and Inverse Quadratic (IQ) functions.

There is no exact solution for this model. Once we make sure the algorithm is correct, then we consider the approximation solution for 100 interpolation points using the MQ function as a reference value, which is shown in Figure 2.

Also, for evaluating our results by using the RBF method, we compare our results with the Chebyshev Spectral Collocation (CSC) method. In this approach, the numerical solution is interpolated using Lagrange polynomials at the roots of the derivative of Chebyshev polynomials. The interested reader is referred to [27] for detailed descriptions of the method. We consider the root-mean-square errors

$$\mathrm{RMSE}:=E(\parallel u-\tilde{u}{{\parallel }_{L^2}^2)}^{\frac{1}{2}},$$

where u can be considered the RBF-reference value (R-value) or the approximation solutions using CSC, and $\tilde{u}$ is the numerical solution obtained by the RBF method. To calculate this error we employ the Monte-carlo method with 500 realizations for all numerical results.

In Figure 2, it is easy to observe the electrostatic potential energy in the molecular domain $D_M$. The potential in the liquid region $D_L$ is small. Additionally, it can be seen the nonlinearly changing near the molecule region. These observations are also visible for subsequent numerical solutions. The numerical solution and error completely depend on the value of the shape parameter c. So at first, we discuss finding an appropriate shape parameter. We survey numerical results for different shape parameters. Here, we present these results for two values of c. For this purpose, we apply the MQ function. In Figure 3 and Figure 4 we have shown the numerical solution for $c=hn$, $c=h{n}^2$ which using R-value estimation have RMSE $0.0169$ and $0.5389$, respectively.

Hardy [28] proposed $c=0.815d$ where $d=1/n{\sum }_{j=1}^n{d}_j$ and $d_j$ is the distance of the jth collocation point from its nearest neighbor which for equispaced node we consider $d=h$. Additionally, we found a meaningful relation between the shape parameter, the number of interpolation points, and d, which we propose in the following formulation to reach a proper shape parameter for our method. This relation is

$$c=d(\sqrt{n}+\frac{1}{i}),i=1,2,\dots$$

(7)

.

lWe gained the least $\mathrm{RMSE}=0.0013$ for the $c=d\sqrt{n}$ using R-value estimation. This nuemical solution is shown in Figure 5.

In

Table 2. Results for various shape parameters using the MQ function with n = 40.

i	$\mathbf{Condition}\mathbf{Number}$	$\mathbf{RMSE}$
1	1.1312 $\times {10}^8$	4.0 $\times {10}^{-2}$
2	1.9033 $\times {10}^7$	3.7 $\times {10}^{-3}$
3	1.0420 $\times {10}^7$	4.6 $\times {10}^{-3}$
4	7.6998 $\times {10}^6$	5.1  $\times {10}^{-3}$
5	6.4686  $\times {10}^6$	2.7  $\times {10}^{-3}$
6	5.6856 $\times {10}^6$	3.2  $\times {10}^{-3}$
7	5.1838 $\times {10}^6$	3.2 $\times {10}^{-3}$
8	4.8844 $\times {10}^6$	3.7 $\times {10}^{-3}$
9	4.6431 $\times {10}^6$	2.0 $\times {10}^{-3}$
10	4.4586 $\times {10}^6$	2.0 $\times {10}^{-3}$
$c_p=d\sqrt{n}$	3.0942 $\times {10}^6$	1.3 $\times {10}^{-3}$

and Figure 6, we show RMSE for the RBF method with shape parameters gained from relation (7) for $i=1,\cdots ,10$. These results are obtained using R-value estimation and the MQ function with n = 40.

In

Table 3. Comparison of RSMEs using two different shape parameters and six different runs with n = 70.

$\mathbf{RBF}$	$\mathbf{MQ}-\mathbf{RMSE}$	$\mathbf{IMQ}-\mathbf{RMSE}$	$\mathbf{IQ}-\mathbf{RMSE}$
	$\mathbf{c}=\mathbf{d}\sqrt{\mathbf{n}}$       $\mathbf{c}=\mathbf{0}.\mathbf{815}\mathbf{d}$	$\mathbf{c}=\mathbf{d}\sqrt{\mathbf{n}}$       $\mathbf{c}=\mathbf{0}.\mathbf{815}\mathbf{d}$	$\mathbf{c}=\mathbf{d}\sqrt{\mathbf{n}}$       $\mathbf{c}=\mathbf{0}.\mathbf{815}\mathbf{d}$
$\mathrm{run}$
1	6.0 $\times {10}^{-4}$             5.6 $\times {10}^{-3}$	7.9 $\times {10}^{-4}$             4.9 $\times {10}^{-3}$	6.0 $\times {10}^{-3}$            1.2 $\times {10}^{-3}$
2	2.1 $\times {10}^{-4}$            9.6 $\times {10}^{-3}$	5.3 $\times {10}^{-4}$             2.1 $\times {10}^{-3}$	3.2 $\times {10}^{-3}$             6.7 $\times {10}^{-4}$
3	3.4 $\times {10}^{-4}$            8.9 $\times {10}^{-3}$	8.4 $\times {10}^{-4}$            4.3 $\times {10}^{-3}$	5.9 $\times {10}^{-3}$            5.2 $\times {10}^{-4}$
4	5.7 $\times {10}^{-4}$            6.0 $\times {10}^{-3}$	7.5 $\times {10}^{-4}$             3.7 $\times {10}^{-3}$	5.1 $\times {10}^{-3}$            7.3 $\times {10}^{-4}$
5	3.7 $\times {10}^{-4}$             7.8 $\times {10}^{-3}$	4.5 $\times {10}^{-4}$             5.0 $\times {10}^{-3}$	2.7 $\times {10}^{-3}$            3.8 $\times {10}^{-3}$
6	5.5 $\times {10}^{-4}$            3.4 $\times {10}^{-3}$	7.4 $\times {10}^{-4}$             2.2 $\times {10}^{-3}$	8.6 $\times {10}^{-3}$            1.4 $\times {10}^{-3}$
Average	4.2 $\times {10}^{-4}$            6.4 $\times {10}^{-3}$	6.8 $\times {10}^{-4}$            3.7 $\times {10}^{-3}$	5.2 $\times {10}^{-3}$             1.3 $\times {10}^{-3}$

, we compare gained results using R-value estimation by the MQ, IMQ, and IQ function for the proposed shape parameter and shape parameter recommended by Hardy for six different runs with n = 70.

According to the results in Table 3, it is clear that the proposed shape parameter provides higher accuracy for the MQ and IMQ functions as the basis function. But using our suggested shape parameter for the IQ function does not yield quite reliable results. After finding the proper shape parameter and fixing it, we have surveyed results for various numbers of interpolation points. In Figure 7 and Figure 8, numerical solutions for $n=50$, and $n=70$ by utilizing the MQ function can be seen, which via R-value estimation have RMSE = 0.0079 and 0.00027, respectively. Actually, by increasing the number of interpolation points, we could reach a more accurate approximation solution.

Additionally, we consider the model contains two molecules with the same assumptions of same as the first example. The implementation of the RBF method by employing the MQ function and using R-value estimation would yield $\mathrm{RMSE}=0.0018$. Figure 9 illustrates the approximation of the numerical solution in this case.

Additionally, we surveyed the numerical results for different amounts of $\kappa \left(x\right)$. We provided results by applying the MQ, IMQ, and IQ functions.

Table 4. Comparison of $\mathrm{RSMEs}$ for various $\kappa \left(x\right)$ with $n=70$.

$\mathbf{RBF}$	$\mathbf{MQ}-\mathbf{RMSE}$	$\mathbf{IMQ}-\mathbf{RMSE}$	$\mathbf{IQ}-\mathbf{RMSE}$	$\mathbf{CSC}-\mathbf{RSME}$
$\kappa \left(x\right)$
0.01	2.7 $\times {10}^{-4}$	2.3 $\times {10}^{-4}$	2.3 $\times {10}^{-3}$	1.0 $\times {10}^{-4}$
10	9.0 $\times {10}^{-4}$	1.2 $\times {10}^{-3}$	5.3 $\times {10}^{-3}$	9.3 $\times {10}^{-5}$
50	1.0 $\times {10}^{-3}$	8.0 $\times {10}^{-4}$	1.9 $\times {10}^{-3}$	4.7 $\times {10}^{-4}$
100	1.1 $\times {10}^{-3}$	4.7 $\times {10}^{-4}$	1.4 $\times {10}^{-3}$	9.6 $\times {10}^{-5}$

and Figure 10 show RMSEs by using the R-valued estimate and SCS method with $n=70$.

As the last case for numerical analysis, we investigated some results using the MQ, IMQ, and IQ functions for the described model including one, two, and five molecular domains ($D_M$). These results are provided in

Table 5. Comparison of RMSEs for various $D_M$ with n = 70.

$\mathbf{RBF}$	$\mathbf{MQ}-\mathbf{RMSE}$	$\mathbf{IMQ}-\mathbf{RMSE}$	$\mathbf{IQ}-\mathbf{RMSE}$	$\mathbf{CSC}-\mathbf{RMSE}$
${\mathrm{D}}_{\mathrm{M}}$
1	2.7 $\times {10}^{-4}$	1.0 $\times {10}^{-3}$	8.2 $\times {10}^{-4}$	4.5 $\times {10}^{-4}$
2	7.5 $\times {10}^{-4}$	3.6 $\times {10}^{-4}$	3.7 $\times {10}^{-3}$	9.2 $\times {10}^{-5}$
5	1.3 $\times {10}^{-4}$	3.7 $\times {10}^{-4}$	4.3 $\times {10}^{-3}$	8.8 $\times {10}^{-5}$

and Figure 11 by using the R-value estimation and CSC method. Figure 12 displays the numerical solution for the model includes five molecular subdomains using the MQ function with n = 70.

The proof of stability for collocation methods is difficult. More specifically for the described model in Section 3 which is non-linear and stochastic, it is even more complicated. However, there are some ways to indicate stability numerically. One of them is as follows: if by changing the shape parameter around a determined value the results do not differ much, it can be concluded that the method is stable for that shape parameter. Due to Table 2, we can observe that whatever the shape parameter would be close to our suggested shape parameter, the condition number of the matrix obtained in (6) would change only insignificantly. Therefore, this is evidence for the stability of the RBF method when employing the suggested shape parameter.

5. Conclusions

In this paper, we have thoroughly demonstrated the efficiency and performance of the RBF method for solving non-linear stochastic elliptic models. The numerical results clearly show that this method provides accurate solutions with ease. Not only is the RBF method efficient, but it is also straightforward to implement, making it more convenient than other methods, such as finite differences or finite elements. Thus, we recommend the RBF method as a powerful tool for solving non-linear stochastic equations, particularly those with complex geometry. Looking ahead, we plan to extend and improve this method for solving biological PDEs, which are often defined on complex domains. With its accuracy and simplicity, we believe the RBF method can be a valuable resource for tackling real-world problems.