1. Introduction. Fully nonlinear elliptic partial differential equations (PDE) arise naturally from differential geometry, optimal mass transportation, stochastic optimal control and other fields of science and engineering. In spite of their wide range of applications, the numerical methods for this type of PDEs is still under development and this is particularly the case if one wants to apply the finite element method (FEM). A major difficulty in their numerical approximation is that, for fully nonlinear PDEs, the correct notion of solution is the so-called viscosity solution, which is based on the maximum principle instead of a variational one. This is reflected in the fact that, in contrast to an extensive literature for linear and quasilinear elliptic PDEs in divergence form, the numerical approximation via finite elements of fully nonlinear PDEs reduces to a few papers; we refer the reader to [14,31] for an overview. The situation is somewhat more satisfactory for finite difference approximations, where convergence to the viscosity solution without rates was studied in the early works [2,24]. Rates of convergence for convex/concave fully nonlinear elliptic equations have been established, e.g. by Krylov [21,22]; Barles and Jakobsen [1] and Debrabant and Jakobsen [12]. However, rates of convergence for the nonconvex/nonconcave case remained an open problem as the techniques used in the convex case could not be generalized to this scenario. For many years the existing results were rather specialized. For instance, [17] considered a one dimensional problem and obtained rates of convergence but its arguments do not extend to more dimensions. A particular nonconvex equation: an obstacle problem for a Hamilton Jacobi Bellman equation was studied in [18], where this particular structure is exploited to obtain convergence rates. This changed when, in 2008, Caffarelli and Souganidis [9] established a rate of convergence for finite difference approximations of elliptic equations of the form

with Dirichlet boundary conditions. This result was extended by Turanova [37] to the case where F is also dependent on x and by Krylov [23] to the Isaacs equation, where he allows the operator F to depend on x, ∇u and u. We finally comment that, to our knowledge, no rates of convergence are available for semi-Lagrangian schemes. Reference [12] provides convergence, but no rates are available for nonconvex operators of order two. It is remarkable, however, to note that explicit rates of convergence for nonlocal operators of order µ ∈ (1, 2) are known. Namely, [3] shows that the rate of convergence is O(h 2-µ µ ). Note that this rate degenerates as the order of the operators approaches two and that the authors of this work state that this is optimal and expected.

In all the works mentioned above, convergence hinges on operator consistency and monotonicity, which are nontrivial properties to be satisfied, especially in the FEM. To overcome this, the recent work [32] introduced a two scale FEM for linear elliptic PDEs in nondivergence form and, on the basis of monotonicity and weak operator consistency, the authors were able to prove convergence of the method with rates. The purpose of this work is to extend these ideas to a particular kind of uniformly elliptic Isaacs equations on convex domains, namely,

where

) is an open, bounded and convex domain, f ∈ C 0,1 ( Ω), the sets A and B are arbitrary finite sets and the matrices A α,β ∈ C 0,1 ( Ω, GL d (R)) are symmetric and uniformly elliptic, in the sense that there are constants λ, Λ ∈ R, with 0 < λ ≤ Λ, such that

We will introduce, following [32], a two-scale FEM, show its convergence to the viscosity solution of (1.1) and provide rates of convergence. The method is based on the approximation of (1.1) proposed by Caffarelli and Silvestre in [8]: we formally rewrite

and we approximate the operator above by the integro-differential operator

where

Hereafter, ϕ is a radially symmetric function with support in the unit ball B 1 of R d , where d ≥ 1 is the dimension, that verifies ˆRd |y| 2 ϕ(y) dy = d and du(x, y) := u(x + y) -2u(x) + u(x -y) (1.4) is the centered second difference operator. The operator

in the sense that if u is a quadratic polynomial, then

see Lemma 2.1.

We discretize (1.1) by using this approximation as follows: Introduce a triangulation T h of Ω and let N h be the set of internal nodes. We introduce a finite element space consisting of piecewise linear functions and denote by {φ z } z∈N h its Lagrange nodal basis. Now, multiply (1.3) by φ z and integrate over Ω. Integrate the first term by parts and apply mass lumping to the second to finally obtain that the finite element approximation, u ε h , satisfies

We show that the FEM is monotone provided that meshes are weakly acute. To show existence and uniqueness, we employ a discrete version of Perron's method, which seems to not have been considered before, especially in the finite element literature. Exploiting monotonicity, and using the notion of weak consistency introduced by Jensen and Smears in [19] we show convergence to the viscosity solution following [2,24]. We also derive rates of convergence for the method. The main difficulty to derive a rate of convergence is to establish a suitable notion of stability for the FEM applied to this fully nonlinear PDE. To address this issue, we resort to the discrete Alexandrov Bakelman Pucci (ABP) estimate of [32], which reads sup

, and, for each node z ∈ N h , we set ω z = ´Ω φ z (x) dx. Note that the nodal contact set is just a finite collection of nodes. With (1.6) we obtain control of the negative part of u ε h . If we consider the concave envelope and corresponding (upper) contact set, we can estimate the positive part. A combination of these bounds yields stability, in the L ∞ -norm, of u ε h in terms of the L d -norm of the right hand side f. Suitable notions of consistency, monotonicity and stability yield rates of convergence. Moreover, by combining them with the regularity estimates for the Isaacs equation of [8], we show that, for some σ ∈ (0, 1],

