# A Narrow-stencil finite difference method for approximating viscosity solutions of fully nonlinear elliptic partial differential equations with applications to Hamilton-Jacobi-

This paper presents a new narrow-stencil finite difference method for approximating the viscosity solution of second order fully nonlinear elliptic partial differential equations including Hamilton-Jacobi-Bellman equations. The proposed finite difference method naturally extends the Lax-Friedrichs method for first order problems to second order problems by introducing a key stabilization and guiding term called a "numerical moment". The numerical moment uses the difference of two (central) Hessian approximations to resolve the potential low-regularity of viscosity solutions. It is proved that the proposed scheme is well posed (i.e, it has a unique solution) and stable in both the l-2 norm and the l-infinity norm. The highlight of the paper is to prove the convergence of the proposed scheme to the viscosity solution of the underlying fully nonlinear second order problem using a novel discrete comparison argument. This paper extends the one-dimensional analogous method of Feng, Kao, and Lewis to the higher-dimensional setting. Numerical tests are presented to gauge the performance of the proposed finite difference methods and to validate the convergence result of the paper.

## Authors

• 11 publications
• 1 publication
03/12/2019

### Solving Nonlinear Parabolic Equations by a Strongly Implicit Finite-Difference Scheme

We discuss the numerical solution of nonlinear parabolic partial differe...
03/17/2021

### Convergent Finite Difference Methods for Fully Nonlinear Elliptic Equations in Three Dimensions

We introduce a generalized finite difference method for solving a large ...
03/25/2021

### A Nitsche Hybrid multiscale method with non-matching grids

We propose a Nitsche method for multiscale partial differential equation...
11/16/2020

### Convergence of Lasserre's hierarchy: the general case

Lasserre's moment-SOS hierarchy consists of approximating instances of t...
10/28/2019

### FD-Net with Auxiliary Time Steps: Fast Prediction of PDEs using Hessian-Free Trust-Region Methods

Discovering the underlying physical behavior of complex systems is a cru...
03/27/2020

### On a generalized Collatz-Wielandt formula and finding saddle-node bifurcations

We introduce the nonlinear generalized Collatz-Wielandt formula λ^*...
02/19/2021

### A convergent finite difference method for computing minimal Lagrangian graphs

We consider the numerical construction of minimal Lagrangian graphs, whi...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

This paper develops and analyzes a narrow-stencil finite difference approximation of the viscosity solution to the following fully nonlinear second order Dirichlet boundary value problem:

 (1a) F[u](x)≡F(D2u,∇u,u,x) =0, ∀x∈Ω, (1b) u(x) =g(x), ∀x∈∂Ω,

where () is a bounded domain and denotes the Hessian matrix of at . The PDE operator is a fully nonlinear second order differential operator in the sense that is nonlinear in (or in at least one of its components). We assume is continuous and is Lipschitz in its first three arguments. Moreover, is assumed to be elliptic and satisfy a comparison principle (defined in Section 2.3). The paper focusses on a particular class of fully nonlinear PDEs called Hamilton-Jacobi-Bellman (HJB) equations which correspond to , where

 (2a) H[u] ≡minθ∈Θ(Lθu−fθ), (2b) Lθu =Aθ(x):D2u+bθ(x)⋅∇u+cθ(x)u.

HJB equations arise in various applications including stochastic optimal control, game theory, and mathematical finance (cf.

[15]).

Numerical fully nonlinear PDEs has experienced some rapid developments in recent years. Various numerical methods have been proposed and tested (cf. [12, 23] and the references therein). In order to guarantee convergence to the viscosity solution of problem (1), most of these works follow the framework laid out by Barles-Souganidis in [1] which requires that a numerical scheme be monotone, consistent, and stable. It follows from a well-known result of [22] that, in general, monotone schemes are necessarily wide-stencil schemes as seen in [9, 10, 26, 29]. Consequently, to design convergent narrow-stencil schemes, one has to abandon the standard monotonicity requirement (in the sense of Barles-Souganidis [1]) and, instead, deal with non-monotone schemes. This in turn prohibits one to directly use the powerful and handy Barles-Souganidis (analysis) framework. On the other hand, it should be pointed out that monotone narrow-stencil schemes could be constructed for problem (1) in special cases. For example, when the matrix is diagonally dominant for all , the methods proposed in [20, 2] become narrow-stencil schemes even though they require wide-stencils for general negative semi-definite matrices . For isotropic coefficient matrices (i.e., diagonal), a monotone linear finite element method was proposed in [18]. Finally, we also mention that when satisfies a Cordes condition, a least square-like interior penalty discontinuous Galerkin (IP-DG) method was proposed in [30]. This method is obviously a narrow-stencil scheme, and it was proved to be convergent even though it is not monotone.

