A Radial Basis Function Method with Improved Accuracy for Fourth Order Boundary Value Problems ()
1. Introduction
This work examines Radial Basis Function (RBF) methods for the solution of fourth order boundary value problems (BVPs) involving the linear Biharmonic equation
$\mathcal{L}u={\Delta}^{2}u=\frac{{\partial}^{4}u}{\partial {x}^{4}}+2\frac{{\partial}^{4}u}{\partial {x}^{2}\partial {y}^{2}}+\frac{{\partial}^{4}u}{\partial {y}^{4}}=f\left(x\mathrm{,}y\right)$
(1)
For the problem to be well posed, two boundary conditions
${\mathcal{B}}_{1}u={g}_{1}\left(x\mathrm{,}y\right)\mathrm{,}\text{\hspace{1em}}{\mathcal{B}}_{2}u={g}_{2}\left(x\mathrm{,}y\right)$
must be applied at each boundary point where ${\mathcal{B}}_{1}$
and ${\mathcal{B}}_{2}$
may be Dirichlet, Neumann, Robin, or Laplace type boundary conditions. The Biharmonic equation is used in numerous applications such as the flex of elastic plates, the flow of particulate suspensions, the flow of molten metals and in the modeling of fluid flow.
Two inherent difficulties exist for fourth order BVPs. First, it is more difficult to accurately approximate fourth order derivatives than it is to approximate lower order derivatives. A variable shape parameter is used to alleviate this issue. Second is the need to apply two boundary conditions at each boundary point. Two approaches have been used previously to enforce the dual boundary conditions. The first is to solve an overdetermined linear system via least squares (LSQ) [1]. The second is to add ghost points (also called fictitious points or ghost centers) in order to apply the second boundary condition [2] [3]. This leads to a square linear system to solve, but there is no theory that dictates how to locate the ghost points. The best locations seem to be problem dependent.
In addition to the RBF method with ghost points or least squares, other approaches to solving higher order BVPs have been used. Some of the approaches include: wavelet methods in references [4] and [5], a differential transformation method in [6], and a decomposition shooting method in [7].
2. RBF Methods for BVPs
RBF methods for BVPs are based on differentiating a RBF interpolant. RBF interpolation uses a set of N distinct points $X=\left\{{x}_{1}^{c}\mathrm{,}\cdots \mathrm{,}{x}_{N}^{c}\right\}$
in ${\mathcal{R}}^{d}$
called centers. No restrictions are placed on the shape of problem domains or on the location of the centers. The RBF interpolant has the form
${\mathcal{I}}_{N}f\left(x\right)={\displaystyle \sum _{k=1}^{N}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{a}_{k}\varphi \left({\Vert x{x}_{k}^{c}\Vert}_{2}\mathrm{,}{\epsilon}_{k}\right)$
(2)
where a is a vector of expansion coefficients. The RBF expansion coefficients are determined by enforcing the interpolation conditions
${\mathcal{I}}_{N}f\left({x}_{k}^{c}\right)=f\left({x}_{k}^{c}\right)\mathrm{,}\text{\hspace{1em}}k=\mathrm{1,2,}\cdots \mathrm{,}N$
(3)
which result in a $N\times N$
linear system
$Ba=f$
(4)
to be solved for the expansion coefficients. The matrix B with entries
${b}_{jk}=\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2}\mathrm{,}{\epsilon}_{k}\right)\mathrm{,}\text{\hspace{1em}}j\mathrm{,}k=\mathrm{1,}\cdots \mathrm{,}N$
(5)
is called the system matrix. The shape parameter ${\epsilon}_{k}$
may be the same at each center, a constant shape parameter, or may be different at each center, a variable shape parameter. Throughout, the inverse multiquadric (IMQ) RBF
$\varphi \left(r\right)=\frac{1}{\sqrt{1+{\epsilon}^{2}{r}^{2}}}$
(6)
is used in all examples. The IMQ is a representative member of the class of global, infinitely differentiable RBFs containing a shape parameter and that interpolate with exponential accuracy under suitable circumstances [8]. With a constant shape parameter, the IMQ system matrix is symmetric positive definite (SPD) and therefore is invertible.
It is well established that the RBF method is most accurate when the system matrix is poorly conditioned [9]. Typically it is most accurate when the condition number of the system matrix B is $\mathcal{O}\left({10}^{16}\right)$
. In some cases, the theoretically SPD system matrix may not be numerically SPD and a Cholesky factorization will fail. However, the linear system can be effectively regularized by the method of diagonal increments (MDI) [10] [11] so that the matrix is numerically SPD and also potentially restore several decimal places of accuracy when compared to the unregularized system.
Derivatives are approximated by differentiating the RBF interpolant as
$\mathcal{L}\left({\mathcal{I}}_{N}f\left(x\right)\right)={\displaystyle \sum _{k=1}^{N}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{a}_{k}\mathcal{L}\varphi \left({\Vert x{x}_{k}^{c}\Vert}_{2}\mathrm{,}{\epsilon}_{k}\right)$
(7)
where $\mathcal{L}$
is a linear differential operator. The operator $\mathcal{L}$
may be a single differential operator or a linear differential operator such as the Laplacian or Biharmonic operator. Evaluating (7) at the centers X can be accomplished by multiplying the expansion coefficients by the evaluation matrix ${H}_{\mathcal{L}}$
with entries
${h}_{jk}=\mathcal{L}\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2},{\epsilon}_{k}\right),\text{\hspace{1em}}j,k=1,\cdots ,N.$
(8)
That is, $\mathcal{L}f\approx {H}_{\mathcal{L}}a$
. Alternatively, derivatives can be approximated by multiplying the function values at the center locations, ${\left\{f\left({x}_{k}^{c}\right)\right\}}_{k=1}^{N}$
, by the differentiation matrix $D={H}_{\mathcal{L}}{B}^{1}$
since
$\mathcal{L}f\approx {H}_{\mathcal{L}}a={H}_{\mathcal{L}}\left({B}^{1}f\right)=\left({H}_{\mathcal{L}}{B}^{1}\right)f\mathrm{.}$
(9)
The solution of fourth order BVPs is approximated by solving a linear system
$Ha=F$
(10)
using LU factorization for the expansion coefficients of the solution and then the solution is obtain by multiplying the expansion coefficients by the system matrix
$u=Ba\mathrm{.}$
(11)
The form and shape of H depends on how the boundary conditions are enforced.
Let ${N}_{I}$
be the number of interior points, ${N}_{B}$
the number of boundary points, ${N}_{G}$
the number of ghost centers, $N={N}_{I}+{N}_{B}$
the number of centers, and $M=N+{N}_{G}$
the number of centers and ghost points. Additionally, let ${x}_{I}^{c}$
be interior centers, ${x}_{B}^{c}$
be boundary centers, and ${x}_{G}^{c}$
be ghost centers.
In the least squares approach, centers are in an array that is ordered as $X=\left[{x}_{I}^{c}\mathrm{;}{x}_{\mathcal{B}}^{c}\right]$
and the linear system $Ha=f$
for the expansion coefficients has the form
$\left[\begin{array}{c}\mathcal{L}\varphi \\ {\mathcal{B}}_{1}\varphi \\ {\mathcal{B}}_{2}\varphi \end{array}\right]a=\left[\begin{array}{c}f\left({x}_{I}^{c}\right)\\ {g}_{1}\left({x}_{B}^{c}\right)\\ {g}_{2}\left({x}_{B}^{c}\right)\end{array}\right]\mathrm{.}$
(12)
The three blocks of H have elements
${\left(\mathcal{L}\varphi \right)}_{jk}=\mathcal{L}\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2},{\epsilon}_{k}\right),\text{\hspace{1em}}j=1,\cdots ,{N}_{I}\text{\hspace{1em}}k=1,\cdots ,N$
${\left({\mathcal{B}}_{1}\varphi \right)}_{jk}={\mathcal{B}}_{1}\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2},{\epsilon}_{k}\right),\text{\hspace{1em}}j={N}_{I}+1,\cdots ,N\text{\hspace{1em}}k=1,\cdots ,N$
${\left({\mathcal{B}}_{2}\varphi \right)}_{\left(j+{N}_{B}\right)k}={\mathcal{B}}_{2}\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2},{\epsilon}_{k}\right),\text{\hspace{1em}}j={N}_{I}+1,\cdots ,N\text{\hspace{1em}}k=1,\cdots ,N.$
In the overdetermined system H is $M\times N$
, a is $N\times 1$
and F is $M\times 1$
. The approximate solution is then $u=Ba$
where the elements of the system matrix B are given by Equation (5).
If two boundary conditions are applied at each boundary point, centers are in an array that is ordered as $X=\left[{x}_{I}^{c}\mathrm{;}{x}_{\mathcal{B}}^{c}\mathrm{;}{x}_{\mathcal{B}}^{c}\right]$
and the linear system $Ha=f$
for the expansion coefficients has the form
$\left[\begin{array}{c}\mathcal{L}\varphi \\ {\mathcal{B}}_{1}\varphi \\ {\mathcal{B}}_{2}\varphi \end{array}\right]a=\left[\begin{array}{c}f\left({x}_{I}^{c}\right)\\ {g}_{1}\left({x}_{B}^{c}\right)\\ {g}_{2}\left({x}_{B}^{c}\right)\end{array}\right]$
(13)
The three blocks of H have elements
${\left(\mathcal{L}\varphi \right)}_{jk}=\mathcal{L}\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2},{\epsilon}_{k}\right),\text{\hspace{1em}}j=1,\cdots ,{N}_{I}\text{\hspace{1em}}k=1,\cdots ,M$
${\left({\mathcal{B}}_{1}\varphi \right)}_{jk}={\mathcal{B}}_{1}\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2},{\epsilon}_{k}\right),\text{\hspace{1em}}j={N}_{I}+1,\cdots ,N\text{\hspace{1em}}k=1,\cdots ,M$
${\left({\mathcal{B}}_{2}\varphi \right)}_{jk}={\mathcal{B}}_{2}\varphi \left({\Vert {x}_{j}^{c}{x}_{k}^{c}\Vert}_{2},{\epsilon}_{k}\right),\text{\hspace{1em}}j=N+1,\cdots ,M\text{\hspace{1em}}k=1,\cdots ,M$
The approximate solution is then $u=B\widehat{a}$
where the elements of the system matrix B are given by Equation (5) and the vector $\widehat{a}$
contains the first N elements of a.
When ghost centers are used, the linear system still has the structure of Equation (13), but instead of the boundary points being repeated, ghost centers are added to the end of the array of centers, that is $X=\left[{x}_{I}^{c}\mathrm{;}{x}_{\mathcal{B}}^{c}\mathrm{;}{x}_{\mathcal{G}}^{c}\right]$
. This changes the last ${N}_{G}$
elements of the right vector in (13) from ${g}_{2}\left({x}_{B}^{c}\right)$
to ${g}_{2}\left({x}_{G}^{c}\right)$
.
The described RBF method for boundary value problems is known as the asymmetric RBF collocation method due to the fact that the matrix H is not symmetric. The method is also known as Kansa’s method [12] in reference to the researcher who first described the method.
Even when employing a constant shape parameter, the evaluation matrix H cannot be shown to always be invertible. In fact, examples have been constructed in which the evaluation matrix is singular [13]. Despite the lack of a firm theoretical underpinning, extensive computational evidence indicates that the matrix H is very rarely singular and the asymmetric method has become well established for steady problems.
Kansa’s method is well documented in the research monographs [8] [14][16] which can be consulted for more information on the method and for more examples of application of Kansa’s method.
3. Variable Shape Parameter
In the majority of implementations of RBF methods a constant shape parameter ${\epsilon}_{k}=\epsilon $
is used. Another possibility is to use a different value of the shape parameter ${\epsilon}_{k}$
at each center k which is called a variable shape parameter strategy. Variable shape strategies result in shape parameters that are different in each column of the system matrix (5) and the evaluation matrix (8). Figure 1 plots several IMQ basis functions with a constant and variable shape. In the left image, the basis functions are all horizontal translations of the same basis function when a constant shape is used. In the right image, a variable shape results in the basis functions having different shapes as is the case in polynomial based methods. One argument for using a variable shape parameter is that it leads to more distinct entries in the RBF matrices which in turn lead to lower condition numbers. A negative consequence of using a variable shape is that the system matrix is no longer symmetric and the standard results of the invertibility of the system matrix for a constant shape parameter no longer apply. However, with Kansa’s method for BVPs this is not an issue as stated in Section 2, the matrix H (11) cannot be shown to always be invertible with a constant shape.
Figure 1. Left: basis functions with the same shape parameter. Right: basis functions with different shape parameters.
Several variable shape parameter strategies have been proposed and can be found in reference [15]. The particular variable shape parameter strategy employed in this work is the random variable shape that is described in [17]. In this strategy, the shape parameter is specified as
${\epsilon}_{k}={\epsilon}_{\mathrm{min}}+\left({\epsilon}_{\mathrm{max}}{\epsilon}_{\mathrm{min}}\right)\times \text{rand}\left(\mathrm{1,}N\right)\mathrm{.}$
(14)
Equation (14) returns N random shape parameters between ${\epsilon}_{\mathrm{min}}$
and ${\epsilon}_{\mathrm{max}}$
. The function rand returns random samples from a uniform distribution over $\left[\mathrm{0,1}\right)$
. Each time the function is called, different random numbers are returned. However, the results can be made reproducible if the rand function is given the same seed before each function call as is illustrated in the follow Python code fragment:
In the numerical results that follow, the random variable shape implementation is compared with constant shape implementations using the average shape parameter
${\epsilon}_{\text{avg}}=\frac{1}{2}\left({\epsilon}_{\mathrm{max}}+{\epsilon}_{\mathrm{min}}\right)$
(15)
as well as ${\epsilon}_{\mathrm{min}}$
and ${\epsilon}_{\mathrm{max}}$
.
In subsequent sections, the random variable shape parameter and constant shape parameter are compared in approximating higher order derivatives and the solution of fourth order BVPs. Additional comparisons with interpolation problems, lower order derivatives, and second order BVPs can be found in reference [17].
4. Approximating Fourth Order Derivatives
Polynomial based pseudospectral methods, such as the Chebyshev pseudospectral (CPS) method, have closed form formulas for the elements of their differentiation matrices (DMs) [18] [19]. Despite being able to accurately specify the elements of its DMs, the CPS method is less accurate for higher order derivatives than it is for lower order derivatives. Formulas for the RBF DM elements do not exist. Instead an illconditioned linear system must be solved to form DMs which makes the problem of accurately calculating higher order derivatives even worse. The random variable shape parameter strategy of the previous section can be used to somewhat alleviate this problem.
To illustrate, the first four derivatives of the infinitely differentiable function
$f\left(x\right)={\text{e}}^{\mathrm{sin}\left(\pi x\right)}$
(16)
are approximated on the interval $\left[\mathrm{1,1}\right]$
at $N=100$
centers located at the ChebyshevGaussLobatto points
${x}_{k}=\mathrm{cos}\left(\frac{k\pi}{N1}\right),\text{\hspace{1em}}k=0,1,\cdots ,N1$
(17)
that are required in the implementation of the CPS method. The results are in Table 1. The RBF method is implemented with a random variable shape ranging between ${\epsilon}_{\mathrm{min}}=2.4$
and ${\epsilon}_{\mathrm{max}}=4$
with ${\epsilon}_{\text{avg}}=3.2$
and with three constant shapes that correspond to ${\epsilon}_{\mathrm{min}}$
, ${\epsilon}_{\mathrm{max}}$
, and ${\epsilon}_{\text{avg}}$
. Moving from the first to fourth derivative, the CPS method loses approximately three accurate decimal places with each increase in derivative order. All three constant shape RBF methods are less than half as accurate as the CPS method and do not have any decimal places of accuracy in approximating the fourth derivative. The RBF method with a random variable shape parameter has approximately four more decimal places of accuracy than any of the three constant shape methods. All four of the RBF methods have system matrices with condition numbers that are either $\mathcal{O}\left({10}^{15}\right)$
or $\mathcal{O}\left({10}^{16}\right)$
.
Table 1. Order $\rho $
derivative accuracy.
$\rho $