which, at best, gives us a O(h 1/3 | log h| 1/3 ) error estimate.

We also discuss how to practically realize the method in question. We study, following [6], a variant of Howard's algorithm to solve the ensuing discrete (nonlinear) systems.

Before proceeding any further, we must also comment that our ideas and methods also apply to a version of the Hamilton Jacobi Bellman equation. Namely, setting #A = 1 or #B = 1 we arrive at a convex or concave operator. In this case much more could be said for this problem, and it is possible that our analysis yields suboptimal error estimates. In particular, the rate of convergence σ, which comes form (2.6) below, might be known. Since our main interest here is to deal with the Isaacs equation, we will not explore this refinement any further.

We must also mention that the linear case, i.e., #A = #B = 1, was already treated in [32]. Our analysis, however, departs from theirs in several respects: • To apply the regularity and approximation results of [8], as opposed to [32], the second difference operator d must be of the form (1.4), and should not be modified near the boundary. This forces us to consider extensions of functions outside the domain Ω, and our analysis must take this into account. • Since our problem is nonlinear existence is not guaranteed by uniqueness, as it is the case in [32]. Thus, although uniqueness is obtained by similar (monotonicity and maximum principle based) arguments, we must develop a discrete Perron's method to assert existence of solutions. • Due to the fact that solutions to the Isaacs equation (1.1) exhibit very low regularity, the convergence arguments that we employ must take this into account.

In particular this means that the relation between the two scales, ε and h, of our method must be modified, and that suitable barrier functions must be constructed. The rest of this paper is organized as follows. In Section 2, we recall the approximation of elliptic problems by integro-differential equations and the main regularity results that follow from it. Section 3 presents our discretization and proves convergence to the viscosity solution. We recall the discrete ABP estimate for finite element methods in Section 4. The discrete ABP estimate allows us to show stability and rates of convergence for our discretization. We discuss some implementation details in Section 5, where we present a convergent iterative scheme.

We will follow standard notation concerning differential operators and function spaces. The relation A B means that there is a nonessential constant C such that A ≤ CB. The value of this constant might change at each occurrence. By A B we mean B A.

Let us review the approximation, proposed by Caffarelli and Silvestre in [8], of fully nonlinear elliptic PDEs by an integro-differential equation and the convergent finite element scheme of [32] for the approximation of elliptic problems in nondivergence form based on this idea. In this work, we exploit this approximation to discretize the fully nonlinear problem (1.1).

Let us begin by fixing some notation. Recall that, owing to (1.2), for all x ∈ Ω the matrices A α,β (x) are uniformly positive definite. Thus, we can define

Let ε > 0 and ϕ be a radially symmetric function with compact support in the unit ball B 1 and such that ´|x| 2 ϕ(x) dx = d. Let Q = 2/λ and notice that, again because of (1.2), for all α ∈ A and β ∈ B we have that if

For this reason, we define

a Qε-neighborhood of ∂Ω and

Given a function w ∈ C(Ω ε ) we define, for each α and β the function

where we denoted by dw(x, y) the second order difference of the function w at the point x in the direction y which is given in (1.4). Notice that, for x ∈ Ω and y ∈ supp ϕ α,β x,ε , we have that x ± y ∈ Ω ε so that (2.3) is well defined. It is easy to check that if w is a quadratic polynomial, then we have dw(x, y) = D 2 w(x) : (y ⊗ y) for all x ∈ Ω.

The integral operator I α,β ε [w] is a consistent approximation of the differential operator A α,βλ 2 I : D 2 w in the following sense [32,8]. Lemma 2.1 (approximation properties of I α,β ε [•]). Let I α,β ε [•] be the integral operator defined by (2.3) and assume that ε > 0 and x ∈ Ω. Then, 1. If p ∈ P 2 , i.e., it is a quadratic polynomial, then

2. Let ε 0 > 0 be fixed. If w ∈ C 2 (Ω ε0 ) then, as ε → 0, we have that

Let us, from now on, assume that the matrices M α,β have a uniform modulus of continuity . In other words, there is a nondecreasing function such that lim t↓0 (t) = (0) = 0 and sup

(2.4)

Under this assumption, we can show that, for every ε > 0, the integral operator maps continuously

, then, for all α ∈ A and β ∈ B and all x, z ∈ Ω

where the hidden constant is independent of α, β, x, z, ε and w.

Proof. By definition of integral operator (2.3) and a change of variable y := εM α,β (x)ỹ, we obtain

By the radial symmetry of ϕ

Thus,

It remains then to estimate the first term. To do so we, again, use that w ∈ C 0,1 (Ω ε ) to obtain

an application of the triangle inequality, together with (2.4) imply

where the hidden constant is uniform in α and β. Gathering the obtained bounds for I and II allows us to conclude.

In [8] Caffarelli and Silvestre proposed an approximation of the Isaacs equation (1.1) by the following integro-differential problem

(2.5)

