On Multiscale RBF Collocation Methods for Solving the Monge–Ampère Equation

Abstract

This paper considers some multiscale radial basis function collocation methods for solving the two-dimensional Monge–Ampère equation with Dirichlet boundary. We discuss and study the performance of the three kinds of multiscale methods. The first method is the cascadic meshfree method, which was proposed by Liu and He (2013). The second method is the stationary multilevel method, which was proposed by Floater and Iske (1996), and is used to solve the fully nonlinear partial differential equation in the paper for the first time. The third is the hierarchical radial basis function method, which is constructed by employing successive refinement scattered data sets and scaled compactly supported radial basis functions with varying support radii. Compared with the first two methods, the hierarchical radial basis function method can not only solve the present problem on a single level with higher accuracy and lower computational cost but also produce highly sparse nonlinear discrete system. These observations are obtained by taking the direct approach of numerical experimentation.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The Monge–Ampère equation with Dirichlet boundary is a fully nonlinear partial differential equation, which is given by$$\begin{array}{}\text{(1)}& \text{det}\left(D^{2}u\left(z\right)\right)=f\left(z\right),\hspace{1em}\text{in\hspace{0.17em}}\Omega ,\text{(2)}& u\left(z\right)=g\left(z\right),\hspace{1em}\text{in\hspace{0.17em}}\partial \Omega ,\end{array}$$where z has n independent variables in a bounded domain $\Omega \subset {\mathrm{ℝ}}^d$ and $D^{2}u$ is the Hessian of the function u. When restricting it to $\Omega \subset {\mathrm{ℝ}}^2$, we can rewrite the equation as$$\begin{array}{}\text{(3)}& Fu=f,\hspace{1em}\text{in\hspace{0.17em}}\Omega ,\text{(4)}& u=g,\hspace{1em}\text{in\hspace{0.17em}}\partial \Omega ,\end{array}$$with a nonlinear map $F:u\mapsto u_{xx}{u}_{yy}-u_{xy}^2$.

The Monge–Ampère equation is originated in geometric surface theory and has been widely applied in dynamic meteorology, elasticity, geometric optics, image processing, conformal geometry, optimal transport, and others. And the numerical solution of the Monge–Ampère equation has been a subject of increasing interest recently. Refer Brenner et al. [1] and Feng and Neilan [2] study for some finite element methods and Oberman’s study [3] for finite difference methods. In order to design a certain type of numerical method which can avoid the tedious mesh generation and the domain integration and can be capable of dealing with any irregular distribution of nodes, some meshfree collocation methods have been researched by Liu and He [4, 5], Böhmer and Schaback [6], and Rashidinia and collaborators [7–9]. However, we are obtaining highly accurate solutions from severely ill-conditioned algebraic systems and high computational cost when using meshfree approximation. Compactly supported radial basis functions can lead to a very well-conditioned sparse system (see [10–12]) but at the cost of a poor approximation accuracy. This is a trade-off principle. That is to say, small support leads to a well-conditioned system but also poor accuracy, while large support yields excellent accuracy at the price of ill-conditioned system. Consequently, to speed up the meshfree collocation for the Monge–Ampère equation, we need to design some efficient multiscale radial basis function collocation methods. Undoubtedly, this is the motivation of the paper.

The remainder of this paper is organized as follows. In next section, we introduce several multiscale RBF collocation methods. Section 3 presents some numerical results, where we will compare these multiscale algorithms through direct numerical observation. And Section 4 closes the paper with a concise summery.

2. Multiscale RBF Collocation Methods

Let $X=\left\{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\right\}$ be a finite point set in $\Omega \subseteq {\mathrm{ℝ}}^2$. We define a commonly used trial discretization parameter:$$\begin{array}{}\text{(5)}& h_{X,\Omega }≔\underset{\mathbf{x}\in \Omega }{\sup }\underset{{\mathbf{x}}_j\in X}{\min }{‖\mathbf{x}-{\mathbf{x}}_j‖}_2,\end{array}$$which can be regarded as the radius of the largest empty ball that can be placed among the data sites ${\mathbf{x}}_j$. Generate nested point sets$$\begin{array}{}\text{(6)}& X_1\subset X_2\subset \cdots \subset X_j\subset \cdots ,\end{array}$$which have trial discretization parameters $h_j=h_{X_j,\Omega }$. Let $X_j^{\ast }$ be newly added point set in $X_j$, then we have $X_1=X_1^{\ast }$ and $X_j=X_{j-1}\cup X_j^{\ast }$ for $j=2,\dots ,$.

Given a kind of the compactly supported radial functions $\Phi :{\mathrm{ℝ}}^2⟶\mathrm{ℝ}$, we can rescale it (in the translation-invariant kernel-based case) as$$\begin{array}{}\text{(7)}& {\Phi }_j\left(\cdot ,y\right)={\epsilon }_j^d\Phi \left({\epsilon }_j\left(\cdot -y\right)\right),\hspace{1em}y\in X_j^{\ast }\text{\hspace{0.17em}or\hspace{0.17em}}y\in X_j.\end{array}$$

In order to make the support radii of ${\Phi }_j$ become smaller and smaller with the addition of new points, we will select$$\begin{array}{}\text{(8)}& {\epsilon }_j=b\frac{1}{h_j},\hspace{1em}\text{with\hspace{0.17em}a\hspace{0.17em}given\hspace{0.17em}constant\hspace{0.17em}}0<b\le 1.\end{array}$$

