A Radial Basis Function Method with Improved Accuracy for Fourth Order Boundary Value Problems

Abstract

Accurately approximating higher order derivatives is an inherently difficult problem. It is shown that a random variable shape parameter strategy can improve the accuracy of approximating higher order derivatives with Radial Basis Function methods. The method is used to solve fourth order boundary value problems. The use and location of ghost points are examined in order to enforce the extra boundary conditions that are necessary to make a fourth-order problem well posed. The use of ghost points versus solving an overdetermined linear system via least squares is studied. For a general fourth-order boundary value problem, the recommended approach is to either use one of two novel sets of ghost centers introduced here or else to use a least squares approach. When using either ghost centers or least squares, the random variable shape parameter strategy results in significantly better accuracy than when a constant shape parameter is used.

Share and Cite:

Sarra, S. , Musgrave, D. , Stone, M. and Powell, J. (2024) A Radial Basis Function Method with Improved Accuracy for Fourth Order Boundary Value Problems. Journal of Applied Mathematics and Physics, 12, 2559-2573. doi: 10.4236/jamp.2024.127151.

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

$ℒ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,y\right)$ (1)

For the problem to be well posed, two boundary conditions

${ℬ}_{1}u={g}_{1}\left(x,y\right),\text{ }{ℬ}_{2}u={g}_{2}\left(x,y\right)$

must be applied at each boundary point where ${ℬ}_{1}$ and ${ℬ}_{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},\cdots ,{x}_{N}^{c}\right\}$ in ${ℛ}^{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

${ℐ}_{N}f\left(x\right)=\sum _{k=1}^{N}\text{ }\text{ }{a}_{k}\varphi \left({‖x-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right)$ (2)

where a is a vector of expansion coefficients. The RBF expansion coefficients are determined by enforcing the interpolation conditions

${ℐ}_{N}f\left({x}_{k}^{c}\right)=f\left({x}_{k}^{c}\right),\text{ }k=1,2,\cdots ,N$ (3)

which result in a $N×N$ linear system

$Ba=f$ (4)

to be solved for the expansion coefficients. The matrix B with entries

${b}_{jk}=\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j,k=1,\cdots ,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

$ℒ\left({ℐ}_{N}f\left(x\right)\right)=\sum _{k=1}^{N}\text{ }\text{ }{a}_{k}ℒ\varphi \left({‖x-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right)$ (7)

where $ℒ$ is a linear differential operator. The operator $ℒ$ 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}_{ℒ}$ with entries

${h}_{jk}=ℒ\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j,k=1,\cdots ,N.$ (8)

That is, $ℒf\approx {H}_{ℒ}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}_{ℒ}{B}^{-1}$ since

$ℒf\approx {H}_{ℒ}a={H}_{ℒ}\left({B}^{-1}f\right)=\left({H}_{ℒ}{B}^{-1}\right)f.$ (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.$ (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};{x}_{ℬ}^{c}\right]$ and the linear system $Ha=f$ for the expansion coefficients has the form

$\left[\begin{array}{c}ℒ\varphi \\ {ℬ}_{1}\varphi \\ {ℬ}_{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].$ (12)

The three blocks of H have elements

${\left(ℒ\varphi \right)}_{jk}=ℒ\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j=1,\cdots ,{N}_{I}\text{ }k=1,\cdots ,N$

${\left({ℬ}_{1}\varphi \right)}_{jk}={ℬ}_{1}\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j={N}_{I}+1,\cdots ,N\text{ }k=1,\cdots ,N$

${\left({ℬ}_{2}\varphi \right)}_{\left(j+{N}_{B}\right)k}={ℬ}_{2}\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j={N}_{I}+1,\cdots ,N\text{ }k=1,\cdots ,N.$

In the overdetermined system H is $M×N$ , a is $N×1$ and F is $M×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};{x}_{ℬ}^{c};{x}_{ℬ}^{c}\right]$ and the linear system $Ha=f$ for the expansion coefficients has the form

