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:
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
HJB equations arise in various applications including stochastic optimal control, game theory, and mathematical finance (cf.).
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  which requires that a numerical scheme be monotone, consistent, and stable. It follows from a well-known result of  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 ) 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 . Finally, we also mention that when satisfies a Cordes condition, a least square-like interior penalty discontinuous Galerkin (IP-DG) method was proposed in . 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 . 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 ). 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 Section5, 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.
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  and  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
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
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 ,
(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
based on the Dirichlet boundary condition (1b), and then define the differential operator by
Then we have
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 , 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. ).
(i) Equation (3) is said to be uniformly elliptic if, for all , there holds
whenever , where and means that is a nonnegative definite matrix.
(ii) Equation (3) is said to be proper elliptic if, for all , there holds
whenever and , where and .
(iii) Equation (3) is said to be degenerate elliptic if, for all , there holds
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. ). 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.
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
We shall assume that is uniformly (in ) negative definite in the sense that there exists constants such that
We also assume there exists a constant such that
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. ). This case will also be considered. However, such problems may not satisfy a comparison principle as discussed in  and .
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
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  and locally Lipschitz. In light of the recent work , 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.
denote the canonical basis vectors for. Define the (first order) forward and backward difference operators by
for a function defined on and
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:
Lastly, we define the “sided” and central gradient operators , , and by
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 :
which in turn leads to the definition of the following four approximations of the Hessian operator :
To construct our numerical methods in the next section, we also need to introduce the following three sets of averaged second order difference operators:
for all . For notation brievity, we also set
Using the above difference operators, we define the following three “centered” approximations of the Hessian operator :
We will also consider the (standard) central approximation of the Hessian operator defined by
We now record some facts that can be easily verified for the various second order operators. First,
Second, a simple computation reveals
Lastly, there holds the relationships
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
for all with . Then, a direct computation yields
relationships that are illustrated in Figure 3. Lastly, there holds
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 , 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