$\epsilon =2.4$

$\epsilon =3.2$

$\epsilon =4$

$2.4\le \epsilon \le 4$

CPS 
1 
7.0e−6 
7.6e−6 
8.2e−6 
5.2e−10 
1.5e−12 
2 
3.1e−3 
4.5e−3 
6.0e−3 
4.3e−7 
5.8e−9 
3 
6.6e−1 
1.3e+0 
2.1e+0 
1.2e−4 
3.8e−6 
4 
8.5e+1 
2.2e+2 
4.5e+2 
1.7e−2 
5.1e−3 
5. Ghost Centers
Ghost centers have been used in a variety of different methods for the numerical solution of differential equations. The points have been used for several reasons. The reasons include eliminating spurious eigenvalues and time stepping difficulties in the pseudospectral method solution of high order time dependent PDEs [20], for applying boundary conditions for fourth order BVPs in the RBF method, and for improving the accuracy of the RBF method for the solution of second order BVPs by adding ghost centers in order to collocate both the PDE and apply boundary conditions at the boundary points [21]. This work focuses solely on the use of ghost centers for applying boundary conditions in fourth order BVPs.
There is no theoretical guidance available that specifies how to locate ghost points in order to obtain the best result. Rather, their use seems more art than science. A fundamental difference exists in using ghost points in polynomial based methods, such as pseudospectral and finite difference methods, and RBF methods since polynomial based methods depend on the locations of centers and RBF methods do not depend on center locations, but only the distances between centers. A large number of strategies to locate the points have been considered in the literature and it is not possible to evaluate each strategy. The following four representative strategies of how to add N_{b} ghost centers have been considered. Figure 2 illustrates the four strategies for a domain that is a unit square.
• Strategy G1 reuses the boundary centers as ghost centers. The ghost centers coincide with the boundary centers (upper left image of Figure 2).
• Strategy G2 uses ghost centers with locations that are small random perturbations of the boundary centers. With strategy G2 the ghost centers may be inside or outside of the domain (upper right image of Figure 2).
• Strategy G3 places ghost centers in an outline of the domain boundary that is the same shape as the boundary. The outlining ghost centers may be close or far from the domain boundary (lower left image of Figure 2). References [1] [2] [22] and [23] have used this strategy.
Figure 2. Interior centers (black), boundary centers (blue), and ghost centers (red). Upper left: strategy (G1). Upper right: strategy (G2). Lower left: strategy (G3). Lower right: strategy (G4).
• Strategy G4 places ghost centers in a circle outside of the domain (lower right image of Figure 2). Reference [3] uses a similar strategy, but locates centers randomly outside a domain and not necessarily in a circle.
While strategies G3 and G4 have been used elsewhere in the literature, we are unaware of strategies G1 and G2 previous being used.
6. Numerical Examples
A fourth order ODE BVP and three problems containing the Biharmonic operator in two space dimensions are considered. In the two dimensional problems, scattered centers are located according to a minimal energy algorithm implemented in the software package from reference [24]. The centers are scattered but cover domains in a fairly uniform manner. Each problem is evaluated using each of the four ghost center strategies and also by using a least squares approach which does not add ghost centers but instead solves an overdetermined system.
6.1. 1d BVP
The first example considered is the ODE BVP
${u}^{\left(4\right)}\left(x\right)u\left(x\right)=4{\text{e}}^{x}$
(18)
on the interval $\Omega =\left[\mathrm{0,1}\right]$
. Both Dirichlet and Neumann boundary conditions are applied at each boundary point as follows:
$u\left(0\right)=1,\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{u}^{\prime}\left(0\right)=2,\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}u\left(1\right)=2e\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{and}\text{\hspace{0.05em}}\text{\hspace{0.17em}}{u}^{\prime}\left(1\right)=3\text{e}\text{.}$
The exact solution is $u\left(x\right)=\left(1+x\right){\text{e}}^{x}$
. On the interior and boundary, $N=50$
nonuniformly spaced centers are distributed according to
${x}_{k}=\frac{1}{2}\frac{\mathrm{arcsin}\left(0.99\mathrm{cos}\left(\pi k/\left(N1\right)\right)\right)}{\mathrm{arcsin}\left(0.99\right)}+\frac{1}{2},\text{\hspace{1em}}k=0,1,\cdots ,N1.$
(19)
The centers are more densely clustered near the boundary points than in the interior of the domain. Ghost centers are located outside the domain a distance $\delta $
from each endpoint.
The numerical results are in Table 2. Taking $\delta =0$
corresponds to ghost point strategy G1 and with a variable shape parameter produces the most accurate results with the solution having eight accurate decimal places. However, with a constant shape, the results are highly inaccurate. The reason why will be explored in the next example. With $\delta =0.01$
(strategy G2), the constant and variable shape approaches all result in about one accurate decimal place in the solution. The second most accurate approach is the LSQ method with a variable shape parameter which has about five more accurate decimal places than the LSQ method with a constant shape.
Table 2. Max errors from the solution of problem (18).