Let $V_j=V_{j-1}+W_j$, with $V_1=W_1$ and $W_j=\text{span}\left\{{\Phi }_j\left(\cdot ,\mathbf{y}\right):\mathbf{y}\in X_j^{\ast }\right\}$ for $j=1,\dots ,$. Then, we can build two types of trial spaces with the form$$\begin{array}{c}\text{(9)}& V_j=\text{span}\left\{\underbrace{{\Phi }_1\left(\mathbf{x},{\mathbf{x}}_1^1\right),\dots ,{\Phi }_1\left(\mathbf{x},{\mathbf{x}}_1^{M_1}\right)},\dots ,\underbrace{{\Phi }_j\left(\mathbf{x},{\mathbf{x}}_j^1\right),\dots ,{\Phi }_j\left(\mathbf{x},{\mathbf{x}}_j^{M_j}\right)}\right\},U_j=\text{span}\left\{\underbrace{{\Phi }_j\left(\mathbf{x},{\mathbf{x}}_1^1\right),\dots ,{\Phi }_j\left(\mathbf{x},{\mathbf{x}}_1^{M_1}\right)},\dots ,\underbrace{{\Phi }_j\left(\mathbf{x},{\mathbf{x}}_j^1\right),\dots ,{\Phi }_j\left(\mathbf{x},{\mathbf{x}}_j^{M_j}\right)}\right\}.\end{array}$$

Here, we have written $X_k^{\ast }=\left\{{\mathbf{x}}_k^1,\dots ,{\mathbf{x}}_k^{M_k}\right\}$ for $k=\mathrm{1,2},\dots ,j$. $V_j$ and $U_j$ have multiscale expressions, which differ from the traditional radial basis trial spaces$$\begin{array}{}\text{(10)}& Q_j=\text{span}\left\{\Phi \left(\mathbf{x},{\mathbf{x}}_1^1\right),\dots ,\Phi \left(\mathbf{x},{\mathbf{x}}_1^{M_1}\right),\dots ,\Phi \left(\mathbf{x},{\mathbf{x}}_j^1\right),\dots ,\Phi \left(\mathbf{x},{\mathbf{x}}_j^{M_j}\right)\right\}.\end{array}$$

Let $M={\displaystyle {\sum }_{k=1}^j{M}_k}$ be the total number of the centers and ${\left\{B\left(\mathbf{x},{\mathbf{x}}_k\right)\right\}}_{k=1}^M$ be basis functions from one of the trial spaces above. Then, with a trial function $v\left(\mathbf{x}\right)={\displaystyle {\sum }_{k=1}^M{c}_{k}B\left(\mathbf{x},{\mathbf{x}}_k\right)}$, the nonsymmetric kernel-based collocation method for equations (3) and (4) will generate the following nonlinear system$$\begin{array}{c}\text{(11)}& Fv\left(\mathbf{y}\right)=f\left(\mathbf{y}\right),\hspace{1em}\mathbf{y}\in Y^I,v\left(\mathbf{y}\right)=g\left(\mathbf{y}\right),\hspace{1em}\mathbf{y}\in Y^B,\end{array}$$for some given interior collocation sites $Y^I$ and boundary collocation sites $Y^B$ and some unknown coefficients $c_k$. We need to compute the derivatives of $B\left(\mathbf{x},{\mathbf{x}}_k\right)$ because the nonlinear operator F acts on the basis functions now. The derivatives can be obtained by the chain rule from Table 1.

Table 1

Derivatives of $\varphi \left(r\right)=\varphi \left(‖\mathbf{x}‖\right)$ with variables $\mathbf{x}=\left(x,y\right)\in {\mathrm{ℝ}}^2$.

	Partial derivatives
First-order	$\left(\partial /\partial x\right)\varphi \left(r\right)=\left(\left(x/r\right)\left(\text{d}/\text{d}r\right)\right)\varphi \left(r\right)$
	$\left(\partial /\partial y\right)\varphi \left(r\right)=\left(\left(y/r\right)\left(\text{d}/\text{d}r\right)\right)\varphi \left(r\right)$
	$\left(\partial /\partial \mathbf{n}\right)\varphi \left(r\right)=\left(\left(⟨\mathbf{x},\mathbf{n}⟩/r\right)\left(\text{d}/\text{d}r\right)\right)\varphi \left(r\right)$
Second-order	$\left({\partial }^2/\partial x^2\right)\varphi \left(r\right)=\left(\left(x^2/r^2\right)\left({\text{d}}^2/\text{d}{r}^2\right)\right)\varphi \left(r\right)+\left(\left(y^2/r^3\right)\left(\text{d}/\text{d}r\right)\right)\varphi \left(r\right)$
	$\left({\partial }^2/\partial y^2\right)\varphi \left(r\right)=\left(\left(y^2/r^2\right)\left({\text{d}}^2/\text{d}{r}^2\right)\right)\varphi \left(r\right)+\left(x^2/r^3\right)\left(\text{d}/\text{d}r\right)\varphi \left(r\right)$
	$\left({\partial }^2/\partial x\partial y\right)\varphi \left(r\right)=\left(\left(xy/r^2\right)\left({\text{d}}^2/\text{d}{r}^2\right)\right)\varphi \left(r\right)-\left(\left(xy/r^3\right)\left(\text{d}/\text{d}r\right)\right)\varphi \left(r\right)$

For example, the Wendland’s $C^6$ function $\varphi \left(r\right)={\left(1-\epsilon r\right)}_{+}^8\left(32{\left(\epsilon r\right)}^3+25{\left(\epsilon r\right)}^2+8\epsilon r+1\right)$ with scaling parameter ε has second-order partial derivatives:$$\begin{array}{c}\text{(12)}& \frac{{\partial }^2}{\partial x^2}\varphi \left(r\right)=22{\epsilon }^2{\left(1-\epsilon r\right)}_{+}^6\left(24{\epsilon }^2{x}^2\left(6\epsilon r+1\right)+\left(16{\left(\epsilon r\right)}^3-9{\left(\epsilon r\right)}^2-6\epsilon r-1\right)\right),\frac{{\partial }^2}{\partial y^2}\varphi \left(r\right)=22{\epsilon }^2{\left(1-\epsilon r\right)}_{+}^6\left(24{\epsilon }^2{y}^2\left(6\epsilon r+1\right)+\left(16{\left(\epsilon r\right)}^3-9{\left(\epsilon r\right)}^2-6\epsilon r-1\right)\right),\\ \frac{{\partial }^2}{\partial x\partial y}\varphi \left(r\right)=528{\epsilon }^{4}xy{\left(1-\epsilon r\right)}_{+}^6\left(6\epsilon r+1\right).\end{array}$$

