2021 UNC Greensboro PDE Conference, Electronic Journal of Differential Equations, Conference 26 (2022), pp. 59–95. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu A NARROW-STENCIL FRAMEWORK FOR CONVERGENT NUMERICAL APPROXIMATIONS OF FULLY NONLINEAR SECOND ORDER PDES XIAOBING FENG, THOMAS LEWIS, KELLIE WARD Abstract. This article develops a unified general framework for designing convergent finite difference and discontinuous Galerkin methods for approxi- mating viscosity and regular solutions of fully nonlinear second order PDEs. Unlike the well-known monotone (finite difference) framework, the proposed new framework allows for the use of narrow stencils and unstructured grids which makes it possible to construct high order methods. The general frame- work is based on the concepts of consistency and g-monotonicity which are both defined in terms of various numerical derivative operators. Specific methods that satisfy the framework are constructed using numerical moments. Ad- missibility, stability, and convergence properties are proved, and numerical experiments are provided along with some computer implementation details. 1. Introduction This article develops a unified general framework for designing convergent narrow- stencil finite difference (FD) and discontinuous Galerkin (DG) methods for approx- imating the viscosity (and regular) solution to the following fully nonlinear second order Dirichlet boundary value problem: F [u](x) ≡ F (D2u,∇u, u,x) = 0, ∀x ∈ Ω, (1.1a) u(x) = g(x), ∀x ∈ ∂Ω, (1.1b) where Ω ⊂ Rd (d = 2, 3) is a bounded domain and D2u(x) denotes the Hessian matrix of u at x. The partial differential equation (PDE) operator F : Sd×d×Rd× R×Ω→ R, where Sd×d ⊂ Rd×d denotes the set of d×d symmetric real matrices, is a fully nonlinear second order differential operator in the sense that F is nonlinear in at least one component of the Hessian D2u. The boundary data g is assumed to be continuous with F Lipschitz with respect to its first three arguments. Moreover, F is assumed to be uniformly and proper elliptic and satisfy a comparison principle (see Section 2 for the definitions). In this paper we focus our attention on two main classes of fully nonlinear second order PDEs, namely, the Monge-Ampère-type and Hamilton-Jacobi-Bellman-type equations (cf. [21, 9]). 2020 Mathematics Subject Classification. 65N06, 65N12. Key words and phrases. Fully nonlinear PDEs; viscosity solutions; Monge-Ampère equation; Hamilton-Jacobi-Bellman equation; narrow-stencil; generalized monotonicity; g-monotonicity; numerical operators; numerical moment. ©2022 This work is licensed under a CC BY 4.0 license. Published August 25, 2022. 59 60 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 Fully nonlinear second order PDEs arise from many scientific and engineering applications such as antenna design, astrophysics, economics, differential geometry, stochastic optimal control, and optimal mass transport; yet, they are a class of PDEs which are difficult to study analytically and even more challenging to ap- proximate numerically. Due to the fully nonlinear structure, there is no general variational (or weak) formulation. As a result, its weak solution concept (called viscosity solutions, see Section 2 for the definition) is complicated and, in partic- ular, very difficult to address numerically. Nevertheless, driven by the need for solving many emerging and intriguing application problems, numerical fully non- linear PDEs has garnered a lot of attention and experienced rapid developments in recent years. See [9, 27] and the references therein for an overview of various numerical methods that have been proposed, analyzed, and tested. To the best of our knowledge, there are only two main approaches in the litera- ture which aim to approximate viscosity solutions. Both approaches have been used successfully to design and analyze (practical) numerical methods for approximating second order fully nonlinear PDEs. The first approach, which was adopted in the overwhelming majority of the existing works, is the Barles-Souganidis’ monotone (wide-stencil) finite difference framework (cf. [1, 29, 7, 10, 22, 28]). We call the sec- ond approach the numerical moment-enhanced g-monotone (narrow-stencil) finite difference and DG framework (cf. [11, 12, 14, 24] and also the original vanishing moment method [17]). It is well-known that monotone (in the sense of Barles- Souganidis [1]) methods are difficult to construct; moreover, they are intrinsically low order, require the use of wide stencils, and yield strongly coupled nonlinear algebraic problems that need to be solved. The narrow-stencil approach aims to sidestep these limitations of the wide-stencil approach so high order methods can be constructed on both structured and non-structured grids. On the other hand, since the narrow-stencil approach abandons the standard monotonicity require- ment, it prevents one to directly use the powerful Barles-Souganidis’ framework for the convergence analysis. Consequently, new machineries and techniques must be developed for the analysis of the proposed narrow-stencil methods which, so far, has only been done on a case-by-case basis in [11, 12, 14, 24]. The primary goal of this paper is to re-examine (with a top-down view) and refine the numerical moment-enhanced g-monotone (narrow-stencil) finite difference and DG approach, which was initiated by us in [11, 12, 14, 24], and to formulate it into a unified framework which is parallel to the Barles-Souganidis’ framework. A far-reaching goal is to provide a blueprint/framework for designing and analyzing practical convergent numerical methods for approximating viscosity (and regular) solutions of fully nonlinear second order PDEs. The remainder of this paper is organized as follows. In Section 2, we first recall the basics of viscosity solution theory and elliptic operators as well as the necessary notations. We then introduce various FD and DG finite element numerical de- rivative operators (cf. [14, 16]) and unify the notations. Those discrete derivative operators are the building blocks for our narrow-stencil framework. In Section 3, we formulate our abstract framework and introduce the key concepts of numerical op- erators, consistency, g-monotonicity, and numerical moments. We then state a few structure conditions/assumptions for numerical operators. It should be noted that all of these concepts and conditions are motivated by and abstractions of similar ones in our earlier works [11, 12, 14, 24]. The numerical moment will play a critical EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 61 role in the specific examples of numerical operators that are presented. In Sections 4–6, we present a complete convergence, admissibility, and stability analysis for the narrow-stencil FD methods proposed in Section 3. Unlike the Barles-Souganidis’ monotone (wide-stencil) framework where admissibility and the `∞-norm stabil- ity of the underlying numerical methods are almost free to obtain (thanks to the monotonicity), our results require entirely new techniques and their proofs given in Sections 4–6 are more technical and involved (as well as much longer). These technical issues are precisely the price to pay for using narrow stencils. Assuming the admissibility and `∞ stability, we first establish the convergence of the numer- ical solution to the viscosity solution of the underlying PDE problem in Section 4. The proof is adapted from the much more detailed version in [14]. We then prove the desired admissibility and `∞ stability in Sections 5 and 6. Our main idea of proving the admissibility is to use the Contractive Mapping Theorem in `2 instead of `∞. The `∞ stability is obtained by a novel numerical embedding technique first introduced in [14]. Finally, in Section 7, we present some numerical experiments to demonstrate the effectiveness of the proposed framework and to address some computer implementation issues. 2. Preliminaries and numerical derivatives 2.1. Notation and definitions. The narrow-stencil framework will rely upon two different partial orderings for matrices. Let Sd×d ⊂ Rd×d denote the set of symmet- ric d×d matrices. We will utilize the convention that A ≥ B if A−B is symmetric nonnegative definite for A,B ∈ Sd×d. We will also introduce the alternative con- vention that A � B if each component of A − B is nonnegative. Note that the partial ordering induced by � does not require symmetric matrices. We also let A : B denote the Frobenius inner product with A : B ≡ ∑d i=1 ∑d j=1 aijbij for all matrices A,B ∈ Rd×d. For a bounded open domain Ω ⊂ Rd, let B(Ω), USC(Ω), and LSC(Ω) denote, respectively, the spaces of bounded, upper semi-continuous, and lower semicontin- uous functions on Ω. For any v ∈ B(Ω), we define v∗(x) := lim sup y→x v(y) and v∗(x) := lim inf y→x v(y). Then, v∗ ∈ USC(Ω) and v∗ ∈ LSC(Ω), and they are called the upper and lower semicontinuous envelopes of v, respectively. Below we use the convention of writing the boundary condition as a discontinuity of the PDE (cf. [1, p.274]). The Dirichlet boundary condition is assumed to hold in the viscosity sense. The following two definitions can be found in [21, 3, 1]. Definition 2.1. Equation (1.1) is said to be proper elliptic if for all (q,x) ∈ Rd×Ω, it holds F (A,q, v,x) ≤ F (B,q, w,x) ∀A,B ∈ Sd×d, A ≥ B, v, w ∈ R, v ≤ w. Definition 2.2. Equation (1.1) is said to be uniformly elliptic if there exists Λ ≥ λ > 0 such that, for all (q, v,x) ∈ Rd × R× Ω, it holds 0 ≥ −λ tr(A−B) ≥ F (A,q, v,x)− F (B,q, v,x) ≥ −Λ tr(A−B) for all A,B ∈ Sd×d with A ≥ B. 62 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 We note that when F (A,q, v,x) is differentiable with respect to the first param- eter, then the proper ellipticity definition is equivalent to requiring that the matrix ∂F ∂A is negative semi-definite and the value ∂F ∂v is nonnegative (cf. [21, p. 441]). If F is also uniformly elliptic, then it holds 0 > −λ|~ξ|2 ≥ ~ξ · ∂F∂A ~ξ ≥ −Λ|~ξ|2 for all ~ξ 6= ~0. Thus, λI ≤ −∂F∂A ≤ ΛI. Definition 2.3. A function u ∈ B(Ω) is called a viscosity subsolution (resp. su- persolution) of (1.1) if, for all ϕ ∈ C2(Ω), if u∗ − ϕ (resp. u∗ − ϕ) has a local maximum (resp. minimum) at x0 ∈ Ω, then we have F∗(D 2ϕ(x0),∇ϕ(x0), u∗(x0),x0) ≤ 0 (resp. F ∗(D2ϕ(x0),∇ϕ(x0), u∗(x0),x0) ≥ 0). The function u is said to be a vis- cosity solution of (1.1) if it is simultaneously a viscosity subsolution and a viscosity supersolution of (1.1). Definition 2.4. Problem (1.1) is said to satisfy a comparison principle if the following statement holds. For any upper semicontinuous function u and lower semicontinuous function v on Ω, if u is a viscosity subsolution and v is a viscosity supersolution of (1.1), then u ≤ v on Ω. Since we assume that F in (1.1) satisfies the comparison principle, we have that the underlying viscosity solution u must be continuous. Furthermore, if F is continuous with respect to x, then by the Lipschitz continuity with respect to D2u and u, we can drop the upper and lower ∗ indices in Definition 2.3. 2.2. Finite difference derivative operators. We introduce several difference operators for approximating first and second order partial derivatives. The narrow- stencil framework will use multiple difference operators to help resolve the underly- ing viscosity solution. The notation and difference operators used are the same as those in [14]. The section ends with a formal result for comparing various discrete second order operators. 2.2.1. Finite difference grids. Assume Ω is a d-rectangle, i.e., Ω = (a1, b1)×(a2, b2)× · · ·×(ad, bd), and let {ei}di=1 denote the canonical basis vectors for Rd. We shall only consider grids that are uniform in each coordinate xi, i = 1, 2, . . . , d. Let Ji(≥ 2) be an integer and hi = bi−ai Ji−1 for i = 1, 2, . . . , d. Define h = (h1, h2, . . . , hd) ∈ Rd, h = maxi=1,2,...,d hi, J = ∏d i=1 Ji, and NJ = {α = (α1, α2, . . . , αd) | 1 ≤ αi ≤ Ji, i = 1, 2, . . . , d}. Then, |NJ| = J . We partition Ω into ∏d i=1 (Ji − 1) sub-d- rectangles with grid points xα = ( a1+(α1−1)h1, a2+(α2−1)h2, . . . , ad+(αd−1)hd ) for each multi-index α ∈ NJ. We call Th = {xα}α∈NJ a mesh (set of nodes) for Ω. We also introduce an extended mesh T ′h which extends Th 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 x such that x = y ± 2hiei for some y ∈ Th ∩ Ω, i ∈ {1, 2, . . . , d}. We set J ′i = Ji + 2 and NJ ′ is defined by replacing Ji by J ′i in the definition of NJ and then removing the extra multi-indices that would correspond to ghost grid points that are not in the set T ′h to ensure |NJ ′| = |T ′h|. 2.2.2. First order finite difference operators. Define the (first order) forward and backward difference operators by δ+ xi,hi v(x) ≡ v(x + hiei)− v(x) hi , δ−xi,hiv(x) ≡ v(x)− v(x− hiei) hi (2.1) EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 63 for a function v defined on Rd. We also consider the central difference operator δxi,hi defined by δxi,hi ≡ 1 2δ + xi,hi + 1 2δ − xi,hi so that δxi,hiv(x) ≡ v(x + hiei)− v(x− hiei) 2hi . The corresponding forward, backward, and central discrete gradient operators are denoted by ∇+ h , ∇−h , and ∇h, respectively. 2.2.3. Second order finite difference operators. 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 ∂2 xixj given by Dµν h,ij ≡ δνxj ,hjδ µ xi,hi for µ, ν ∈ {+,−}, which in turn leads to the definition of the following four approximations of the Hessian operator D2 ≡ [∂2 xixj ]: Dµν h ≡ [ Dµν h,ij ]d i,j=1 for µ, ν ∈ {+,−}. To analyze our narrow-stencil framework in the next section, we also need to introduce the following two sets of averaged second order difference operators: δ̂2 xi,xj ;hi,hj ≡ 1 2 ( D+− h,ij +D−+ h,ij ) = 1 2 ( δ−xj ,hjδ + xi,hi + δ+ xj ,hj δ−xi,hi ) , δ̃2 xi,xj ;hi,hj ≡ 1 2 ( D−−h,ij +D++ h,ij ) = 1 2 ( δ+ xj ,hj δ+ xi,hi + δ−xj ,hjδ − xi,hi ) , (2.2) for all i, j = 1, 2, . . . , d. Note that, while the various components of Dµν h may not be self-adjoint operators, the discrete second order partial derivative operators δ̂2 xi,xj ;hi,hj and δ̃2 xi,xj ;hi,hj defined by (2.2) are self-adjoint as can be verified by the component forms in [14]. Using the above difference operators, we define the following two “centered” approximations of the Hessian operator D2 ≡ [∂2 xixj ]: D̂2 h ≡ [ δ̂2 xi,xj ;hi,hj ]d i,j=1 , D̃2 h ≡ [ δ̃2 xi,xj ;hi,hj ]d i,j=1 . (2.3) We will also consider the average central approximation of the Hessian operator D 2 h defined by D 2 h ≡ 1 2 D̂2 h + 1 2 D̃2 h. (2.4) A visual representation of the local stencils for the various discrete Hessians D̃2 h, D̂2 h, and D 2 h can be found in Figure 1. For notation brevity, we set δ2 xi,hi ≡ δ̂2 xi,hi ≡ δ̂2 xi,xi;hi,hi and δ 2 xi,hi ≡ δ 2 xi,xj ;hi,hj . Then δ2 xi,hiv(x) = v(x− hiei)− 2v(x) + v(x + hiei) h2 i , δ 2 xi,hi = δxi,hiδxi,hi = δ2 xi,2hi for all i = 1, 2, . . . , d. Lastly, we denote the discrete Laplacian operator by ∆h ≡∑d i=1 δ 2 xi,hi . 64 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 D̂2 h D̃2 h D 2 h uxx uyy uxy Figure 1. Illustration of the local stencils for the three discrete Hessian operators D̂2 h, D̃2 h, and D 2 h. 2.2.4. Properties of second order finite difference operators. We now derive a rela- tionship between D̂2 h and D̃2 h that will be essential to the analysis of our narrow- stencil methods. The following expands upon observations in [14]. Let V be a grid function defined over T ′h. Suppose V = 0 over Th ∩ ∂Ω and V satisfies ∆hVα = 0 for all xα ∈ Sh, where Sh ⊂ Th ∩ ∂Ω is defined by Sh ≡ { xα ∈ Th ∩ ∂Ω : xα + hiei ∈ Th ∩ Ω or xα − hiei ∈ Th ∩ Ω (2.5) for some i ∈ {1, 2, . . . , d} } . Since V = 0 over Th ∩ ∂Ω, it also holds δ2 xi,hi Vα = 0 for all xα ∈ Sh and i ∈ {1, 2, . . . , d} such that either xα−ei ∈ Th∩Ω or xα+ei ∈ Th∩Ω. Note that the additional boundary information is needed to incorporate the ghost points associated with the extended grid T ′h. Let D̂ij,0 and D̃ij,0 denote the matrix representations of [D̂2 h]ij and [D̃2 h]ij , re- spectively, restricted to grid functions defined over Th ∩ Ω with the boundary as- sumptions for V built in. Notationally, the zero subscript is used to denote the boundary conditions. Then D̂ii,0 is symmetric positive definite for all i = 1, 2, . . . , d. Lastly, let δ2 xi,hi;0 denote the central difference operator δ2 xi,hi with the zero Dirich- let boundary data assumption. It is easy to check that (cf. [14]) it holds δ̃2 xi,xj ;hi,hjVα − δ̂ 2 xi,xj ;hi,hjVα = hihj 2 δ2 xi,hiδ 2 xj ,hjVα for all i, j = 1, 2, . . . , d and xα ∈ Th ∩ Ω. Suppose i 6= j. If xα±ej ∈ Th ∩ ∂Ω, then δ2 xi,hi Vα±ej = 0 using the fact that xi is orthogonal to xj and V = 0 over Th ∩ ∂Ω. Then δ2 xi,hi;0δ 2 xj ,hj ;0Vα = δ2 xi,hi;0δ 2 xj ,hjVα = δ2 xi,hiδ 2 xj ,hjVα, and, by a simple computation, δ2 xi,hiδ 2 xj ,hjVα = δ2 xj ,hjδ 2 xi,hiVα = δ2 xj ,hj ;0δ 2 xi,hi;0Vα for all xα ∈ Th ∩ Ω. Thus, D̃ij,0 − D̂ij,0 = hihj 2 D̂ii,0D̂jj,0 = hihj 2 D̂jj,0D̂ii,0 for all i 6= j, and it follows that the matrix is symmetric positive definite. Choose i ∈ {1, 2, . . . , d}. Observe that, since δ2 xi,hi Vα = 0 over Sh, it holds δ2 xi,hi δ2 xi,hi Vα = δ2 xi,hi;0 δ2 xi,hi Vα. Thus, by the Dirichlet boundary condition, we have δ2 xi,hiδ 2 xi,hiVα = δ2 xi,hi;0δ 2 xi,hi;0Vα, EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 65 and it follows that D̃ii,0 − D̂ii,0 = h2 i 2 D̂ii,0D̂ii,0. Therefore, the matrix D̃ii,0 − D̂ii,0 is symmetric positive definite. We have proved the following lemma, where the result for Dij,0 ≡ 1 2D̂ij,0+ 1 2D̃ij,0 is an immediate consequence: Lemma 2.5. The matrix D̃ij,0 − D̂ij,0 is symmetric positive definite for all i, j ∈ {1, 2, . . . , d}. Furthermore, it holds −D̂ij,0 > −Dij,0 > −D̃ij,0. 2.3. Discontinuous Galerkin finite element derivative operators. We note that the above FD numerical derivative operators are only defined on uniform Cartesian grids. In order to extend them to arbitrary grids, and, in particular, to triangular/tetrahedral grids, we utilize finite element DG numerical derivatives which were first introduced in [16] (also see [12]). Below we recall their definitions and some useful properties, but we shall use new notations which are consistent with the above FD discrete derivative operators. 2.3.1. DG mesh and space notations. Let Ω be a polygonal domain, and let Th denote a locally quasi-uniform shape-regular partition of the domain Ω with h ≡ maxK∈Th(diamK). We introduce the broken L2-space, broken H1-space, and bro- ken C0-space defined by L2(Th) ≡ ∏ K∈Th L2(K), H1(Th) ≡ ∏ K∈Th H1(K), C0(Th) ≡ ∏ K∈Th C0(K) and the broken L2-inner product (·, ·)Th defined by (v, w)Th ≡ ∑ K∈Th ∫ K vw dx ∀v, w ∈ L2(Th). Let EIh denote the set of all interior faces/edges of Th, EBh denote the set of all boundary faces/edges of Th, and Eh ≡ EIh ∪ EBh . Then, for a set Sh ⊂ Eh, we define the broken L2-space by L2(Sh) ≡ ∏ e∈Sh L 2(e) and the broken L2-inner product over Sh by 〈v, w〉Sh ≡ ∑ e∈Sh ∫ e v w ds ∀v, w ∈ L2(Sh). For a fixed integer r ≥ 0, we define the standard DG finite element space V h ⊂ H1(Th) ⊂ L2(Ω) by V h ≡ ∏ K∈Th Pr(K), where Pr(K) denotes the set of all polynomials on K with degree not exceeding r. For K,K ′ ∈ Th, let e = ∂K ∩ ∂K ′ ∈ EIh, and let vK , vK′ denote the restrictions of v to K,K ′, respectively. Without a loss of generality, we assume that the global labeling number of K is smaller than that of K ′ and define the following (standard) jump and average notations: [v] ∣∣ e ≡ vK ∣∣ e − vK′ ∣∣ e , {v} ∣∣ e ≡ vK ∣∣ e + vK′ ∣∣ e 2 (2.6) for any v ∈ Hm(Th). We also define ne ≡ nK = −nK′ as the normal vector to e. Based on the formulation in [12], we will extend the jump and average operators to the boundary of the domain in a nonstandard way. Let K ∈ Th such that e ⊂ 66 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 ∂K ∈ EBh , and define ne as the unit outward normal for the underlying boundary simplex. To unify notation, we impose the convention that the set exterior to the domain Ω has a global labeling number of 0 with the indexing starting at 1 for the “first” label for the simplices in Th. We also assume v ∈ L2(Th) ∩ C0(Th) ⊂ L2(Ω) can be extended to be in ( L2(Th)×L2(Ωc) ) ∩ ( C0(Th)×C0(Ωc) ) , where Ωc = Rd\Ω. Then, we can use the same convention as for interior edges and define [v] ∣∣ e ≡ vΩc ∣∣ e − vK ∣∣ e , {v} ∣∣ e ≡ vΩc ∣∣ e + vK ∣∣ e 2 , (2.7) where vΩc , vK denote the restrictions of v to Ωc,K, respectively. Below we will specify how to choose values for vΩc and how to interpret vK in order to naturally impose a boundary condition. Using the jump and average operators for Eh, we define the labelling-dependent trace operators T±i (v) : Eh → R for each i = 1, 2, . . . , d for a given function v ∈ Hm(Th) by T±i (v) ≡ { v } ∓ 1 2 sgn(n(i) e ) [ v ] where sgn(y) =  1 if y > 0, −1 if y < 0, 0 if y = 0 (2.8) for all y ∈ R and n (i) e denoting the i-th component of ne. Note that the exact trace values for v along the boundary still need to be specified for the jump and average operators. Let K ∈ Th such that e ⊂ ∂K ∈ EBh . Suppose we have Dirichlet boundary data for the given function v, denoted by g. Then, for r ≥ 1, we use the convention that vΩc ∣∣ e = vK ∣∣ e = g ∣∣ e so that T±i (v) = g. The convention yields the standard interpretation as introduced in [16]. If r = 0, we use the convention that vΩc ∣∣ e = g ∣∣ e and vK ∣∣ e is given by the interior limit for v. Thus, T±i (v) is given by either g or the interior limit for v depending on the choice for ± and the sign of the i-th component of the unit normal vector. Such a nonstandard approach allows for weighting degrees of freedom associated with the Dirichlet data against degrees of freedom associated with the value on the interior of K. Since r = 0 implies only one degree of freedom is available on K, such a weighting is essential to not overly emphasize the boundary condition. Notationally, we write T±,gi to denote the natural enforcement of the Dirichlet boundary data g. When no boundary data is explicitly given and r ≥ 1, we use the convention that vΩc ∣∣ e = vK ∣∣ e for vK given by the interior limit for v. If r = 0, we always let vK ∣∣ e be given by the interior limit for v and have to assign values for vΩc ∣∣ e appropriately. Typically we use the convention that vΩc ∣∣ e = vK ∣∣ e or T+ i (v) = T−i (v) with the understanding that vK ∣∣ e is given by the interior limit for v. More explicit values for vΩc ∣∣ e can be assigned based on context as in [12]. Notationally, the use of T±i denotes the lack of given Dirichlet data. Remark 2.6. (a) The trace operators T±i and T±,gi are nonstandard in that their values depend on the individual components of the edge normal ne. The standard definition used for LDG assigns a single-value (called a numerical flux) based on the edge normal vector as a whole. EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 67 (b) A labelling-independent definition can also be used so that T±i can be associated with the upwind or downwind direction with respect to the xi axis. The conventions are equivalent on a uniform Cartesian mesh using the natural ordering. See [25] for more details. 2.3.2. First order DG derivative operators. The main idea in [16] for defining DG derivative operators is to use the following local integration by parts formula for a given function v ∈ H1(Th) ∩ C0(Ω):∫ K vxiϕdx = ∫ ∂K vϕ(xI)ni ds− ∫ K vϕxi dx, i = 1, 2, . . . , d, K ∈ Th, (2.9) with test functions ϕ chosen from the DG space V h, ϕ(xI) denoting the limit from the interior of K, and ni denoting the i-th component of the unit outward normal vector for K. Thus, the DG (partial in xi) derivative intends to approximate the weak partial derivative vxi for all v ∈ H1(Th). To this end, the trace value v|∂K must be appropriately chosen/defined when v is not continuous. We define DG first order partial derivative operators ∂±xi,h for i = 1, 2, . . . , d as follows: for u ∈ H1(Th),∫ K ∂±xi,huφ dx ≡ ∫ ∂K T±i (u)φ(xI)n (i) K ds− ∫ K uφxi dx ∀φ ∈ V h (2.10) for all K ∈ Th. Notice that the “forward/backward” DG first order derivative oper- ators ∂±xi,h are different if the values of T±i (u) are different due to a discontinuity in u. It is easy to check that ∂±xi,h coincides with the FD operators δ±xi,hi on Cartesian grids when using the natural ordering (cf. [16]). Hence, the forward/backward DG derivative operators are indeed generalizations of the forward/backward difference operators to unstructured grids. When boundary trace data g is known, we de- fine the DG first order partial derivative operators ∂±,gxi,h that naturally enforce the boundary data by∫ K ∂±,gxi,h uφ dx = ∫ ∂K T±,gi (u)φ(xI)n (i) K ds− ∫ K uφxi dx ∀φ ∈ V h (2.11) for all K ∈ Th and i = 1, 2, . . . , d. Using the DG first order partial derivative operators as building blocks, we can define various central DG first order derivative operators and corresponding DG finite element gradient operators. Let i ∈ {1, 2, . . . , d}. Then, we define ∂xi,h ≡ 1 2 ( ∂+ xi,h + ∂−xi,h ) as a generalization of the central difference operator δxi,hi . If boundary data g is given, then we define the following two central DG first derivative operators that naturally enforce the boundary condition: ∂gxi,h ≡ 1 2 ( ∂+,g xi,h + ∂−,gxi,h ) , ∂ g xi,h ≡ 1 2 ( ∂gxi,h + ∂xi,h ) . The first operator ∂gxi,h naturally generalizes the form of the central difference op- erator δxi,hi when acting on a grid function with known boundary values, and the operator ∂ g xi,h generalizes the central difference operator δxi,hi in the sense that both correspond to antisymmetric matrices when vectorized with g = 0 (see [15] for the motivation for ∂ g xi,h). In general, we use ∂gxi,h when r = 0 and ∂ g xi,h when 68 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 r ≥ 1 to naturally enforce boundary conditions while appropriately weighting inte- rior degrees of freedom versus fixed boundary data. Corresponding finite element gradient operators ∇±h , ∇h, ∇±,gh , ∇gh, and ∇gh are naturally defined by letting all components be given by the appropriate DG first order partial derivative operator. For example, ∇±h ≡ ( ∂±x1,h , ∂±x2,h , . . . , ∂±xd,h )T . 2.3.3. Second order DG derivative operators. Similar to the finite difference (and to the classical calculus) construction, using first order DG derivative operators as the building blocks, we can easily define their high order extensions. Below we only define the second order operators. We also only define the operators that will be used directly in our framework for approximating fully nonlinear elliptic equations. We will only consider the case when boundary conditions correspond to Dirichlet boundary data. More information about Neumann boundary data can be found in [8, 16]. We first define the following one-sided second order DG partial derivatives: ∂µνxixj ,h ≡ ∂ ν xj ,h∂ µ xi,h , ∂µν,gxixj ,h ≡ ∂νxj ,h∂ µ,g xi,h µ, ν ∈ {+,−}, (2.12) where g corresponds to given Dirichlet boundary data. Then we define the eight “sided” d× d matrix-valued DG Hessian operators Dµν h ≡ [ ∂µνxixj ,h ]d i,j=1 , Dµν,g h ≡ [ ∂µν,gxixj ,h ]d i,j=1 , µ, ν ∈ {+,−}, (2.13) and the six central d× d matrix-valued DG Hessian operators D̂2 h ≡ 1 2 ( D+− h +D−+ h ) , D̂2,g h ≡ 1 2 ( D+−,g h +D−+,g h ) , (2.14a) D̃2 h ≡ 1 2 ( D++ h +D−−h ) , D̃2,g h ≡ 1 2 ( D++,g h +D−−,gh ) , (2.14b) D 2 h ≡ 1 2 ( D̂2 h + D̃2 h ) , D 2,g h ≡ 1 2 ( D̂2,g h + D̃2,g h ) (2.14c) that can be used when assuming the underlying method has reduced form as intro- duced below. Remark 2.7. It can be shown ([16]) that the above second order DG operators coincide with their corresponding FD operators on Cartesian grids. Moreover, it is easy to see that all of the DG operators defined above can be applied to any piecewise “nice” functions on Th including those in V h. 3. A narrow-stencil and g-monotone numerical framework In this section we formulate a general framework for both FD and DG methods that can be used to approximate fully nonlinear elliptic boundary value problems using narrow-stencil methods. We first introduce the ideas using FD methods. We then provide examples and extend the ideas to DG methods. 3.1. A narrow-stencil FD framework. The narrow-stencil FD schemes that we consider will all correspond to seeking a grid function Uα : NJ ′ → R such that F̂ [Uα,xα] = 0 for xα ∈ Th ∩ Ω, (3.1a) Uα = g(xα) for xα ∈ Th ∩ ∂Ω, (3.1b) ∆hUα = 0 for xα ∈ Sh ⊂ Th ∩ ∂Ω (3.1c) EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 69 for all α ∈ NJ, where F̂ [Uα,xα] ≡ F̂ ( D−−h Uα, D −+ h Uα, D +− h Uα, D ++ h Uα,∇hUα, Uα,xα ) . (3.2) Since the schemes only depend upon the discrete Hessian operators Dµν h and the discrete gradient operator ∇h, they are inherently narrow-stencil. The multiple Hessian operators are used to avoid the directional resolution approach used for monotone schemes that often lead to the use of wide-stencils. We also note that the auxiliary boundary condition (3.1c) needed to define the ghost points that arise when calculating D±±h Uα for nodes xα near the boundary could be generalized to setting ∆hUα = h(xα) for some bounded function h. A well chosen h can increase the accuracy of the underlying scheme by removing any boundary layer error associated with the auxiliary boundary condition. Such an h can be chosen using a refining process by solving various iterations of (3.1) with increasingly better chosen functions h based on the previous iteration. The main goal for this paper is to define sufficient conditions that F̂ can satisfy in order to guarantee the scheme (3.1) is admissible and convergent. The following definitions are adapted from the 1D definitions presented in [11]. Definition 3.1. (i) A function F̂ : ( Rd×d )4 × Rd × R× Ω→ R is called a numerical operator. (ii) A numerical operator F̂ is said to be consistent (with the differential oper- ator F ) if F̂ satisfies lim inf Pµν→P ;µ,ν=−,+ q→v,λ→v,ξ→x F̂ (P−−, P−+, P+−, P++,q, λ, ξ) ≥ F∗(P,v, v,x), lim sup Pµν→P ;µ,ν=−,+ q→v,λ→v,ξ→x F̂ (P−−, P−+, P+−, P++,q, λ, ξ) ≤ F∗(P,v, v,x), where F∗ and F ∗ denote, respectively, the lower and upper semi-continuous envelopes of F . (iii) A numerical operator F̂ is said to be generalized-monotone or g-monotone if F̂ (A++, B+−, B−+, A−−,q, v,x) ≤ F̂ (B++, A+−, A−+, B−−,q, w,x) for all Aµν , Bµν ∈ Rd×d; q ∈ Rd; v, w ∈ R; x ∈ Ω such that Bµν � Aµν and w ≥ v for all µ, ν ∈ {+,−}. A numerical operator F̂ is said to be uniformly g-monotone if there exists a constant κ∗ > 0 such that F̂ (P++, P+−, P−+, P−−,q, v,x) is increasing in P++, P−−, and v at a rate bounded below by κ∗ and decreasing in P+− and P−+ at a rate bounded above by −κ∗ using the partial ordering imposed by �. (iv) A numerical operator F̂ can be written in reduced form if there exists a function Ĝ : ( Rd×d )2 × R× Ω→ R such that F̂ (P−−, P−+, P+−, P++,v, v,x) = Ĝ(P̃ , P̂ ,v, v,x) for P̃ ≡ 1 2 (P−− + P++) and P̂ ≡ 1 2 (P−+ + P+−) for all P−−, P−+, P+−, P++ ∈ Rd×d; v ∈ Rd; v ∈ R; and x ∈ Ω. Remark 3.2. 70 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 (a) When F and F̂ are continuous, the definition of consistency can be simpli- fied to F̂ (P, P, P, P,v, v,x) = F (P,v, v,x) for all P ∈ Rd×d, v ∈ Rd, v ∈ R, and x ∈ Ω. (b) When F̂ is differentiable, g-monotonicity can be defined by requiring that the matrices ∂F̂ ∂P−− and ∂F̂ ∂P++ have all nonnegative entries, the matri- ces ∂F̂ ∂P−+ and ∂F ∂P+− have all nonpositive entries, and ∂F̂ ∂v is nonnegative. In other words, F̂ (↑, ↓, ↓, ↑, ·, ↑, ·). For a uniformly g-monotone numeri- cal operator, ∂F̂ ∂P−− � κ∗1d×d, ∂F̂ ∂P++ � κ∗1d×d, and ∂F̂ ∂v ≥ κ∗ while − ∂F̂ ∂P+− � κ∗1d×d and − ∂F̂ ∂P−+ � κ∗1d×d, where 1d×d denotes the matrix with all components equal to 1. (c) The g-monotonicity approach for narrow-stencil methods uses a component partial ordering instead of the SPD partial ordering for symmetric matrices. The approach also directly compares high order differences instead of look- ing directly at function values as in the standard monotonicity approach of Barles and Souganidis making it easier to design g-monotone schemes for a wide class of problems. The consistency of F̂ with F will allow the scheme to also take advantage of the SPD partial ordering associated with a proper elliptic operator. (d) To simplify notation, we will assume F̂ can be written in reduced form and write F̂ (P̃ , P̂ ,v, v,x) instead of introducing the new function Ĝ. A key tool for designing the g-monotone numerical operators in Section 3.2 is the introduction of a numerical moment as defined in [14]: Definition 3.3. Let A : RJ × Th → Rd×d and V be a given grid function. The discrete operator M : RJ → R defined by M [V,xα] ≡ A ( Vα,xα ) : ( D̃2 hVα − D̂2 hVα ) for all xα ∈ Th ∩ Ω is called a numerical moment operator. 3.2. Examples of g-monotone FD methods. We now introduce particular ex- amples of g-monotone FD methods that fulfill the structure assumptions of the narrow-stencil framework. The first method is the Lax-Friedrichs-like method pro- posed in [14] that uses both a numerical moment and a numerical viscosity (where the g-monotone definition could be extended for multiple discrete gradient argu- ments). The (general) method is defined by F̂ [Uα,xα] ≡ F ( D 2 hUα,∇hUα, Uα,xα ) +A(Uα,xα) : ( D̃2 hUα − D̂2 hUα ) − ~β(Uα,xα) · ( ∇+ hUα −∇ − hUα ) , where ~β : RJ × Th → Rd is a vector-valued function and −~β(Uα,xα) · ( ∇+ hUα − ∇−hUα ) is called a numerical viscosity. Note that the method is (globally) g- monotone and consistent using the framework above for the particular choices ~β = ~0 and A = σ1d×d for the constant σ > K/2, where K denotes the global Lipschitz constant of F with respect to the Hessian argument and 1d×d denotes the matrix with all entries equal to one. The choices ~β = ~0 and A = σ1d×d for σ > K/2 were the focus in the admissibility, stability, and convergence analysis in [14]. EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 71 In Section 7 we test the performance of the consistent FD method corresponding to the choice F̂γ,σ[Uα,xα] = F ( D 2 hUα,∇hUα, Uα,xα ) + (Mα + γId×d + σ1d×d) : ( D̃2 hUα − D̂2 hUα ) (3.3) for σ ≥ 0 and γ + σ ≥ 0, where [Mα]ij ≡ 1 2 ∣∣ ∂F ∂Pij ∣∣ (D2 hUα,∇hUα,Uα,xα) for all i, j = 1, 2, . . . , d using the convention that F = F (P,v, v,x) for P ∈ Rd×d, v ∈ Rd, v ∈ R, and x ∈ Ω. By construction, the method is only locally g-monotone in the sense that the linearization of the method is g-monotone. Furthermore, for σ > 0, the method is only locally uniformly g-monotone. The admissibility proof in Section 5 when applied to F̂γ,σ holds for γ ≥ 0 and σ = 0 and γ ≥ −σ with σ ≥ 0 if the problem is uniformly elliptic. If the operator F is only degenerate elliptic but globally Lipschitz, the proof would hold for either γ > 0 or σ > 0. Consider the linear problem A : D2u = f . Then, the method F̂0,0 is equivalent to F̂0,0[Uα,xα] = F ( D2 hUα,xα ) = A : D2 hUα − f(xα) for the discrete Hessian D2 h defined by [ D2 hUα ] ij = { δ̃2 xi,xj ;hi,hj Uα if aij(xα) > 0, δ̂2 xi,xj ;hi,hj Uα if aij(xα) ≤ 0 (3.4) for all i, j ∈ {1, 2, . . . , d} and xα ∈ Th ∩ Ω. Thus, the method takes the form of an upwinding-type method where, instead of matching the choice of the discrete partial derivative approximation to the advection field, we match the choice of the discrete second-order partial derivative approximation to the sign of the corresponding dif- fusion coefficient in A. For −A symmetric nonnegative definite, we would have aii ≤ 0 for all i = 1, 2, . . . , d. Consequently, when γ = σ = 0, the method would only have a nine-point stencil instead of a 13-point stencil in two-dimensions, and the auxiliary boundary condition ∆hUα = 0 would not be required. Similarly, the choice γ = −σ would also only have a nine-point stencil and would not require the auxiliary boundary condition. We would further reduce the stencil to only seven points if the diffusion coefficient a12 has a fixed sign. Thus, the methods based on choosing γ = σ = 0 or γ = −σ are of particular interest since the choices γ = −σ or γ ≥ 0 with σ = 0 represent limiting choices for enforcing the g-monotonicity of a numerical operator F̂ while remaining consistent with the PDE operator F . 3.3. A narrow-stencil DG framework. We can also naturally formulate narrow- stencil and g-monotone DG methods by seeking a piecewise polynomial function uh ∈ V h such that( F̂ [uh], ϕh ) Th + γB ∑ e∈EB 1 he 〈 uh − g, ϕh 〉 e + γI ∑ e∈EI 1 he 〈 [uh], [ϕh] 〉 e = 0 (3.5) for all ϕh ∈ V h, where γB , γI ≥ 0 and F̂ [uh] = F̂ ( D−−,gh uh, D −+,g h uh, D +−,g h uh, D ++,g h uh,∇ g huh, uh, · ) 72 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 is the same as the numerical operator used for FD methods but is now evaluated using DG derivatives. Notationally, he denotes the diameter of e and uh, ϕh are evaluated using interior limits over e ∈ EB when γB > 0. When r = 0, we use ∇gh instead of ∇gh to approximate the gradient operator since the corresponding trace operators naturally weight exterior limits versus interior limits when assuming only the exterior limit corresponds to g. We also set γB = γI = 0 when r = 0 to ensure consistency with the DG method and the underlying FD method on uniform Cartesian grids. Observe that the term ( F̂ [uh], ϕh ) Th corresponds to projecting the numerical operator F̂ [uh] into the discrete space V h using an L2 projection. For a quasi- uniform mesh, the penalty terms can be controlled using the uniform ellipticity assumption for F and the properties of the DWDG method for approximating Poisson’s equation derived in [25]. Consequently, we can choose γB = γI = 0 even when r ≥ 1. As such, the formulation for DG methods requires projecting the FD formulation into the discrete space and optionally adding penalization. We note that the auxiliary boundary condition is not required for r ≥ 1 based on the definitions of the DG derivative operators and, in particular, the way in which the boundary trace operators are defined. For r = 0, explicit rules for defining the exterior values for the boundary trace operators are provided in [12], and they are consistent with the Dirichlet boundary data and the auxiliary boundary condition (3.1c). Letting q±i = ∂±,gxi,h uh, the difficulty addressed in [12] is how to define q±i ∣∣ Ωc , which can be thought of as defining ghost points for the partial derivative with respect to xi when ni 6= 0. We refer the reader to [12] for the complete formulation when r = 0. Remark 3.4. (a) By construction, the proposed DG methods can be considered “narrow- stencil.” (b) When utilizing ∇gh instead of ∇gh, the DG method (3.5) is equivalent to the nonstandard LDG methods in [12] written in a compact form using the DG finite element calculus. For the unified framework we utilize ∇gh to more closely mimic the antisymmetric property of the FD operator δxi,hi as inspired by [15] where DG methods were formulated for approximating stationary Hamilton-Jacobi equations. (c) The DG method (3.5) is equivalent to the FD method (3.1) when Th is a uniform Cartesian mesh and the natural ordering is used. As such, all of the analytical results in Sections 4, 5, and 6 can be extended to (3.5) in this special case while, in general, the DG approach formally allows for higher degree bases and more general meshes. 4. Convergence analysis In this section we prove that consistent, stable, and g-monotone methods con- verge to the underlying viscosity solution of (1.1). Similar to the proof in [11], the result will assume the numerical operator can be written in reduced form. The theorem requires the methods defined by (3.1) are admissible in the sense that a solution Uα exists. The admissibility of the schemes will be verified in Section 5. We also use the definition in [14] that defines a piecewise constant extension uh for EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 73 a given grid function U ∈ S(T ′h) by uh(x) ≡ Uα, x ∈ Bα (4.1) for all α ∈ N′J , where Bα ≡ ∏ i=1,2,...,d ( xα − hi 2 ei,xα + hi 2 ei ] for all x ∈ Ω′ ≡ ∪α∈N′JBα ⊃ Ω. For transparency, we will only explicitly consider operators F that have the form F [u](x) = F ( D2u, u,x ) in (1.1a). We note that the proof can be readily extended to the more general case F [u](x) = F ( D2u,∇u, u,x ) using the techniques in [14]. Indeed, in the proof below, we would have ∇hk uhk (zk) → ∇ϕ(x0) in Case (i) which exploits the con- sistency of the scheme and (4.6) could be rewritten as F∗ ( D̂2 hk uhk (xk),∇hk uhk (xk), uhk (xk),xk ) − F∗ ( D2ϕ(x0),∇ϕ(x0), ϕ(x0),xk ) ≥ −λδ2 x`,h (k`) ` uhk (xk)− λ |ϕx`,x`(x0)| −K (|uhk (xk)|+ |ϕ(x0)|) −K d∑ i=1 ( |[∇+ hk uhk (xk)]i|+ |[∇−hk uhk (xk)]i|+ |ϕxi(x0)| ) −K d∑ i=1 d∑ j=1, (i,j)6=(`,`) ( |[D̂2 hk uhk (xk)]ij |+ |ϕxixj (x0)| ) so that the fact f` ( h ([k0]`) ` ) d∑ i=1 ( |[∇+ hk uhk (xk)]i|+ |[∇−hk uhk (xk)]i|+ |ϕxi(x0)| ) → 0 for min k sufficiently large and k` →∞ can be exploited in Case (ii). Theorem 4.1. Suppose the operator F in (1.1) is proper and uniformly elliptic with λ > 0, g is continuous on ∂Ω, F is Lipschitz continuous with respect to its first two arguments, and (1.1) satisfies the comparison principle. Suppose F̂ is consistent, is uniformly g-monotone, and can be written in reduced form, and suppose that F̂ is Lipschitz continuous with respect to its first three arguments when written in reduced form. Let U ∈ S(T ′h) be the solution to the scheme (3.1), and let uh denote the piecewise constant extension of U defined by (4.1). If (3.1) is admissible and `∞-norm stable, then uh converges to u locally uniformly as h→ 0+. Proof. The following is a sketch of the proof that highlights the differences from the complete convergence proof for the Lax-Friedrich’s-like method given in [14]. Step 1. Since the underlying FD scheme is assumed to be `∞-norm stable, there exists a constant C > 0 such that ‖uh‖L∞(Ω) ≤ C independent of h. Define the upper and lower semicontinuous functions u and u by u(x) ≡ lim sup h→0+, ξ→x uh(ξ), u(x) ≡ lim inf h→0+, ξ→x uh(ξ), where the limits are understood as multi-limits. We show u is a viscosity subsolution of (1.1). The proof that u is a viscosity supersolution of (1.1) is analogous. By the comparison principle, we have u = u, and it follows that u = u = u is the viscosity solution of (1.1). 74 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 Let ϕ ∈ C2(Ω) be a quadratic polynomial such that u − ϕ takes a strict local maximum at x0 ∈ Ω with u(x0) = ϕ(x0). Then there exists a ball, Br0(x0) ⊂ Rd, centered at x0 with radius r0 > 0 (in the `∞ metric) such that u(x)− ϕ(x) < u(x0)− ϕ(x0) = 0 ∀x ∈ ( Br0(x0) ∩ Ω ) \ {x0}. (4.2) Suppose x0 ∈ Ω. We show that F∗ ( D2ϕ(x0), ϕ(x0),x0 ) ≤ 0 (4.3) based on various cases determined by the regularity of u at x0. Note that if x0 ∈ ∂Ω, then, by the argument in [14], u can be shown to satisfy the boundary condition (1.1b) in the viscosity sense. By the definition of u and (4.2), there exists (maximizing) sequences {hk}, {xk}, and {zk} and a constant K0 > 0 such that hk → 0+, (4.4a) xk → x0 with xk ∈ Thk , (4.4b) uhk (xk)→ u(x0), (4.4c) zk → x0 with |x(ki) i − z(ki) i | ≤ 1 2 h (ki) i and uhk (zk) = uhk (xk), (4.4d) uhk (z)− ϕ(z) is locally maximized at z = zk for all min k ≥ K0. (4.4e) Let H(k) ∈ Rd×d be defined by H(k) = D2u(k)(zk) using the convention in [14] to define the local interpolation functions u(k). Case (i): {H(k)} has a uniformly bounded subsequence. In this case, there exists a symmetric matrix H ∈ Rd×d and a subsequence (not relabeled) such that H(k) → H, D̃2 hk uhk (zk) → H, and D̂2 hk uhk (zk) → H with D2ϕ − H symmetric positive semidefinite (see [14] for details). Thus, 0 = lim mink→∞ F̂ [uhk ,xk] = lim mink→∞ F̂ [uhk , zk] = lim mink→∞ F̂ ( D̃2 hk uhk (zk), D̂2 hk uhk (zk), uhk (zk), zk ) ≥ F∗ (H,ϕ(x0),x0) ≥ F∗ ( D2ϕ(x0), ϕ(x0),x0 ) by the consistency of the scheme and the ellipticity of F . Case (ii): {H(k)} does not have a uniformly bounded subsequence and there is no set of local interpolation functions ũhk such that the sequence H̃(k) has a bounded subsequence (see [14] for the definition of ũhk ). If such functions ũhk exist, then the argument in Case (i) can be easily updated to show F∗ ( D2ϕ(x0), ϕ(x0),x0 ) ≤ 0. There exists a pair of indices (i, j) such that the sequence [D̃2 hk uhk (zk)]ij or [D̂2 hk uhk (zk)]ij does not have a bounded subsequence. Thus, by (4.4e), there exists an index ` ∈ {i, j} and a subsequence such that δ2 ξji ,hk uhk (zk) → −∞, δ2 ηji ,hk uhk (zk)→ −∞, δ2 x`,2h (k`) ` uhk (zk)→ −∞, or δ2 x`,h (k`) ` uhk (zk)→ −∞ using the notation in [14] to rewrite the components of [D̃2 hk uhk (zk)]ij and [D̂2 hk uhk (zk)]ij in terms of central difference operators. Since F̂ is uniformly g-monotone, F̂ must be EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 75 uniformly increasing with respect to δ2 ξji ,hk uhk (zk) → −∞, δ2 ηji ,hk uhk (zk) → −∞, or δ2 x`,2h (k`) ` uhk (zk) → −∞. By sending ki → ∞ and kj → ∞ while ensuring min k is sufficiently large, if it holds δ2 ξji ,hk uhk (zk) → −∞, δ2 ηji ,hk uhk (zk) → −∞, or δ2 x`,2h (k`) ` uhk (zk) → −∞, then there must hold δ2 x`,h (k`) ` uhk (zk) → −∞ to en- sure F̂ [uhk ,xk] = 0 for all k. Therefore, there exists an index ` such that the sequence δ2 x`,h (k`) ` uhk (zk) = δ2 x`,h (k`) ` uhk (xk) does not have a bounded subsequence as min k→∞. Choose sequences {hk}, {xk} that maximize the rate at which δ2 x`,h (k`) ` uhk (xk)→ −∞. By the definition of the scheme, we have 0 = F̂ [uhk ,xk] = F̂ ( D̃2 hk uhk (xk), D̂2 hk uhk (xk), uhk (xk),xk ) = F∗ ( D2ϕ(x0), ϕ(x0),xk ) + F∗ ( D̂2 hk uhk (xk), uhk (xk),xk ) − F∗ ( D2ϕ(x0), ϕ(x0),xk ) + F̂ ( D̂2 hk uhk (xk), D̂2 hk uhk (xk), uhk (xk),xk ) − F∗ ( D̂2 hk uhk (xk), uhk (xk),xk ) + F̂ ( D̃2 hk uhk (xk), D̂2 hk uhk (xk), uhk (xk),xk ) − F̂ ( D̂2 hk uhk (xk), D̂2 hk uhk (xk), uhk (xk),xk ) . (4.5) Then, by the mean value theorem, the Lipschitz continuity of F , and the uniform and proper ellipticity of F , there exists a constant K ≥ 0 such that F∗ ( D̂2 hk uhk (xk), uhk (xk),xk ) − F∗ ( D2ϕ(x0), ϕ(x0),xk ) ≥ −λ ( δ2 x`,h (k`) ` uhk (xk)− ϕx`,x`(x0) ) −K (|uhk (xk)|+ |ϕ(x0)|) −K d∑ i=1 d∑ j=1, (i,j)6=(`,`) ( |[D̂2 hk uhk (xk)]ij |+ |ϕxixj (x0)| ) ≥ −λδ2 x`,h (k`) ` uhk (xk)− λ |ϕx`,x`(x0)| −K (|uhk (xk)|+ |ϕ(x0)|) −K d∑ i=1 d∑ j=1 (i,j) 6=(`,`) ( |[D̂2 hk uhk (xk)]ij |+ |ϕxixj (x0)| ) . (4.6) Using the consistency of the scheme, F̂ ( D̂2 hk uhk (xk), D̂2 hk uhk (xk), uhk (xk),xk ) − F∗ ( D̂2 hk uhk (xk), uhk (xk),xk ) ≥ 0. (4.7) Lastly, by the mean value theorem, the Lipschitz continuity of F̂ , and the g- monotonicity of F̂ , there exists a constant K̂ ≥ 0 and a sequence 0 ≤ ak ≤ K 76 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 such that F̂ ( D̃2 hk uhk (xk), D̂2 hk uhk (xk), uhk (xk),xk ) − F̂ ( D̂2 hk uhk (xk), D̂2 hk uhk (xk), uhk (xk),xk ) ≥ ak [ D̃2 hk uhk (xk)− D̂2 hk uhk (xk) ] `` − K̂ d∑ i=1 d∑ j=1, (i,j)6=(`,`) ( |[D̃2 hk uhk (xk)]ij |+ |[D̂2 hk uhk (xk)]ij | ) . (4.8) Plugging (4.6), (4.7), and (4.8) into (4.5), it follows that 0 ≥ F∗ ( D2ϕ(x0), ϕ(x0),xk ) − λδ2 x`,h (k`) ` uhk (xk) + ak [ D̃2 hk uhk (xk)− D̂2 hk uhk (xk) ] `` − λ|ϕx`,x`(x0)| −K (|uhk (xk)|+ |ϕ(x0)|) −K d∑ i=1 d∑ j=1, (i,j) 6=(`,`) ( |[D̂2 hk uhk (xk)]ij |+ |ϕxixj (x0)| ) − K̂ d∑ i=1 d∑ j=1, (i,j) 6=(`,`) ( |[D̃2 hk uhk (xk)]ij |+ |[D̂2 hk uhk (xk)]ij | ) . (4.9) Choose the corresponding optimal function f` (defined in [14]) and subsequences such that lim mink→∞ f` ( h (k`) ` ) δ2 xk` ,h (k`) ` uhk (xk) = −C` (4.10) for some constant C` > 0. Then, by [14], it holds lim inf mink→∞ f` ( h (k`) ` )[ D̃2 hk uhk (xk)− D̂2 hk uhk (xk) ] `` ≥ 0 implying f` ( h (k`) ` ) ak[D̃2 hk uhk (xk)− D̂2 hk uhk (xk)]`` ≥ − λ 4 C` (4.11) for all min k sufficiently large. Furthermore, by combining the observations in Sub- cases iia, iib, and iic in the proof of Theorem 6.1 in [14], there exists subsequences such that, for min k sufficiently large and k` >> max kj for all j 6= `, λ 4 C` ≥ f` ( h (k`) ` )( λ|ϕx`,x`(x0)|+K (|uhk (xk)|+ |ϕ(x0)|) ) +Kf` ( h (k`) ` ) d∑ i=1 d∑ j=1 (i,j)6=(`,`) ( |[D̂2 hk uhk (xk)]ij |+ |ϕxixj (x0)| ) + K̂f` ( h (k`) ` ) d∑ i=1 d∑ j=1 (i,j)6=(`,`) ( |[D̃2 hk uhk (xk)]ij |+ |[D̂2 hk uhk (xk)]ij | ) . (4.12) (Note that the primary difficulty in showing (4.12) is controlling the contributions of [D̃2 hk uhk (xk)]`,j and [D̂2 hk uhk (xk)]`,j for j 6= ` due to the 1 h (k`) ` factor when approximating mixed derivatives.) Thus, scaling (4.9) by f` ( h (k`) ` ) and plugging EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 77 in (4.10), (4.11), and (4.12), there exists an index k0 with min k0 sufficiently large such that 0 ≥ f` ( h ([k0]`) ` ) F∗ ( D2ϕ(x0), ϕ(x0),xk0 ) + 3λ 4 C` − λ 4 C` − λ 4 C` > f` ( h ([k0]`) ` ) F∗ ( D2ϕ(x0), ϕ(x0),xk0 ) . (4.13) The bound 0 ≥ F∗[ϕ](x0) follows since f` ( h ([k0]`) ` ) > 0 and xk → x0. Hence, (4.3) has been verified. The remainder of the proof is identical to Steps 4-6 in [14], and the result follows. � Remark 4.2. (a) g-monotonicity allowed us to identify a sufficiently positive term when λ > 0. Consequently, we could strongly exploit the uniformly elliptic structure of the PDE operator F . (b) Theorem 4.1 is proved under the assumption that the numerical scheme is admissible and `∞-norm stable. The remainder of the paper verifies sufficient conditions under which the assumptions hold. 5. Admissibility analysis The goal of this section is to show that the proposed narrow stencil scheme (3.1) has a unique solution whenever the numerical operator F̂ is consistent, is g-monotone, and can be written in reduced form. For transparency, we will only consider operators F that have the form F [u](x) = F ( D2u, u,x ) in (1.1a). The proofs can be adapted for the more general case using the techniques in [15] by exploiting the fact that the matrix representation of ∇h is anti-symmetric. The idea for proving the well-posedness is to equivalently reformulate the pro- posed scheme as a fixed point problem and to prove the mapping is contractive in the `2-norm. To this end, let S(T ′h) denote the space of all grid functions on T ′h, and introduce the mapping Mρ : S(T ′h)→ S(T ′h) defined by Û ≡MρU, (5.1) where the grid function Û ∈ S(T ′h) is defined by Ûα = Uα − ρF̂ [Uα,xα], if xα ∈ Th ∩ Ω, (5.2a) Ûα = g(xα), if xα ∈ Th ∩ ∂Ω, (5.2b) ∆hÛα = 0, if xα ∈ Sh (5.2c) for ρ > 0 an undetermined constant. Clearly, the iteration defined in (5.2) is the standard forward Euler method with pseudo time-step ρ complemented with a boundary condition consistent with (3.1). To show Mρ is a contraction, we will linearize the operator via the mean value theorem. As a preliminary result in Section 5.1, we will first consider a simple case when F is linear with constant-valued coefficients and a simple scheme based on D 2 h is used to discretize the Hessian. The general case will be considered in Section 5.2. We do require one additional structure assumption on the numerical operator F̂ to assist in the admissibility and stability proofs. The condition will ensure that the method based on using multiple Hessian operators is compatible with the 78 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 uniform ellipticity property of the corresponding PDE problem. We first motivate the property before defining it. Note that the Lax-Friedrich’s-like method in [14] and the examples in Section 3.2 satisfy the additional structure assumption. Suppose that F = F (P, v, x) is uniformly elliptic and differentiable with respect to its first two arguments and F̂ = F̂ (P̃ , P̂ , v, x) is differentiable with respect to its first three arguments. Then, if F̂ is consistent with F , it holds ∂F ∂P = ∂ ∂P F (P, v, x) = ∂ ∂P F̂ (P, P, v, x) = ∂F̂ ∂P̃ + ∂F̂ ∂P̂ , (5.3a) ∂F ∂v = ∂ ∂v F (P, v, x) = ∂ ∂v F̂ (P, P, v, x) = ∂F̂ ∂v (5.3b) for all P ∈ Rd×d, v ∈ R, and x ∈ Ω. Let A,B ∈ Rd×d, and suppose P = 1 2A+ 1 2B. Then, for the Lax-Friedrich’s-like method in [14] where F̂ (A,B, v, x) = F (1 2 A+ 1 2 B, v, x ) + γ1d×d : (A−B) for γ sufficiently large, it holds ∂ ∂P̃ F̂ (A,B, v, x) + ∂ ∂P̂ F̂ (A,B, v, x) = 1 2 ∂ ∂P F (P, v, x) + γ1d×d + 1 2 ∂ ∂P F (P, v, x)− γ1d×d = ∂ ∂P F (P, v, x) ≤ −λI, and for the g-monotone method F̂ (A,B, v, x) = F (B, v, x) + γ1d×d : (A−B), it holds ∂ ∂P̃ F̂ (A,B, v, x) + ∂ ∂P̂ F̂ (A,B, v, x) = γ1d×d + ∂ ∂P F (B, v, x)− γ1d×d = ∂ ∂P F (B, v, x) ≤ −λI. Note that ∂ ∂P F (B, v, x) may not equal ∂ ∂P F (P, v, x) for P 6= B. However, the same uniform ellipticity bound holds. Similarly, for the g-monotone method F̂ (A,B, v, x) = F (A, v, x) + γ1d×d : (A−B), it holds ∂ ∂P̃ F̂ (A,B, v, x) + ∂ ∂P̂ F̂ (A,B, v, x) = ∂ ∂P F (A, v, x) + γ1d×d − γ1d×d = ∂ ∂P F (A, v, x) ≤ −λI, and again the same uniform ellipticity bound holds. Thus, we assume the following compatibility condition when proving the admissibility and stability of our proposed narrow-stencil schemes. Definition 5.1. Suppose F = F (P, v,x) is proper elliptic and differentiable with respect to its first two arguments with ∂ ∂P F (A,w,x) ≤ −λI and ∂ ∂vF (A,w,x) ≥ κ0 for all A ∈ Sd×d, w ∈ R, and x ∈ Ω. Suppose the numerical operator F̂ is consistent with F and differentiable with respect to its first three arguments. The numerical EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 79 operator F̂ = F̂ (P̃ , P̂ , v,x) is elliptic compatible if there exists a constant c > 0 independent of h such that ∂ ∂P̃ F̂ (A,B,w,x) + ∂ ∂P̂ F̂ (A,B,w,x) ≤ −cλI and ∂ ∂v F̂ (A,B,w,x) ≥ cκ0 for all A,B ∈ Sd×d, w ∈ R, and x ∈ Ω. 5.1. Admissibility of a simple method for linear, constant coefficient PDEs. Consider the linear elliptic boundary value problem L[u] ≡ −A : D2u = − d∑ i=1 d∑ j=1 aijuxixj = f in Ω, (5.4a) u = g on ∂Ω, (5.4b) where A is constant-valued and symmetric positive definite, and consider the simple FD scheme corresponding to finding a grid function Uα : NJ ′ → R such that LhUα ≡ −A : D 2 hUα = f(xα) for xα ∈ Th ∩ Ω, (5.5a) Uα = g(xα) for xα ∈ Th ∩ ∂Ω, (5.5b) ∆hUα = 0 for xα ∈ Sh ⊂ Th ∩ ∂Ω (5.5c) for all α ∈ NJ. We show that (5.5) is equivalent to solving a linear system L~U = b with the matrix L symmetric positive definite. Let λ0 > 0 denote the smallest eigenvalue of A. Define A0 ≡ A−λ0I. Then, A0 is symmetric nonnegative definite. Thus, there exists an eigenvalue decomposition A0 = QΛQT = ∑d k=1 λkqkq T k , where λk ≥ 0 and {qk} forms an orthonormal basis for Rd. Define constants b (k) i ∈ R by qk = ∑d i=1 b (k) i ei, and observe that qkq T k = d∑ i=1 ( b (k) i )2 eie T i + d∑ i=1 d∑ j=1, j 6=i b (k) i b (k) j eie T j . Thus, A0 = d∑ k=1 λk [ d∑ i=1 ( b (k) i )2 eie T i + d∑ i=1 d∑ j=1, j 6=i b (k) i b (i) j eie T j ] = d∑ i=1 [ d∑ k=1 λk ( b (k) i )2] eie T i + d∑ i=1 d∑ j=1, j 6=i [ d∑ k=1 λkb (k) i b (k) j ] eie T j , and it follows that [A0]ii = aii − λ0 = d∑ k=1 λk ( b (k) i )2 , [A0]ij = aij = d∑ k=1 λkb (k) i b (k) j for all i, j = 1, 2, . . . , d with j 6= i. Let J0 = |Th ∩ Ω|, Di ∈ RJ0×J0 denote the matrix representation of δxi,hi with the boundary condition (5.5b) and M ∈ RJ0×J0 denote the matrix representation of −∆2h and the boundary conditions (5.5b) and (5.5c). Then (Di) T = −Di and 80 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 there exists diagonal matrices Bi ∈ RJ0×J0 for i = 1, 2, . . . , d with nonnegative components such that M = − d∑ i=1 (DiDi −Bi) ≥ − d∑ i=1 DiDi = d∑ i=1 (Di) T Di. The positive components of Bi correspond to nodes xα near the boundary and increase the coefficients for Uα. Such corrections are needed in the matrix form to account for the values of δxi,hiUα′ when xα′ ∈ Th ∩ ∂Ω and for the implementation of the auxiliary boundary condition (5.5c). Indeed, suppose xα − hiei ∈ ∂Ω. Then δ2 xi,2hiUα = δxi,hiδxi,hiUα = 1 2hi δxi,hi (Uα+ei − Uα−ei) = 1 2hi δxi,hiUα+ei − 1 2hi δxi,hiUα−ei = 1 2hi δxi,hiUα+ei − 1 4h2 i Uα + 1 4h2 i g(xα − hiei). However, when computing using the matrix representation, Di treats the boundary value Uα−ei as a known value in its representation and removes it when calculating DiU . Consequently, the second application of Di does not act on the boundary node leading to a smaller coefficient for the adjacent interior node involved in the calculation of δxi,hiUα−ei . Thus, M would contain the contribution 1 4h2 i to the coefficient for Uα while −DiDi would not. We can see that (5.5c) ensures that the ghost value Uα±2ei satisfies 1 2h2 i Uα±2ei = − 1 2h2 i Uα + 1 h2 i g(xα±ei)− d∑ j=1, j 6=i δ2 xj ,hjg(xα±ei) for xα±ei ∈ Th ∩ ∂Ω, where 1 2h2 i Uα±2ei is directly involved in the computation of δ 2 xi,hiUα. Lastly, note that correction terms are not needed when considering the relationship of δxi,hiδxj ,hj to DiDj for i 6= j since the computation of δxi,hiUα′ would only include boundary nodes whenever xα ± hjej ∈ Th ∩ ∂Ω. Observe that −A0 : D 2 h = − d∑ i=1 d∑ j=1 [A0]ijδxi,hiδxj ,hj = − d∑ i=1 (aii − λ0)δxi,hiδxi,hi − d∑ i=1 d∑ j=1, j 6=i aijδxi,hiδxj ,hj = − d∑ i=1 d∑ k=1 λk ( b (k) i )2 δxi,hiδxi,hi − d∑ i=1 d∑ j=1, j 6=i [ d∑ k=1 λkb (k) i b (k) j ] δxi,hiδxj ,hj = − d∑ k=1 λk [ d∑ i=1 ( b (k) i )2 δxi,hiδxi,hi + d∑ i=1 d∑ j=1, j 6=i b (k) i b (k) j δxi,hiδxj ,hj ] = − d∑ k=1 λk ([ d∑ i=1 b (k) i δxi,hi ][ d∑ i=1 b (k) i δxi,hi ]) , EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 81 and it follows that 0J0×J0 ≤ d∑ k=1 λk ([ d∑ i=1 b (k) i Di ]T [ d∑ i=1 b (k) i Di ]) = d∑ i=1 d∑ k=1 λk ( b (k) i )2 DT i Di + d∑ i=1 d∑ j=1, j 6=i [ d∑ k=1 λkb (k) i b (k) j ] DT i Dj = d∑ i=1 (aii − λ0)DT i Di + d∑ i=1 d∑ j=1, j 6=i aijD T i Dj = d∑ i=1 d∑ j=1 [A0]ijD T i Dj . Using the fact that DT i = −Di, L ≡ − d∑ i=1 d∑ j=1 aijDiDj + d∑ i=1 aiiBi (5.6) = −λ0 d∑ i=1 ( D2 i −Bi ) − d∑ i=1 d∑ j=1 [A0]ijDiDj + d∑ i=1 (aii − λ0)Bi (5.7) = λ0M + d∑ i=1 d∑ j=1 [A0]ijD T i Dj + d∑ i=1 (aii − λ0)Bi (5.8) ≥ λ0M > 0J0×J0 . (5.9) Therefore, (5.5) has a unique solution and the matrix representation L~U = b yields a symmetric positive definite matrix L. 5.2. Admissibility for fully nonlinear PDEs. To show that the mapping Mρ has a unique fixed point in S(T ′h), we first establish a lemma that specifies conditions under which Mρ is a contraction in `2. The proof will assume F is differentiable; however, the assumption is for ease of notation and the proof can be extended for F Lipschitz but not differentiable. The result will utilize the following result found in [14]: Lemma 5.2. Let B,F ∈ RJ×J such that B is symmetric nonnegative definite and F is symmetric positive definite. Define R ∈ RJ×J such that R is upper triangular and F = R∗R. Then ‖σI − FB‖2 ≤ σ for all positive constants σ such that σI > RBR∗. Lemma 5.3. Suppose the operator F in (1.1) is proper and uniformly elliptic, dif- ferentiable, and Lipschitz continuous with respect to its first two arguments. Suppose F̂ is consistent, g-monotone, can be written in reduced form, and is differentiable with respect to its first three arguments. Choose U, V ∈ S(T ′h) that satisfy the boundary conditions (5.5b) and (5.5c), and let Û =MρU and V̂ =MρV for Mρ defined by (5.1) and (5.2). Then, for F̂ elliptic compatible, it holds ‖Û − V̂ ‖`2(Th) ≤ ( 1− ρcλκ 4 − ρck0 4 ) ‖U − V ‖`2(Th) 82 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 for all ρ > 0 sufficiently small, where c is from the definition of elliptic compatibility, 4 > ρcλκ+ρck0, ∂F ∂D2u ≥ λI, ∂F∂u ≥ k0, and M ≥ κI for M the matrix representation of −∆2h. Proof. Let W ≡ V − U and Ŵ ≡ V̂ − Û . Then, by the boundary conditions, Ŵα = Wα = 0 for all xα ∈ Th∩∂Ω and ∆hŴα = ∆hWα = 0 for all xα ∈ Sh. Thus, by the mean value theorem for F̂ = F̂ (P̃ , P̂ , v,x), it holds Ŵα = Wα − ρ ( F̂ [Vα,xα]− F̂ [Uα,xα] ) = Wα − ρ ( F̂ (D̃2 hVα, D̂ 2 hVα, Vα,xα)− F̂ (D̃2 hUα, D̂ 2 hUα, Uα,xα) ) = Wα − ρF̂ ( D 2 hVα + 1 2 D̃2 hVα − 1 2 D̂2 hVα, D 2 hVα − 1 2 D̃2 hVα + 1 2 D̂2 hVα, Vα,xα ) + ρF̂ ( D 2 hUα + 1 2 D̃2 hUα − 1 2 D̂2 hUα, D 2 hUα − 1 2 D̃2 hUα + 1 2 D̂2 hUα, Uα,xα ) = ( 1− ρ∂F̂ ∂v ) Wα − ρ d∑ i=1 d∑ j=1 ( ∂F̂ ∂P̃ij + ∂F̂ ∂P̂ij ) δ 2 xi,xj ;hi,hjWα − ρ1 2 d∑ i=1 d∑ j=1 ( ∂F̂ ∂P̃ij − ∂F̂ ∂P̂ij ) δ̃2 xi,xj ;hi,hjWα + ρ 1 2 d∑ i=1 d∑ j=1 ( ∂F̂ ∂P̃ij − ∂F̂ ∂P̂ij ) δ̂2 xi,xj ;hi,hjWα (5.10) for all xα ∈ Th ∩ Ω. Let J0 = |Th∩Ω| and Ŵ,W ∈ RJ0 denote the vectorization of the grid functions Ŵ and W restricted to Th∩Ω, respectively. We introduce several matrix operators in RJ0×J0 that act on W and correspond to the FD operators with the boundary data naturally incorporated directly into the definition of the matrix. Then, using the notation of Section 2.2.4, D̃ij,0 and D̂ij,0 denote the matrix representations of δ̃2 xi,xj ;hi,hj and δ̂2 xi,xj ;hi,hj , respectively. Similarly, Dij,0 denotes the matrix repre- sentations of δ 2 xi,xj ;hi,hj . We also let F̂0 denote the diagonal matrix corresponding to the nodal values of ∂F̂ ∂v , F̃ij denote the diagonal matrix corresponding to ∂F̂ ∂P̃ij , and F̂ij denote the diagonal matrix corresponding to ∂F̂ ∂P̂ij . Then F̂0, F̃ij , and −F̂ij are all nonnegative definite, and we have (5.10) becomes Ŵ = (I − ρF̂0)W − ρ d∑ i=1 d∑ j=1 ( F̃ij + F̂ij ) Dij,0W − ρ1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij ) D̃ij,0W + ρ 1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij ) D̂ij,0W. EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 83 Letting M denote the matrix corresponding to −∆2h, it follows that Ŵ = (I − ρF̂0)W − ρ1 2 cλMW − ρ [ d∑ i=1 d∑ j=1 ( F̃ij + F̂ij ) Dij,0 − 1 2 cλM ] W − ρ1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij )( D̃ij,0 −Dij,0 ) W + ρ 1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij )( D̂ij,0 −Dij,0 ) W. (5.11) Choose xα ∈ Th ∩ Ω. Let A(xα) ∈ Rd×d be defined by Aij(xα) = ∂F̂ ∂P̃ij + ∂F̂ ∂P̂ij . Then, −A(xα) ≥ cλI by the elliptic compatible condition. Observe that( A(xα)− 1 2 cλI ) : D 2 hWα = d∑ i=1 d∑ j=1 ( ∂F̂ ∂P̃ij + ∂F̂ ∂P̂ij ) δ 2 xi,xj ;hi,hjWα − 1 2 cλ d∑ i=1 δ 2 xi,hiWα = d∑ i=1 d∑ j=1 ( ∂F̂ ∂P̃ij + ∂F̂ ∂P̂ij ) δ 2 xi,xj ;hi,hjWα − 1 2 cλ∆2hWα. Then, the term [ ∑d i=1 ∑d j=1(F̃ij + F̂ij)Dij,0 − 1 2cλM ]W in (5.11) is the matrix representation of ( A(xα)− 1 2cλI ) : D 2 hWα. We next utilize the frozen coefficient technique. Define the matrices Lα ∈ RJ0×J0 as the matrix representations of ( A(xα) − 1 2cλI ) : D 2 h for all xα ∈ Th ∩ Ω. Then, by Section 5.1, we have Lα is symmetric positive definite for all multi-indices α since the coefficient matrix is constant-valued for each fixed value of α. Define Ek ∈ RJ0×J0 by [Ek]ij = 1 only if i = j = k and 0 otherwise. Notationally, we let α(k) be the multi-index corresponding to the single-index k. Then, (5.11) can be rewritten as Ŵ = (I − ρF̂0)W − ρ1 2 cλMW − ρ J0∑ k=1 EkLα(k)W − ρ1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij )( D̃ij,0 −Dij,0 ) W − ρ1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij )( Dij,0 − D̂ij,0 ) W ≡ GW for the iteration matrix G ∈ RJ0×J0 defined by G ≡ (I − ρF̂0)− ρ1 2 cλM − ρ J0∑ k=1 EkLα(k) − ρ1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij )( D̃ij,0 −Dij,0 ) 84 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 − ρ1 2 d∑ i=1 d∑ j=1 ( F̃ij − F̂ij )( Dij,0 − D̂ij,0 ) . Choose ε > 0. Let N = 2 + J0 + 2d2. Observe that Ek + εI and F̃ij − F̂ij + εI are symmetric positive definite for all k and all i, j. Choose K such that 0J0×J0 ≤∑J0 k=1 Lα(k) ≤ KI, 0J0×J0 ≤ D̃ij,0 −Dij,0 ≤ KI, and 0J0×J0 ≤ Dij,0 − D̂ij,0 ≤ KI for all i, j, where the first bound uses the fact that Lα(k) is symmetric positive definite for all indices k and the other bounds follow from Lemma 2.5. Then, by Lemma 5.2, it holds ‖G‖2 ≤ 1 2 + ‖ 1 2N I − ρF̂0‖2 + ‖ 1 2N I − ρ1 2 cλM‖2 + J0∑ k=1 ‖ 1 2N I − ρ(Ek + εI)Lα(k)‖2 + ρε‖ J0∑ k=1 Lα(k)‖2 + 1 2 d∑ i=1 d∑ j=1 ‖ 1 N I − ρ ( F̃ij − F̂ij + εI )( D̃ij,0 −Dij,0 ) ‖2 + 1 2 d∑ i=1 d∑ j=1 ‖ 1 N I − ρ ( F̃ij − F̂ij + εI )( Dij,0 − D̂ij,0 ) ‖2 + ρε 1 2 d∑ i=1 d∑ j=1 ‖D̃ij,0 −Dij,0‖2 + ρε 1 2 d∑ i=1 d∑ j=1 ‖Dij,0 − D̂ij,0‖2 ≤ 1 2 + 1 2N − ρck0 + 1 2N − ρ1 2 cλκ+ J0∑ k=1 1 2N + d2 1 N + ρε(d2 + 1)K < 1− ρ3 4 ck0 − ρ 1 4 cλκ for all ρ > 0 sufficiently small and ε < c k0+λκ 4(d2+1)K chosen independently of ρ. The bound ‖Ŵ‖`2(Th∩Ω) ≤ ( 1− ρcλκ 4 − ρck0 4 ) ‖W‖`2(Th∩Ω) follows since ‖Ŵ‖2 ≤ ‖G‖2‖W‖2, and the bound over Th follows since Ŵα = Wα = 0 over Th ∩ ∂Ω. The proof is complete. � As an immediate corollary to Lemma 5.3, we have the following well-posedness result by the contractive mapping theorem. Theorem 5.4. Suppose the operator F in (1.1) is proper and uniformly elliptic, differentiable, and Lipschitz continuous with respect to its first two arguments. Sup- pose F̂ is consistent, is g-monotone, can be written in reduced form, is differentiable with respect to its first three arguments, and is elliptic compatible. The scheme (3.1) for approximating problem (1.1) has a unique solution. Remark 5.5. We emphasize that the properties of F and F̂ guarantee a unique solution whenever k0 > 0 or λ > 0. EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 85 6. Stability analysis For transparency and consistency with the results in Section 5, assume the op- erator F in (1.1a) has the form F [u](x) = F ( D2u, u,x ) . Theorem 6.1. Suppose the operator F in (1.1) is proper and uniformly elliptic and is Lipschitz continuous and differentiable with respect to its first two arguments with λ > 0 or k0 > 0. Suppose F̂ is consistent, is g-monotone, can be written in reduced form, is Lipschitz continuous and differentiable with respect to its first three arguments, and is elliptic compatible. Then the solution U to the scheme (3.1) for approximating problem (1.1) is `2-norm stable in the sense that( ∏ i=1,2,...,d hi1/2 ) ‖U‖`2(Th∩Ω) ≤ C, where C is a positive h-independent constant which depends on Ω, the lower (proper) ellipticity constants k0 and λ, F , and g. Proof. Define the function v ∈ C0(Ω) ∩H2(Ω) to be the solution to −∆v = 0 in Ω, (6.1a) v = g on ∂Ω, (6.1b) and define V : T ′h → R by Vα = v(xα) for all xα ∈ Th ∩ Ω and introduce the ghost points so that ∆hVα = 0 for all xα ∈ Sh. Then, by the mean value theorem, there exists a linear operator Lh such that F̂ [Uα,xα]− F̂ [Vα,xα] = Lh[Uα − Vα]. Furthermore, since U is a solution to (3.1), F̂ [Uα,xα]− F̂ [Vα,xα] = 0− F̂ [Vα,xα] = −F̂ [Vα,xα]. Thus, Uα − Vα is a solution to Lh[Uα − Vα,xα] + F̂ [Vα,xα] = 0 for xα ∈ Th ∩ Ω, (6.2a) Uα − Vα = 0 for xα ∈ Th ∩ ∂Ω, (6.2b) ∆h(Uα − Vα) = 0 for xα ∈ Sh ⊂ Th ∩ ∂Ω, (6.2c) where Lh has the form Lh ≡ Ah : D̃2 h −Bh : D̂2 h + ch for Ah, Bh ∈ Sd×d with all nonnegative components and ch ∈ R with ch ≥ 0 by the g-monotonicity of F̂ . Furthermore, Lh is consistent with the linear elliptic boundary value problem Lw ≡ −A : D2w + cw = −F [v], in Ω, (6.3a) u = g, on ∂Ω (6.3b) with the matrix A symmetric positive definite with λI ≤ A ≤ ΛI and c ≥ k0, where Ah − Bh ≤ −c0λI and ch ≥ c0k0 for c0 the constant based on the elliptic compatibility of F̂ . Applying Theorem 5.4, Uα−Vα is the unique solution to (6.2). Furthermore, by Lemma 5.3 and using the technique in [14, Theorem 4.3], it holds ‖U − V ‖`2(Th∩Ω) ≤ 4 c0λκ+ c0k0 ‖Lh[0] + F̂ [Vα,xα]‖`2(Th∩Ω) 86 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 = 4 c0λκ+ c0k0 ‖F̂ [Vα,xα]‖`2(Th∩Ω). Therefore, ‖U‖`2(Th∩Ω) ≤ ‖V ‖`2(Th∩Ω) + 4 c0λκ+ c0k0 ‖F̂ [Vα,xα]‖`2(Th∩Ω) and the result follows by the properties of v. The proof is complete. � Extending the techniques above and following the proofs in Section 5 of [14], we can also prove the following results. Lemma 6.2. Let i ∈ {1, 2, . . . , d}. Suppose the operator F in (1.1) is proper and uniformly elliptic, differentiable, and Lipschitz continuous with respect to its first two arguments. Suppose F̂ is consistent, g-monotone, can be written in reduced form, is differentiable with respect to its first three arguments, and is elliptic com- patible. Choose U, V ∈ S(T ′h) that satisfy the boundary conditions (5.5b) and (5.5c), and let Û =Mρ,2U and V̂ =Mρ,2V , where Ŵ ∈ S(T ′h) is defined by Ŵ =Mρ,2W for some W ∈ S(T ′h) if −δ2 xi,2hiŴα = −δ2 xi,2hiWα − ρF̂ [Wα,xα] ∀xα ∈ Th ∩ Ω, Ŵα = g(xα) ∀xα ∈ Th ∩ ∂Ω, ∆hŴα = 0 ∀xα ∈ Sh. Then it holds ‖δ2 xi,2hi(Û − V̂ )‖`2(Th∩Ω) ≤ ( 1− ρcλκ 4 ) ‖δ2 xi,2hi(U − V )‖`2(Th∩Ω) for all ρ > 0 sufficiently small, where 4 > ρcλκ, ∂F ∂D2u ≥ λI, ∂F ∂u ≥ k0, and M ≥ κI for M the matrix representation of −∆2h. Theorem 6.3. Suppose the operator F in (1.1) is proper and uniformly elliptic, and is Lipschitz continuous and differentiable with respect to its first two arguments with λ > 0. Suppose F̂ is consistent, is g-monotone, can be written in reduced form, is Lipschitz continuous and differentiable with respect to its first three arguments, and is elliptic compatible. Then the solution U to the scheme (3.1) for approximating problem (1.1) satisfies( ∏ i=1,2,...,d hi1/2 ) ‖δ2 xi,2hiU‖`2(Th∩Ω) ≤ C for all i = 1, 2, . . . , d, where C is a positive h-independent constant which depends on Ω, the lower (proper) ellipticity constant λ, F , and g. Lastly, by [14], combining Theorems 6.1 and 6.3 using a novel numerical embed- ding technique yields the following `∞ stability result. Theorem 6.4. Under the assumptions of Lemma 6.2, the numerical solution U is stable in the `∞-norm for d ≤ 3; that is, U satisfies ‖U‖`∞(Th) ≤ C for d ≤ 3, where C is a positive constant independent of h. EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 87 7. Numerical experiments In this section we test the performance of the FD method based on the numerical operator F̂γ,σ defined by (3.3) for σ ≥ 0 and γ + σ ≥ 0. We consider test problems based on choosing F to be linear and uniformly elliptic with non-divergence form, the Hamilton-Jacobi-Bellman operator, the Monge-Ampère operator, and the oper- ator coming from the equation of prescribed Gauss curvature. We will see that for the linear problem, F̂0,0 yields a nonsingular sparse matrix. However, the nonlinear PDE problems considered are degenerate, and, consequently, do not conform to all of the assumptions used in the convergence, admissibility, and stability analysis. Observationally, the nonlinear solver ‘fsolve’ in MATLAB has trouble in the degen- erate cases finding a zero when using the zero function as the initial guess and σ and γ are small. For the Hamilton-Jacobi-Bellman problem, we can successfully solve the limiting cases by forming a sequence of approximations for decreasing values of γ starting with an initial large value for γ when initializing fsolve with the zero function. For the Monge-Ampère operator and in the equation of prescribed Gauss curvature, we need to form the sequence of approximations with σ large to more strictly impose g-monotonicity with respect to all of the components of the discrete Hessian operators. Once fsolve is in a neighborhood of the solution we can success- fully find a zero for all σ ≥ 0 and γ ≥ −σ. For the Monge-Ampére problem and the prescribed Gauss curvature problem, we choose an exact solution for which the Monge-Ampère operator is locally uniformly elliptic. It is worth noting that fsolve does not find a false solution but instead always reported no solution found when staring with a poor initial guess. This is in direct contrast to the experiment in [9] that found false solutions to the Monge-Ampère problem when using the standard nine-point finite difference formula for approximating the discrete Hessian. Additional numerical tests for g-monotone FD methods can be found in [24, 11, 13] and tests for the corresponding DG scheme can be found in [24, 12]. 7.1. Test 1: Linear F with non-aligned grids. In this test, we consider the linear uniformly elliptic problem F [u] = −A : D2u−f for a discontinuous coefficient matrix A and uniform grids chosen such that no monotone finite difference method exists (see [26]). The matrix A(x) will be uniformly symmetric positive definite, and f will be a uniformly bounded function. We form a non-singular linear system LU = F that is solved using MATLAB’s backslash command. The matrix L is formed using sparse storage. Let Ω = (−1, 1)2. We form a sequence of uniform grids with Nx = Ny = 10, 40, 80, 120, 180, 240, 300. Let hx = hy = 2 301 based on the finest mesh, define the directions v(i) ∈ R2 by v(1) = [ 1 hy 2 ] , v(2) = [ hx 5hy ] , v(3) = [ 10hx hy ] , v(4) = [ hx 2 2 ] , and define the unit length vectors q (i) 1 = 1 ‖v(i)‖2 v(i) for all i = 1, 2, 3, 4. Note that, by construction, the directions v(1) and v(4) do not align with grid points in any of the meshes while v(2) aligns with nodes at least 5 layers away and v(3) aligns with nodes at least 10 layers away. For each i = 1, 2, 3, 4, choose unit length vectors q (i) 2 that are orthogonal to q (i) 1 . Define the orthogonal matrices Qi = [q (i) 1 ,q (i) 2 ] for each 88 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 i = 1, 2, 3, 4, and define the diagonal matrix Λ(x, y) ∈ R2×2 by Λ(x, y) = [ 2− sin (e5x) cos (e−3y) 0 0 2− sign ( cos (6πx) sin (6πy) )] for all (x, y) ∈ Ω. Then, Λ11 ∈ [1, 3] oscillates several times and Λ22 ∈ {2, 3, 4} has several discontinuities over Ω. Finally, we define A(x, y) =  Q1 Λ(x, y)QT1 , for x ≥ 0, y ≥ 0, Q2 Λ(x, y)QT2 , for x < 0, y ≥ 0, Q3 Λ(x, y)QT3 , for x < 0, y < 0, Q4 Λ(x, y)QT4 , otherwise so that A is strictly symmetric positive definite and discontinuous. No monotone method exists for this problem on the specified grids due to the choices for Q1 and Q4, and Q2 and Q3 would lead to wider-stencils. We consider the problem −A : D2u = f for two different solutions u(1)(x, y) = sin (π 2 (x+ y)2 ) , u(2)(x, y) = x3 18 ( 3 log(x2)− 11 ) + ( y − 1 2 )8/3√|x+ 1 5 |5 so that u(1) ∈ C∞(Ω) and u(2) ∈ C2(Ω) \ C3(Ω). The source function f and boundary data g are chosen so that the solution is given by either u(1) or u(2). Using F̂0,0, we observe optimal second order rates when approximating u(1) and expected deteriorated rates when approximating u(2) in Table 1. Table 1. Approximation results for Test 1 using F̂0,0. u(1) u(2) h `∞ Error Order `∞ Error Order 1.82e-01 5.77e-02 5.21e-03 4.88e-02 3.86e-03 2.06 8.29e-04 1.40 2.47e-02 9.79e-04 2.01 3.26e-04 1.37 1.65e-02 4.40e-04 1.99 1.86e-04 1.40 1.10e-02 1.96e-04 2.00 1.05e-04 1.41 8.30e-03 1.11e-04 2.00 7.02e-05 1.42 6.64e-03 7.10e-05 2.00 5.11e-05 1.43 7.2. Test 2: Hamilton-Jacobi-Bellman equations. This example is adapted from [32], and it was considered in [14] using the Lax-Friedrich’s-like method. Let Λ = [0, π/3]× SO(2), where SO(2) is the set of 2× 2 rotation matrices and define σθ by σθ ≡ RT [ 1 sin(φ) 0 cos(φ) ] , θ = (φ,R) ∈ Λ. Consider the Hamilton-Jacobi-Bellman equation F [u] ≡ inf θ∈Λ ( Lθu− fθ ) = 0, EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 89 Lθu = −Aθ(x) : D2u+ bθ(x) · ∇u+ cθ(x)u with Ω = (0, 1)2, Aθ = 1 2σ θ ( σθ )T , βθ = ~0, cθ = π2, fθ = √ 3 sin2(φ/π2) + g for g chosen independent of θ, and Dirichlet boundary data chosen such that the exact solution is given by u(x, y) = exy sin(πx) sin(πy). The optimal controls vary significantly throughout the domain and the corresponding diffusion coefficient is not diagonally dominant in parts of Ω. We approximate u using F̂γ,0 for γ = 1000, 10, 1, 0 in Table 2. The case γ = 1000 appears to be in a pre-asymptotic regime with the error dominated by the numerical moment term. Otherwise, as expected, we see that the methods are around second order accuracy with the approximations increasing in accuracy as γ decreases. In the implementation of F̂γ,σ, we do not calculate derivatives of F to define Mα. Instead, we define F̂ θγ,σ by F̂ θγ,σ[Uα] ≡ Aθ(xα) : D2,θ h Uα + cθ(xα)Uα − fθ(xα) + (γId×d + σ1d×d) : ( D̃2 hUα − D̂2 hUα ) for D2,θ h Uα defined analogously to (3.4) for each value of α and θ. We then solve the optimization problem infθ∈Λ F̂ θ γ,σ[Uα] = 0. When choosing γ = 0, the solver does not successfully find a solution even for coarse meshes. The solver also takes several iterations for γ = 1 sometimes struggling to find a solution for the finer meshes. Thus, we first find the solution corresponding to γ = 1000, and then we sequentially find the solutions correspond- ing to γ = 100, 10, 1, 0 by using the solution for the next largest value of γ as an initial guess. This iterative technique allows fsolve to successfully find a solution for each value of γ and often allows fsolve to significantly decrease the number of it- erations needed to converge. We also can form a better initial guess by first solving the problem corresponding to a small finite number of controls θ before optimizing over Λ. 7.3. Test 3: Monge-Ampère equation. Consider the Monge-Ampère problem F [u] = −det(D2u) + f = −uxxuyy + |uxy|2 + f = 0 over Ω = (0, 1)2. The problem has a unique convex viscosity solution whenever f ≥ 0. We choose the source term f and boundary function g such that the exact solution is u(x, y) = exp ( x2+y2 2 ) . The matrix Mα is easily found since | ∂F∂Pij | = |Pk`| for (k, `) = (2, 2) if (i, j) = (1, 1), (k, `) = (1, 1) if (i, j) = (2, 2), and (k, `) = (1, 2) if (i, j) = (1, 2) or (i, j) = (2, 1). This problem is degenerate and the uniformity of the g-monotonicity of F̂γ,σ[U ] strongly depends upon U similarly to how the ellipticity of F [u] strongly depends upon the convexity of u. Consequently, we now form initial guesses for fsolve by solving a sequence of problems based on a decreasing sequence of values for σ starting with σ large. As such, we are using a stronger form of the numerical moment to overcome the conditional ellipticity and potential degeneracy of the problem. Rates of convergence for the various tests can be found in Table 3 where we consider the method F̂γ,0 for γ = 1000, 100, 10, 1, 0 and Table 4 where we consider the method F̂−σ,σ for σ = 1000, 100, 10, 1. We observe optimal / near optimal rates 90 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 Table 2. Approximation results for Test 2 using F̂γ,0 with γ = 1000, 100, 10, 1, 0. γ = 1000 γ = 100 γ = 10 h `∞ Error Order `∞ Error Order `∞ Error Order 1.57e-01 1.30e+00 1.17e+00 6.08e-01 9.43e-02 1.27e+00 0.05 9.83e-01 0.35 3.34e-01 1.18 6.15e-02 1.22e+00 0.09 7.51e-01 0.63 1.73e-01 1.54 4.56e-02 1.15e+00 0.20 5.67e-01 0.94 1.01e-01 1.80 3.63e-02 1.08e+00 0.29 4.34e-01 1.16 6.47e-02 1.95 2.89e-02 9.78e-01 0.43 3.19e-01 1.35 4.05e-02 2.05 2.21e-02 8.32e-01 0.61 2.12e-01 1.53 2.30e-02 2.12 γ = 1 γ = 0 h `∞ Error Order `∞ Error Order 1.57e-01 1.32e-01 2.35e-02 9.43e-02 5.08e-02 1.87 1.03e-02 1.61 6.15e-02 2.15e-02 2.01 4.96e-03 1.71 4.56e-02 1.17e-02 2.03 2.92e-03 1.77 3.63e-02 7.43e-03 1.99 1.94e-03 1.79 2.89e-02 4.74e-03 1.97 1.28e-03 1.81 2.21e-02 2.81e-03 1.96 7.92e-04 1.80 of convergence as γ decreases in Table 3 and as σ decreases in Table 4. For γ and σ large the rates appear to be suboptimal but improving towards a rate of 2 as h decreases. We also note that the method F̂−σ,σ appears more accurate than the analogous method F̂γ,0 when γ = σ. Such a relationship is expected based on the sensitivity of the Monge-Ampère problem to the auxiliary boundary condition. Indeed, setting ∆hUα = 0 along the boundary and observing that g must be convex along the tangential direction implies U must be concave along the normal direction which brings a qualitative error into the interior of the domain. We can see this directly in Figure 2 where we plot the approximation for varying h, γ, and σ values. For γ > 0 on a coarse mesh we can see that the convexity of the approximation is incorrect near the boundary and that for γ large this forces the curvature to be incorrect throughout the interior. Consequently, for γ large, the method does not appear to enforce the convexity of the underlying viscosity solution. Instead, based on the solver’s performance and based on the tests in [24], the numerical moment for γ large appears to minimize the number of times the function changes convexity over the domain by penalizing discontinuities in the second derivative and steering the approximation towards the correct viscosity solution as h → 0. We also note that, in contrast, choosing σ large and minimizing the effect of the auxiliary boundary condition by setting γ = −σ appears to completely eliminate the convexity issue. Another way to decrease the convexity issue when γ > 0 is to set only δ2 xi,hi Uα = 0 along the normal direction instead of setting ∆hUα = 0 or to choose a more appropriate positive value for the auxiliary boundary condition consistent with the convex nature of the viscosity solution. EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 91 Table 3. Approximation results for Test 3 using F̂γ,0 with γ = 1000, 100, 10, 1, 0. γ = 1000 γ = 100 γ = 10 h `∞ Error Order `∞ Error Order `∞ Error Order 2.83e-01 5.59e-01 5.49e-01 4.02e-01 1.29e-01 5.86e-01 -0.06 5.14e-01 0.08 9.65e-02 1.81 6.15e-02 5.64e-01 0.05 2.22e-01 1.14 2.09e-02 2.08 3.01e-02 4.22e-01 0.40 5.07e-02 2.06 5.23e-03 1.94 1.99e-02 2.32e-01 1.46 2.18e-02 2.04 2.39e-03 1.90 1.49e-02 1.31e-01 1.97 1.22e-02 2.00 1.38e-03 1.89 1.19e-02 8.13e-02 2.11 7.95e-03 1.90 8.98e-04 1.91 γ = 1 γ = 0 h `∞ Error Order `∞ Error Order 2.83e-01 4.19e-02 1.57e-02 1.29e-01 9.30e-03 1.91 3.41e-03 1.94 6.15e-02 2.31e-03 1.89 7.87e-04 1.99 3.01e-02 5.87e-04 1.91 1.88e-04 2.01 1.99e-02 2.64e-04 1.94 8.21e-05 2.01 1.49e-02 1.50e-04 1.93 4.58e-05 2.00 1.19e-02 9.70e-05 1.94 2.91e-05 2.00 Table 4. Approximation results for Test 3 using F̂γ,σ for σ = 1000, 100, 10, 1 and γ = −σ. σ = 1000 σ = 100 σ = 10 σ = 1 h `∞ Error Order `∞ Error Order `∞ Error Order `∞ Error Order 2.83e-01 3.65e-02 3.36e-02 1.95e-02 1.47e-02 1.29e-01 3.56e-02 0.03 2.49e-02 0.38 6.22e-03 1.45 3.05e-03 1.99 6.15e-02 3.07e-02 0.20 1.19e-02 1.00 1.63e-03 1.81 7.00e-04 2.00 3.01e-02 1.97e-02 0.62 3.83e-03 1.59 4.07e-04 1.95 1.67e-04 2.00 1.99e-02 1.24e-02 1.13 1.79e-03 1.84 1.79e-04 1.98 7.31e-05 2.00 1.49e-02 8.14e-03 1.44 1.02e-03 1.92 1.00e-04 1.99 4.08e-05 2.00 1.19e-02 5.66e-03 1.62 6.59e-04 1.95 6.41e-05 1.99 2.60e-05 2.00 7.4. Test 4: Prescribed Gauss curvature equations. This example is adapted from [19]. Let K > 0. The equation of prescribed Gauss curvature corresponds to the choice F [u] = − det(D2u) (1 + |∇u|2) d+2 2 +Kf = 0. The problem is based on the Monge-Ampère operator, and it is known that the problem when f = 1 has a unique convex viscosity solution for each K ∈ [0,K∗) for some positive constant K∗. We let K = 0.1 and Ω = (0, 1)2, and we choose f 92 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 Figure 2. Plots of various approximations for Test 3 using F̂γ,σ. The first row corresponds to h = 1.29e-01 with σ = 0 and γ = 1000, 100, 10 from left to right. The second row corresponds to h = 3.01e-02 with σ = 0 and γ = 1000, 100, 10 from left to right. The third row corresponds to h = 3.01e-02 with σ = 0 and γ = 1, 0 from left to right. The first plot on the fourth row corresponds to h = 3.01e-02 with σ = 1000 and γ = −1000. The second plot on the fourth row corresponds to h = 1.19e-02 with γ = σ = 0. and g such that the exact solution is given by u(x, y) = e x2+y2 2 . The matrix Mα is easily found since | ∂F∂Pij | = |Pk`|/(1 + |v|2)2 for (k, `) = (2, 2) if (i, j) = (1, 1), (k, `) = (1, 1) if (i, j) = (2, 2), and (k, `) = (1, 2) if (i, j) = (1, 2) or (i, j) = (2, 1). Rates of convergence for the various tests can be found in Table 5 where we consider the method F̂γ,0 for γ = 1000, 100, 10, 1, 0 and Table 6 where we consider the method F̂−σ,σ for σ = 1000, 100, 10, 1. Overall we observe similar behavior as Test 3 but less accuracy and lower rates when γ > 0 or σ > 0. The method F̂0,0 does exhibit an optimal convergence rate of 2 and is the most accurate of all of the methods tested. Since this problem involves the gradient operator, any monotone method that directly approximates the gradient would in general be limited to only EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 93 first order accuracy until the mesh is fine enough to use the local uniform ellipticity to enforce the monotonicity with respect to ∇±hUα. Table 5. Approximation results for Test 4 using F̂γ,0 with γ = 1000, 100, 10, 1, 0. γ = 1000 γ = 100 γ = 10 h `∞ Error Order `∞ Error Order `∞ Error Order 2.83e-01 1.30e+00 5.57e-01 5.33e-01 1.29e-01 1.27e+00 0.05 5.80e-01 -0.05 4.05e-01 0.35 6.15e-02 1.22e+00 0.09 5.24e-01 0.14 1.67e-01 1.20 3.01e-02 1.15e+00 0.20 2.79e-01 0.88 8.64e-02 0.93 1.99e-02 1.08e+00 0.29 1.72e-01 1.17 5.76e-02 0.98 1.49e-02 9.78e-01 0.43 1.31e-01 0.96 4.24e-02 1.06 1.19e-02 8.32e-01 0.61 1.06e-01 0.91 3.29e-02 1.12 γ = 1 γ = 0 h `∞ Error Order `∞ Error Order 2.83e-01 2.56e-01 2.19e-02 1.29e-01 1.08e-01 1.10 3.89e-03 2.19 6.15e-02 5.52e-02 0.91 8.20e-04 2.11 3.01e-02 2.51e-02 1.10 1.90e-04 2.05 1.99e-02 1.49e-02 1.26 8.25e-05 2.02 1.49e-02 1.01e-02 1.35 4.59e-05 2.01 1.19e-02 7.30e-03 1.43 2.92e-05 2.01 Table 6. Approximation results for Test 4 using F̂γ,σ for σ = 1000, 100, 10, 1 and γ = −σ. σ = 1000 σ = 100 σ = 10 σ = 1 h `∞ Error Order `∞ Error Order `∞ Error Order `∞ Error Order 2.83e-01 3.68e-02 3.63e-02 3.20e-02 1.81e-02 1.29e-01 3.70e-02 -0.01 3.40e-02 0.08 2.11e-02 0.53 6.83e-03 1.24 6.15e-02 3.59e-02 0.04 2.70e-02 0.31 1.02e-02 0.98 2.09e-03 1.60 3.01e-02 3.19e-02 0.17 1.61e-02 0.72 3.84e-03 1.38 5.75e-04 1.81 1.99e-02 2.73e-02 0.38 1.04e-02 1.06 1.96e-03 1.62 2.61e-04 1.92 1.49e-02 2.31e-02 0.57 7.22e-03 1.26 1.18e-03 1.74 1.47e-04 1.96 1.19e-02 1.96e-02 0.73 5.28e-03 1.39 7.88e-04 1.81 9.45e-05 1.97 Acknowledgements. X. Feng was partially supported by the NSF grants DMS- 1620168 and DMS-2012414. T. Lewis and K. Ward were partially supported by the NSF grant DMS-2111059. 94 X. FENG, T. LEWIS, K. WARD EJDE-2022/CONF/26 References [1] G. Barles, P. E. Souganidis; Convergence of approximation schemes for fully nonlinear second order equations, Asymptotic Anal., 4 (1991), 271–283. [2] F. Bonnans, H. Zidani; Consistency of generalized finite difference schemes for the stochastic HJB equation, SIAM J. Numer. Anal., 41 (2003), 1008–1021. [3] L. A. Caffarelli, X. Cabré; Fully nonlinear elliptic equations, Vol. 43 of American Mathemat- ical Society Colloquium Publications, AMS, Providence, RI, 1995. [4] M. G. Crandall, H. Ishii, P.-L. Lions; User’s guide to viscosity solutions of second order partial differential equations, Bull. Amer. Math. Soc., 27 (1992), 1–67. [5] M. G. Crandall, P.-L. Lions; Viscosity solutions of Hamilton-Jacobi equations, Trans. Amer. Math. Soc., 277 (1983), 1–42. [6] M. G. Crandall, L. C. Evans, P.-L. Lions; Some properties of viscosity solutions of Hamilton- Jacobi equations, Trans. Amer. Math. Soc., 282 (1984), 487–502. [7] K. Debrabant, E. Jakobsen; Semi-Lagrangian schemes for linear and fully non-linear diffu- sion equations, Math. Comp., 82 (2013), 1433–1462. [8] W. Feng, T. Lewis, ,S. Wise; Discontinuous Galerkin derivative operators with applications to second order elliptic problems and stability, Mathematical Meth. in App. Sciences, 38(18) (2015), 5160–5182. [9] X. Feng, R. Glowinski, M. Neilan; Recent developments in numerical methods for fully non- linear second order partial differential equations, SIAM Rev., 55 (2013), 205–267. [10] X. Feng, M. Jensen; Convergent semi-Lagrangian methods for the Monge-Ampère equation on unstructured grids, SIAM J. Numer. Anal., 55 (2017), 691–712. [11] X. Feng, C. Kao, T. Lewis; Convergent FD methods for one-dimensional fully nonlinear second order partial differential equations, J. Comput. Appl. Math., 254 (2013), 81–98. [12] X. Feng, T. Lewis; Nonstandard local discontinuous Galerkin methods for fully nonlinear second order elliptic and parabolic equations in high dimensions, J. Scient. Comput., 77 (2018), 1534–1565. [13] X. Feng, T. Lewis; A Narrow-stencil finite difference method for approximating viscosity so- lutions of fully nonlinear elliptic partial differential equations with applications to Hamilton- Jacobi-Bellman equations, http://arxiv.org/abs/1907.10204 (2019). [14] X. Feng, T. Lewis; A Narrow-stencil finite difference method for approximating viscosity solutions of Hamilton-Jacobi-Bellman equations, SIAM J. Numer. Anal., 59 (2021), 886– 924. [15] X. Feng, T. Lewis, A. Rapp; Dual-Wind Discontinuous Galerkin Methods for Stationary Hamilton-Jacobi Equations and Regularized Hamilton-Jacobi Equations, Commun. Appl. Math. Comput., 4 (2022), 563–596. [16] X. Feng, T. Lewis, M. Neilan; Discontinuous Galerkin finite element differential calculus and applications to numerical solutions of linear and nonlinear partial differential equations, J. Comput. Appl. Math., 299 (2016), 68–91. [17] X. Feng, M. Neilan; Vanishing moment method and moment solutions for second order fully nonlinear partial differential equations, J. Scient. Comput., 38 (2008), 74–98. [18] X. Feng, M. Neilan; Mixed finite element methods for the fully nonlinear Monge-Ampére equation based on the vanishing moment method, SIAM J. Numer. Anal., 47 (2009), 1226– 1250. [19] X. Feng, M. Neilan; Finite element approximations of general fully nonlinear second order elliptic partial differential equations based on the vanishing moment method, Computers & Mathematics with Applications, 68 (2014), 2182–2204. [20] W. H. Fleming, H. M. Soner; Controlled Markov Process and Viscosity Solutions, Springer, New York, 2006. [21] D. Gilbarg, N. S. Trudinger; Elliptic Partial Differential Equations of Second Order, Classics in Mathematics, Springer-Verlag, Berlin, 2001. [22] M. Jensen, I. Smears; On the convergence of finite element methods for Hamilton-Jacobi- Bellman equations, SIAM J. Numer. Anal., 51 (2013), 37–162. [23] H. Kushner, P. G. Dupuis; Numerical methods for stochastic optimal control problems in continuous time, Volume 24 of Applications of Mathematics, Springer, New York, 1992. EJDE-2018/CONF/26 A NARROW-STENCIL NUMERICAL FRAMEWORK 95 [24] T. Lewis, Finite Difference and Discontinuous Galerkin Finite Element Methods for Fully Nonlinear Second Order Partial Differential Equations, Ph.D. Thesis, University of Ten- nessee, 2013. [25] T. Lewis, M. Neilan; Convergence analysis of a symmetric dual-wind discontinuous Galerkin method, J. Sci. Comput., 59 (2014), 602-625. [26] T. S. Motzkin, W. Wasow; On the approximation of linear elliptic elliptic differential equa- tions by difference equations with positive coefficients, J. Math. Phys., 31 (1953), 253–259. [27] M. Neilan, A.J. Salgado, W. Zhang; Numerical analysis of strongly nonlinear PDEs, Acta Numerica, 26 (2020), 137–303. [28] R. H. Nochetto, D. Ntogakas, W. Zhang; Two-scale method for the Monge-Ampére equation: convergence rates, IMA J. Numer. Anal., 39(3) (2019), 1085–1109. [29] A. Oberman, Wide stencil finite difference schemes for the elliptic Monge-Ampère equation and functions of the eigenvalues of the Hessian, Discrt. Cont. Dynam. Syst. series B, 10 (2008), 221–238. [30] M. Safonov, Nonuniqueness for second-order elliptic equations with measurable coefficients, SIAM J. Numeri. Anal., 30(4) (1999), 379–395. [31] A. J. Salgado, W. Zhang; Finite element approximation of the Isaacs equation, ESAIM: Math. Model. Numer. Anal., 53(2) (2019), 351–374. [32] I. Smears, E. Süli; Discontinuous Galerkin finite element approximation of Hamilton-Jacobi- Bellman equations with Cordes coefficients, SIAM J. Numer. Anal., 52 (2014), 993–1016. Xiaobing Feng Department of Mathematics, The University of Tennessee, Knoxville, TN 37996, USA Email address: xfeng@utk.edu Thomas Lewis Department of Mathematics and Statistics, The University of North Carolina at Greens- boro, Greensboro, NC 27412, USA Email address: tllewis3@uncg.edu Kellie Ward Department of Mathematics and Statistics, The University of North Carolina at Greens- boro, Greensboro, NC 27412, USA Email address: kmward7@uncg.edu 1. Introduction 2. Preliminaries and numerical derivatives 2.1. Notation and definitions 2.2. Finite difference derivative operators 2.3. Discontinuous Galerkin finite element derivative operators 3. A narrow-stencil and g-monotone numerical framework 3.1. A narrow-stencil FD framework 3.2. Examples of g-monotone FD methods 3.3. A narrow-stencil DG framework 4. Convergence analysis 5. Admissibility analysis 5.1. Admissibility of a simple method for linear, constant coefficient PDEs 5.2. Admissibility for fully nonlinear PDEs 6. Stability analysis 7. Numerical experiments 7.1. Test 1: Linear F with non-aligned grids 7.2. Test 2: Hamilton-Jacobi-Bellman equations 7.3. Test 3: Monge-Ampère equation 7.4. Test 4: Prescribed Gauss curvature equations Acknowledgements References