Stability and convergence of second order backward differentiation schemes for parabolic Hamilton-Jacobi-Bellman equations

We study a second order BDF (Backward Differentiation Formula) scheme for the numerical approximation of parabolic HJB (Hamilton-Jacobi-Bellman) equations. The scheme under consideration is implicit, non-monotone, and second order accurate in time and space. The lack of monotonicity prevents the use of well-known convergence results for solutions in the viscosity sense. In this work, we establish rigorous stability results in a general nonlinear setting as well as convergence results for some particular cases with additional regularity assumptions. While most results are presented for one-dimensional, linear parabolic and non-linear HJB equations, some results are also extended to multiple dimensions and to Isaacs equations. Numerical tests are included to validate the method.

1. Introduction. This paper provides stability and convergence results for a type of implicit finite difference scheme for the approximation of nonlinear parabolic equations using backward differentiation formulae (BDF).
In particular, we consider Hamilton-Jacobi-Bellman (HJB) equations of the following form: is a second order differential operator. Here, (Σ) ij is symmetric non-negative definite for all arguments. Linear parabolic equations, corresponding to the case |Λ| = 1, are a special case for which more comprehensive results are obtained in the paper. It is well known that for nonlinear, possibly degenerate equations the appropriate notion of solutions to be considered is that of viscosity solutions [8]. We assume throughout the whole paper the well-posedness of the problem, namely the existence and uniqueness of a solution in the viscosity sense.
Under such weak assumptions, convergence of numerical schemes can only be guaranteed if they satisfy certain monotonicity properties, in addition to the more standard consistency and stability conditions for linear equations [2]. This in turn reduces the obtainable consistency order to 1 in the general case [11].
On the other hand, in many cases -especially in non-degenerate ones -solutions exhibit higher regularity and are amenable to higher order approximations. The existence of classical solutions and their regularity properties under a strict ellipticity condition have been investigated, for instance, in [14,10].
The higher order of convergence in both space and time of discontinuous Galerkin approximations is demonstrated theoretically and empirically in [17] for sufficiently regular solutions under a Cordes condition for the diffusion matrix, a measure of the ellipticity. More recently, it was shown empirically in [6] that schemes based on first derivative approximations in time and space based on a second order backward differentiation formula (see, e.g., [19], Section 12.11, for the definition of BDF schemes for ODEs) have good convergence properties. In particular, in a non-degenerate controlled diffusion example therein where the second order, non-monotone Crank-Nicolson scheme fails to converge, the (also non-monotone) BDF2 scheme shows second order convergence.
For constant coefficient parabolic PDEs, the L 2 -stability and smoothing properties of the BDF scheme are a direct consequence of the strong A-stability of the scheme. Moreover, [3] shows that for the multi-dimensional heat equation the BDF time stepping solution and its first numerical derivative are stable in the maximum norm. The technique, which is strongly based on estimates for the resolvent of the discrete Laplacian, do not easily extend to variable coefficients or the nonlinear case.
A more general linear parabolic setting is considered in [4], where second order convergence is shown for variable timestep using energy techniques. This result is extended to a semi-linear example in [9]; the application to incompressible Navier-Stokes equations has been analysed in [13]. In [5], a closely reltated BDF scheme is studied for a diffusion problem with an obstacle term (which includes the American option problem in mathematical finance).
The scheme we propose is constructed by using a second order BDF approximation for the first derivatives in both time and space. Combining this with the standard three-point central finite difference for the second spatial derivative in one dimension, the scheme is second order consistent by construction.
For this scheme, we establish new stability results in the H 1 -and L 2 -norms (see Theorems 4 and 6, respectively) for linear parabolic PDEs and their nonlinear HJB counterpart. These generalize some results of [4,9,5] to more general non-linear situations. From this analysis we deduce error bounds for classical smooth and piecewise smooth solutions (see Theorems 16 and 18). Extensions of the results to Isaacs equations and the two-dimensional case are also given.
The outline of the paper is as follows. In Section 2, we define some specific BDF schemes and state the main results concerning well-posedness and stability in discrete H 1 -or L 2 -norms. In Sections 3 and 4 we prove the main results and give an extension from HJB to Isaacs equations. In Section 5, we give further stability results in the discrete L 2 -norm, which are weaker in the sense that they hold only for uncontrolled Lipschitz regulary diffusion coefficients, but stronger in the sense that they allow for degenerate diffusion and can be extended to two dimensions. In Section 6, we deduce error estimates from the stability results and from the truncation error of the scheme for sufficiently regular solutions. Section 7 studies carefully two numerical examples, the Eikonal equation and a second order equation with controlled diffusion. Section 8 concludes. An appendix contains a proof of the existence of solutions for our schemes.
2. Definition of the scheme and main result. We focus in the first instance on the one-dimensional equation It is known (see Theorem A.1 in [1]) that with the following assumptions: -Λ is a compact set, there exists a unique bounded continuous viscosity solution of (2). We will make individual assumptions for each result as we go along, but in general assume a unique and continuous solution (e.g. to define the classical truncation error).
2.1. The BDF2 scheme. For the approximation in the x variable, we will consider the PDE on a truncated domain Ω : Let N ∈ N * ≡ N\{0} the number of time steps, τ := T /N the time step size, and t n = nτ , n = 1, . . . , N . Let I ∈ N * the number of interior mesh points, and define a uniform mesh (x i ) 1≤i≤I with mesh size h by Hereafter, we denote by u a numerical approximation of v, the solution of (1), i.e.
For each time step t k , the unkowns are the values u k i for i = 1, . . . , I. Standard Dirichlet boundary conditions use the knowledge of the values at the boundary, v(t, x min ) and v(t, x max ). Here, as a consequence of the size of the stencil for the spatial BDF2 scheme below, we will assume that values at the two left-and rightmost mesh points are given, that is, v(t, x j ) for j ∈ {−1, 0} as well as j ∈ {I +1, I +2} are known (corresponding to the values at the points (x −1 , x 0 , x I+1 , x I+2 ) ≡ (x min − h, x min , x max , x max + h)). 1 We then consider the following scheme, for k ≥ 2, i ∈ I, where we denote as usual by [u] k i the numerical solution excluding at (t k , x i ), and (the usual second order approximation of v xx ), b + := max(b, 0) and b − := max(−b, 0) denote the positive and negative part of b, respectively, and where a second order leftor right-sided BDF approximation is used for the first derivative in space: Note in particular the implicit form of the scheme (3). The existence of a unique solution of this nonlinear implicit scheme will be addressed later on. We will also define the numerical Hamiltonian associated with the scheme: As discussed above, the scheme is completed by the following boundary conditions: Since (3) is a two-step scheme, for the first time step k = 1, i ∈ I, we use a backward Euler step, is given by the initial condition (2b).
Remark 1. As the backward Euler step is only used once, it does not affect the overall second order of the scheme.
Remark 2. Most of our results also apply to the scheme obtained by replacing the BDF approximation (4) of the drift term by a centred finite difference approximation: However, numerical tests (see Section 7.1) show that the BDF upwind approximation as in (4) has a better behaviour in some extreme cases where the diffusion vanishes. We shall give a rigorous stability estimate for the BDF scheme in the linear case even for possibly vanishing diffusion (Section 5.2).