2.1. Cascadic MeshFree Method

Once a convex solution be chosen, equations (3) and (4) can be rewritten in the equivalent form:$$\begin{array}{c}\text{(13)}& △u=\sqrt{2f\left(x,y\right)+u_{xx}^2+u_{yy}^2+2u_{xy}^2},\hspace{1em}\text{in\hspace{0.17em}}\Omega ,u=g\left(x,y\right),\hspace{1em}\text{in\hspace{0.17em}}\partial \Omega .\end{array}$$

This is still a fully nonlinear partial differential equation although the nonlinear terms are moved to the right side. Benamou et al. [13] have considered an iterative finite difference method for the above equation. But it will become very difficult to construct the approximation of $u_{xx},u_{yy}$, and $u_{xy}$ on unstructured mesh. An efficient iterative algorithm based on kernel-based meshfree discretization was proposed in [14]. The cascadic meshfree method is as follows.

By observing the above calculation process, it is easy to find that algorithm 1 avoids tedious interpolation and makes higher derivative functions $u_{xx},u_{yy}$, and $u_{xy}$ be easy to compute. Compared with the traditional cascadic multigrid method [15], the cascadic meshfree method requires neither domain nor boundary discretization; i.e., it is truly meshfree. It is flexible to deal with physical domain with curved boundary.

Algorithm 1: Cascadic meshfree method.

 

Input: right-hand sides f and $g$; number of levels n.

 

Output: numerical solution $u_n\in U_n$.

 

Set $f_0=2f$.

 

for $j=1,\dots ,n$ do

 

 (1) Determine $s_j\in U_j$ to $f_{j-1}$ and $g$ with

 

    $F{s}_j\left(\mathbf{y}\right)=\sqrt{f_{j-1}\left(\mathbf{y}\right)},\mathbf{y}\in Y_j^I,$

 

    $s_j\left(\mathbf{y}\right)=g\left(\mathbf{y}\right),\mathbf{y}\in Y_j^B.$

 

 (2) Update the solution and the right side:

 

    $u_j=s_j,$

 

    $f_j=f_0+{\left(u_j\right)}_{xx}^2+{\left(u_j\right)}_{yy}^2+2{\left(u_j\right)}_{xy}^2.$

 

end for

2.2. Stationary Multilevel Method

With a stationary multiscale algorithm [16], the condition number of the discrete matrix can be relatively small, and the computation can be performed in $O\left(N\right)$ operations.

In this algorithm, the present problems (3) and (4) are 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 be 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. And the final numerical solution is the sum of all of approximation. For interpolation problems, the linear convergence order has been proved in Sobolev spaces on the sphere in [17] and on bounded domains in [18]. And the applications of this algorithm in solving linear partial differential equations on spheres were proposed in [19] and on bounded domains were proposed in [20–22]. A detailed discussion on the stationary multilevel algorithm is reviewed in the latest papers [23, 24]. Algorithm 2 is the specific form for nonlinear situation.

Algorithm 2: Stationary multilevel RBF collocation.

 

Input: right-hand sides f and $g$; number of levels n.

 

Output: numerical solution $u_n\in U=U_1+\cdots +U_n$.

 

Set $u_0=0$, $f_0=f$, and $g_0=g$.

 

for $j=1,\dots ,n$ do

 

 (1) Determine $s_j\in U_j$ to $f_{j-1}$ and $g_{j-1}$ with

 

    $F{s}_j\left(\mathbf{y}\right)=f_{j-1}\left(\mathbf{y}\right),\mathbf{y}\in Y_j^I,$

 

    $s_j\left(\mathbf{y}\right)=g_{j-1}\left(\mathbf{y}\right),\mathbf{y}\in Y_j^B.$

 

 (2) Update the solution and the residuals:

 

    $u_j=u_{j-1}+s_j,$

 

    $f_j=f_{j-1}-F{s}_j,$

 

    $g_j=g_{j-1}-s_j.$

 

end for

2.3. Hierarchical Radial Basis Function Method

In this subsection, we consider a very simple algorithm for solving Monge–Ampère equations (3) and (4). If we choose to use the trial spaces $V_j$ to construct kernel-based collocation, then we will obtain a single-level Algorithm 3.

Algorithm 3: Hierarchical radial basis function method.

 

Input: right-hand sides f and $g$; Number of levels n.

 

Output: numerical solution $u_n\in V_n$.

 

for $j=1,\dots ,n$ do

 

 (1) Determine $s_j\in V_j$ to f and $g$ with

 

    $F{s}_j\left(\mathbf{y}\right)=f\left(\mathbf{y}\right),\mathbf{y}\in Y_j^I,$

 

    $s_j\left(\mathbf{y}\right)=g\left(\mathbf{y}\right),\mathbf{y}\in Y_j^B.$

 

 (2) Update the solution:

 

    $u_j=s_j.$

 

end for

