A Multiscale RBF Collocation Method for the Numerical Solution of Partial Differential Equations

Abstract

In this paper, we derive and discuss the hierarchical radial basis functions method for the approximation to Sobolev functions and the collocation to well-posed linear partial differential equations. Similar to multilevel splitting of finite element spaces, the hierarchical radial basis functions are constructed by employing successive refinement scattered data sets and scaled compactly supported radial basis functions with varying support radii. Compared with the compactly supported radial basis functions approximation and stationary multilevel approximation, the new method can not only solve the present problem on a single level with higher accuracy and lower computational cost, but also produce a highly sparse discrete algebraic system. These observations are obtained by taking the direct approach of numerical experimentation.

1. Introduction

In this paper, we introduce the hierarchical radial basis functions method for the approximation to Sobolev functions and the collocation to well-posed linear partial differential equations. Hierarchical basis may be a common concept in finite element books, which is formed by decomposing finite element spaces (see an earlier paper by Yserentant, in 1986 [1]). Hierarchical radial basis functions network is a concept in neural networks, which is an approximating neural model and is a self-organizing (by growing) multiscale version of a radial basis function network. In this paper, hierarchical radial basis functions (H-RBFs) refer to some basis functions that are constructed by employing hierarchical data structure and scaled compactly supported radial basis functions. Below, we can explain why hierarchical radial basis functions are necessary.

As the first kernel-based collocation method was introduced by Kansa [2], radial basis functions have been successfully applied to solving various partial differential equations numerically. However, unfortunately, we are obtaining highly accurate solutions from severely ill-conditioned linear systems and high computational cost when using radial basis functions collocation methods. This is a well-known uncertainty or trade-off principle [3]. To overcome this problem, different strategies have been suggested, see a summary of existing methods in the recent monograph [4]. There are about three kinds of popular methods to deal with this problem.

The first method is finding optimal shape parameter $\epsilon$ (which is related to the scattered centers distribution, and usually is inversely proportional to mesh norm). We have become accustomed to scaling one of radial functions through multiplying the independent variable $\parallel ·\parallel$ by a shape parameter $\epsilon$ in the practical approximation. Obviously, a smaller value of $\epsilon$ causes the function to become flatter, whereas increasing $\epsilon$ leads to a more peaked radial function, and therefore localizes its influence. The choice of $\epsilon$ has a profound influence on both the approximation accuracy and numerical stability of the solution to interpolation problem. The general observation was that a large $\epsilon$ leads to a very well-conditioned linear system but also a poor approximation rate, whereas a smaller $\epsilon$ yields excellent approximation at the price of a badly conditioned system. A number of strategies for choosing a “good” value of $\epsilon$ have been suggested, such as the power function as indicator, the Cross Validation algorithm and the Contour–Padé algorithm. Some discussion of the existing parametrization schemes is provided in books [4,5]. Madych [6] indicated that the interpolation error goes to zero as $\epsilon \to 0$, this does not seem to be true in practice. Beyond that, there is no case known where the error and the sensitivity are both reasonably small.

The second method is the utilization of compactly supported radial basis functions (CSRBFs). Since 1995, a series of compactly supported radial basis functions have emerged gradually, such as Wendland’s function [7], Wu’s function [8], Buhmann’s function [9], and others [5]. The compact support automatically ensures the strict positive definiteness of CSRBFs. However, an additional trade-off principle remains, now depending on the choice of the support size. Specifically, a small support leads to a well-conditioned system but also poor approximation accuracy, whereas a larger support yields excellent accuracy at the price of ill-conditioned system.

The third method is stationary multilevel interpolation. With a stationary multiscale algorithm [10], the condition number of the discrete matrix can be relatively small, and computation can be performed in $O(N)$ operations. In this method, the present problem (interpolation or numerical PDEs) is solved first on the coarsest level by one of the compactly supported radial basis functions with a larger support (usually scaling the size of the support with the fill distance). Then, the residual can be formed, and then computed on the next finer level by the same compactly supported radial basis function but with a smaller support. This process can be repeated and be stopped on the finest level. The final approximation is the sum of all of interpolants. For interpolation problems, the linear convergence order has been proved in Sobelev spaces on the sphere in [11], and on bounded domains in [12]. Applications of this algorithm in solving PDEs on spheres are proposed in [13], on bounded domains were proposed in [14,15,16]. However, a series of global interpolation or approximation problems must be solved on different levels, although two kinds of local algorithms have been derived in [17,18] recently.

This paper considers solving the present problem on a single level. At the current scattered data set the trial function is represented as linear combination of hierarchical radial basis functions. The approximation method can produce a sparse discrete algebraic system, because hierarchical radial basis functions are derived from CSRBFs with different support radii. Compared with compactly supported radial basis functions approximation [7,8] and stationary multilevel approximation [11,12,13,14,15,16,17,18], the new method can solve the present problem on a single level with higher accuracy and lower computational cost. The effectiveness of H-RBFs collocation method will be conformed by several numerical observations.

2. H-RBFs Trial Spaces

In this section, we build H-RBFs trial spaces and assemble the discrete spaces with appropriate norms. In particular, we will prove two norm equivalence theorems and describe a function spaces commuting diagram. To avoid writing constants repeatedly, we shall use notations ≲ and ≅. Short notation $x\lesssim y$ means $x\le C_{1}y$, and $x\cong y$ means $C_{1}x\le y\le C_{2}x$, and $C_1$ and $C_2$ are positive constants.