The primary goal of this paper is to develop a convergent narrow-stencil finite difference method for problem (1) in the general case when is uniformly elliptic. Thus, we only require to be uniformly negative definite as long as the underlying PDE satisfies the comparison principle. The method to be introduced in this paper can be regarded as a high dimensional extension of the one-dimensional finite difference method developed and analyzed by the authors in [14]. We note that this extension is far from trivial in terms of both method design and convergence analysis because of the additional difficulties caused by the second order mixed derivatives on the off-diagonal entries of the Hessian matrix for . Indeed, constructing “correct” finite difference discretizations for those mixed derivatives and controlling their “bad” effect in the convergence analysis poses some significant new challenges. Recall that the main reason to use wide-stencils in the works cited above is exactly to “correctly” approximate the second order mixed derivatives. Hence, the new challenges which must be dealt with in this paper can be regarded as the price we need to pay for not using wide-stencils.

To achieve our goal, several key ideas are introduced and utilized in this paper. First, we introduce multiple discrete approximations of and to locally capture the behavior of the underlying function when it has low regularity. The reason for using such multiple approximations is to compensate for the low resolution of a fixed grid, especially when using a narrow-stencil finite difference approximation. Second, as explained above, using narrow-stencils forces us to abandon the hope of constructing monotone finite difference schemes (again, in the sense of Barles-Souganidis [1]). However, in order to obtain a convergent scheme, the anticipated finite difference method must have some “good” structures. Our key motivating idea is to impose such a structure, called generalized monotonicity, on our finite difference scheme. We shall demonstrate that besides allowing for the use of narrow-stencils, generalized monotonicity is a more natural concept/structure than the usual monotonicity concept/structure, and it is much easier to verify. Third, to obtain practical generalized monotone finite difference methods, two important questions to answer are: how to build such a scheme and what is the primary mechanism that makes the method work? Inspired by the vanishing moment method developed in [11, 13], we propose a higher order stabilization term/mechanism, called a numerical moment, as a means to address these questions. We will see that the stabilization term/mechanism is the key for us to overcome the lack of monotonicity in the admissibility, stability, and convergence analysis.

The remainder of this paper is organized as follows. In Section 2 we introduce the basics of viscosity solution theory and elliptic operators as well as the corresponding notation. We also define our primary finite difference operators that will be used throughout the paper and list some identities for those difference operators. In Section 3 we first formulate our finite difference method and then prove several properties of the scheme. Section 4 is devoted to studying the well-posedness of the scheme. A contraction mapping technique based on a pseudo-time explicit Euler scheme is employed to prove the existence and uniqueness for the proposed finite difference scheme and to derive an

-norm stability estimate. In Section

5, an -norm stability estimate for the second-order central difference quotients of the numerical solution is first derived using a Sobolev iteration technique. Then an -norm stability estimate is derived based on combining both estimates. The results in Sections 4 and 5 rely upon some auxiliary linear algebra results that can be found in Appendix A. In Section 6 we present a detailed convergence analysis for the proposed finite difference method. In Section 7 we provide some numerical experiments that validate and complement the theory. The paper is concluded in Section 8 with a brief summary and some final remarks.

## 2 Preliminaries

We begin by recalling the viscosity solution concept for problem (1) as documented in [6, 8, 5] and defining basic difference operators that will be used to construct our finite difference method. We define the relevant function space notation in Section 2.1, state the definition of viscosity solutions in Section 2.2, and recall the definition of ellipticity and the comparison principle in Section 2.3. The comparison principle ensures uniqueness of solutions to (1). The specific assumptions for HJB equations will be given in Section 2.4. Lastly, in Section 2.5, we define our basic difference operators which will serve as the building blocks for constructing various discrete gradient and Hessian operators. The discrete gradient and Hessian operators will be used to define our narrow-stencil finite difference method in Section 3.