$\epsilon =2.2$

$\epsilon =2.4$

$\epsilon =2.6$

$2.2\le \epsilon \le 2.6$

$\delta =0$

4.0e+9 
1.1e+7 
5.0e+6 
1.32e−9 
$\delta =0.01$

1.9e−2 
1.9e−2 
1.8e−2 
1.8e−2 
$\delta =1$

2.7 
3.3 
3.8 
2.4 
LSQ 
2.3e−3 
6.8e−3 
4.9e−3 
5.2e−8 
6.2. Linear Biharmonic Equation
Equation (1) is considered on a domain which is the unit square and with the forcing term
$\begin{array}{c}f\left(x,y\right)=16{\pi}^{4}\left(\mathrm{cos}\left(2\pi x\right)1\right)\mathrm{cos}\left(2\pi y\right)+16{\pi}^{4}\left(\mathrm{cos}\left(2\pi y\right)1\right)\mathrm{cos}\left(2\pi x\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}+32{\pi}^{4}\mathrm{cos}\left(2\pi x\right)\mathrm{cos}\left(2\pi y\right).\end{array}$
Dirichlet
$u\left(x\mathrm{,}y\right)=\left(1\mathrm{cos}\left(2\pi x\right)\right)\left(1\mathrm{cos}\left(2\pi y\right)\right)$
and Neumann conditions
$\frac{\partial u}{\partial n}=0$
are enforced. The exact solution is
$u\left(x\mathrm{,}y\right)=\left(1\mathrm{cos}\left(2\pi x\right)\right)\left(1\mathrm{cos}\left(2\pi y\right)\right).$
Example center and ghost center locations for the problem are in Figure 2. The numerical results were obtained with ${N}_{I}=541$
, ${N}_{b}=84$
and $M=709$
. Ghost center strategy G4 used a circle with radius $R=0.4$
and center $\left(\mathrm{1.5,1.5}\right)$
. The numerical results are in Table 3.
Table 3. Max errors from the solution of the linear Biharmonic equation on the unit square.