2.1. H-RBFs Trial Spaces

Let $X=\{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\}$ be a finite point set in $\Omega \subseteq {\mathbb{R}}^d$. We define a commonly used trial discretization parameter:

$h_{X,\Omega }:=\underset{\mathbf{x}\in \Omega }{sup}\underset{{\mathbf{x}}_j\in X}{min}{\parallel \mathbf{x}-{\mathbf{x}}_j\parallel }_2,$

which can be regarded as the radius of the largest empty ball that can be placed among the data sites ${\mathbf{x}}_j$. To build H-RBFs trial spaces, we need the nested point sets

$X_1\subset X_2\subset \cdots \subset X_j\subset \cdots$

which have trial discretization parameters $h_j=h_{X_j,\Omega }$. Let $X_j^{*}$ be newly added point set in $X_j$, then we have $X_1=X_1^{*}$ and $X_j=X_{j-1}\bigcup X_j^{*}$ for $j=2,\dots$.

Given a kind of the compactly supported radial functions $\Phi :{\mathbb{R}}^d\to \mathbb{R}$, we can rescale it (by a scaling parameter ${\epsilon }_j>0$, and in the translation-invariant kernel-based case) as

(1)${\Phi }_j(·,\mathbf{y})={\epsilon }_j^d\Phi ({\epsilon }_j(·-\mathbf{y})),\mathbf{y}\in X_j^{*}.$

To make the support radii of ${\Phi }_j$ become smaller and smaller with the addition of new points, we select

(2)${\epsilon }_j=b\frac{1}{h_j},\mathrm{with}\mathrm{a}\mathrm{given}\mathrm{constant}0<b\le 1.$

Let $V_j=V_{j-1}\oplus W_j$, with $V_1=W_1$ and $W_j=\mathrm{span}\{{\Phi }_j(·,\mathbf{y}):\mathbf{y}\in X_j^{*}\}$ for $j=1,\dots$. Then we can build the H-RBFs trial spaces of the form

$V_j=\mathrm{span}\{\underbrace{{\Phi }_1(\mathbf{x},{\mathbf{x}}_1^1),\cdots ,{\Phi }_1(\mathbf{x},{\mathbf{x}}_1^{M_1})},\cdots ,\underbrace{{\Phi }_j(\mathbf{x},{\mathbf{x}}_j^1),\cdots ,{\Phi }_j(\mathbf{x},{\mathbf{x}}_j^{M_j})}\},$

which differs from the RBFs approximation spaces

$U_j=\mathrm{span}\{\underbrace{\Phi (\mathbf{x},{\mathbf{x}}_1^1),\cdots ,\Phi (\mathbf{x},{\mathbf{x}}_1^{M_1})},\cdots ,\underbrace{\Phi (\mathbf{x},{\mathbf{x}}_j^1),\cdots ,\Phi (\mathbf{x},{\mathbf{x}}_j^{M_j})}\}.$

Here, we have written $X_k^{*}=\{{\mathbf{x}}_k^1,\cdots ,{\mathbf{x}}_k^{M_k}\}$ for $k=1,2,\dots ,j$.

2.2. Norms of the Discrete Spaces $V_j$

The following is a relation between several function spaces.

$W^{m,2}(\Omega )\underset{\begin{array}{c}(3)\\ \mathit{When}\widehat{\Phi }\mathit{decays}\\ \mathit{only}\mathit{algebraically}\\ \mathit{with}m>\frac{\mathit{d}}{\mathit{2}}\end{array}}{⟵}{\mathcal{N}}_{\Phi }(\Omega )\underset{\begin{array}{c}(2)\\ \mathit{When}\Phi \\ \mathit{is}\mathit{strictly}\\ \mathit{positive}\mathit{definite}\end{array}}{⟵}{H}_{\Phi }(\Omega )\underset{\begin{array}{c}(1)\\ \mathit{When}\mathit{point}\\ \mathit{evaluation}\mathit{functionals}\\ \mathit{are}\mathit{continuous}\end{array}}{⟵}H(\Omega )$

Adding a restrictive condition (1), the Hilbert space $H(\Omega )$ becomes the reproducing kernel Hilbert space $H_{\Phi }(\Omega )$ with reproducing $\Phi :\Omega \times \Omega \to \mathbb{R}$. ${\mathcal{N}}_{\Phi }(\Omega )$ is a reproducing kernel Hilbert space with a strictly positive definite kernel $\Phi$. ${\mathcal{N}}_{\Phi }(\Omega )$ is the native space of $\Phi$, which contains all functions of the form $g={\sum }_{j=1}^N{c}_j\Phi (·,{\mathbf{x}}_j)$ provided ${\mathbf{x}}_j\in \Omega$, with

(3)${\parallel g\parallel }_{{\mathcal{N}}_{\Phi }(\Omega )}^2=\sum _{j=1}^N\sum _{k=1}^N{c}_j{c}_k\Phi ({\mathbf{x}}_j,{\mathbf{x}}_k).$

When $\Omega ={\mathbb{R}}^d$, we get a characterization of ${\mathcal{N}}_{\Phi }({\mathbb{R}}^d)$ in terms of the Fourier transforms

(4)${\parallel g\parallel }_{{\mathcal{N}}_{\Phi }({\mathbb{R}}^d)}^2={\int }_{{\mathbb{R}}^d}\frac{|\widehat{g}{(\omega )|}^2}{\widehat{\Phi }(\omega )}d\omega ,g\in L^2({\mathbb{R}}^d).$

