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

Authors: Zhiyong Liu; Qiuyan Xu
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 1. Introduction

The Monge–Ampère equation with Dirichlet boundary is a fully nonlinear partial differential equation, which is given by(1)detD2uz=fz,inΩ,(2)uz=gz,in∂Ω,where z has n independent variables in a bounded domain Ω⊂ℝd and D2u is the Hessian of the function u. When restricting it to Ω⊂ℝ2, we can rewrite the equation as(3)Fu=f,inΩ,(4)u=g,in∂Ω,with a nonlinear map F:u↦uxxuyy−uxy2.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. Section3 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

LetX=x1,…,xN be a finite point set in Ω⊆ℝ2. We define a commonly used trial discretization parameter:(5)hX,Ω≔supx∈Ωminxj∈Xx−xj2,which can be regarded as the radius of the largest empty ball that can be placed among the data sites xj. Generate nested point sets(6)X1⊂X2⊂⋯⊂Xj⊂⋯,which have trial discretization parameters hj=hXj,Ω. Let Xj∗ be newly added point set in Xj, then we have X1=X1∗ and Xj=Xj−1∪Xj∗ for j=2,…,.Given a kind of the compactly supported radial functionsΦ:ℝ2⟶ℝ, we can rescale it (in the translation-invariant kernel-based case) as(7)Φj⋅,y=εjdΦεj⋅−y,y∈Xj∗ory∈Xj.In order to make the support radii ofΦj become smaller and smaller with the addition of new points, we will select(8)εj=b1hj,with a given constant0<b≤1.LetVj=Vj−1+Wj, with V1=W1 and Wj=spanΦj⋅,y:y∈Xj∗ for j=1,…,. Then, we can build two types of trial spaces with the form(9)Vj=spanΦ1x,x11,…,Φ1x,x1M1︸,…,Φjx,xj1,…,Φjx,xjMj︸,Uj=spanΦjx,x11,…,Φjx,x1M1︸,…,Φjx,xj1,…,Φjx,xjMj︸.Here, we have writtenXk∗=xk1,…,xkMk for k=1,2,…,j. Vj and Uj have multiscale expressions, which differ from the traditional radial basis trial spaces(10)Qj=spanΦx,x11,…,Φx,x1M1,…,Φx,xj1,…,Φx,xjMj.LetM=∑k=1jMk be the total number of the centers and Bx,xkk=1M be basis functions from one of the trial spaces above. Then, with a trial function vx=∑k=1MckBx,xk, the nonsymmetric kernel-based collocation method for equations (3) and (4) will generate the following nonlinear system(11)Fvy=fy,y∈YI,vy=gy,y∈YB,for some given interior collocation sites YI and boundary collocation sites YB and some unknown coefficients ck. We need to compute the derivatives of Bx,xk 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ϕr=ϕx with variables x=x,y∈ℝ2.
Partial derivatives
First-order
∂
/
∂
x
ϕ
r
=
x
/
r
d
/
d
r
ϕ
r
∂
/
∂
y
ϕ
r
=
y
/
r
d
/
d
r
ϕ
r
∂
/
∂
n
ϕ
r
=
x
,
n
/
r
d
/
d
r
ϕ
r
Second-order
∂
2
/
∂
x
2
ϕ
r
=
x
2
/
r
2
d
2
/
d
r
2
ϕ
r
+
y
2
/
r
3
d
/
d
r
ϕ
r
∂
2
/
∂
y
2
ϕ
r
=
y
2
/
r
2
d
2
/
d
r
2
ϕ
r
+
x
2
/
r
3
d
/
d
r
ϕ
r
∂
2
/
∂
x
∂
y
ϕ
r
=
x
y
/
r
2
d
2
/
d
r
2
ϕ
r
−
x
y
/
r
3
d
/
d
r
ϕ
rFor example, the Wendland’sC6 function ϕr=1−εr+832εr3+25εr2+8εr+1 with scaling parameter ε has second-order partial derivatives:(12)∂2∂x2ϕr=22ε21−εr+624ε2x26εr+1+16εr3−9εr2−6εr−1,∂2∂y2ϕr=22ε21−εr+624ε2y26εr+1+16εr3−9εr2−6εr−1,∂2∂x∂yϕr=528ε4xy1−εr+66εr+1.


### 2.1. Cascadic MeshFree Method