Notice that, for x ∈ Ω, the definition of F ε [u ε ](x) requires values of u ε in Ω ε . Thus in (2.5) the boundary condition has been replaced by a so-called volume constraint over (∂Ω) ε . To define this volume constraint, the function that defines the boundary condition in (1.1) (in this case g ≡ 0) has been extended in a way that its extension belongs to C 0,1 ((∂Ω) ε ).

The approximation property of the integral operator I α,β ε [•] given in Lemma 2.1, allows us to relate u ε and u, which solve (2.5) and (1.1), respectively, as follows [8].

Proposition 2.3 (approximation of the Isaacs equation). Let Ω ⊂ R d be convex and bounded, and f ∈ C 0,s ( Ω) for some s ∈ [0, 1). There exists a unique function

where U is any open set compactly contained in Ω and the hidden constant depends on the distance between U and ∂Ω. In addition, there is a σ > 0 such that

where u is the (unique) viscosity solution to (1.1).

An important consequence of Proposition 2.3 is that, in the case when the coefficient matrices A α,β are independent of x, we have the following regularity result for the function u, see Theorem 4.8 of [8] and the comments at the beginning of Section 4, for the C 1,s estimate and Theorem 5.2 and the comment at the beginning of Section 5 for the C 0,1 estimate up to the boundary.

Corollary 2.4 (regularity of u). Assume that, for all α ∈ A and

The importance of this result lies in the fact that F is not convex nor concave and, for that reason, the maximal regularity we can assert for u is C 1,s (Ω), for some s > 0. We refer the reader, for instance, to the works by Nadirashvili and Vlȃdut [29,30,28], who have constructed viscosity solutions to nonconvex fully nonlinear elliptic equations whose Hessian is not bounded. This is in sharp contrast with, for instance, the uniformly elliptic Hamilton Jacobi Bellman equation where, under suitable assumptions on Ω and the "coefficients" of the equation, it can be shown that the solution belongs to C 2,s (Ω) for some s > 0. Let us, finally, comment that the rate of convergence given in (2.6) cannot be improved; see the last paragraph of [8].

In this paper, based on the idea of an integro-differential approximation, we propose a finite element method for Isaacs equation (1.1).

The approximation properties of the integral operator in Lemma 2.1 and 2.2 are essential to ensure the consistency of the method. Let us briefly discuss here the use of numerical quadrature rules that preserve these properties. For each fixed α ∈ A and β ∈ B, we aim to design a quadrature and choice of ϕ to approximate I α,β ε [w]. The change of variable y = εM α,β (x)ζ reveals that we must approximate

This shows that, to preserve the consistency and convergence properties of the integral operator, we must require the quadrature to satisfy the following three conditions:

• For all j = 1, . . . , J, the quadrature points q j must satisfy q j ∈ supp(ϕ).

• The weights must be positive: ω j > 0 for all j = 1, . . . , J.

• The quadrature must be exact for all quadratic polynomials, i.e., for a given kernel function ϕ with ´Rd |y| 2 ϕ(y) dy = d we must have that, if p ∈ P 2 , then

The construction of quadrature formulas of this form which are of maximal degree of precision can be found in classical references like [27], see also [4,5,33]. For instance, for dimension d = 2, given ϕ(y) = 4/π in the unit ball and ϕ(y) = 0 outside the unit ball, reference [36] shows that the following collection of weights and nodes has the requisite properties:

where the weights q j = (ρ j , θ j ) are given in polar coordinates. Indeed, since ϕ(q j ) = 4/π, this formula is exact for quadratic polynomials and it preserves the consistency in Lemma 2.1 and Lemma 2.2.

Here we describe the scheme we use to approximate (1.1) and show that it converges to the viscosity solution. To keep technicalities to a minimum, in what follows we assume that Ω is a convex polytope.

Let {T h } h>0 be a quasi-uniform, in the classical sense for finite elements [13,11,7], family of triangulations of size h of the domain Ω. We denote by N h and N ∂ h the collection of interior and boundary nodes of T h , respectively. We define the finite element space

where we denoted by P 1 the space of polynomials of degree one.

We also need to introduce a triangulation of Ω ε . Since this may not be a polytope we proceed in a standard way, namely, we introduce a triangulation T ε h of the polytopal domain Ω T ε ⊂ Ω ε in such a way that:

and, as Ω and Ω ε are convex, we have that Ω ⊂ Ω T ε ⊂ Ω ε and we can extend continuously functions in these spaces from Ω T ε to Ω ε by a constant in the direction normal to ∂Ω T ε . We denote by {φ z } z∈N h ∪N ε h the Lagrange nodal basis of V ε h . Notice that, by restriction, we can consider V ε h to be continuously embedded in V h . Below we will require T h to be weakly acute, i.e., for all z 1 , z 2 ∈ N h with z 1 = z 2 we must have

It is known, see [10,Section III.20] and [31,Lemma 3.42], that this is a sufficient condition for the discrete Laplacian, defined below, to be a monotone operator. Let z ∈ N h , we set

For w h ∈ V h and z ∈ N h , we define the discrete Laplacian of w h at z by

Finally, for z ∈ N h , we set

With this notation at hand, we now define our scheme. We seek a function u ε h ∈ V 0 h such that, for every z ∈ N h , satisfies

3.2. Galerkin projection. In the analysis that follows we will make repeated use of the Galerkin projection and its properties. The Galerkin projection is the map