$\epsilon =1.8$

$\epsilon =2.0$

$\epsilon =2.2$

$1.8\le \epsilon \le 2.2$

G1 
1.2e+3 
2.9e+4 
3.4e+2 
3.4e−6 
G2, $\delta =0.01$

3.4e−4 
5.5e−5 
4.1e−5 
1.9e−6 
G3, $\delta =0.1$

1.3e+3 
4.1e+3 
1.9e+4 
1.4e−6 
G3, $\delta =0.01$

3.6e+6 
1.8e+7 
6.7e+5 
3.0e−5 
G4 
2.2e−4 
4.9e−3 
2.1e−2 
1.8e−3 
LSQ 
2.5e−3 
2.3e−3 
2.6e−3 
2.2e−4 
The least accurate constant shape approaches are with G1 and G3 with both $\delta =0.1$
and $\delta =0.01$
. In the right image of Figure 3, for a large range of shape parameter, the condition number of the matrix H is in excess of 10^{19}. Such large condition numbers are difficult to calculate accurately in double precision floating point arithmetic and indicate that the matrix is essentially numerically singular. Note how much better the matrices are conditioned with a variable shape in the left image of the figure. This improvement in conditioning results in the significant increase in accuracy of the three approaches when a variable shape is used. The most accurate results are with a variable shape parameter and with ghost center strategies G1, G2 with $\delta =0.01$
, and G3 with $\delta =0.1$
.
Figure 3. Condition number of H with the four ghost point strategies. Left: variable shape. Right: constant shape.
6.3. A Complexly Shaped Domain
This example considers a Biharmonic type equation
${\Delta}^{2}u+xyu+2y\mathrm{sin}\left(x\right)\frac{\partial u}{\partial x}y\mathrm{cos}\left(x\right)\frac{\partial u}{\partial y}=f\left(x,y\right)$
(20)
where
$\begin{array}{c}f\left(x,y\right)=xy\left(x\mathrm{cos}\left(y\right)+y\mathrm{sin}\left(x\right)\right)+x\mathrm{cos}\left(y\right)+y\left(x\mathrm{sin}\left(y\right)\mathrm{sin}\left(x\right)\right)\mathrm{cos}\left(x\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}+2y\left(y\mathrm{cos}\left(x\right)+\mathrm{cos}\left(y\right)\right)\mathrm{sin}\left(x\right)+y\mathrm{sin}\left(x\right)\end{array}$
(21)
in a complexly shaped domain. The domain is bounded by the curve
$\partial \Omega =\left\{\left(x,y\right)x=\rho \mathrm{cos}\theta \text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}y=\rho \mathrm{sin}\theta \right\}$
where $0\le \theta \le 2\pi $
and $\rho =\sqrt{\mathrm{cos}\left(2\theta \right)+\sqrt{1.1{\mathrm{sin}}^{2}\left(2\theta \right)}}$
(Figure 4). Both Dirichlet boundary conditions
Figure 4. Left: domain shape and scattered centers for problem (20). Right: pointwise errors from variable shape parameter and ghost center G1 approach.
$u\left(x\mathrm{,}y\right)=y\mathrm{sin}\left(x\right)+x\mathrm{cos}\left(y\right)$
and Laplace type boundary conditions
$\Delta u\left(x\mathrm{,}y\right)=x\mathrm{cos}\left(y\right)y\mathrm{sin}\left(x\right)$
are applied. The exact solution is
$u\left(x\mathrm{,}y\right)=y\mathrm{sin}\left(x\right)+x\mathrm{cos}\left(y\right)\mathrm{.}$
(22)
The numerical results in Table 4 were obtained with ${N}_{I}=343$
, ${N}_{b}=77$
and $M=497$
. This is a relatively sparse discretization of the domain with such a small number of centers. Ghost point strategy G4 was implemented with a circle of radius $R=0.5$
and center (5, 5). Strategy G3 was easily implemented for the simply shaped square domain in the previous example but it is significantly more difficult to effectively implement G3 with a more complexly shaped domain and thus its results are omitted.
Table 4. Max errors for problem (20).