Once a convex solution be chosen, equations (3) and (4) can be rewritten in the equivalent form:(13)△u=2fx,y+uxx2+uyy2+2uxy2,inΩ,u=gx,y,in∂Ω.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 uxx,uyy, and uxy 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 algorithm1 avoids tedious interpolation and makes higher derivative functions uxx,uyy, and uxy 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 un∈Un.
Setf0=2f.
for j=1,…,n do
(1) Determinesj∈Uj to fj−1 and g with
Fsjy=fj−1y,y∈YjI,
sjy=gy,y∈YjB.
(2) Update the solution and the right side:
uj=sj,
fj=f0+ujxx2+ujyy2+2ujxy2.
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 ON 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 un∈U=U1+⋯+Un.
Setu0=0, f0=f, and g0=g.
for j=1,…,n do
(1) Determinesj∈Uj to fj−1 and gj−1 with
Fsjy=fj−1y,y∈YjI,
sjy=gj−1y,y∈YjB.
(2) Update the solution and the residuals:
uj=uj−1+sj,
fj=fj−1−Fsj,
gj=gj−1−sj.
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 Vj 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 un∈Vn.
for j=1,…,n do
(1) Determinesj∈Vj to f and g with
Fsjy=fy,y∈YjI,
sjy=gy,y∈YjB.
(2) Update the solution:
uj=sj.
end forAlgorithm3 is called the hierarchical radial basis function method because it is similar to multilevel splitting of finite element spaces [25]. In fact, Vj 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 Bj on the j level, which has the form(14)Bj=Φ1y1,x11Φ1y1,x12⋯Φjy1,xjMjΦjy1,xjMjΦ1y2,x12⋯Φjy2,xjMj⋮⋮⋱⋮Φ1yNI,x11Φ1yNI,x12⋯ΦjyNI,xjMjΦ1yNI+1,x11Φ1yNI+1,x12⋯ΦjyNI+1,xjMj⋮⋮⋱⋮Φ1yN,x11Φ1yN,x12⋯ΦjyN,xjMjN×M.Here,y1,…,yNI are interior collocation data and yNI+1,…,yN are boundary collocation data in Ω¯⊂ℝ2. Let Bjxx,Bjyy, and Bjxy be some matrices formed by the derivatives of the corresponding functions in Bj, respectively. Then, the approximate values of u, uxx,uyy, and uxy on the data y1,…,yN can be represented by the column vectors as(15)u˜=Bjc,u˜xx=Bjxxc,u˜yy=Bjyyc,u˜xy=Bjxyc,with c=c1,…,cMT, respectively. And the nonlinear discrete system in Algorithm 3 can be written as(16)u˜xxku˜yyk−u˜xyk2−fyk=0,fork=1,…,NI,u˜k−gyk=0,fork=NI+1,…,N.In the next numerical experiments, the linear discrete system will be solved by the least square method, which is a “∖” 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.


## 2.1. Cascadic MeshFree Method

Once a convex solution be chosen, equations (3) and (4) can be rewritten in the equivalent form:(13)△u=2fx,y+uxx2+uyy2+2uxy2,inΩ,u=gx,y,in∂Ω.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 uxx,uyy, and uxy 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 algorithm1 avoids tedious interpolation and makes higher derivative functions uxx,uyy, and uxy 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 un∈Un.
Setf0=2f.
for j=1,…,n do
(1) Determinesj∈Uj to fj−1 and g with
Fsjy=fj−1y,y∈YjI,
sjy=gy,y∈YjB.
(2) Update the solution and the right side:
uj=sj,
fj=f0+ujxx2+ujyy2+2ujxy2.
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 ON 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 un∈U=U1+⋯+Un.
Setu0=0, f0=f, and g0=g.
for j=1,…,n do
(1) Determinesj∈Uj to fj−1 and gj−1 with
Fsjy=fj−1y,y∈YjI,
sjy=gj−1y,y∈YjB.
(2) Update the solution and the residuals:
uj=uj−1+sj,
fj=fj−1−Fsj,
gj=gj−1−sj.
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 Vj 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 un∈Vn.
for j=1,…,n do
(1) Determinesj∈Vj to f and g with
Fsjy=fy,y∈YjI,
sjy=gy,y∈YjB.
(2) Update the solution:
uj=sj.
end forAlgorithm3 is called the hierarchical radial basis function method because it is similar to multilevel splitting of finite element spaces [25]. In fact, Vj 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 Bj on the j level, which has the form(14)Bj=Φ1y1,x11Φ1y1,x12⋯Φjy1,xjMjΦjy1,xjMjΦ1y2,x12⋯Φjy2,xjMj⋮⋮⋱⋮Φ1yNI,x11Φ1yNI,x12⋯ΦjyNI,xjMjΦ1yNI+1,x11Φ1yNI+1,x12⋯ΦjyNI+1,xjMj⋮⋮⋱⋮Φ1yN,x11Φ1yN,x12⋯ΦjyN,xjMjN×M.Here,y1,…,yNI are interior collocation data and yNI+1,…,yN are boundary collocation data in Ω¯⊂ℝ2. Let Bjxx,Bjyy, and Bjxy be some matrices formed by the derivatives of the corresponding functions in Bj, respectively. Then, the approximate values of u, uxx,uyy, and uxy on the data y1,…,yN can be represented by the column vectors as(15)u˜=Bjc,u˜xx=Bjxxc,u˜yy=Bjyyc,u˜xy=Bjxyc,with c=c1,…,cMT, respectively. And the nonlinear discrete system in Algorithm 3 can be written as(16)u˜xxku˜yyk−u˜xyk2−fyk=0,fork=1,…,NI,u˜k−gyk=0,fork=NI+1,…,N.In the next numerical experiments, the linear discrete system will be solved by the least square method, which is a “∖” 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 Algorithms1–3. We consider the multiscale collocation methods for Monge–Ampère equations (3) and (4) with the true solution:(17)u=expx2+y22.Before the nonlinear solver starts working in Algorithms2 and 3, a reasonable choice for the initial data is solving the Laplace equation △u=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 points9,25,81,289,and1089 in the interior of the regular domain Ω=0,12 and create additional 4N−1 equally spaced collocation points for the boundary conditions on each level. So, the boundary data are 8,16,32,64,and128 successively. In our experiments, Wendland’s C6 functions carried the parameters εj used to construct trial spaces Vj and Uj. M=40×40 equally spaced evaluation points are used to compute(18)RMS−error=1M∑k=1Muξk−ujξk2=1Mu−uj2.Table2 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 equalTT2: trial data and testing data are equal only on boundary, while additional collocation points are created in the internal domainTT3: trial data only include interior points, while collocation is implemented in internal domain and the boundaryTable 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 + 64Letεj=2j−1ε 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ε=0.25
Rate
TT2ε=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.352537Table 4
Stationary multilevel RBF collocation (square domain).
TT1ε=0.25
Rate
TT2ε=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.079134Table 5
Hierarchical radial basis function method (square domain).
TT3ε=0.5
Rate
TT3ε=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.328189By TT1 strategy and with an initial valueε=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 2f+u1xx2+u1yy2+2u1xy2 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 Table4 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 ε=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 εj=2j−1ε, 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 Table5, 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 ε=0.5, the convergence rate of this method is close to 1.5. With an initial parameter ε=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.Figure 1
Initial scattered centers and its refinement.
(a)
(b)In Figure2, 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 2
Errors of Algorithms1 to 3 for different testing strategies and initial ε.