that is, the discrete Laplacian of the Galerkin projection at a node is a weighted average of the Laplacian of the original function around this node.

Combining the well known near optimal approximation properties [34,25], in the L ∞ -norm, of the Galerkin projection with a standard interpolation estimate over

Owing to (3.5), for every z ∈ N h we obtain

By definition (2.3) of the integral operator I α,β ε [•], we deduce that for all nodes z ∈ N h we have

uniformly in A and B.

Based on the notions of monotonicity and weak consistency advanced in [2,19,32], we show existence and uniqueness of solutions to (3.3). Moreover, the family {u ε h } h>0,ε>0 of solutions to (3.3) converges to u, the unique viscosity solution of (1.1), provided the mesh size h and ε satisfy a suitable relation.

We begin with two elementary properties, namely monotonicity and continuity, of the inf-sup operator.

Lemma 3.1 (monotonicity and continuity of inf-sup). Let {X α,β : α ∈ A, β ∈ B} ⊂ R and {Y α,β : α ∈ A, β ∈ B} ⊂ R be two families parametrized by two index sets A and B. If for every fixed α ∈ A and β ∈ B we have that

Moreover, if there is a constant C such that for all α and β we have

Proof. If X α,β ≤ Y α,β for all fixed α and β, taking supremum on Y α,β with respect to β first, we have

Taking infimum on sup β∈B X α,β with respect to α, we get

This proves the first inequality (3.7).

To prove the second inequality, we only need to note that if X α,β ≤ Y α,β + C for all α and β, then by (3.7)

This proves the second inequality (3.8).

We now establish the monotonicity and weak consistency of the proposed scheme. Lemma 3.2 (monotonicity of the numerical scheme). Assume that for every discretization parameter h > 0 the mesh T h is weakly acute. Then the family of operators

h and equality holds at some node z ∈ N h , then

Proof. Let v h , w h ∈ V ε h be as indicated. Since T h is weakly acute, the operator ∆ h is monotone, i.e.,

It remains then to deal with the inf-sup over the approximating integral operators. Since, over T ε h , we have that v h ≤ w h and equality holds at z ∈ N h , we can conclude that, whenever y ∈ B Qε , the second order differences verify dv h (z, y) ≤ dw h (z, y). Therefore, for fixed α ∈ A and β ∈ B, we get

Using the monotonicity of the inf-sup operator (3.7), we obtain

Combining both inequalities, we then deduce that

This completes the proof.

As a consequence, we obtain a discrete maximum principle for our numerical scheme.

Corollary 3.3 (discrete maximum principle). Assume that T h is weakly acute.

and, for all z ∈ N h , we have

We argue by contradiction by assuming that v h -w h attains a positive value. We then necessarily have that there is z 0 ∈ N h for which v h (z 0 ) -w h (z 0 ) > 0.

We now set p(x) = |x| 2 -C, where the constant C is so large that, for all x ∈ Ω ε , we have p(x) ≤ 0. Let p h ∈ V ε h be the Lagrange interpolant of p and, for δ > 0, we define e h = v h -w h + δp h ∈ V ε h . By construction we have that e h ≤ 0 in N ε h . In addition, for δ > 0 sufficiently small, e h (z 0 ) > 0. Consequently, there is z 1 ∈ N h where e h attains a, necessarily positive, maximum. At this node we must have, for all α ∈ A and β ∈ B that

The weak acuteness of the mesh T h allows us to conclude also that

Adding these two inequalities yields that, for all α ∈ A and β ∈ B we have

Notice that since p is convex and p h is its piecewise linear interpolant we have that, for all y ∈ B Qε , dp h (z 1 , y) ≥ dp(z 1 , y). This, in conjunction with Lemma 2.1, implies that

Observe now that the linear function ψ

where, to get rid of the constant term -C, we used the partition of unity property of the Lagrange basis functions {φ z } z∈N h ∪N ∂ h . Finally, we recall that the weak acuteness of T h implies that, for z = z 1 , ´Ω ∇φ z (x)•∇φ z1 (x) dx ≤ 0 to obtain ∆ h p h (z 1 ) ≥ 0.

(3.11)

Inserting (3.10) and (3.11) in (3.9) yields

Applying Lemma 3.1 to this inequality we deduce that

which is a contradiction. The analysis of our scheme is based on a careful estimation of the consistency error, which we define as the difference between the result of applying the operator to the Galerkin projection of the solution to (2.5) and the approximate right hand side, i.e., F ε h [G h u ε ](z) -f z . The critical step in this analysis is to understand what happens with the integral operators. With this in mind, and for future use, for a function w ∈ C(Ω ε ) and z ∈ N h we denote

The following result bounds R h,ε [w](z) uniformly in z ∈ N h when w ∈ C 0,1 (Ω ε ). Lemma 3.4 (Galerkin projection vs. integral operator). Let w ∈ C 0,1 (Ω ε ). If the coefficient matrices are such that (2.4) holds then, for every z ∈ N h , we have

where the hidden constant is independent of h, ε and w.

Proof. Let z ∈ N h and consider, for x ∈ supp φ z ,

Estimate (3.6) immediately yields

On the other hand, Lemma 2.2 implies