Algorithm 3 is called the hierarchical radial basis function method because it is similar to multilevel splitting of finite element spaces [25]. In fact, $V_j$ are obtained by splitting the reproducing Hilbert space (or native space of the corresponding strictly positive define radial basic function). We consider the Matlab implementation of Algorithm 3. Let the radial basis function matrix be $B_j$ on the j level, which has the form$$\begin{array}{}\text{(14)}& B_j={\left[\begin{array}{cccc}{\Phi }_1\left({\mathbf{y}}_1,{\mathbf{x}}_1^1\right)& {\Phi }_1\left({\mathbf{y}}_1,{\mathbf{x}}_1^2\right)& \cdots & {\Phi }_j\left({\mathbf{y}}_1,{\mathbf{x}}_j^{M_j}\right)\\ {\Phi }_j\left({\mathbf{y}}_1,{\mathbf{x}}_j^{M_j}\right)& {\Phi }_1\left({\mathbf{y}}_2,{\mathbf{x}}_1^2\right)& \cdots & {\Phi }_j\left({\mathbf{y}}_2,{\mathbf{x}}_j^{M_j}\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\Phi }_1\left({\mathbf{y}}_{N_I},{\mathbf{x}}_1^1\right)& {\Phi }_1\left({\mathbf{y}}_{N_I},{\mathbf{x}}_1^2\right)& \cdots & {\Phi }_j\left({\mathbf{y}}_{N_I},{\mathbf{x}}_j^{M_j}\right)\\ {\Phi }_1\left({\mathbf{y}}_{N_I+1},{\mathbf{x}}_1^1\right)& {\Phi }_1\left({\mathbf{y}}_{N_I+1},{\mathbf{x}}_1^2\right)& \cdots & {\Phi }_j\left({\mathbf{y}}_{N_I+1},{\mathbf{x}}_j^{M_j}\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\Phi }_1\left({\mathbf{y}}_N,{\mathbf{x}}_1^1\right)& {\Phi }_1\left({\mathbf{y}}_N,{\mathbf{x}}_1^2\right)& \cdots & {\Phi }_j\left({\mathbf{y}}_N,{\mathbf{x}}_j^{M_j}\right)\end{array}\right]}_{N\times M}.\end{array}$$

Here, ${\mathbf{y}}_1,\dots ,{\mathbf{y}}_{N_I}$ are interior collocation data and ${\mathbf{y}}_{N_I+1},\dots ,{\mathbf{y}}_N$ are boundary collocation data in $\overline{\Omega }\subset {ℝ}^2$. Let ${\left(B_j\right)}_{xx},{\left(B_j\right)}_{yy}$, and ${\left(B_j\right)}_{xy}$ be some matrices formed by the derivatives of the corresponding functions in $B_j$, respectively. Then, the approximate values of u, $u_{xx},u_{yy}$, and $u_{xy}$ on the data $\left\{{\mathbf{y}}_1,\dots ,{\mathbf{y}}_N\right\}$ can be represented by the column vectors as$$\begin{array}{c}\text{(15)}& \tilde{\mathbf{u}}=B_j\mathbf{c},{\tilde{\mathbf{u}}}_{\mathbf{xx}}={\left(B_j\right)}_{xx}\mathbf{c},\\ {\tilde{\mathbf{u}}}_{\mathbf{yy}}={\left(B_j\right)}_{yy}\mathbf{c},\\ {\tilde{\mathbf{u}}}_{\mathbf{xy}}={\left(B_j\right)}_{xy}\mathbf{c},\end{array}$$with $\mathbf{c}={\left[c_1,\dots ,c_M\right]}^T$, respectively. And the nonlinear discrete system in Algorithm 3 can be written as$$\begin{array}{c}\text{(16)}& {\left({\tilde{\mathbf{u}}}_{\mathbf{xx}}\right)}_k{\left({\tilde{\mathbf{u}}}_{\mathbf{yy}}\right)}_k-{\left({\tilde{\mathbf{u}}}_{\mathbf{xy}}\right)}_k^2-f\left({\mathbf{y}}_k\right)=0,\hspace{1em}\text{for\hspace{0.17em}}k=1,\dots ,N_I,{\left(\tilde{\mathbf{u}}\right)}_k-g\left({\mathbf{y}}_k\right)=0,\hspace{1em}\text{for\hspace{0.17em}}k=N_I+1,\dots ,N.\end{array}$$

In the next numerical experiments, the linear discrete system will be solved by the least square method, which is a “$\setminus$” operation in Matlab. And the Matlab function lsqnonlin will be used for solving nonlinear discrete system. Computations were performed on a laptop with 2.4 GHz Intel Core i7 processor, using Matlab running in the Windows 7 operating system.

3. Numerical Examples

In this section, we discuss the implementation of Algorithms 1–3. We consider the multiscale collocation methods for Monge–Ampère equations (3) and (4) with the true solution:$$\begin{array}{}\text{(17)}& u=\exp \left(\frac{\left(x^2+y^2\right)}{2}\right).\end{array}$$

Before the nonlinear solver starts working in Algorithms 2 and 3, a reasonable choice for the initial data is solving the Laplace equation $△u=\sqrt{2f}$. This method of selecting initial values is natural, which has been used in [13, 14].

3.1. Square Domain

First, we generate a nested Halton points $\mathrm{9,25,81,289},\text{\hspace{0.17em}and\hspace{0.17em}}1089$ in the interior of the regular domain $\Omega ={\left[\mathrm{0,1}\right]}^2$ and create additional $4\left(\sqrt{N}-1\right)$ equally spaced collocation points for the boundary conditions on each level. So, the boundary data are $\mathrm{8,16,32,64},\text{\hspace{0.17em}and\hspace{0.17em}}128$ successively. In our experiments, Wendland’s $C^6$ functions carried the parameters ${\epsilon }_j$ used to construct trial spaces $V_j$ and $U_j$. $M=40\times 40$ equally spaced evaluation points are used to compute$$\begin{array}{}\text{(18)}& \text{RMS}-\text{error}=\sqrt{\frac{1}{M}{\displaystyle \sum }_{k=1}^M{\left[u\left({\xi }_k\right)-u_j\left({\xi }_k\right)\right]}^2}=\frac{1}{\sqrt{M}}{‖u-u_j‖}_2.\end{array}$$

Table 2 shows three kinds of strategies for choosing trial data and testing data. A detailed explanation is as follows:

 