Definitions and main results.
In the remainder of this paper, we prove various stability and convergence results for the scheme (3). We state in this section the first main well-posedness and stability results.
Let u denote the solution of (3) and let v be the solution of (1). The error associated with the scheme is then defined by For any function φ we will also use the notation φ k i := φ(t k , x i ) as well as φ k := (φ k i ) 1≤i≤I and [φ] k i := (φ m j ) (j,m) =(i,k) , and the error vector at time t k is defined by The consistency error will be denoted by E k (φ) := (E k i (φ)) 1≤i≤I ∈ R I and is defined in the classical way as follows, for any smooth enough function φ: By extension, for the exact solution v of (1), we will simply define Note that (9) is well-defined for any continuous function.
In particular for the scheme (3) it is clear that we have second order consistency in space and time, that is, for sufficiently regular data φ. Throughout the paper, A will denote the finite difference matrix associated to the second order derivative, i.e.
Let x, y A := x, Ay . Then we consider the A-norm defined as follows: (with the convention in (12) that x 0 = x I+1 = 0). Hence, √ h|x| A approximates the H 1 semi-norm in Ω. Similarly, we will consider later the standard Euclidean norm defined by x 2 := x, x , such that √ h x approximates the L 2 -norm. Our first result concerns the solvability of the numerical scheme S (τ,h) (seen as an equation for u k , with [u] k given) and is the following.
The scheme is hence well-defined even if σ vanishes. A uniform ellipticity condition for σ will be needed for proving the H 1 stability of the scheme.
We provide a relaxation of the ellipticity condition for stability in the Euclidean norm in Section 5.2.
Our main stability result is the following.
Theorem 4. Assume (A1), (A2), as well as the CFL condition (13). Then there exists a constant C ≥ 0 (independent of τ and h) and τ 0 > 0 such that, for any τ ≤ τ 0 , The proof of Theorem 4 will be the subject of Section 4.
Remark 5. As a consequence of the stability result and under further mild regularity assumptions on the boundary data, we can deduce that the scheme (3) is A-norm bounded: where the constant C depends only on T and on the data but not on τ and h.
The analysis of the controlled case is made complicated by the fact that even if the solution to (2) is classical and the supremum is attained for each x and t (and similarly for each i and k in (3)), we cannot make any assumptions on the regularity of this optimal control as a function of x and t (or i and k, respectively).
In certain circumstances, the previous bound holds with the A-norm replaced by the Euclidean norm. In particular, we consider the following assumption: Assumption (A3). The diffusion coefficient is independent of the control, i.e. σ ≡ σ(t, x) and there exists L ≥ 0 such that Theorem 6. Assume (A1), (A2), (A3), as well as the CFL condition (13). Then there exists C ≥ 0 (independent of τ and h) and τ 0 > 0 such that, for any τ ≤ τ 0 , As a consequence, error estimates will be obtained under the main assumptions (A1), (A2) and (A3) or under some specific assumptions, see Sections 5 and 6.
The extension of the presented results to other type of nonlinear operators (inf, sup inf or inf sup) and corresponding equations will also be discussed.
3. Proof of Theorem 3 (well-posedness of the scheme). The scheme (3) at time t k (for k ≥ 2) can be written in the following form: a ∈ R I and M k a ∈ R I×I with the following non-zero entries: . For k = 1, the terms are different but the form (and analysis) is similar. The fact that (M a ) i,i±2 are nonnegative breaks the monotonicity of the scheme and makes the analysis more difficult. We will use the following lemma, whose proof is given in appendix A: Asssume that Λ is some set, (q a ) a∈Λ is a family of vectors in R I , (M a ) a∈Λ is a family of matrices in R I×I such that: Then there exists a unique solution X in R n of Remark 8. For a fixed a ∈ Λ, we have Proof of Theorem 3. We are going to prove properties (i) and (ii) in Lemma 7. Condition (i) is immediately verified, and we turn to proving (ii). We have (omitting the dependency on k and a in σ, b ± , r) and By the CFL condition (13), there exists > 0 such that and therefore Then by using Taking τ small enough such that for instance 2 + τ r ≥ 4 , and since b(.) and σ(.) are bounded functions (by (A1)), we obtain the bound Since the last bound is a constant < 1, we can apply Lemma 7 to obtain the existence and uniqueness of the solution of the BDF2 scheme.