where we used that |x -z| h. We now invoke Lemma 3.1 with

which immediately implies (3.13).

The monotonicity property given in Lemma 3.2 ensures that, for every ε > 0 and h > 0, the scheme (3.3) has a unique solution. The main idea behind the the existence and uniqueness of the numerical solution is a discrete variant of Perron's method [15, §6.1]. Theorem 3.5 (existence and uniqueness). Let the family of meshes {T h } h>0 be weakly acute. For every h > 0 and ε > 0 the finite element scheme (3.3) has a unique solution.

Proof. We begin by proving existence in several steps. First, we define the set of discrete super-solutions

We claim that the set S h is nonempty. Set

where C > 0 is a constant to be chosen later. Let F = min{0, minΩ f} -δ h,ε and let R > 0 be so large that, for some ξ ∈ Ω the ball ξ

, which is a nonnegative quadratic polynomial in ξ + B R , p(x) ≥ 0 on (∂Ω) ε , and D 2 p = λ -1 F I ≤ 0. Owing to the uniform ellipticity condition (1.2), for each α and β, we obtain

We now claim that G h p ∈ S h . To see this, first notice that for every z ∈ N ε h we have G h p(z) ≥ 0. Now, for z ∈ N h , consider

Using the consistency of I α,β ε [•] for quadratics, see Lemma 2.1, we conclude that

Lemma 3.4 and the fact that p C 0,1 (Ωε) f L ∞ (Ω) , show that we can choose C > 0 sufficiently large so that R h,ε [p](z) -δ h,ε ≤ 0 for all z ∈ N h . Therefore,

Thus, S h is nonempty. We now show that the minimum between two super-solutions is a super-solution.

Since, for every z ∈ N ε h , we have v h (z) ≥ 0 and w h (z) ≥ 0, then

The monotonicity result of Lemma 3.2 implies then that

Since for every z ∈ N h we either have that (v ∧w) h (z) = v h (z) or (v ∧w) h (z) = w h (z), we conclude that (v ∧ w) h ∈ S h . Finally, we show that the smallest super-solution is a solution. Let u h ∈ V ε h be defined by

The reasoning given above shows that u h ∈ S h . We claim that u h is a solution, for if that is not the case, then either: Case I: There is a node z 0 ∈ N ε h such that u h (z 0 ) > 0. Define, for δ > 0, v h = u h -δφ z0 and notice that for δ sufficiently small, we have v h (z 0 ) ≥ 0. Consequently, this function verifies

In addition, for all z ∈ N h , we have v h (z ) = u h (z ). The monotonicity of the scheme shown in Lemma 3.2 implies then that

so that v h ∈ S h . However, this contradicts the fact that u h is the smallest supersolution of (3.3). Consequently, u h ∈ V 0 h . Case II: There exists an interior node z 0 ∈ N h such that

We, again, define v h = u h -δφ z0 . Obviously v h (z 0 ) ≤ u h (z 0 ) and, v h (z) = u h (z) for all z = z 0 . Invoking Lemma 3.2 we see that

In addition, since φ z0 attains its absolute strict maximum at z 0 we must have

This, in conjunction with Lemma 3.1 yields that, for δ sufficiently small,

In this case, we also obtain a contradiction. This proves the existence of a solution.

Uniqueness immediately follows from Corollary 3.3. Remark 3.6 (alternative proof of existence). An alternative proof of existence of solutions will be given in Theorem 5.2 via the convergence of Howard's algorithm.

The discrete Alexandrov Bakelman Pucci estimate. The next step towards proving convergence of the discrete solutions to the viscosity solution is to show that solutions of (3.3) are bounded uniformly with respect to h and ε. To achieve this, we will employ the discrete ABP maximum principle of [32,Theorem 5.1], which will also be our main tool to derive a rate of convergence. As Lemma 3.2 shows, weak acuteness of the mesh is sufficient to ensure the monotonicity of the numerical scheme. However, the discrete ABP maximum principle requires a slightly stronger assumption on the mesh, which we state below. We say that the mesh T h is weakly acute with respect to faces if for every face F and z 1 , z 2 ∈ N h with z 1 = z 2 we have ˆωF ∇φ z1 (x) • ∇φ z2 (x) dx ≤ 0, (3.16) where ω F = ∪{K ± ∈ T h : F ⊂ K ± }. This condition must be compared with weak acuteness as expressed in (3.1). In two dimensions (d = 2), condition (3.16) is equivalent to (3.1) and is valid if and only if the sum of two angles opposite to a face is not greater than π. However, condition (3.16) is stronger than weak acuteness for d > 2. We refer the reader to [32,Section 3.3] for further details and references on this condition.

To state the discrete ABP estimate, we also introduce the following notation. Let B R be a ball of radius R which contains the domain Ω. We define the convex envelope of a function where we denote by v - h the negative part of v h in Ω and v - h := 0 in B R \ Ω. We also define the (lower) nodal contact set

Lemma 3.7 (discrete ABP estimate). Let the family of meshes {T h } h>0 be quasiuniform and satisfy (3.16).