TT1: trial data and testing data are equal

 

TT2: trial data and testing data are equal only on boundary, while additional collocation points are created in the internal domain

 

TT3: trial data only include interior points, while collocation is implemented in internal domain and the boundary

Table 2

Three kinds of strategies for arranging trial and collocation data.

TT1		TT2		TT3
Trial	Test	Trial	Test	Trial	Test
9 + 8	9 + 8	9 + 8	25 + 8	9	9 + 8
25 + 16	25 + 16	25 + 16	81 + 16	25	25 + 16
81 + 32	81 + 32	81 + 32	289 + 32	81	81 + 32
289 + 64	289 + 64	289 + 64	1089 + 64	289	289 + 64

Let ${\epsilon }_j=2^{j-1}\epsilon$ be scaled parameter defined by (8) on the j-level. Then, for different values of ε, the numerical results (including RMS errors and convergence rates) are displayed in Tables 3 to 5. It should be noted that Algorithms 1 and 2 are terrible and nonconvergent if the centers just being placed in the interior of the domain. This is the reason for the lack of corresponding TT3 experiments in Tables 3 to 4. But in Algorithm 3, we let the centers only include Halton data while the collocation sites contain the centers and additional data on the boundary. This is because the boundary data are not nested.

Table 3

Cascadic meshfree method (square domain).

TT1$\epsilon =0.25$	Rate	TT2$\epsilon =0.25$	Rate
6.093110e − 02		3.075884e − 01
1.499562e − 02	2.022638	6.131756e − 01	−0.995299
6.510642e − 03	1.203669	1.470053e + 00	−1.261496
9.034329e − 03	−0.472618	1.876969e + 00	−0.352537

Table 4

Stationary multilevel RBF collocation (square domain).

TT1$\epsilon =0.25$	Rate	TT2$\epsilon =0.25$	Rate
1.760832e − 02		8.699399e − 01
5.185207e − 02	−1.558144	2.448284e − 01	1.829145
1.910704e − 01	−1.881631	1.454338e − 01	0.751409
2.965791e − 01	−0.634312	1.536339e − 01	−0.079134

Table 5

Hierarchical radial basis function method (square domain).

TT3$\epsilon =0.5$	Rate	TT3$\epsilon =0.25$	Rate
2.187937e + 00		1.967965e + 00
7.361101e − 01	1.571578	8.919338e − 02	4.463624
1.622483e − 01	2.181718	8.846755e − 03	3.333716
8.991480e − 02	0.851573	2.202111e − 04	5.328189

By TT1 strategy and with an initial value $\epsilon =0.25$, the cascadic meshfree method is convergent although it seems that the convergence ceases at a later stage. But when implementing the algorithm by TT2 strategy with the same initial value, the cascadic meshfree method is nonconvergent. This is not contradictory to the linear and nonlinear theory of Schaback [26] and Böhmer and Schaback [6, 27]. The theory holds that a testing with many more degrees of freedom is necessary because the small trial space must be tested on a fine-grained space discretization. But TT2 strategy does not bring a good approximation on the beginning level. This makes $\sqrt{2f+{\left(u_1\right)}_{xx}^2+{\left(u_1\right)}_{yy}^2+2{\left(u_1\right)}_{xy}^2}$ be no longer an ideal right side on the next level.

The same set of experiments as for the cascadic meshfree method is displayed in Table 4 for the stationary multilevel RBF collocation method. Different from the numerical results of the stationary multilevel method for solving linear elliptic problems (see Table 41.4 in Fasshauer’s book [28]), Algorithm 2 with an initial $\epsilon =0.25$ is nonconvergent. It seems that the TT2 strategy is more numerically stable than the TT1 strategy, but with limited accuracy. Maybe a more reasonable parameter scheme should not be ${\epsilon }_j=2^{j-1}\epsilon$, when using the stationary multilevel RBF collocation method to solve fully nonlinear partial differential equation. On lack of the theoretical guidance, we have not obtained some satisfactory numerical results.

In Table 5, we list RMS errors and convergence rate for the hierarchical radial basis function collocation method. We observe that this method has an ideal convergence behavior. Even to a relatively large initial parameter $\epsilon =0.5$, the convergence rate of this method is close to 1.5. With an initial parameter $\epsilon =0.25$, the RMS error on 4 level is remarkably reduced to ${10}^{-4}$. Algorithm 3 is a single-level collocation method but has a multiscale trial space. So there is no relationship between the former numerical solution and the latter one. We do not obtain an ideal numerical result on 1 level where only 9 centers are arranged in the interior of the domain. In addition, these 9 centers are far from the boundary. But this situation has been significantly improved after the addition of the added 16 centers. Refer Figure 1 for the initial scattered centers and its refinement. Clearly, compared with the first two algorithms, Algorithm 3 can solve the model problem with a higher accuracy and lower computational cost.

[figures omitted; refer to PDF]

In Figure 2, we show plots of the RMS errors of the above three algorithms. It shows the exponential decay of the error by the hierarchical radial basis function collocation.

[figure omitted; refer to PDF]

3.2. Sector Domain

In this subsection, we consider numerically solving the problems (3) and (4) in a sector domain (see Figure 3) but with the same true solution as in previous example. We first generate a set of successive Halton points: $N=\mathrm{9,25,81,289,1089}$. We place some equidistant points on the rectangular boundary and the curved boundary, respectively, and make the total number of boundary points be still equal to $4\left(\sqrt{N}-1\right)$. The first $M=100$ interior points are used to compute the RMS error, and the same ${\epsilon }_j=2^{j-1}\epsilon$ is used for j-level experiments. As in example 1, Wendland’s $C^6$ functions carried the parameters ${\epsilon }_j$ used to construct trial spaces $V_j$ and $U_j$. The numerical results are displayed in Tables 6–8 and Figure 4.