$\epsilon =1.5$

$\epsilon =1.6$

$\epsilon =1.7$

$1.5\le \epsilon \le 1.7$

G1 
2.4e+4 
1.5e+4 
1.4e+3 
5.3e−4 
G2, $\delta =0.001$

4.5e−4 
4.9e−4 
8.6e−4 
2.2e−3 
G4 
1.6e−1 
1.1 
9.8e−2 
3.5e−1 
LSQ 
1.5e−1 
2.4e−1 
3.2e−1 
3.7e−3 
With a variable shape parameter, method G1 (right image of Figure 4) has three accurate decimal places followed by G2 and LSQ with 2 accurate decimal places. With a constant shape, G2 with very small perturbations of the boundary points with $\delta =0.001$
was most accurate. In addition to the results in Table 4, similar or even better accuracy, can be obtained with other variable shape parameter ranges. For example, if the LSQ approach is taken with the same center locations that produced the results in Table 4 and with the shape parameter range $0.9\le \epsilon \le 1.9$
, the maximum error is 6.3e−4.
6.4. Steady Incompressible Fluid Flow past a Cylinder
This example considers the Biharmonic Equation (1) on the semiannular domain in Figure 5. The forcing term is $f\left(x\mathrm{,}y\right)=0$
. On the two linear portions of the boundary from (−5, 0) to (−1, 0) and from (1, 0) to (5, 0), a Dirichlet BC, $u=0$
, and a BC $\frac{{\partial}^{2}u}{\partial {y}^{2}}=0$
are enforced. On the outer semicircular portion of the boundary a tangential BC, $\frac{\partial u}{\partial s}=0$
, and a Laplace type BC, $\Delta u=0$
, are specified. On the inner semicircular portion of, the boundary a Dirichlet BC, $u=0$
, and a Neumann BC, $\frac{\partial u}{\partial n}=0$
, are applied. Reference [22] derives the exact solution for the problem which is
$u\left(x\mathrm{,}y\right)=\frac{1}{5}\left(\frac{{k}_{1}y}{{x}^{2}+{y}^{2}}+{k}_{2}y+{k}_{3}y\mathrm{ln}\left(\sqrt{{x}^{2}+{y}^{2}}\right)+{k}_{4}y\left({x}^{2}+{y}^{2}\right)\right)$
where
${k}_{1}=\frac{1225}{1776+2500\mathrm{ln}\left(5\right)},\text{\hspace{1em}}{k}_{2}=\frac{12}{\frac{444}{25}25\mathrm{ln}\left(5\right)}$
and
${k}_{3}=\frac{1}{\frac{444}{625}+\mathrm{ln}\left(5\right)},\text{\hspace{1em}}{k}_{4}=\frac{25}{2\left(8881250\mathrm{ln}\left(5\right)\right)}.$
Figure 5. Left: Domain and scattered centers for fluid flow past a cylinder problem. Right: pointwise errors from variable shape parameter and ghost center G2 approach.
The PDE models creeping fluid flow through a periodic arrangement of cylinders [25].
The numerical results in Table 5 in the relatively large domain were obtained with ${N}_{I}=2325$
, ${N}_{b}=175$
and $M=2675$
. Ghost point strategy G4 is implemented with a circle of radius $R=1$
that is centered at (20, 20). In this example the variable shape LSQ approach is most accurate followed by method G2 with $\delta =0.05$
with a variable shape parameter (right image of Figure 5).
Table 5. Max errors from the cylinder flow problem.