, where the hidden constant is independent of h. This inequality gives us an estimate on the negative part of v h , while a bound for its positive part can be derived in the same fashion by considering a concave envelope and the corresponding (upper) contact set. Combining these bounds yields that the L ∞ -norm of v h is controlled by a discrete L d -like norm of {f z } z∈N h . We make this idea rigorous below.

Proposition 3.8 (uniform boundedness). Let the family of meshes {T h } h>0 be quasi-uniform and satisfy (3.16). Then, the family {u ε h ∈ V 0 h } h>0 of solutions of (3.3) is uniformly bounded with respect to h and ε, i.e.,

where the hidden constant is independent of h and ε.

Proof. Notice that, by construction u ε h = 0 on N ∂ h and, consequently (u ε h ) -= 0 only in the interior of Ω. Now, by definition of the convex envelope (3.17) we have that

is convex, we obtain that, for every z ∈ C - h (u ε h ) and y ∈ B Qε we have du ε (z, y) ≥ dΓ(u ε h )(z, y) ≥ 0. Consequently,

In conclusion, for every z ∈ C - h (u ε h ), the scheme (3.3) reduces to

The previous inequality, together with the discrete ABP estimate of Lemma 3.7 then yields sup

where the last inequality follows from the shape regularity of T h . We have obtained a lower bound for u ε h . By considering the positive part, and the corresponding concave envelope, we can obtain an upper bound on u ε h and this yields its uniform boundedness.

Combining the discrete maximum principle of Corollary 3.3 and the consistency estimate encoded in Lemma 3.4 it is now possible, following the guidelines established in [2,24], to prove convergence of our scheme (3.3). We begin by recalling the essence of the convergence results established in [2,24]:

If (1.1) has a strong comparison principle, then uniformly bounded solutions of a consistent, monotone, stable scheme converge to the viscosity solution of (1.1). To be able to assert consistency, we must guarantee that the right hand side of (3.13) tends to zero as we refine the mesh and let ε → 0. This imposes a relation between the parameter ε and the mesh size h of the form

which from now on we assume. At this point the reason why we call our scheme a two-scale one becomes evident. The fine scale is given by the mesh size h and provides a scale for the discretization of functions. On the other hand, the coarse scale ε can be understood as a stencil width that, as is well known [26,20], must be large enough to be able to guarantee monotonicity of a scheme. In this regard, our method bears resemblance to semi-Lagrangian schemes [12] in that functions are approximated with piecewise linears, but their derivatives are approximated at a coarser scale.

To establish convergence, we begin by constructing boundary barrier functions, which will be essential to show that the boundary conditions are attained in a classical sense.

Lemma 3.9 (barrier functions). Assume that Ω is convex and x 0 ∈ ∂Ω. Let h and ε be small enough and satisfy (3.19). With these assumptions, for every positive constant E > 0 there exist finite element functions p ± x0,h ∈ V ε h such that: 1. for all x ∈ (∂Ω) ε we have

2. for all z ∈ N h they satisfy

and, 3. for all

where the hidden constant depends only on the domain Ω.

Proof. Let x 0 ∈ ∂Ω ε be the closest point in ∂Ω ε to x 0 . Without loss of generality we may assume that x 0 is the origin. Now, since Ω is convex, it is not difficult to see that Ω ε is convex as well. Consequently, it has a supporting hyperplane at x 0 which, again without losing generality, we may assume is given by x d = 0. Finally, since Ω ε is bounded, there is L > 0 such that

With this notation at hand we define the functions

and the discrete barriers are

Let us now show that these functions satisfy the claimed properties. 1. To show this we recall that, if z ∈ N ε h ,

Similarly, we have

For ε and h small enough, assumption (3.19)

≤ -λE and we can choose ε and h so that

It remains to establish the growth bound (3.22), but this is an easy consequence of the fact that p ± x0 (0) = 0 and estimate (3.5) for the Galerkin projection. The functions p ± x0,h have thus been constructed. We are now in position to establish convergence of our scheme. Theorem 3.10 (convergence). Assume that f ∈ C 0,1 ( Ω), that the matrices {A α,β : α ∈ A, β ∈ B} are uniformly elliptic (1.2) and satisfy

where M α,β = (A α,βλ 2 I) 1/2 . For h > 0 and ε > 0 let u ε h ∈ V 0 h be the solution to the numerical scheme (3.3), obtained over a mesh T h that is weakly acute with respect to faces, and belongs to a quasi-uniform family. Assume also that (3.19) holds. For

Then the upper semi-continuous function u is a viscosity subsolution of (1.1) with u ≤ 0 on ∂Ω and the lower semi-continuous function u is a viscosity supersolution of (1.1) with u ≥ 0 on ∂Ω. Proof. We begin by remarking that, in light of the uniform boundedness shown in Proposition 3.8, the functions u and u are bounded.

Let us now show that u is a subsolution. Let p be a quadratic polynomial such that u -p has a local maximum at x 0 ∈ Ω. We need to show that

Without loss of generality we can assume that this is a strict maximum.

Let us now show that there are {z h ∈ N h } h>0 such that u ε h -G h p attains its maximum at z h and, as h → 0, we have z h → x 0 . If that is not the case, for any {z h } h>0 , there is a subsequence {z h k } ⊂ {z h } h>0 that converges to some y 0 = x 0 . By definition (3.23) we have