Here

$\widehat{\Phi }(\mathit{\omega })={(2\pi )}^{-\frac{d}{2}}{\int }_{{\mathbb{R}}^d}\Phi (\mathbf{x})e^{-i{\mathbf{x}}^T\mathit{\omega }}d\mathbf{x}.$

In fact, there exist such radial functions whose Fourier transforms decay only algebraically, satisfying

(5)$\widehat{\Phi }(\mathit{\omega })\cong {\left(1+{\parallel \mathit{\omega }\parallel }_2^2\right)}^{-\tau },\mathit{\omega }\in {\mathbb{R}}^d.$

Then, scaled radial functions ${\Phi }_j(\mathbf{x})={\epsilon }_j^d\Phi ({\epsilon }_j\mathbf{x})$ have Fourier transforms satisfying

(6)$\widehat{{\Phi }_j}(\mathit{\omega })\cong {\left(1+\frac{1}{{\epsilon }_j^2}{\parallel \mathit{\omega }\parallel }_2^2\right)}^{-\tau },\mathit{\omega }\in {\mathbb{R}}^d.$

To clarify the relation between Sobolev spaces and native spaces, we cite the following extension operator theorems.

Theorem 1.

(Sobolev extension operator, see Section 5.17 in [19].)

Suppose $\Omega \subseteq {\mathbb{R}}^d$ is open and has a Lipschitz boundary. Let $\tau \ge 0$. Then. for all $g\in H^{\tau }(\Omega )$, there exists a linear operator, $E:H^{\tau }(\Omega )\to H^{\tau }({\mathbb{R}}^d)$, such that

(1) ${Eg|}_{\Omega }=g$,

(2) ${\parallel Eg\parallel }_{H^{\tau }({\mathbb{R}}^d)}\lesssim {\parallel g\parallel }_{H^{\tau }(\Omega )}$.

Theorem 2.

(Native extension operator, see Section 10.7 in [20].)

Suppose Φ is a strictly positive definite kernel. Let ${\Omega }_1\subseteq {\Omega }_2\subseteq {\mathbb{R}}^d$. Then, each function $g\in {\mathcal{N}}_{\Phi }({\Omega }_1)$ has a natural extension to a function, $\tilde{E}g\in {\mathcal{N}}_{\Phi }({\Omega }_2)$, such that

(7) $\parallel \tilde{E}{g\parallel }_{{\mathcal{N}}_{\Phi }({\Omega }_2)}={\parallel g\parallel }_{{\mathcal{N}}_{\Phi }({\Omega }_1)}.$

The main difference between these two theorems is that the extensions for functions from Sobolev spaces impose a restriction on regions $\Omega$, while extensions from native spaces are adequate for more general regions.

Then, we have the first of the following norm equivalence theorems.

Theorem 3.

Suppose $\Phi \in L^1({\mathbb{R}}^d)\cap C({\mathbb{R}}^d)$ is a strictly positive definite kernel satisfying (5) with $\tau >\frac{d}{2}$, and ${\Phi }_j$ is defined by (1) with ${\epsilon }_j\ge 1$. Then we have $H^{\tau }({\mathbb{R}}^d)={\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)$, and for every $g\in H^{\tau }({\mathbb{R}}^d)$, we have

${\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)}\lesssim {\parallel g\parallel }_{H^{\tau }({\mathbb{R}}^d)}\lesssim {\epsilon }_j^{\tau }{\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)}.$

Proof. 

Using norm definition of Sobolev space $H^{\tau }({\mathbb{R}}^d)$, definition (4), and inequalities (5) and (6), we have

$\begin{array}{ccc}\hfill {\parallel g\parallel }_{H^{\tau }({\mathbb{R}}^d)}^2& =& {\int }_{{\mathbb{R}}^d}{\left|\widehat{g}(\mathit{\omega })\right|}^2{\left(1+{\parallel \mathit{\omega }\parallel }_2^2\right)}^{\tau }d\mathit{\omega }\hfill \\ \hfill & =& {\epsilon }_j^{2\tau }{\int }_{{\mathbb{R}}^d}{\left|\widehat{g}(\mathit{\omega })\right|}^2{\left(\frac{1}{{\epsilon }_j^2}+\frac{1}{{\epsilon }_j^2}{\parallel \mathit{\omega }\parallel }_2^2\right)}^{\tau }d\mathit{\omega }\hfill \\ \hfill & \le & {\epsilon }_j^{2\tau }{\int }_{{\mathbb{R}}^d}{\left|\widehat{g}(\mathit{\omega })\right|}^2{\left(1+\frac{1}{{\epsilon }_j^2}{\parallel \mathit{\omega }\parallel }_2^2\right)}^{\tau }d\mathit{\omega }\hfill \\ \hfill & \lesssim & {\epsilon }_j^{2\tau }{\int }_{{\mathbb{R}}^d}\frac{|\widehat{g}{(\mathit{\omega })|}^2}{\widehat{{\Phi }_j}(\mathit{\omega })}d\mathit{\omega }\hfill \\ \hfill & =& {\epsilon }_j^{2\tau }{\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)}^2.\hfill \end{array}$

By similar scaling method, the lower bound follows from