### 2.1 Function spaces and notation

Standard function and space notation as in [4] and [16] will be adopted in this paper. For example, for a bounded open domain , , , and are used to denote, respectively, the spaces of bounded, upper semicontinuous, and lower semicontinuous functions on . For any , we define

 v∗(x)≡limsupy→xv(y)andv∗(x)≡liminfy→xv(y).

Then, and , and they are called the upper and lower semicontinuous envelopes of , respectively. We also let denote the set of symmetric real-valued matrices.

In presenting the relevant viscosity solution concept for fully nonlinear second order PDEs, we shall use to denote a general fully nonlinear second order operator. More precisely, . The general second order fully nonlinear PDE problem is to seek a locally bounded function such that is a viscosity solution of

 (3) G(D2u,∇u,u,x) =0in ¯¯¯¯Ω,

where we have adopted the convention of writing the boundary condition as a discontinuity of the PDE (cf. [1, p.274]). We shall treat the Dirichlet boundary condition explicitly in the formulation of our finite difference method for (3).

### 2.2 Viscosity solutions

Due to the full nonlinearity of the differential operator in (3), the standard weak solution concept based upon lowering the order of derivatives on the solution function using integration by parts is not applicable. Thus, a different solution concept is necessary for fully nonlinear PDEs. Among several available weak solution concepts for fully nonlinear PDEs, the best known one is the viscosity solution notion whose definition is stated below.

Let denote the second order operator in (3).

(i) A locally bounded function is called a viscosity subsolution of (3) if , when has a local maximum at ,

(ii) A locally bounded function is called a viscosity supersolution of (3) if , when has a local minimum at ,

(iii) A locally bounded function is called a viscosity solution of (3) if is both a viscosity subsolution and a viscosity supersolution of (3).

###### Remark 2.1

(a) Applying the above definition to the Dirichlet boundary value problem (1), we require the viscosity solution satisfies the Dirichlet boundary condition (1b) in the viscosity sense. Using the notation of [5, Section 7], this can be described as follows. First, define the boundary operator by

 (4) B(A,q,v,x)≡v(x)−g(x)