Proof of Theorem 4 (stability in the A-norm)
. The proof consists of three main steps: first, we show a "linear" recursion for the error (Lemma 9); second, we pass from such a recursion for the error in vector form to a scalar recursion (Lemma 10); finally, we show the stability estimate from this scalar recursion (Lemma 11).
4.1. Treatment of the nonlinearity. First, we have the following: Lemma 9. Let u be the solution of scheme (3) and v the solution of equation (2).
Proof. By definition of the consistency error (9), one has (for k ≥ 2, The scheme simply reads , the following recursion is obtained for the error in R I : For simplicity of presentation, we first consider the case when b and r vanish, i.e. b(.) ≡ 0 and r(.) ≡ 0. In this case, To simplify the presentation, we will assume that σ and are continuous functions of a so that the supremum is attained. 2 For each given k, i, let thenā k i ∈ Λ denote an optimal control in (26). In the same way, letb k i denote an optimal control for H[v k ] i . By using the optimality ofā k i , it holds The general case is obtained easily by considering sequences of -optimal controls and letting → 0, such that (30) below still holds for a suitably definedσ 2 ,b + ,b − ,r. and, in the same way, Therefore, combining (27) and (28), In particular, we can write Isaacs equations. The same technique used above to deal with the nonlinear operator applies also to Isaacs equations, i.e. equations of the following form: To simplify the presentation, let us consider again b, r ≡ 0, and now also ≡ 0. By analogous definitions and reasoning to above, we get (25), where, for φ = u, v, Let (ā k i ,b k i ) ∈ Λ 1 × Λ 2 denote an optimal control in (32). 3 One has Analogously, one can prove )). At this point, it is sufficient to take forâ k i ∈ Λ 1 andb k i ∈ Λ 2 optimal controls in (33) and (34), respectively, to be able to write Proof. For simplicity of presentation we will assume that b has constant positive sign. The terms coming from the negative part of b can be treated in a similar way. We remark that for E ∈ R I , −D 2 E = AE, where A is the finite difference matrix defined in (11). By (22), we get the following: where We form the scalar product of (36) with AE k . By using the identity 2 a − b, a A = |a| 2 A + |a − b| 2 A − |b| 2 A , one has: where we have also used |a + b| 2 ≤ 2|a| 2 + 2|b| 2 . From (σ 2 ) k i ≥ η > 0 for all k, i: where · denotes the canonical Euclidean norm in R I . In order to estimate the drift component, let us introduce the notation with the convention that E i = 0 for all indices i which are not in I. It holds: By using the boundedness of the drift term, and δE k , For the last term, using the boundedness of r and the Cauchy-Schwarz inequality, Therefore, putting (38), (40) and (41) together, Easy calculus shows that the minimal eigenvalue of A is λ min (A) = 4 h 2 sin 2 ( πh 2 ) ≥ 4. Hence X, AX ≥ 4 X, X and therefore X ≤ 1 2 |X| A . In the same way, we have also |X| A ≤ 1 2 AX . Hence it holds Then, combining (37) and (44), we obtain the desired inequality with C := 2 η C 2 1 .