$\begin{array}{ccc}\hfill {\parallel g\parallel }_{H^{\tau }({\mathbb{R}}^d)}^2& =& {\int }_{{\mathbb{R}}^d}{\left|\widehat{g}(\mathit{\omega })\right|}^2{\left(1+{\parallel \mathit{\omega }\parallel }_2^2\right)}^{\tau }d\mathit{\omega }\hfill \\ \hfill & \ge & {\int }_{{\mathbb{R}}^d}{\left|\widehat{g}(\mathit{\omega })\right|}^2{\left(1+\frac{1}{{\epsilon }_j^2}{\parallel \mathit{\omega }\parallel }_2^2\right)}^{\tau }d\mathit{\omega }\hfill \\ \hfill & \gtrsim & {\int }_{{\mathbb{R}}^d}\frac{|\widehat{g}{(\mathit{\omega })|}^2}{\widehat{{\Phi }_j}(\mathit{\omega })}d\mathit{\omega }\hfill \\ \hfill & =& {\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)}^2.\hfill \end{array}$

 □

A similar norm equivalence theorem with inverse of ${\epsilon }_j$ as scaling parameter can be found in [12]. If further assumptions are made on boundary of $\Omega$, the results of Theorem 3 will be still true for ${\mathcal{N}}_{{\Phi }_j}(\Omega ),\Omega \subseteq {\mathbb{R}}^d$. Then, we have another norm equivalence theorem.

Theorem 4.

Suppose $\Phi \in L^1({\mathbb{R}}^d)\cap C({\mathbb{R}}^d)$ is a strictly positive definite kernel satisfying (5) with $\tau >\frac{d}{2}$, and ${\Phi }_j$ is defined by (1) with ${\epsilon }_j\ge 1$. If $\Omega \subseteq {\mathbb{R}}^d$ is open and has a Lipschitz boundary, then $H^{\tau }(\Omega )={\mathcal{N}}_{{\Phi }_j}(\Omega )$ with equivalent norms, and, for every $g\in H^{\tau }(\Omega )$, we have

${\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}(\Omega )}\lesssim {\parallel g\parallel }_{H^{\tau }(\Omega )}\lesssim {\epsilon }_j^{\tau }{\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}(\Omega )}.$

Proof. 

Every $g\in {\mathcal{N}}_{{\Phi }_j}(\Omega )$ has an extension $\tilde{E}g\in {\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)$. By Theorem 3, $\tilde{E}g\in {\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)=H^{\tau }({\mathbb{R}}^d)$. Therefore, we have $g=\tilde{E}{g|}_{\Omega }\in H^{\tau }(\Omega )$ and

${\parallel g\parallel }_{H^{\tau }(\Omega )}\le \parallel \tilde{E}{g\parallel }_{H^{\tau }({\mathbb{R}}^d)}\lesssim {\epsilon }_j^{\tau }\parallel \tilde{E}{g\parallel }_{{\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)}\lesssim {\epsilon }_j^{\tau }{\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}(\Omega )}.$

In addition, every $g\in H^{\tau }(\Omega )$ has an extension $Eg\in H^{\tau }({\mathbb{R}}^d)={\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)$. Thus, ${g=Eg|}_{\Omega }\in {\mathcal{N}}_{{\Phi }_j}(\Omega )$ and

${\parallel g\parallel }_{{\mathcal{N}}_{{\Phi }_j}(\Omega )}\le {\parallel Eg\parallel }_{{\mathcal{N}}_{{\Phi }_j}({\mathbb{R}}^d)}\lesssim {\parallel Eg\parallel }_{H^{\tau }({\mathbb{R}}^d)}\lesssim {\parallel g\parallel }_{H^{\tau }(\Omega )}.$

 □

We can understand these two extension operator theorems (Theorems 1 and 2) and two norm equivalence theorems (Theorems 3 and 4) using a diagram. More concretely, under the following conditions,

• (c1) $\Omega \subseteq {\mathbb{R}}^d$ is open and has a Lipschitz boundary

• (c2) ${\Phi }_j(\mathbf{x})={\epsilon }_j^d\Phi ({\epsilon }_j\mathbf{x})$ with ${\epsilon }_j\ge 1$

• (c3) $\Phi \in L^1({\mathbb{R}}^d)\cap C({\mathbb{R}}^d)$

• (c4) $\widehat{\Phi }(\omega )\cong {\left(1+{\parallel \omega \parallel }_2^2\right)}^{-\tau }$, $\tau >\frac{d}{2}$,

we have the following commuting diagram and corollary.

Corollary 1.

Under conditions (c1)–(c4), $V_j\subseteq {\mathcal{N}}_{\Phi }(\Omega )$, and we can assemble $V_j$ with ${\mathcal{N}}_{\Phi }(\Omega )$ norm.

3. Interpolation via H-RBFs

In this section we discuss the scattered data interpolation with hierarchical radial basis functions.

Given a target function $f:\Omega (\subseteq {\mathbb{R}}^d)\to \mathbb{R}$, now we can find a interpolant $f_j\in V_j$ of the form

(8)$f_j(\mathbf{x})=\sum _{i=1}^j\sum _{k=1}^{M_i}{c}_i^k{\Phi }_i(\mathbf{x},{\mathbf{x}}_i^k).$

The coefficients $c_i^k$ are found by enforcing the interpolation conditions

(9)$f_j(\mathbf{x})=\sum _{i=1}^j\sum _{k=1}^{M_i}{c}_i^k{\Phi }_i(\mathbf{x},{\mathbf{x}}_i^k)=f(\mathbf{x}),\mathrm{for}\mathrm{every}\mathbf{x}\in Y_j,$