[figures omitted; refer to PDF]

Table 6

Cascadic meshfree method (sector domain).

TT1$\epsilon =0.25$	Rate	TT2$\epsilon =0.25$	Rate
1.221624e − 02		3.332465e − 02
1.853373e − 03	2.720575	7.971312e − 02	−1.258228
1.380550e − 03	0.424910	3.077875e − 01	−1.949046
2.619970e − 03	−0.924307	4.410556e − 01	−0.519026

Table 7

Stationary multilevel RBF collocation (sector domain).

TT1$\epsilon =0.25$	Rate	TT2$\epsilon =0.25$	Rate
1.192816e − 01		4.395848e − 01
3.150129e − 01	−1.401039	1.605481e − 01	1.453136
3.808680e − 01	−0.273880	9.773937e − 02	0.715994
4.758944e − 01	−0.321350	8.697644e − 02	0.168315

Table 8

Hierarchical radial basis function method (sector domain).

TT3$\epsilon =0.5$	Rate	TT3$\epsilon =0.25$	Rate
1.980441e + 00		4.540010e − 01
4.899895e − 01	2.014999	1.935033e − 02	4.552266
9.722933e − 02	2.333287	5.096518e − 04	5.246702
3.016228e − 02	1.688646	5.384921e − 05	3.242515

[figure omitted; refer to PDF]

From numerical results, we observe that TT2 testing is not an ideal strategy for the cascadic meshfree method. And given a relatively small initial ε stationary multilevel method will be convergent although it becomes very slow at a later stage. But the performance of hierarchical radial basis function method is different from first two kinds of algorithms. With an initial parameter $\epsilon =0.25$, the RMS error on the 4-level is reduced to ${10}^{-5}$. Figure 4 clearly shows the superiority of the hierarchical radial basis function method.

4. Conclusions

The paper considered three kinds of multiscale RBF collocation methods for solving the fully nonlinear Monge–Ampère equation. By a specific example, we found that the hierarchical radial basis function collocation method has higher accuracy and lower computational complexity. However, a convergence proof for the hierarchical radial basis function collocation method is missing. This will depend on the approximation of trial spaces $V_j$, new inverse inequality, and sampling theorem. Based on the convergence theory of Böhmer and Schaback [6, 27], we will consider the following theoretical aspects in future work [29]:

(1)

Trial Approximation. We will consider the approximation of spaces $V_j$ under some Sobolev norms. That is to say, with maps $I_j:U⟶V_j$, we will prove error bounds:

$$\begin{array}{}\text{(19)}& {‖u-I_j\left(u\right)‖}_{W^{m,\mu }\left(\Omega \right)}\le C{r}_j^{\alpha }{‖u‖}_U,\hspace{1em}\text{for\hspace{0.17em}all\hspace{0.17em}}u\in U,\end{array}$$

 

where $r_j$ is the trial discretization parameter, α is the order of approximation, and U usually is the reproducing Hilbert space which satisfies additional regularity requirement.

(2)

Stability of Testing Strategy. We will prove the boundedness and coerciveness of nonlinear testing operator. This will depend on some Poincaré-type inequality and Bernstein-type inequality. The detailed discussion will be given in [29].

Acknowledgments

This paper used some Halton datasets and drew lessons from partial codes from Fasshauer’s book [28]. The authors are grateful to [28] for its free CD. The research of the first author was partially supported by the Natural Science Foundations of China (no. 11501313), the Natural Science Foundations of Ningxia Province (no. 2019AAC02001), 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 (no. TJGC2019012).

Appendix

Appendix

Following are the Matlab codes of the hierarchical radial basis function method for Monge–Ampère equation. The DistanceMatrixCSRBF function and 2D Halton data can be cited from [28]:

(1)

% Main routine

(2)

rbf = @(e,r) r. ^∧ 8. ∗ (66 ∗ spones(r)–…

(3)

 154 ∗ r + 121 ∗ r. ^∧ 2–32 ∗ r. ^∧ 3);

(4)

Lrbf = @(e,r) 44 ∗ e ^∧ 2 ∗ r. ^∧ 6 ∗ (84 ∗ spones(r)−…

(5)

 264 ∗ r + 267 ∗ r. ^∧ 2–88 ∗ r. ^∧ 3);

(6)

% exact solution and its Laplacian

(7)

u = @(x,y) exp((x. ^∧ 2 + y. ^∧ 2)/2);

(8)

Fu = @(x,y) (1 + x. ^∧ 2 + y. ^∧ 2). ∗ exp(x. ^∧ 2 + y. ^∧ 2);

(9)

K = 4; neval = 40; gridtype = ’h’;

(10)

N = (2 ^∧ K + 1) ^∧ 2;

(11)

name = sprintf(’Data2D_%d%s’ ,N,gridtype);

(12)

load(name);

(13)

xdata = dsites(:,1); ydata = dsites(:,2);

(14)

ep = 0.5 ∗ 2. ^∧ [ 0:K–1];

(15)

% create neval–by–neval equally spaced points

(16)

% in the unit square

(17)

grid = linspace(0,1,neval); [xe,ye] = meshgrid(grid);

(18)

epoints = [xe(:) ye(:)];

(19)

% compute exact solution

(20)

exact = u(epoints(:,1),epoints(:,2));

(21)

N1 = [0 9 25 81 289 1089 4225];

(22)

for k = 1:K

(23)

 N2 = (2 ^∧ {k} + 1) ^∧ 2;

(24)

 intdata = [xdata(1:N2) ydata(1:N2)];

(25)

 % add boundary points

(26)

 sn = sqrt(N2); bdylin = linspace(0,1,sn)’;

(27)

 bdy0 = zeros(sn–1,1); bdy1 = ones(sn–1,1);