4.4.
A universal stability lemma. In the following, it is assumed that | · | is any vectorial norm. We will use the result for the canonical Euclidean norm | · | ≡ · and the A-norm | · | ≡ | · | A .
In order to prove the following Lemma 11, we will exploit properties of the matrix in particular the fact that (M τ ) −1 ≥ 0 for τ small enough (which we prove).

Stability in the Euclidean norm.
The fundamental stability result given by Lemma 11 applies to any vectorial norm. In this section, we discuss some special cases where (46) can be obtained for the Euclidean norm | · | = · .
We first prove the stability result for this norm under the extra assumption (A3), i.e., the control may appear except in the diffusion term, which must also be Lipschitz continuous in the following proof.

Proof of Theorem 6 (stability in the Euclidean norm) .
We consider the scalar product of (36) directly with E k (instead of AE k previously used), again in the situation where b ≥ 0 to simplify the argument. We obtain: As in Section 4.3, we have We now focus on bounding the other terms on the left-hand side of (53). By using the Lipschitz continuity of σ 2 one has Therefore, by the Cauchy-Schwarz inequality, one obtains where δE k is defined by (39). Moreover, for the first order term one has where for the last equality we have used that δ 2 E k ≤ δE k . Putting together estimates (55) and (56), using the fact that E k , R k E k ≥ − r ∞ E k 2 , we get where we have denoted C 1 := L + 4 b ∞ and have used again the Cauchy-Schwarz inequality. Hence, together with (54), this gives (46) with | · | = · and the constant C := 4( 4η + r ∞ ). By using Lemma 11, this concludes the proof of Theorem 6. 2 5.2. Linear equation with degenerate diffusion term. The next result concerns the case of a possibly degenerate diffusion term. It will require more restrictive assumptions on the drift and diffusion terms, and we shall assume that there is no control here. Indeed, in this case, one cannot count on the positive term coming from the non-degenerate diffusion which, in the proof of Theorem 6, is used to compensate the negative correction terms coming from the drift term. This leads us to consider the following assumptions: Assumption (A4). r is bounded. The drift and diffusion coefficients are independent of the control, i.e. b ≡ b(t, x) and σ ≡ σ(t, x), and there exist L 1 , L 2 ≥ 0 such that, for all t, x, h: (The last condition is equivalent to (σ 2 ) xx ≥ −L 2 in the differentiable case.) Proposition 12. Let assumption (A4) be satisfied. Then (46) holds for |·| = · .
Proof. We consider again the scalar recursion (53). For any vector E = (E i ) 1≤i≤I (with E j = 0 for j ∈ {−1, 0, I + 1, I + 2}), it holds: Hence, by the semi-concavity assumption (58) on σ 2 , Now we focus on a lower bound for E k , F k BE k . Let We assume again b i ≥ 0 for all i to simplify the presentation. The case where b i ≤ 0 for some i is similar. Then, the following bound holds: . Then, by the Lipschitz continuity of b(.) and the bound y k By combining the bounds (59) and (60), we obtain Therefore, inequality (46) is obtained with C := 4( L2 4 + 3L 1 + r ∞ ), which leads to the desired stability estimate.