where we use $Y_j$ as testing data on $j$-th level. This may lead to an unsymmetric or even nonsquare discrete system. The linear system will be solved by least square method (it is a Matlab ’\’ operation in our experiments). Computations were performed on a laptop with 2.4 GHz Intel Core i7 processor, using MATLAB running in the Windows 7 operating system.

First, we generate the nested spaced Halton points $N=$ 9, 25, 81, 289, 1089, and 4225 in the interior of the domain $\Omega ={[0,1]}^2$ (blue points in Figure 1). In Figure 1, we display six data sets used in the experiments. We choose all of blue points as the whole centers for the basis functions, and use Wendland’s compactly supported function $\varphi (r)={(1-r)}_{+}^6(35r^2+18r+3)$ to interpolate 2D Franke’s function. $M=40\times 40$ equally spaced evaluation points are used to compute RMS-error:

$\mathrm{RMS}-\mathrm{error}=\sqrt{\frac{1}{M}\sum _{k=1}^M{[f({\xi }_k)-f_j({\xi }_k)]}^2}=\frac{1}{\sqrt{M}}{\parallel f-f_j\parallel }_2.$

We now present four sets of interpolation experiments: nonstationary CSRBFs interpolation ($\epsilon$ fixed), stationary CSRBFs interpolation (with an initial $\epsilon$, then double its value for every successive level), multilevel interpolation, and H-RBFs interpolation. For the first three methods, we can let testing data be taken to be the same as the centers, because the discrete matrices produced by these methods are nonsingular under this choice. However, a nonsquare testing is necessary for H-RBFs interpolation method. To do this simply, we use $N+4(\sqrt{N}-1)$ Halton points to test (9) on each level. We list numerical results (including RMS-error, convergence rates, and total CPU time) in Table 1, Table 2, Table 3 and Table 4.

In the nonstationary case (Table 1), we have convergence, although it is not obvious what the rate might be. However, the computation requires lost of time because the matrices become increasingly dense. The stationary CSRBFs interpolation is also numerically stable, but there will be essentially no convergence (see Table 2). Table 3 is the numerical results for multilevel interpolation. The linear convergence behavior for Sobolev functions interpolation has been proved by Wendland in [12]. The corresponding 3D multilevel experiment can be found in Fasshauer’s book [5], and 2D uniformly data interpolation experiment is in [12]. From Table 3, we observe that the convergence seems cease at a later stage when given a initial $\epsilon =0.5$. Compared with nonstationary CSRBFs interpolation, multilevel method saves much time but with limited accuracy.

The corresponding RMS-error and observed convergence rates for the H-RBFs method is listed in Table 4. Several observations can be made by looking at Table 4. The H-RBFs interpolation method is numerically stable and has relatively small errors, even for $\epsilon =0.5$ case. The RMS-error on 4225 level has been remarkably reduced to ${10}^{-6}$. Compared with CSRBFs and multilevel interpolation, the present H-RBF method can doubtlessly solve problem (9) with a higher accuracy and lower computational cost. The sparsity behavior of H-RBFs interpolation matrices are displayed in Figure 2.

4. Collocation via H-RBFs

In this section, we discuss the implementation of the H-RBFs collocation method for linear partial differential equations.

4.1. Example 1

We consider the following Poisson problem with Dirichlet boundary conditions,

(10)$\begin{array}{ccc}\hfill {\nabla }^{2}u(x,y)& =& -\frac{5}{4}{\pi }^{2}sin(\pi x)cos(\frac{\pi y}{2}),(x,y)\in \Omega ={[0,1]}^2,\hfill \\ \hfill u(x,y)& =& sin(\pi x),(x,y)\in {\Gamma }_1,\hfill \\ \hfill u(x,y)& =& 0,(x,y)\in {\Gamma }_2,\hfill \end{array}$

where ${\Gamma }_1=\{(x,y):0\le x\le 1,y=0\}$ and ${\Gamma }_2=\partial \Omega \setminus {\Gamma }_1$. The exact solution of the problem (10) is given by

$u(x,y)=sin(\pi x)cos(\frac{\pi y}{2}).$

We use the unsymmetric Kansa method to solve problem (10). At this time, a nonsquare testing is necessary because the small trial space must be tested on a fine-grained space discretization according to Schaback’s theory [21,22]. As always, the interior collocation data are taken to be the same as the centers. We create additional $4(\sqrt{N}-1)$ equally spaced collocation points for the boundary conditions on each level (red points in Figure 1). That is, on each level, the centers only include blue points, whereas the collocation sites contain blue points in the interior of the domain and red points on the boundary. Consequently, the collocation matrix becomes nonsquare, because the test side has more degrees of freedom than trial side. In our experiments, Wendland’s $C^6$ function $\varphi (r)={(1-r)}_{+}^8(32r^3+25r^2+8r+1)$ is used to construct trial spaces. We list RMS-error and total CPU time (last row in the Tables) for CSRBFs collocation method, multilevel collocation method and H-RBFs collocation method in Table 5, Table 6, Table 7, Table 8, Table 9 and Table 10, respectively.

Table 5 is the nonstationary CSRBFs collocation with fixed scaling parameters of $\epsilon =0.5$ and $\epsilon =0.25$. Table 6 is the stationary setting with two different initial parameter values, $0.5$ and $0.25$, and then doubles its value for every successive level. We note that the nonconvergence behavior can be observed for the stationary collocation method. Similar observations are reported in detail in Fasshauer’s book [5], where the boundary centers were selected to coincide with the boundary collocation points.