### 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=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 4N−1. The first M=100 interior points are used to compute the RMS error, and the same εj=2j−1ε is used for j-level experiments. As in example 1, Wendland’s C6 functions carried the parameters εj used to construct trial spaces Vj and Uj. The numerical results are displayed in Tables 6–8 and Figure 4.Figure 3
Scattered centers and its refinement in sector domain.
(a)
(b)Table 6
Cascadic meshfree method (sector domain).
TT1ε=0.25
Rate
TT2ε=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.519026Table 7
Stationary multilevel RBF collocation (sector domain).
TT1ε=0.25
Rate
TT2ε=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.168315Table 8
Hierarchical radial basis function method (sector domain).
TT3ε=0.5
Rate
TT3ε=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.242515Figure 4
Errors of Algorithms1 to 3 for different testing strategies and initial ε for sector example.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 ε=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.


## 3.1. Square Domain

First, we generate a nested Halton points9,25,81,289,and1089 in the interior of the regular domain Ω=0,12 and create additional 4N−1 equally spaced collocation points for the boundary conditions on each level. So, the boundary data are 8,16,32,64,and128 successively. In our experiments, Wendland’s C6 functions carried the parameters εj used to construct trial spaces Vj and Uj. M=40×40 equally spaced evaluation points are used to compute(18)RMS−error=1M∑k=1Muξk−ujξk2=1Mu−uj2.Table2 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 equalTT2: trial data and testing data are equal only on boundary, while additional collocation points are created in the internal domainTT3: trial data only include interior points, while collocation is implemented in internal domain and the boundaryTable 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 + 64Letεj=2j−1ε 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ε=0.25
Rate
TT2ε=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.352537Table 4
Stationary multilevel RBF collocation (square domain).
TT1ε=0.25
Rate
TT2ε=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.079134Table 5
Hierarchical radial basis function method (square domain).
TT3ε=0.5
Rate
TT3ε=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.328189By TT1 strategy and with an initial valueε=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 2f+u1xx2+u1yy2+2u1xy2 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 Table4 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 ε=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 εj=2j−1ε, 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 Table5, 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 ε=0.5, the convergence rate of this method is close to 1.5. With an initial parameter ε=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.Figure 1
Initial scattered centers and its refinement.
(a)
(b)In Figure2, 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 2
Errors of Algorithms1 to 3 for different testing strategies and initial ε.


## 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=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 4N−1. The first M=100 interior points are used to compute the RMS error, and the same εj=2j−1ε is used for j-level experiments. As in example 1, Wendland’s C6 functions carried the parameters εj used to construct trial spaces Vj and Uj. The numerical results are displayed in Tables 6–8 and Figure 4.Figure 3
Scattered centers and its refinement in sector domain.
(a)
(b)Table 6
Cascadic meshfree method (sector domain).
TT1ε=0.25
Rate
TT2ε=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.519026Table 7
Stationary multilevel RBF collocation (sector domain).
TT1ε=0.25
Rate
TT2ε=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.168315Table 8
Hierarchical radial basis function method (sector domain).
TT3ε=0.5
Rate
TT3ε=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.242515Figure 4
Errors of Algorithms1 to 3 for different testing strategies and initial ε for sector example.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 ε=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 spacesVj, 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 Vj under some Sobolev norms. That is to say, with maps Ij:U⟶Vj, we will prove error bounds:(19)u−IjuWm,μΩ≤CrjαuU,for allu∈U,whererj 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].



---

*Source: 1748037-2020-03-19.xml*