$\epsilon =0.5$

$\epsilon =1.0$

$\epsilon =1.5$

$0.5\le \epsilon \le 1.5$

G1 
7.9e+14 
4.3e+14 
1.2e+15 
1.6e−2 
G2, $\delta =0.05$

3.2e−2 
2.2e−2 
1.1 
5.9e−4 
G4 
2.2e−1 
3.1e−2 
9.5e−2 
1.5e−2 
LSQ 
9.1e−4 
1.9e−2 
3.3e−1 
1.1e−4 
As in the previous examples, the shape parameters used many not be optimal for all approaches. For example, if $0.5\le \epsilon \le 1.4$
is used the LSQ error is 8.2e−5.
Reference [2] solves the same problem using the RBF method with the multiquadric RBF [15] using a ghost center strategy somewhat similar to the G3 strategy described here excepted that the ghost centers are not the same distance from the domain boundary around the entire boundary. The best reported accuracy that is realized in reference [2] is about one accurate decimal place. In reference [2], the location of the ghost points are displayed in an image but a formula for their exact location is not given. Thus, we were not able to fairly compare to their results.
7. Conclusions
The numerical solution of fourth order boundary value problems is challenging due to the inherent difficultly in accurately approximating fourth order derivatives and due to the necessity of applying two boundary conditions at each boundary point. It has been demonstrated that a random variable shape parameter strategy can allow the Radial Basis Function method to realize several more decimal places of accuracy in approximating any order derivative when compared to the more commonly used constant shape parameter approach.
Two approaches have been used to apply the dual boundary conditions that are necessary to make a fourth order BVP wellposed: a least squares solution of an overdetermined linear system and the use of ghost centers. Ghost centers have been successfully used by multiple researchers but a theoretical underpinning that states where the centers should be located does not exist. The optimal locations seem largely problem dependent.
In all numerical examples, the random variable shape parameter strategy performed better than using a constant shape parameter regardless of whether least squares or a ghost center methodology was used. The variable shape least squares approach was not the most accurate on the majority of examples, but it produced accurate results on each of the problems. Instead of the location of ghost points being problem and boundary condition dependent, strategies G1 and G2 with a variable random shape parameter provided good results in all example problems and were easily implemented for any shaped domain.
The least squares approach and ghost center strategies G1 and G2 with a random variable shape parameter are the recommended approaches. Each can easily be implemented for a fourthorder boundary value problem regardless of the type of boundary conditions or the shape of the domain. All results in this manuscript are reproducible via scripts in the folder/papers/JAMP2024/ in version 2.1 of the Python Radial Basis Function Toolbox [26].