Table 7 and Table 8 are numerical results for multilevel collocation method. Several observations can be made by looking at these two tables. We observe that multilevel collocation method is nonconvergent with some relatively large initial values, such as $0.5$, $0.25$, and $0.125$. With an initial parameter $0.0625$, the method has convergence behavior. However, it seems that the convergence ceases at a later stage. But the multilevel collocation method will have linear convergence if we putting additional centers on boundary. The corresponding numerical experiments have been shown in book [5]. Generally, the bandwidth of the collocation matrix must be allowed to increase slowly from one level to the next (namely, the support radii go slower to zero than the mesh norms) when solving linear partial differential equations by multilevel collocation method. In fact, these phenomena have been noted in Fasshauer’s numerical observations in [23] and theoretical explanations in [16]. Through above experiments (Table 7 and Table 8), the conclusion that we need to add is that the multilevel collocation method is terrible if the centers just being placed in the interior of the domain.

In Table 9 and Table 10, we list RMS-error and CPU time for the hierarchical radial basis functions collocation method. We observe that H-RBFs method has an ideal convergence behavior. Even to a relatively large initial parameter $\epsilon =0.5$, the convergence rate of H-RBFs method is close to 2. With an initial parameter $\epsilon =0.0625$, RMS-error on 4225 level is remarkably reduced to ${10}^{-8}$, whereas the CPU consumption time is only ~$37.6$ s. Compared with CSRBFs and multilevel collocation method, the present H-RBFs method can solve the model problem with a higher accuracy and lower computational cost. The sparsity behavior of H-RBFs collocation matrices are displayed in Figure 3.

4.2. Example 2

In this subsection, we consider the following Helmholtz test problem,

(11)$\begin{array}{ccc}\hfill -{\nabla }^{2}u(x,y)+u(x,y)& =& cos(\pi x)cos(\pi y),(x,y)\in \Omega ={[-1,1]}^2,\hfill \\ \hfill \frac{\partial }{\partial \mathbf{n}}u(x,y)& =& 0,\mathrm{on}\partial \Omega ,\hfill \end{array}$

where $\mathbf{n}$ denotes the unit outer normal vector. It is easy to verify that the exact solution for the problem (11) is given by

$u(x,y)=\frac{cos(\pi x)cos(\pi y)}{2{\pi }^2+1}.$

For this example, we use the same centers data, collocation sites, and Wendland’s $C^6$ function as in Section 4.1. The RMS-error, total CPU time, and observed convergence rates are listed in Table 11 and Table 12, and the sparsity behavior of the H-RBF collocation matrices of this problem are displayed in Figure 4. We observe that it seems that the convergence ceases at a later stage when using a relatively large initial value $\epsilon =0.5$. However, the H-RBF method has ideal convergence behavior and lower computational cost, as in example 1. We also display the plots of the absolute error in Figure 5.

5. Conclusions

To handle the long-standing trade-off principle in the RBF collocation method, we derived hierarchical radial basis functions (H-RBFs) collocation method in the paper. Based on a nested scattered point sets, H-RBF trial spaces are constructed using scaled compactly supported radial basis functions with varying support radii. Several numerical observations were demonstrated in the paper. The experiments showed that H-RBFs collocation method can retain high accuracy with a lower computational cost, when compared with existing CSRBFs collocation method and multilevel RBFs collocation method.

There are many possibilities for enhancement of this method:

• (1). A convergence proof for H-RBFs collocation method will depend on the approximation of H-RBFs trial spaces, new inverse inequality (a frequently used inequality has been given in [24] for RBFs case), and sampling theorem.

• (2). This method can be used for solving well-posed nonlinear partial differential equations, and the convergence analysis of hierarchical radial basis function collocation method for nonlinear discretization will depend on B$\ddot{\mathrm{o}}$hmer/Schaback theory [25].

Author Contributions

Conceptualization and Formal analysis, Z.L.; Methodology and Writing–original draft preparation, Q.X.

Funding

The research of the first author was partially supported by the Natural Science Foundations of Ningxia Province (No.2019AAC02001), the Natural Science Foundations of China (No.11501313), the Project funded by China Postdoctoral Science Foundation (No.2017M621343), and the Third Batch of Ningxia Youth Talents Supporting Program (No.TJGC2018037). Research of the second author was partially supported by the Natural Science Foundations of Ningxia Province (No.NZ2018AAC03026), and the Fourth Batch of Ningxia Youth Talents Supporting Program.

Acknowledgments

The authors are extremely grateful to Zongmin Wu from Fudan University, and would like to thank two unknown reviewers who made valuable comments on an earlier version of this paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Figures and Tables

Figure 1. Six scattered data sets used in experiments: N interior Halton points (blue) and 4(N−1) equally spaced boundary collocation points (red).

Figure 2. Sparsity behavior of H-RBF interpolation matrices at level 1, 2, 3, 4, 5, and 6.

Figure 3. Possion problem: Sparsity behavior of discrete matrices by H-RBFs collocation at levels 1, 2, 3, 4, 5, and 6.

Figure 4. Helmholtz problem: Sparsity behavior of discrete matrices by H-RBFs collocation at levels 1, 2, 3, 4, 5, and 6.

Figure 5. Final error graphs on level 6. Top left to bottom right: ε = 0.5, 0.25, 0.125, and 0.0625.

Table 1

Compactly supported radial basis function (CSRBF) interpolation: nonstationary.