(28)

 bdydata = [bdylin(1:end–1) bdy0;…

(29)

  bdy1 bdylin(1:end–1); …

(30)

  flipud(bdylin(2:end)) bdy1;…

(31)

  bdy0 flipud(bdylin(2:end))];

(32)

 rhs1 = Fu(intdata(:,1),intdata(:,2));

(33)

 rhs2 = u(bdydata(:,1),bdydata(:,2));

(34)

 % solve linear system on k–level

(35)

 c{k} = linfh(k,N1,ep,xdata, ydata,…

(36)

  intdata, bdydata,rhs1,rhs2);

(37)

 % solve nonlinear system on k–level

(38)

 coef{k} = lsqnonlin(@(c)nonfh(c,k,N1,ep,xdata,…

(39)

  ydata, intdata,bdydata,rhs1,rhs2),c{k});

(40)

 % assemble evaluation matrix

(41)

 for j = 1:k

(42)

  xctrs{j} = [xdata(N1(j) + 1:N1(j + 1))…

(43)

   ydata(N1(j) + 1:N1(j + 1))];

(44)

  DM_eval = DistanceMatrixCSRBF(epoints,…

(45)

   xctrs{j},ep(j));

(46)

  EMK{j} = rbf(ep(j),DM_eval);

(47)

 end

(48)

 EM = EMK{1};

(49)

 if(k > 1)

(50)

  for j = 2:k

(51)

   EM = [EM EMK{j}];

(52)

  end

(53)

 end

(54)

 % evaluate RBF approximation

(55)

  Pf = EM ∗ coef{k};

(56)

 % compute maximum error on evaluation grid

(57)

  rms_err = norm(Pf–exact)/neval;

(58)

  if (k > 1)

(59)

   rms_rate = log(rms_err_old/rms_err)/log(2);

(60)

  end

(61)

 rms_err_old = rms_err;

(62)

end

(63)

%——————————–

(64)

% linear solution function

(65)

function c = linfh(k,N1,ep,xdata, ydata,idata,…

(66)

  bdata,rhs1,rhs2)

(67)

 rbf = @(e,r) r. ^∧ 8. ∗ (66 ∗ spones(r)–154 ∗ r + 121 ∗ …

(68)

  r.2 ^∧ –32 ∗ r.3 ^∧);

(69)

 dxxrbf = @(e,r,dx) 528 ∗ e ^∧ 4 ∗ dx. ^∧ 2. ∗ r. ^∧ 6. ∗ (7 ∗ spones(r)…

(70)

  –6 ∗ r)–22 ∗ e ^∧ 2 ∗ r. ^∧ 7. ∗ (24 ∗ spones(r)…

(71)

  –39 ∗ r + 16 ∗ r. ^∧ 2);

(72)

 dyyrbf = @(e,r,dy) 528 ∗ e ^∧ 4 ∗ dy. ^∧ 2. ∗ r. ^∧ 6. ∗ (7 ∗ spones(r)…

(73)

  –6 ∗ r)–22 ∗ e ^∧ 2 ∗ r. ^∧ 7. ∗ (24 ∗ spones(r)…

(74)

  –39 ∗ r + 16 ∗ r. ^∧ 2);

(75)

 dxyrbf = @(e,r,dx,dy) 528 ∗ e ^∧ 4 ∗ dx. ∗ dy. ∗ r. ^∧ 6. ∗ …

(76)

  (7 ∗ spones(r)–6 ∗ r);

(77)

 Lrbf = @(e,r) 44 ∗ e ^∧ 2 ∗ r. ^∧ 6. ∗ (84 ∗ spones(r)–264 ∗ r + …

(78)

  267 ∗ r. ^∧ 2–88 ∗ r. ^∧ 3);

(79)

 % build trial data

(80)

 for j = 1:k

(81)

  xctrs{j} = [xdata(N1(j) + 1:N1(j + 1))…

(82)

   ydata(N1(j) + 1:N1(j + 1))];

(83)

  cent = xctrs{j};

(84)

  D1 = DistanceMatrixCSRBF(idata, cent, ep(j));

(85)

  D2 = DistanceMatrixCSRBF(bdata, cent, ep(j));

(86)

  CMK{j} = Lrbf(ep(j),D1);

(87)

  BMK{j} = rbf(ep(j),D2);

(88)

 end

(89)

 CM = CMK{1};

(90)

 BM = BMK{1};

(91)

 if(k$>$1)

(92)

  for j = 2:k

(93)

   CM = [CM CMK{j}];

(94)

   BM = [BM BMK{j}];

(95)

  end

(96)

 end

(97)

 LCM = [CM; BM];

(98)

 Lrhs = [sqrt(2 ∗ rhs1); rhs2];

(99)

 c = LCM\Lrhs;

References

[1] S. C. Brenner, T. Gudi, M. Neilan, L.-Y. Sung, "C 0 penalty methods for the fully nonlinear Monge-Ampère equation," Mathematics of Computation, vol. 80 no. 276,DOI: 10.1090/s0025-5718-2011-02487-7, 2011.

[2] X. Feng, M. Neilan, "Mixed finite element methods for the fully nonlinear monge-ampère equation based on the vanishing moment method," SIAM Journal on Numerical Analysis, vol. 47 no. 2, pp. 1226-1250, DOI: 10.1137/070710378, 2009.

[3] A. M. Oberman, "Wide stencil finite difference schemes for the elliptic Monge-Ampère equations and functions of the eigenvalues of the Hessian," Discrete and Continuous Dynamical Systems—Series B, vol. 10 no. 1, pp. 221-238, 2008.

[4] Z. Liu, Y. He, "Solving the elliptic Monge-Ampère equation by Kansa’s method," Engineering Analysis with Boundary Elements, vol. 37 no. 1, pp. 84-88, DOI: 10.1016/j.enganabound.2012.09.004, 2013.