On the other hand, since

Passing to the limit, we have

which contradicts that u -p attains a strict maximum at x 0 .

In what follows we consider this sequence of nodes and, for simplicity, we suppress the subindex h. As in Corollary 3.3, the fact that u ε h -G h p attains its maximum at z yields

which combined with Lemma 3.1 allows us to conclude that

Since f is continuous f z → f(x 0 ) and so it remains to study the behavior of the right hand side in this inequality.

Repeating the computations in the proof of Theorem 3.5 that led to (3.14) allows us to obtain

where R h,ε [p](z) is defined in (3.12). Consequently, if we are able to show that, for every z ∈ N h , R h,ε [p](z) → 0 as h → 0, ε → 0 we would get that

which is what we need to prove. The bound on R h,ε [p](z), obtained in Lemma 3.4, and the choice of ε allow us to conclude. We now show that if x 0 ∈ ∂Ω, then we must have u (x 0 ) ≤ 0. To do so, let ε and h be small enough and set E = 2 λ f L ∞ (Ω) in the barrier function p - x0,h of Lemma 3.9. Property (3.21) yields that, for every z ∈ N h

Moreover, by (3.20), we have p - x0,h ≥ 0 = u ε h in (∂Ω) ε . Consequently, invoking the discrete maximum principle of Corollary 3.3 we obtain

where the upper bound follows from (3.22). Letting ε → 0 and x → x 0 we conclude that u (x 0 ) ≤ 0.

Finally, the fact that u is a super-solution can be shown in a similar fashion. For brevity, we skip the details.

Corollary 3.11 (convergence). Let the family of triangulations {T h } h>0 be quasi-uniform and assume that, for every h > 0, T h is weakly acute with respect to faces (3.16), and that, as ε, h → 0, we have (3.19). In this setting the family {u ε h } h>0,ε>0 of solutions to (3.3) converges pointwise to u, the viscosity solution of the uniformly elliptic Isaacs equation (1.1).

Proof. Using the notation of Theorem 3.10, we observe that u ≥ 0 ≥ u on ∂Ω and that u and u are viscosity sub and super solutions, respectively. By a comparison principle for (1.1); see [31, Theorem 2.52(b)], we then must have

On the other hand, by definition (3.23), u ≤ u . Both inequalities readily imply that u = u = u, the viscosity solution of (1.1). This proves convergence of u ε h to u. Remark 3.12 (nonhomogeneous boundary conditions). All the considerations given above can be extended to the case of nonhomogeneous boundary conditions, i.e., when in (1.1) the boundary condition reads u = g. For that we need to assume that g = tr ∂Ω G for some G ∈ C 2 ( Ω), and a suitable extension of g to (∂Ω) ε must be provided in the second equation of (2.5). As described in [8], this can be achieved by using the so-called sup-convolution:

which given the smoothness of g and Ω satisfies g

Let us now briefly indicate the modifications that need to be made to our arguments: • Existence: To show existence we need to redefine the set of super-solutions. We now require that v h (z) ≥ G h g(z) for z ∈ N ε h . We need to show that this set is nonempty, and for that it suffices to show G h p ∈ S h with

• Uniform boundedness: We consider the function w h = u ε h + M where M is sufficiently large so that w h ≥ 0 on N ∂ h . This is possible because g is bounded. Notice that, for all z ∈ N h ,

and by ABP

• Barrier functions: When proving convergence, we must be able to assert that the boundary conditions are attained in a classical sense. For that it suffices to show that u (x) ≥ g(x) for x ∈ ∂Ω and this is achieved by constructing suitable barrier functions (cf. Lemma 3.9). In this case, the barrier functions will be defined as

We leave the details of this program to the reader.

Let us now study the rate of convergence of u ε h to u, under the smoothness assumptions of Corollary 2.4. We will achieve this by comparing the solution of the discrete scheme (3.3) with the solution of the integrodifferential approximation (2.5). In light of the technique that we are adopting, we must immediately note that since the approximation results of Proposition 2.3 and, as a consequence, the regularity obtained in Corollary 2.4 only apply in the case of constant coefficients, A α,β (x) = A α,β , in what follows we must assume this. In this setting (2.4) is trivially satisfied with ≡ 0. It is possible that the results of [37] can be used to extend Proposition 2.3 to variable coefficients and if that is the case, our results will immediately follow as well.

The main technical tool we will employ to obtain rates of convergence in the L ∞norm will be the discrete ABP estimate of Lemma 3.7. Thus, we must require that the mesh is weakly acute with respect to faces, as defined in (3.16). We finally remark that since the relation (3.19) is required for convergence of our method we will study the rate of convergence under this assumption.

Subtract from this equation the scheme (3.3) to obtain, for all z ∈ N h ,

Notice, finally, that G h u ε -u ε h ≡ 0 over N ε h . 4.2. Rate of convergence. With the error equation (4.1) at hand, we now readily obtain an error estimate. This is the content of the next result.

Theorem 4.1 (rate of convergence). Let u ∈ C 1,s (Ω) ∩ C 0,1 ( Ω) be the viscosity solution of the Isaacs equation (1.1) and u ε h ∈ V ε h be the solution of (3.3). Choose ε so that (3.19) holds with ≡ 0. Then, there is σ > 0 such that