$\left[\begin{array}{c}ℒ\varphi \\ {ℬ}_{1}\varphi \\ {ℬ}_{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(ℒ\varphi \right)}_{jk}=ℒ\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j=1,\cdots ,{N}_{I}\text{ }k=1,\cdots ,M$

${\left({ℬ}_{1}\varphi \right)}_{jk}={ℬ}_{1}\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j={N}_{I}+1,\cdots ,N\text{ }k=1,\cdots ,M$

${\left({ℬ}_{2}\varphi \right)}_{jk}={ℬ}_{2}\varphi \left({‖{x}_{j}^{c}-{x}_{k}^{c}‖}_{2},{\epsilon }_{k}\right),\text{ }j=N+1,\cdots ,M\text{ }k=1,\cdots ,M$

The approximate solution is then $u=B\stackrel{^}{a}$ where the elements of the system matrix B are given by Equation (5) and the vector $\stackrel{^}{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};{x}_{ℬ}^{c};{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)×\text{rand}\left(1,N\right).$ (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[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 ill-conditioned 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[-1,1\right]$ at $N=100$ centers located at the Chebyshev-Gauss-Lobatto points

${x}_{k}=-\mathrm{cos}\left(\frac{k\pi }{N-1}\right),\text{ }k=0,1,\cdots ,N-1$ (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 Nb 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[0,1\right]$ . Both Dirichlet and Neumann boundary conditions are applied at each boundary point as follows:

$u\left(0\right)=1,\text{ }\text{\hspace{0.17em}}\text{ }{u}^{\prime }\left(0\right)=2,\text{ }\text{ }\text{\hspace{0.17em}}u\left(1\right)=2e\text{ }\text{\hspace{0.17em}}\text{ }\text{and}\text{ }\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$ non-uniformly spaced centers are distributed according to

${x}_{k}=\frac{1}{2}\frac{\mathrm{arcsin}\left(-0.99\mathrm{cos}\left(\pi k/\left(N-1\right)\right)\right)}{\mathrm{arcsin}\left(0.99\right)}+\frac{1}{2},\text{ }k=0,1,\cdots ,N-1.$ (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{ }\text{ }+32{\pi }^{4}\mathrm{cos}\left(2\pi x\right)\mathrm{cos}\left(2\pi y\right).\end{array}$

Dirichlet

$u\left(x,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,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(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 1019. 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{ }+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,y\right)=y\mathrm{sin}\left(x\right)+x\mathrm{cos}\left(y\right)$

and Laplace type boundary conditions

$\Delta u\left(x,y\right)=-x\mathrm{cos}\left(y\right)-y\mathrm{sin}\left(x\right)$

are applied. The exact solution is

$u\left(x,y\right)=y\mathrm{sin}\left(x\right)+x\mathrm{cos}\left(y\right).$ (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 semi-annular domain in Figure 5. The forcing term is $f\left(x,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 semi-circular 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,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{ }{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{ }{k}_{4}=\frac{25}{2\left(888-1250\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 well-posed: 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 fourth-order 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].

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

 [1] Yao, G., Tsai, C.H. and Chen, W. (2010) The Comparison of Three Meshless Methods Using Radial Basis Functions for Solving Fourth-Order Partial Differential Equations. Engineering Analysis with Boundary Elements, 34, 625-631. https://doi.org/10.1016/j.enganabound.2010.03.004 [2] Mardanov, R.F. and Zaripov, S.K. (2016) Solution of Stokes Flow Problem Using Biharmonic Equation Formulation and Multiquadrics Method. Lobachevskii Journal of Mathematics, 37, 268-273. https://doi.org/10.1134/s1995080216030161 [3] Liu, C. and Ku, C. (2022) A Simplified Radial Basis Function Method with Exterior Fictitious Sources for Elliptic Boundary Value Problems. Mathematics, 10, Article 1622. https://doi.org/10.3390/math10101622 [4] Haq, F. and Ali, A. (2011) Numerical Solution of Fourth Order Boundary-Value Problems Using HAAR Wavelets. Applied Mathematical Sciences, 5, 3131-3146. [5] Arifeen, S.U., Haq, S., Ghafoor, A., Ullah, A., Kumam, P. and Chaipanya, P. (2021) Numerical Solutions of Higher Order Boundary Value Problems via Wavelet Approach. Advances in Difference Equations, 2021, Article No. 347. https://doi.org/10.1186/s13662-021-03495-6 [6] B. Ogunrinde, R. and M. Ojo, O. (2018) Application of Differential Transformation Method to Boundary Value Problems of Order Seven and Eight. American Journal of Computational Mathematics, 8, 269-278. https://doi.org/10.4236/ajcm.2018.83022 [7] Al-Zaid, N., Alzahrani, K., Bakodah, H. and Al-Mazmumy, M. (2023) Efficient Decomposition Shooting Method for Solving Third-Order Boundary Value Problems. International Journal of Modern Nonlinear Theory and Application, 12, 81-98. https://doi.org/10.4236/ijmnta.2023.123006 [8] Buhmann, M.D. (2003). Radial Basis Functions. Cambridge University Press. https://doi.org/10.1017/cbo9780511543241 [9] Schaback, R. (1995) Error Estimates and Condition Numbers for Radial Basis Function Interpolation. Advances in Computational Mathematics, 3, 251-264. https://doi.org/10.1007/bf02432002 [10] Sarra, S.A. (2014) Regularized Symmetric Positive Definite Matrix Factorizations for Linear Systems Arising from RBF Interpolation and Differentiation. Engineering Analysis with Boundary Elements, 44, 76-86. https://doi.org/10.1016/j.enganabound.2014.04.019 [11] Sarra, S.A. (2023) Local Radial Basis Function Methods: Comparison, Improvements, and Implementation. Journal of Applied Mathematics and Physics, 11, 3867-3886. https://doi.org/10.4236/jamp.2023.1112245 [12] Kansa, E.J. (1990) Multiquadrics—A Scattered Data Approximation Scheme with Applications to Computational Fluid-Dynamics—II Solutions to Parabolic, Hyperbolic and Elliptic Partial Differential Equations. Computers & Mathematics with Applications, 19, 147-161. https://doi.org/10.1016/0898-1221(90)90271-k [13] Hon, Y.C. and Schaback, R. (2001) On Unsymmetric Collocation by Radial Basis Functions. Applied Mathematics and Computation, 119, 177-186. https://doi.org/10.1016/s0096-3003(99)00255-6 [14] Fasshauer, G.E. (2007) Meshfree Approximation Methods with Matlab. World Scientific. https://doi.org/10.1142/6437 [15] Sarra, S.A. and Kansa, E.J. (2009) Multiquadric Radial Basis Function Approximation Methods for the Numerical Solution of Partial Differential Equations. Advances in Computational Mechanics, 2, 1-206. [16] Wendland, H. (2004). Scattered Data Approximation. Cambridge University Press. https://doi.org/10.1017/cbo9780511617539 [17] Sarra, S.A. and Sturgill, D. (2009) A Random Variable Shape Parameter Strategy for Radial Basis Function Approximation Methods. Engineering Analysis with Boundary Elements, 33, 1239-1245. https://doi.org/10.1016/j.enganabound.2009.07.003 [18] Hesthaven, J.S., Gottlieb, S. and Gottlieb, D. (2007) Spectral Methods for Time-Dependent Problems. Cambridge University Press. https://doi.org/10.1017/cbo9780511618352 [19] Welfert, B.D. (1997) Generation of Pseudospectral Differentiation Matrices I. SIAM Journal on Numerical Analysis, 34, 1640-1657. https://doi.org/10.1137/s0036142993295545 [20] Fornberg, B. (2006) A Pseudospectral Fictitious Point Method for High Order Initial-Boundary Value Problems. SIAM Journal on Scientific Computing, 28, 1716-1729. https://doi.org/10.1137/040611252 [21] Fedoseyev, A.I., Friedman, M.J. and Kansa, E.J. (2002) Improved Multiquadric Method for Elliptic Partial Differential Equations via PDE Collocation on the Boundary. Computers & Mathematics with Applications, 43, 439-455. https://doi.org/10.1016/s0898-1221(01)00297-8 [22] Kuwabara, S. (1959) The Forces Experienced by Randomly Distributed Parallel Circular Cylinders or Spheres in a Viscous Flow at Small Reynolds Numbers. Journal of the Physical Society of Japan, 14, 527-532. https://doi.org/10.1143/jpsj.14.527 [23] Golbabai, A. and Rabiei, H. (2012) A Meshfree Method Based on Radial Basis Functions for the Eigenvalues of Transient Stokes Equations. Engineering Analysis with Boundary Elements, 36, 1555-1559. https://doi.org/10.1016/j.enganabound.2012.04.001 [24] Hines, T. (2016) Python Package Containing Tools for Radial Basis Function (RBF) Applications. https://github.com/treverhines/RBF [25] Marshall, H., Sahraoui, M. and Kaviany, M. (1994) An Improved Analytic Solution for Analysis of Particle Trajectories in Fibrous, Two-Dimensional Filters. Physics of Fluids, 6, 507-520. https://doi.org/10.1063/1.868346 [26] Sarra, S.A. (2016) The Radial Basis Function Toolkit. http://www.scottsarra.org/rbf/rbf.html