[5] Z. Liu, Y. He, "An iterative meshfree method for the elliptic monge-ampère equation in 2D," Numerical Methods for Partial Differential Equations, vol. 30 no. 5, pp. 1507-1517, DOI: 10.1002/num.21849, 2014.

[6] K. Böhmer, R. Schaback, "A meshfree method for solving the Monge-Ampère equation," Numerical Algorithms, vol. 82 no. 2, pp. 539-551, DOI: 10.1007/s11075-018-0612-1, 2019.

[7] J. Rashidinia, G. E. Fasshauer, M. Khasi, "A stable method for the evaluation of Gaussian radial basis function solutions of interpolation and collocation problems," Computers & Mathematics with Applications, vol. 72 no. 1, pp. 178-193, DOI: 10.1016/j.camwa.2016.04.048, 2016.

[8] J. Rashidinia, M. Khasi, G. E. Fasshauer, "A stable Gaussian radial basis function method for solving nonlinear unsteady convection-diffusion-reaction equations," Computers & Mathematics with Applications, vol. 75 no. 5, pp. 1831-1850, DOI: 10.1016/j.camwa.2017.12.007, 2018.

[9] J. Rashidinia, M. Khasi, "Stable Gaussian radial basis function method for solving Helmholtz equations," Computational Methods for Differential Equations, vol. 7, pp. 138-151, 2019.

[10] H. Wendland, "Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree," Advances in Computational Mathematics, vol. 4 no. 1, pp. 389-396, DOI: 10.1007/bf02123482, 1995.

[11] Z. Wu, "Compactly supported positive definite radial functions," Advances in Computational Mathematics, vol. 4 no. 1, pp. 283-292, DOI: 10.1007/bf03177517, 1995.

[12] Q. Xu, Z. Liu, "Scattered data interpolation and approximation with truncated exponential radial basis function," Mathematics, vol. 7 no. 11,DOI: 10.3390/math7111101, 2019.

[13] J.-D. Benamou, B. D. Froese, A. M. Oberman, "Two numerical methods for the elliptic Monge-Ampère equation," ESAIM: Mathematical Modelling and Numerical Analysis, vol. 44 no. 4, pp. 737-758, DOI: 10.1051/m2an/2010017, 2010.

[14] Z. Liu, Y. He, "Cascadic meshfree method for the elliptic Monge-Ampère equation," Engineering Analysis with Boundary Elements, vol. 37 no. 7-8, pp. 990-996, DOI: 10.1016/j.enganabound.2013.04.005, 2013.

[15] F. A. Bornemann, P. Deuflhard, "The cascadic multigrid method for elliptic problems," Numerische Mathematik, vol. 75 no. 2, pp. 135-152, DOI: 10.1007/s002110050234, 1996.

[16] M. S. Floater, A. Iske, "Multistep scattered data interpolation using compactly supported radial basis functions," Journal of Computational and Applied Mathematics, vol. 73 no. 1-2, pp. 65-78, DOI: 10.1016/0377-0427(96)00035-0, 1996.

[17] Q. T. L. Gia, I. H. Sloan, H. Wendland, "Multiscale analysis in Sobolev spaces on the sphere," SIAM Journal on Numerical Analysis, vol. 48 no. 6, pp. 2065-2090, DOI: 10.1137/090774550, 2010.

[18] H. Wendland, "Multiscale analysis in Sobolev spaces on bounded domains," Numerische Mathematik, vol. 116 no. 3, pp. 493-517, DOI: 10.1007/s00211-010-0313-8, 2010.

[19] Q. T. Le Gia, I. H. Sloan, H. Wendland, "Multiscale RBF collocation for solving PDEs on spheres," Numerische Mathematik, vol. 121 no. 1, pp. 99-125, DOI: 10.1007/s00211-011-0428-6, 2012.

[20] A. Chernih, Q. T. Le Gia, "Multiscale methods with compactly supported radial basis functions for Galerkin approximation of elliptic PDEs," IMA Journal of Numerical Analysis, vol. 34 no. 2, pp. 569-591, DOI: 10.1093/imanum/drt004, 2014.

[21] A. Chernih, Q. T. Le Gia, "Multiscale methods with compactly supported radial basis functions for the Stokes problem on bounded domains," Advances in Computational Mathematics, vol. 42 no. 5, pp. 1187-1208, DOI: 10.1007/s10444-016-9458-z, 2016.

[22] P. Farrell, H. Wendland, "RBF multiscale collocation for second order elliptic boundary value problems," SIAM Journal on Numerical Analysis, vol. 51 no. 4, pp. 2403-2425, DOI: 10.1137/120898383, 2013.

[23] H. Wendland, "Multiscale radial basis functions," Frames and Other Bases in Abstract and Function Spaces–Novel Methods in Harmonic Analysis, vol. 1, pp. 265-299, 2017.

[24] H. Wendland, "Solving partial differential equations with multiscale radial basis functions," Contemporary Computational Mathematics–A Celebration of the 80th Birthday of Ian Sloan, pp. 1191-1213, DOI: 10.1007/978-3-319-72456-0_55, 2018.

[25] H. Yserentant, "On the multi-level splitting of finite element spaces," Numerische Mathematik, vol. 49 no. 4, pp. 379-412, DOI: 10.1007/bf01389538, 1986.

[26] R. Schaback, "Unsymmetric meshless methods for operator equations," Numerische Mathematik, vol. 114 no. 4, pp. 629-651, DOI: 10.1007/s00211-009-0265-z, 2010.

[27] K. Böhmer, R. Schaback, "A nonlinear discretization theory," Journal of Computational and Applied Mathematics, vol. 254, pp. 204-219, DOI: 10.1016/j.cam.2013.03.029, 2013.

[28] G. E. Fasshauer, Meshfree Approximation Methods with MATLAB, 2007.

[29] Z. Liu, Q. Xu, "Convergence of hierarchical radial basis function collocation method," , . In preparation