Centers	$\mathit{\epsilon }=0.5$	Rate	$\mathit{\epsilon }=0.25$	Rate
9	1.113967 $\times {10}^{-1}$		1.537637 $\times {10}^{-1}$
25	6.950489 $\times {10}^{-2}$	0.680520	8.362841 $\times {10}^{-2}$	0.878650
81	7.783806 $\times {10}^{-3}$	3.158567	8.430085 $\times {10}^{-3}$	3.310374
289	4.849611 $\times {10}^{-4}$	4.004535	3.210030 $\times {10}^{-4}$	4.714889
1089	2.253331 $\times {10}^{-5}$	4.427738	2.115728 $\times {10}^{-5}$	3.923360
4225	1.128459 $\times {10}^{-6}$	4.319634	3.572115 $\times {10}^{-6}$	2.566304
	4.42e+01(s)		4.42e+01(s)

Table 2

CSRBF interpolation: stationary.

Centers	$\mathit{\epsilon }=0.5$	Rate	$\mathit{\epsilon }=0.25$	Rate
9	1.113967 $\times {10}^{-1}$		1.537637 $\times {10}^{-1}$
25	4.205301 $\times {10}^{-2}$	1.405425	6.950489 $\times {10}^{-2}$	1.145528
81	2.288906 $\times {10}^{-2}$	0.877551	9.036469 $\times {10}^{-3}$	2.943283
289	2.330558 $\times {10}^{-2}$	−0.026018	5.143627 $\times {10}^{-3}$	0.812973
1089	1.986107 $\times {10}^{-2}$	0.230732	4.590088 $\times {10}^{-3}$	0.164264
4225	2.010532 $\times {10}^{-2}$	−0.017633	4.395273 $\times {10}^{-3}$	0.062569
	1.46e+00(s)		3.16e+00(s)

Table 3

Multilevel interpolation.

Centers	Initialε = 0.5	Rate	Initialε = 0.25	Rate
9	1.113967 $\times {10}^{-1}$		1.537637 $\times {10}^{-1}$
25	5.321739 $\times {10}^{-2}$	1.065737	7.271406 $\times {10}^{-2}$	1.080408
81	8.346277 $\times {10}^{-3}$	2.672693	5.457674 $\times {10}^{-3}$	3.735876
289	3.142086 $\times {10}^{-3}$	1.409410	9.231828 $\times {10}^{-4}$	2.563598
1089	1.646392 $\times {10}^{-3}$	0.932414	2.752895 $\times {10}^{-4}$	1.745667
4225	9.882529 $\times {10}^{-4}$	0.736356	1.053762 $\times {10}^{-4}$	1.385401
	2.08e+00(s)		5.53e+00(s)

Table 4

Hierarchical radial basis function (H-RBF) interpolation.

Centers	Initialε = 0.5	Rate	Initialε = 0.25	Rate
9	8.715764 $\times {10}^{-2}$		9.975067 $\times {10}^{-2}$
25	4.380040 $\times {10}^{-2}$	0.992683	5.449353 $\times {10}^{-2}$	0.872242
81	8.765413 $\times {10}^{-3}$	2.321050	2.426089 $\times {10}^{-2}$	1.167453
289	1.146738 $\times {10}^{-3}$	2.934287	8.471348 $\times {10}^{-4}$	4.839897
1089	8.960631 $\times {10}^{-5}$	3.677791	8.592933 $\times {10}^{-6}$	6.623297
4225	9.886575 $\times {10}^{-6}$	3.180058	1.835490 $\times {10}^{-6}$	2.226986
	4.42e+00(s)		9.97e+00(s)

Table 5

CSRBFs collocation: nonstationary.

Centers	$\mathit{\epsilon }=0.5$	Rate	$\mathit{\epsilon }=0.25$	Rate
9	5.154990 $\times {10}^{-1}$		2.312114 $\times {10}^{-1}$
25	2.523773 $\times {10}^{-1}$	1.030387	2.501198 $\times {10}^{-2}$	3.208521
81	5.129004 $\times {10}^{-2}$	2.298832	8.196424 $\times {10}^{-3}$	1.609552
289	1.169028 $\times {10}^{-2}$	2.133369	1.976536 $\times {10}^{-3}$	2.052021
1089	9.400882 $\times {10}^{-4}$	3.636370	2.150290 $\times {10}^{-4}$	3.200371
4225	7.877153 $\times {10}^{-5}$	3.577050	4.789840 $\times {10}^{-5}$	2.166482
	5.61e+01(s)		5.59e+01(s)

Table 6

CSRBFs collocation: stationary.

Centers	Initialε = 0.5	Rate	Initialε = 0.25	Rate
9	5.154990 $\times {10}^{-1}$		2.312114 $\times {10}^{-1}$
25	1.880555 $\times {10}^{-1}$	1.454811	2.523773 $\times {10}^{-1}$	−0.126370
81	1.845714 $\times {10}^{-1}$	0.026979	2.415339 $\times {10}^{-1}$	0.063357
289	3.232111 $\times {10}^{-1}$	−0.808298	1.619474 $\times {10}^{-1}$	0.576701
1089	4.385499 $\times {10}^{-1}$	−0.440264	2.008842 $\times {10}^{-1}$	−0.310839
4225	4.782622 $\times {10}^{-1}$	−0.125061	2.480243 $\times {10}^{-1}$	−0.304117
	2.11e+00(s)		6.14e+00(s)

Table 7

Multilevel collocation: with initial 0.5 and 0.25.