based on the Dirichlet boundary condition (1b), and then define the differential operator by

 G(D2u,∇u,u,x)≡{F[u](x),if x∈Ω,B[u](x),if x∈∂Ω.

Then we have

 G∗(A,q,v,x)={F(A,q,v,x),if x∈Ω,min{F(A,q,v,x),B(A,q,v,x)},if x∈∂Ω

and

 G∗(A,q,v,x)={F(A,q,v,x),if x∈Ω,max{F(A,q,v,x),B(A,q,v,x)},if x∈∂Ω.

Using the boundary operator , we have explicitly written the boundary condition as a discontinuity of the PDE operator .

(b) If is continuous, then another way to enforce the Dirichlet boundary condition is in the pointwise sense. As explained in [5], the enforcing of the Dirichlet boundary condition in the viscosity sense and in the pointwise sense may not be equivalent in general.

The above definitions can be informally interpreted as follows. Without a loss of generality, we may assume whenever has a local maximum at or whenever has a local minimum at . Then, is a viscosity solution of (3) if for all smooth functions such that “touches” the graph of from above at , we have , and for all smooth functions such that “touches” the graph of from below at , we have . A geometric interpretation of the definition of viscosity solutions can be found in Figure 1. Since the definition is based on locally moving derivatives from the viscosity solution to a smooth test function, the concept of viscosity solutions relies upon a “differentiation-by-parts” approach. This is in contrast to weak solution theory where derivatives are moved to a test function using an “integration-by-parts” approach, which, in general, is not a local operation.

### 2.3 Ellipticity and the comparison principle

An advantage of the viscosity solution concept for fully nonlinear problems is that, due to its entirely local structure, it is flexible enough to seek solutions in the space of bounded functions. To ensure the existence and uniqueness of viscosity solutions, additional structure conditions on the differential operator are often necessary. The most common such conditions are an ellipticity requirement and a comparison principle. The following definitions of ellipticity and the comparison principle are standard (cf. [5]).

(i) Equation (3) is said to be uniformly elliptic if, for all , there holds

 λtr (A−B)≤G(B,p,v,x)−G(A,p,w,x)≤Λtr (A−B)

whenever , where and means that is a nonnegative definite matrix.

(ii) Equation (3) is said to be proper elliptic if, for all , there holds

 G(A,p,v,x)≤G(B,p,w,x)

whenever and , where and .

(iii) Equation (3) is said to be degenerate elliptic if, for all , there holds

 G(A,p,v,x)≤G(B,p,v,x)

whenever , where .

Problem (3) is said to satisfy a comparison principle if the following statement holds. For any upper semicontinuous function and lower semicontinuous function on , if is a viscosity subsolution and is a viscosity supersolution of (3), then on .

The comparison principle is sometimes called a strong uniqueness property due to the fact it immediately yields the uniqueness of viscosity solutions for problem (3) (cf. [1]). When is differentiable, the degenerate ellipticity condition can also be defined by requiring that the matrix is negative semi-definite while the proper ellipticity condition additionally requires the number is nonnegative (cf. [16, p. 441]). We then clearly see that the ellipticity concept for nonlinear problems generalizes the notion for linear elliptic equations.

###### Remark 2.2

Since we have assumed that satisfies the comparison principle, it follows that that problem (3) has a unique viscosity solution .

### 2.4 Hamilton-Jacobi-Bellman equations

As mentioned in Section 1, a specific class of fully nonlinear PDEs that fits the structure requirements in this paper are (stationary) Hamilton-Jacobi-Bellman (HJB) equations described by (2).

We shall assume that is uniformly (in ) negative definite in the sense that there exists constants such that

 (5) −Λ|ξ|2≤Aθ(x)ξ⋅ξ≤−λ|ξ|2∀ξ∈Rd,x∈Ω,θ∈Θ.

We also assume there exists a constant such that

 ∥Aθ∥L∞(Ω),∥bθ∥L∞(Ω),∥cθ∥L∞(Ω)≤K

and for some constant for all and . Thus, the equation is proper elliptic and satisfies the comparison principle. We note that for stochastic optimal control applications, we have (cf. [15]). This case will also be considered. However, such problems may not satisfy a comparison principle as discussed in [24] and [28].

Under the above assumptions on the coefficients and the control set, problem (1)–(2) has a unique solution and is globally Lipschitz continuous. Thus, is differentiable almost everywhere and the weak derivatives are bounded by the constant . Since is closed and bounded, there actually exists a function such that

 H[u]=minθ∈Θ(Lθu−fθ)=Aθ∗:D2u+bθ∗⋅∇u+cθ∗u=0in Ω.

Thus, , , and , and we have is differentiable everywhere with bounded first order derivatives.

The analysis assumes a global Lipschitz condition for ease of presentation; however, locally Lipschitz should be sufficient. The global Lipschitz assumption excludes Monge-Ampère-type equations which are only conditionally degenerate elliptic [4] and locally Lipschitz. In light of the recent work [10], we know that the Monge-Ampère equation has an equivalent HJB reformulation. Much of the work of this paper can be extended to the Monge-Ampère equation via its HJB reformulation.

### 2.5 Difference operators

We introduce several difference operators for approximating first and second order partial derivatives. The multiple difference operators will be used to help resolve the low regularity of viscosity solutions and will play a key role in motivating and defining our new narrow-stencil FD methods.

Assume is a -rectangle, i.e., . We shall only consider grids that are uniform in each coordinate , . Let be an integer and for . Define , , , and . Then, . We partition into sub--rectangles with grid points for each multi-index . We call a mesh (set of nodes) for . We also introduce an extended mesh which extends by a collection of ghost grid points that are at most one layer exterior to in each coordinate direction. In particular, we choose ghost grid points such that for some , . We set and is defined by replacing by in the definition of and then removing the extra multi-indices that would correspond to ghost grid points that are not in the set to ensure .

#### 2.5.1 First order difference operators

We define the standard forward and backward difference operators as well as the central difference operator for approximating first order derivatives including the gradient operator. The forward and backward difference operators will serve as the building blocks for constructing all of the first and second order difference operators in this paper.

Let

denote the canonical basis vectors for

. Define the (first order) forward and backward difference operators by

 (6) δ+xi,hiv(x)≡v(x+hiei)−v(x)hi,δ−xi,hiv(x)≡v(x)−v(x−hiei)hi

for a function defined on and

 δ+xi,hiVα≡Vα+ei−Vαhi,δ−xi,hiVα≡Vα−Vα−eihi

for a grid function defined on the grid . Note that “ghost-values” may need to be introduced in order for the above difference operators to be well-defined on the boundary of . We also define the following central difference operator:

 (7) δxi,hi≡12(δ+xi,hi+δ−xi,hi)for i=1,2,⋯,d.

Lastly, we define the “sided” and central gradient operators , , and by

 (8) ∇±h≡[δ±x1,h1,δ±x2,h2,⋯,δ±xd,hd]T,∇h≡[δx1,h1,δx2,h2,⋯,δxd,hd]T.

#### 2.5.2 Second order difference operators

We now introduce a number of difference operators that approximate second order derivatives including the Hessian operator. Using the forward and backward difference operators introduced in the previous subsection, we have the following four possible approximations of the second order differential operator :

 Dμνh,ij≡δνxj,hjδμxi,hifor μ,ν∈{+,−},

which in turn leads to the definition of the following four approximations of the Hessian operator :

 (9) Dμνh≡[Dμνh,ij]di,j=1for μ,ν∈{+,−}.

To construct our numerical methods in the next section, we also need to introduce the following three sets of averaged second order difference operators:

 (10a) ˆδ2xi,xj;hi,hj ≡12(D+−h,ij+D−+h,ij)=12(δ−xj,hjδ+xi,hi+δ+xj,hjδ−xi,hi), (10b) ˜δ2xi,xj;hi,hj ≡12(D−−h,ij+D++h,ij)=12(δ+xj,hjδ+xi,hi+δ−xj,hjδ−xi,hi), (10c) ¯¯¯δ2xi,xj;hi,hj ≡12(ˆδ2xi,xj;hi,hj+˜δ2xi,xj;hi,hj)=14(D++h,ij+D+−h,ij+D−+h,ij+D−−h,ij) =14(δ+xj,hjδ+xi,hi+δ−xj,hjδ+xi,hi+δ+xj,hjδ−xi,hi+δ−xj,hjδ−xi,hi)

for all . For notation brievity, we also set

 δ2xi,hi≡ˆδ2xi,hi≡ˆδ2xi,xi;hi,hi,˜δ2xi,hi≡˜δ2xi,xi;hi,hi,¯¯¯δ2xi,hi≡¯¯¯δ2xi,xi;hi,hi.

Using the above difference operators, we define the following three “centered” approximations of the Hessian operator :

 (11) ˆD2h=[ˆδ2xi,xj;hi,hj]di,j=1,˜D2h=[˜δ2xi,xj;hi,hj]di,j=1,¯¯¯¯¯D2h=[¯¯¯δ2xi,xj;hi,hj]di,j=1.

We will also consider the (standard) central approximation of the Hessian operator defined by

 (12) [D2h]ij≡⎧⎨⎩δ2xi,hiif i=j,¯¯¯δ2xi,xj;hi,hj,if i≠j.

We now record some facts that can be easily verified for the various second order operators. First,

 2ˆD2h=D+−h+D−+h,2˜D2h=D−−h+D++h,2¯¯¯¯¯D2h=ˆD2h+˜D2h.

Second, a simple computation reveals

 (13a) ˆδ2xi,xj;hi,hjv(x) =−v(x+hiei−hjej)−v(x+hiei)−[v(x−hjej)−v(x)]2hihj −v(x−hiei+hjej)−v(x−hiei)−[v(x+hjej)−v(x)]2hihj, (13b) ˜δ2xi,xj;hi,hjv(x) =v(x+hiei+hjej)−v(x+hiei)−[v(x+hjej)−v(x)]2hihj +v(x−hiei−hjej)−v(x−hiei)−[v(x−hjej)−v(x)]2hihj, (13c) ¯¯¯δ2xi,xj;hi,hjv(x) =v(x+hiei+hjej)+v(x−hiei−hjej)4hihj −v(x−hiei+hjej)+v(x+hiei−hjej)4hihj.

Lastly, there holds the relationships

 (14a) δ2xi,hi =δ±xi,hiδ∓xi,hi, (14b) δ2xi,2hi =¯¯¯δ2xi,hi=δxi,hiδxi,hi.

We have is the standard (second order) three-point central difference operator for approximating second order derivatives in one-dimension, is a lower-resolution (second-order) three-point central difference operator for approximating second order derivatives in one-dimension, is a (second order) five-point central difference operator for approximating second order derivatives in one-dimension, is the standard (second order) central difference operator for approximating second order mixed derivatives, and and are two alternative (second order) central difference operators for approximating second order mixed derivatives. Observe that in two dimensions, corresponds to a 7-point stencil, and correspond to 9-point stencils, and corresponds to an 11-point stencil, as seen in Figure 2. The operators and use nearest neighbors in the stencil (including diagonal directions) while and use neighbors two steps away in each Cartesian direction.

We lastly introduce the (diagonal) variables and so that we can define the central difference operators and by

 (15a) δ2ξji,hv(x) ≡v(x−hiei+hjej)−2v(x)+v(x+hiei−hjej)h2i+h2j, (15b) δ2ηji,hv(x) ≡v(x+hiei+hjej)−2v(x)+v(x−hiei−hjej)h2i+h2j

for all with . Then, a direct computation yields

 (16a) [¯¯¯¯¯D2hv]ij =[D2hv]ij=h2i+h2j4hihj(δ2ηji,hv−δ2ξji,hv), (16b) [ˆD2hv]ij =hi2hjδ2xi,hiv+hj2hiδ2xj,hjv−h2i+h2j2hihjδ2ξji,hv, (16c) [˜D2hv]ij =−hi2hjδ2xi,hiv−hj2hiδ2xj,hjv+h2i+h2j2hihjδ2ηji,hv,

relationships that are illustrated in Figure 3. Lastly, there holds

 (17) [˜D2hVα]ii−[ˆD2hVα]ii =˜δ2xi,hiVα−δ2xi,hiVα =2(δ2xi,2hiVα−δ2xi,hiVα) =12(δ2xi,hiVα−ei−2δ2xi,hiVα+δ2xi,hiVα+ei)

for all .

## 3 A narrow-stencil finite difference method

We now propose a new finite difference method for approximating viscosity solutions of the fully nonlinear second order problem (1). In the construction of such a narrow-stencil scheme, our main idea is to design a numerical operator that utilizes the various numerical Hessians defined in the previous section so that a relaxed (or generalized) monotonicity condition is satisfied. The crux of the formulation is to introduce the concept of numerical moments which can be considered analogous to the concept of numerical viscosities used in the Crandall and Lions’ finite difference framework for (first order) Hamilton-Jacobi equations. Below, we first define our narrow-stencil finite difference method, and we then describe the idea and details of numerical moments. We also prove a generalized monotonicity result for our proposed numerical operator.

### 3.1 Formulation of the narrow-stencil finite difference method

In Section 2.5.2 we defined four basic sided discrete Hessian operators for . These four operators are the building blocks for constructing our numerical operator that approximates the differential operator . However, inspired by our earlier 1-D work [14], we construct the numerical operator so that it only depends explicitly upon the central difference operators and with regard to the Hessian argument in . Recall that denotes the extended mesh of that includes ghost mesh points and is the corresponding extended index set of (see Section 2.5).

Our specific narrow-stencil finite difference method is defined by seeking a grid function such that for all

 (18a) ˆF[Uα,xα] =0 for xα∈Th∩Ω, (18b) Uα =g(xα) for xα∈Th∩∂Ω,

where

 (19) ˆF[Uα,xα] ≡˜F(D−−hUα,D−+hUα,D+−hUα,D++hUα,∇+hUα,∇−hUα,Uα,xα) ≡ˆF(ˆD2hUα,˜D2hUα,∇+hUα,∇−hUα,Uα,xα) ≡F(¯¯¯¯¯D2hUα,∇hUα,Uα,