Extension to a two-dimensional case.
Under suitable assumptions, the result of Theorem 6 can be extended to multi-dimensional equations. The nonlinearity can be treated exactly as in Section 4.1 (or 4.2), so that we can focus on the linear for a positive definite matrix Σ and a drift vector b. For simplicity, we furthermore consider the two-dimensional case d = 2, with r, ≡ 0, and omit the dependence of the coefficients on the time variable, then with where σ 1 , σ 2 ≥ 0 and ρ ∈ [−1, 1] is the correlation parameter, the equation reads The computational domain is given by Ω := (x min , x max ) × (y min , y max ). We introduce the discretization in space defined by the steps h x , h y > 0 and we denote by G (hx,hy) the associated mesh. In what follows, given any function φ of (x, y) ∈ Ω, we will denote φ ij = φ(x i , y j ) for (i, j) ∈ I := I 1 × I 2 , where I 1 = {1, . . . , I 1 }, I 2 = {1, . . . , I 2 }.
Assuming that ρ ≥ 0 everywhere (the case when ρ ≤ 0 is similar), we consider a 7-point stencil for the second order derivatives (see [12,Section 5.1.4] and the BDF approximation of the first order derivatives The scheme is therefore defined, for k ≥ 2, by A straightforward calculation shows that the second order term also reads The scheme is completed with the following boundary conditions: For simplicity, assume h x = h y =: h. We consider the following assumptions: We then have the following result. The proof is similar to the one of Theorem 6, using (62) with α ij , β ij , γ ij ≥ 0 by assumption (A2'), and is therefore omitted.
Proposition 13. Let assumptions (A1'),(A2') and (A3') be satisfied. Then the stability estimate (47) holds for | · | = · . Remark 14. (i) If h x = h y and for instance h y = Ch x for some C ≥ 1, (A2') has to hold with σ 2 replaced by σ 2 /C as a result of the scaling properties of the scheme.
(ii) Observe that assumption (A2') is equivalent to requiring strong diagonal dominance of the covariance matrix.
(iii) When the strong diagonal dominance of the matrix Σ is not guaranteed, one can consider the generalized finite difference scheme in [7]. However, determining the precise set of assumptions on the coefficients needed to apply the previous arguments does not seem easy from the construction in [7].
6. Error estimates. In this section, we give detailed error estimates for the implicit BDF2 scheme (3). We consider the following rescaled norms on R I : corresponding to discrete approximations of L 2 (Ω)-and H 1 (Ω) norms, respectively. Both these norms will be used in the forthcoming numerical section.
In addition, we define the following semi-norm on some interval I = (a, b): For a given open subset Ω * T of (0, T )×Ω, we define C k, (Ω * T ) as the set of functions v : Ω * T → R which admit continuous derivatives ( ∂ i v ∂t i ) 0≤i≤k and ( ∂ j v ∂x j ) 0≤j≤ on Ω * T . We also denote by C k, b (Ω * T ) the subset of functions with bounded derivatives on Ω * T . Assumption (A5). v ∈ C 1,2 ((0, T ) × Ω) and for some constant C ≥ 0: Remark 15. By results in [10] and [14], assumption (A5) is satisfied for sufficiently smooth data and given a uniform ellipticity condition.
where C is a constant which depends on the derivatives of v of order 3 and 4 in t and x, respectively. (ii) If (A5) holds for some δ ∈ (0, 1], then the numerical solution u of (3), (5) converges to v in the L 2 -norm with for some constant C (possibly different from the one in (A5)).
Proof. We first prove (ii). By Taylor expansion, we can write for instance, for some Similarly, using the higher spatial regularity, there exists a constant The result (ii) now follows directly by inserting the obtained truncation error into the stability estimate of Theorem 6. For the proof of (i) (smooth case), expansion up to order 3 and 4 gives the truncation error of higher order for k ≥ 2, and we use the fact that the error from the first backward Euler step is bounded by E 1 ≤ Cτ (τ + h 2 ); in particular, we use that (E 1 − E 0 )/τ + (∆ 1 A + F 1 B + R 1 )E 1 = −E 1 , with E 1 ≤ C(τ + h 2 ), E 0 = 0 and the bound is otherwise similar and simpler than that for k ≥ 2.
The previous arguments can also be used to derive error estimates for piecewise smooth solutions. In this case, we will need to limit the number of non-regular points that may appear in the exact solution (assumption (A6)(i) is similar to [5]).
Such a situation will be illustrated in the numerical example of Section 7.2.
Remark 19. (i) Similar results can be derived for errors in the A-norm, however derivatives of one order higher are required due to the derivative in the definition of the norm. (ii) The estimates in Theorem 16 are not always sharp, as symmetries and the smoothing behaviour of the scheme can result in higher order convergence. We discuss such special cases for Examples 1 and 2 in Section 7, Remarks 21 and 22, respectively. (iii) These error estimates can be compared with [5], where an error bound of order h 1/2 was obtained for diffusion problems with an obstacle term, under the main assumption that v xx is a.e. bounded with a finite number of singularities (instead of (A5)) . In the present context it seems natural to assume the Hölder regularity of u t and u xx coming from the ellipticity assumption (see Remark 15).
7. Numerical tests. We now compare the performance of the BDF2 scheme with other second order finite difference schemes on two examples. 7.1. Test 1: Eikonal equation. The first example is based on a deterministic control problem (σ ≡ 0) and motivates the choice of the BDF2 approximation for the drift term in (4), compared to the more classical centered scheme (7). We consider with v 0 (x) = max(0, 1 − x 2 ) 4 and T = 0.2. The initial datum is shown in Figure 1 (dashed line). The exact solution is Remark 20. The Eikonal equation can be written as v t +max a∈{−1,1} (av x ) = 0 in HJB form. Note that our theoretical analysis does not cover this example, however, since in the degenerate case assumption (A4) is required, which is not satisfied here. In Figure 1, we show the results obtained at the terminal time T = 0.2 using schemes (3)-(7) (left) and (3)-(4) (right) with τ /h = 0.5. We numerically observe that the centered approximation generates undesirable oscillations, whereas the BDF2 scheme is stable.
As stated in Theorem 3, in case of a degenerate diffusion, a CFL condition of the form τ ≤ Ch has to be satisfied for well-posedness of the BDF2 scheme. Table 1 shows numerical convergence of order 2 in both time and space, although the solution is globally only Lipschitz.
Remark 21. The full convergence order here is due to the particular symmetry of the solution. To confirm this, we report in Table 2 the results obtained for the same equation with initial data v(0, x) = − max(0, 1 − x 2 ) 4 (see also Figure 2). In this case, there is no such symmetry around the two singular points and as a result the full convergence order is lost: the scheme is globally only of order 1 in the H 1 norm and roughly 1.5 in the L 2 and L ∞ norm.  Table 2
In spite of the apparent simplicity of the equation under consideration, in [16] an example of non-convergence of the Crank-Nicolson scheme is given for a similar optimal control problem. The BDF2 scheme, in contrast, has shown good performance for that same problem in [6]. Figure 3 (top row) shows the initial data and the value function at terminal time computed using the BDF2 scheme. The error and convergence rate in different norms are reported in Table 3. Here an accurate numerical solution computed by an implicit Euler scheme (in order to ensure convergence) is used for comparison.
Taking τ ∼ h the BDF2 scheme gives clear second order convergence, see Table 3. This is not the case for CN as shown in Table 4. The CN scheme also exhibits some instability in the second order derivative for high CFL number, i.e. τ /h, see Figure 3 (this is analogous to the finding in [16]). One can verify that for a small CFL number, i.e. τ ∼ h 2 , the CN scheme shows second order of convergence.
Remark 22. In this example, due to the strict ellipticity, Assumption (A5) is guaranteed for some δ > 0 (see Remark 15). Then Theorem 16 gives convergence   with order δ. Furthermore, Fig. 3, bottom row, suggests Hölder continuity of u xx in x, which is expected by virtue of the control being piecewise constant. Therefore, we conjecture that Assumption (A6) is satisfied, such that Theorem 18 would give the higher order 1/2 + δ. In the test, in fact the full order 2 is observed (see Table 3).

Conclusion.
We have proved the well-posedness and stability in L 2 and H 1 norms of a second order BDF scheme for HJB equations with enough regularity of the coefficients. The significance of the results is that this was achieved for a second   (6) with I + 1 = 20 × 2 9 , N = 2 22 is used. order (and hence) non-monotone scheme. For smooth or piecewise smooth solutions, as is often the case, one can use the recursion we derived to bound the error of the numerical solution in terms of the truncation error of the scheme. The latter depends on the regularity of the solution and has to be estimated for individual examples.
The numerical tests demonstrate convergence at least as good as predicted by the theoretical results, and often better, due to symmetries of the solution or smoothing properties of the equation and the scheme. This is in contrast to some alternative second order schemes, such as the central spatial difference in the case of a first order equation, or the Crank-Nicolson time stepping scheme for a second order equation, which can show poor or no convergence.
Appendix A. Proof of Lemma 7.
In order to prove the existence and uniqueness of a solution to (21), we consider a fixed-point approach. The initial problem (21) can be written as follows: sup a∈Λ (L a X − (q a − U a X)) = 0, where L a and U a are two matrices such that M a ≡ L a + U a . We consider in particular L a to be the lower triangular part of M a including the diagonal terms, (L a ) ij := (M a ) ij 1 i≥j , and U a the remaining upper triangular part, (U a ) ij := (M a ) ij 1 i<j . For a given vector c ∈ R I , let g(c) := X denote the (unique) solution of the following simplified problem: Therefore, solving (64) amounts to solving g(X) = X. By elementary compuations one can show that g is δ-Lipschitz for the . ∞ norm, with δ := sup a (L a ) −1 U a ∞ . For a diagonally dominant matrix, the following classical estimate holds (this is related to the Gauss-Seidel relaxation method; see for instance, Th. 8.2.12 in [18]). By using the assumptions on the matrices M a , we have δ < 1. Hence, g is a contraction mapping on R I and therefore we obtain the existence and uniqueness of a solution of (64) as desired. 2