Centers	Initialε = 0.5	Rate	Initialε = 0.25	Rate
9	5.154990 $\times {10}^{-1}$		2.312114 $\times {10}^{-1}$
25	5.196034 $\times {10}^{-1}$	−0.011441	1.827791 $\times {10}^{-1}$	0.339111
81	5.185346 $\times {10}^{-1}$	0.002970	1.886535 $\times {10}^{-1}$	−0.045637
289	5.185903 $\times {10}^{-1}$	−0.000155	1.890747 $\times {10}^{-1}$	−0.003218
1089	5.185993 $\times {10}^{-1}$	−0.000025	1.891628 $\times {10}^{-1}$	−0.000672
4225	5.186014 $\times {10}^{-1}$	−0.000006	1.891758 $\times {10}^{-1}$	−0.000099
	3.19e+00(s)		9.47e+00(s)

Table 8

Multilevel collocation: with initial 0.125 and 0.0625.

Centers	Initialε = 0.125	Rate	Initialε = 0.0625	Rate
9	3.345335 $\times {10}^{-1}$		2.877511 $\times {10}^{-1}$
25	1.189874 $\times {10}^{-1}$	1.491342	4.878530 $\times {10}^{-2}$	2.560303
81	1.238221 $\times {10}^{-1}$	−0.057460	2.944292 $\times {10}^{-2}$	0.728526
289	1.166663 $\times {10}^{-1}$	0.085881	1.706502 $\times {10}^{-2}$	0.786879
1089	1.158203 $\times {10}^{-1}$	0.010500	7.585827 $\times {10}^{-3}$	1.169664
4225	1.159396 $\times {10}^{-1}$	−0.001485	6.469642 $\times {10}^{-3}$	0.229621
	2.03e+01(s)		4.21e+01(s)

Table 9

H-RBFs collocation of example 1: with initial 0.5 and 0.25.

Centers	Initialε = 0.5	Rate	Initialε = 0.25	Rate
9	5.154990 $\times {10}^{-1}$		2.312114 $\times {10}^{-1}$
25	1.200069 $\times {10}^{-1}$	2.102853	4.838983 $\times {10}^{-2}$	2.256436
81	2.112357 $\times {10}^{-2}$	2.506192	3.451755 $\times {10}^{-3}$	3.809302
289	7.398191 $\times {10}^{-3}$	1.513609	1.075014 $\times {10}^{-4}$	5.004903
1089	2.440784 $\times {10}^{-3}$	1.599828	5.835322 $\times {10}^{-6}$	4.203399
4225	5.606710 $\times {10}^{-4}$	2.122118	2.751613 $\times {10}^{-7}$	4.406463
	4.07e+00(s)		9.69e+00(s)

Table 10

H-RBFs collocation of example 1: with initial 0.125 and 0.0625.

Centers	Initialε = 0.125	Rate	Initialε = 0.0625	Rate
9	3.345335 $\times {10}^{-1}$		2.877511 $\times {10}^{-1}$
25	1.060521 $\times {10}^{-2}$	4.979305	6.205372 $\times {10}^{-3}$	5.535160
81	1.898035 $\times {10}^{-3}$	2.482195	1.655041 $\times {10}^{-3}$	1.906651
289	2.565084 $\times {10}^{-5}$	6.209357	2.628040 $\times {10}^{-5}$	5.976736
1089	6.264679 $\times {10}^{-7}$	5.355622	4.950454 $\times {10}^{-7}$	5.730282
4225	2.553269 $\times {10}^{-8}$	4.616823	1.102677 $\times {10}^{-8}$	5.488478
	1.97e+01(s)		3.76e+01(s)

Table 11

H-RBFs collocation of example 2: with initial 0.5 and 0.25.

Centers	Initialε = 0.5	Rate	Initialε = 0.25	Rate
9	3.117164 $\times {10}^{-2}$		1.719666 $\times {10}^{-1}$
25	2.035525 $\times {10}^{-2}$	0.614833	8.951139 $\times {10}^{-2}$	0.941985
81	1.837566 $\times {10}^{-2}$	0.147605	8.238622 $\times {10}^{-3}$	3.441596
289	1.749727 $\times {10}^{-2}$	0.070666	3.459194 $\times {10}^{-3}$	1.251967
1089	1.744016 $\times {10}^{-2}$	0.004717	1.833905 $\times {10}^{-3}$	0.915517
4225	1.787914 $\times {10}^{-2}$	−0.035864	1.024535 $\times {10}^{-4}$	4.161878
	1.94e+00(s)		5.48e+00(s)

Table 12

H-RBFs collocation of example 2: with initial 0.125 and 0.0625.

Centers	Initialε = 0.125	Rate	Initialε = 0.0625	Rate
9	5.473827 $\times {10}^{-1}$		3.289428 $\times {10}^{-1}$
25	1.411404 $\times {10}^{-2}$	5.277347	6.743427 $\times {10}^{-2}$	2.286283
81	2.332565 $\times {10}^{-3}$	2.597141	2.429136 $\times {10}^{-3}$	4.794967
289	4.059524 $\times {10}^{-4}$	2.522535	2.783964 $\times {10}^{-4}$	3.125231
1089	1.776398 $\times {10}^{-5}$	4.514284	1.239397 $\times {10}^{-6}$	7.811358
4225	5.102425 $\times {10}^{-7}$	5.121628	3.532572 $\times {10}^{-8}$	5.132776
	1.25e+01(s)		2.62e+01(s)