where the hidden constant depends on λ and Λ, but is independent of ε and h. Proof. We write

and we examine each term separately. By Proposition 2.3, there is σ > 0 such that

Estimate (3.5) on the Galerkin projection immediately yields

where, using Proposition 2.3, we bounded u ε C 0,1 ( Ω) by f C 0,1 ( Ω) . Let us denote e h = G h u ε -u ε h ∈ V 0 h . Using the convex envelope Γ(e h ) of the error we can, as in the proof of Proposition 3.8, conclude that for every α ∈ A and β ∈ B we have

so that using the monotonicity property of the inf-sup operator (3.7), we have

In conclusion, at the nodal contact set C -(e h ) equation (4.1) reduces to

An application of the discrete ABP estimate of Lemma 3.7 then yields

sup

This yields a lower bound for e h . An upper bound can be derived in a similar fashion by considering -e h . Hence, we infer

Using the fact that u ε ∈ C 0,1 ( Ω) uniformly in ε we can invoke the bounds on R h,ε [u ε ](z) obtained in Lemma 3.4 which, recalling that ≡ 0, can be combined with the estimate of Proposition 2.3 to yield max

Notice now that (3.19) implies that h| log h|ε -2 → 0 as h → 0 and ε → 0. Therefore, we obtain

Combining the estimates of these three steps yields the result. Notice that in the best case scenario, that is σ = 1, we would obtain a rate of convergence of order O(h 1/3 | log h| 1/3 ).

Remark 4.3 (explicit rate of convergence). The rate of convergence in Theorem 4.1 is given rather implicitly. It seems that this is a recurring feature in the literature; see for instance the main result in [23].

5. Implementation details. Let us discuss how to obtain a solution to the nonlinear problem that scheme (3.3) entails. The main difficulty in devising a convergent algorithm for the solution of (3.3) is the fact that, due to the inf-sup operations, the underlying operator is neither convex nor concave. This is in sharp contrast with, for instance, the interpretation of Howard's algorithm [16] as a semi-smooth Newton method described in [19,35] for the Hamilton-Jacobi-Bellman equation since, as it is shown in [6,Remark 5.3] such a method may not converge.

On the other hand, [6, Section 5] presents a convergent generalization of Howard's algorithm for max-min problems, which we readily adapt here. We will present an algorithm which requires the solution, at every iteration step, of a Hamilton Jacobi Bellman equation, which can be realized via a semi-smooth Newton method. We will also comment on an algorithm with inexact solves.

For a given w h ∈ V 0 h and z ∈ N h define α(w h , z) ∈ A as the element that infimizes the supremum of the integral operators when applied to w h at the point z, that is Our algorithm can then be described as follows:

• Initialization: Choose w -1 h ∈ V 0 h and set α 0 = α(w -1 h ). • Iteration: For k ≥ 0 find w k h ∈ V 0 h that solves

and set α k+1 = α(w k h ).

(5.3)

• Convergence test: If F ε h [w k h ](z) = f z for all z ∈ N h stop. Notice that (5.3) is equivalent to

The analysis of the algorithm (5.2)-(5.3) relies on the following properties of the operators F α h,ε . Lemma 5.1 (monotonicity and comparison). For every α ∈ A #N h the operator F α h,ε is monotone and satisfies a comparison principle, i.e., if v h , w h ∈ V 0 h are such that

Proof. The proof repeats the arguments of Lemma 3.2 and Corollary 3.3. For brevity, we skip details.

The comparison principle will allow us to obtain convergence. Theorem 5.2 (convergence). The sequence {w k h } k≥0 ⊂ V 0 h obtained by algorithm (5.2)-(5.3) converges in a finite number of steps to u ε h ∈ V 0 h , solution of (3.3).

Proof. We proceed in two steps: First, we show that the sequence is monotone, i.e., for every z ∈ N h , w k h (z) ≥ w k+1 h (z). By construction:

Subtracting f z from this inequality we realize that

which by the comparison principle established in Lemma 5.1 implies w k h ≥ w k+1 h over N h .

Next, since the set A is finite, there are at most (#A) #N h different variables and there must be two indices κ > for which (5.3) yields α κ = α . This implies that w κ h = w h and by monotonicity w k h = w h for all k ≥ l. But then, by uniqueness w h = u h . This concludes the proof. Notice that (5.2) requires the exact solution of a discrete version of a Hamilton Jacobi Bellman problem, which can be achieved by Howard's algorithm. We also propose a scheme with inexact solves in this step:

• Initialization: Choose w -1 h ∈ V 0 h and set α 0 = α(w -1 h ).

Set α k+1 = α(w k h ).

(5.5)

• Convergence test: If F ε h [w k h ](z) = f z for all z ∈ N h stop. The convergence of this algorithm follows mutatis mutandis the proof of Theorem 5.4 of [6].

Theorem 5.3 (convergence with inexact solves). Assume that the sequence of errors {η k } k∈N ∈ 1 . Then the sequence {w k h } k≥0 ⊂ V 0 h , obtained by algorithm (5.4)-(5.5), converges to u ε h ∈ V 0 h , solution of (3.3).