Stability and convergence of second order backward differentiation schemes for parabolic Hamilton-Jacobi-Bellman equations
SSTABILITY AND CONVERGENCE OF SECOND ORDERBACKWARD DIFFERENTIATION SCHEMES FOR PARABOLICHAMILTON-JACOBI-BELLMAN EQUATIONS
OLIVIER BOKANOWSKI ∗ , ATHENA PICARELLI † , AND
CHRISTOPH REISINGER † Abstract.
We study a second order BDF (Backward Differentiation Formula) scheme for thenumerical approximation of parabolic HJB (Hamilton-Jacobi-Bellman) equations. The scheme underconsideration is implicit, non-monotone, and second order accurate in time and space. The lack ofmonotonicity prevents the use of well-known convergence results for solutions in the viscosity sense. Inthis work, we establish rigorous stability results in a general nonlinear setting as well as convergenceresults for some particular cases with additional regularity assumptions. While most results arepresented for one-dimensional, linear parabolic and non-linear HJB equations, some results are alsoextended to multiple dimensions and to Isaacs equations. Numerical tests are included to validatethe method.
1. Introduction.
This paper provides stability and convergence results for atype of implicit finite difference scheme for the approximation of nonlinear parabolicequations using backward differentiation formulae (BDF).In particular, we consider Hamilton-Jacobi-Bellman (HJB) equations of the fol-lowing form: v t ( t, x ) + sup a ∈ Λ (cid:110) L a [ v ]( t, x ) + r ( t, x, a ) v + (cid:96) ( t, x, a ) (cid:111) = 0 , (1)where ( t, x ) ∈ [0 , T ] × R d , Λ ⊂ R m is a compact set and L a [ v ]( t, x ) = −
12 tr[Σ( t, x, a ) D x v ( t, x )] + b ( t, x, a ) D x v ( t, x ) is a second order differential operator. Here, (Σ) ij is symmetric non-negative definitefor all arguments. Linear parabolic equations, corresponding to the case | Λ | = 1 , area special case for which more comprehensive results are obtained in the paper.It is well known that for nonlinear, possibly degenerate equations the appropriatenotion of solutions to be considered is that of viscosity solutions [8]. We assumethroughout the whole paper the well-posedness of the problem, namely the existenceand uniqueness of a solution in the viscosity sense.Under such weak assumptions, convergence of numerical schemes can only beguaranteed if they satisfy certain monotonicity properties, in addition to the morestandard consistency and stability conditions for linear equations [2]. This in turnreduces the obtainable consistency order to 1 in the general case [11].On the other hand, in many cases – especially in non-degenerate ones – solutionsexhibit higher regularity and are amenable to higher order approximations. Theexistence of classical solutions and their regularity properties under a strict ellipticitycondition have been investigated, for instance, in [14, 10].The higher order of convergence in both space and time of discontinuous Galerkinapproximations is demonstrated theoretically and empirically in [17] for sufficientlyregular solutions under a Cordes condition for the diffusion matrix, a measure of the ∗ Laboratoire J.-L. Lions, Université Pierre et Marie Curie 75252 Paris Cedex 05, France, andUFR de Mathématiques, Site Chevaleret, Université Paris-Diderot, 75205 Paris Cedex, France([email protected]). † Mathematical Institute, University of Oxford, Andrew Wiles Building, Woodstock Rd, OxfordOX2 6GG, UK ([email protected], [email protected])1 a r X i v : . [ m a t h . NA ] F e b O. BOKANOWSKI, A. PICARELLI, C. REISINGER ellipticity. More recently, it was shown empirically in [6] that schemes based on firstderivative approximations in time and space based on a second order backward differ-entiation formula (see, e.g., [19], Section 12.11, for the definition of BDF schemesfor ODEs) have good convergence properties. In particular, in a non-degeneratecontrolled diffusion example therein where the second order, non-monotone Crank-Nicolson scheme fails to converge, the (also non-monotone) BDF2 scheme shows sec-ond order convergence.For constant coefficient parabolic PDEs, the L -stability and smoothing prop-erties of the BDF scheme are a direct consequence of the strong A-stability of thescheme. Moreover, [3] shows that for the multi-dimensional heat equation the BDFtime stepping solution and its first numerical derivative are stable in the maximumnorm. The technique, which is strongly based on estimates for the resolvent of thediscrete Laplacian, do not easily extend to variable coefficients or the nonlinear case.A more general linear parabolic setting is considered in [4], where second orderconvergence is shown for variable timestep using energy techniques. This result isextended 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 isstudied for a diffusion problem with an obstacle term (which includes the Americanoption problem in mathematical finance).The scheme we propose is constructed by using a second order BDF approximationfor the first derivatives in both time and space. Combining this with the standardthree-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 - and L -norms (seeTheorems 4 and 6, respectively) for linear parabolic PDEs and their nonlinear HJBcounterpart. These generalize some results of [4, 9, 5] to more general non-linearsituations. From this analysis we deduce error bounds for classical smooth and piece-wise smooth solutions (see Theorems 16 and 18). Extensions of the results to Isaacsequations and the two-dimensional case are also given.The outline of the paper is as follows. In Section 2, we define some specific BDFschemes and state the main results concerning well-posedness and stability in discrete H - or L -norms. In Sections 3 and 4 we prove the main results and give an extensionfrom HJB to Isaacs equations. In Section 5, we give further stability results in thediscrete L -norm, which are weaker in the sense that they hold only for uncontrolledLipschitz regulary diffusion coefficients, but stronger in the sense that they allow fordegenerate diffusion and can be extended to two dimensions. In Section 6, we deduceerror estimates from the stability results and from the truncation error of the schemefor sufficiently regular solutions. Section 7 studies carefully two numerical examples,the Eikonal equation and a second order equation with controlled diffusion. Section 8concludes. 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 instanceon the one-dimensional equation v t + sup a ∈ Λ (cid:18) − σ ( t, x, a ) v xx + b ( t, x, a ) v x + r ( t, x, a ) v + (cid:96) ( t, x, a ) (cid:19) = 0 ,t ∈ [0 , T ] , x ∈ R , (2a) v (0 , x ) = v ( x ) x ∈ R . (2b)It is known (see Theorem A.1 in [1]) that with the following assumptions:– Λ is a compact set, DF2 SCHEMES FOR HJB EQUATIONS
3– for some C > the functions φ ≡ σ, b, r, (cid:96) : [0 , T ] × R × Λ → R and v : R → R satisfy for any t, s ∈ [0 , T ] , x, y ∈ R , a ∈ Λ | v ( x ) | + | φ ( t, x, a ) | ≤ C , | v ( t, x ) − v ( s, y ) | + | φ ( t, x, a ) − φ ( s, y, a ) | ≤ C ( | x − y | + | t − s | / ) , there exists a unique bounded continuous viscosity solution of (2).We will make individual assumptions for each result as we go along, but in generalassume a unique and continuous solution (e.g. to define the classical truncation error). For the approximation in the x variable, we willconsider the PDE on a truncated domain Ω := ( x min , x max ) , where x min < x max .Let N ∈ N ∗ ≡ N \{ } 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 auniform mesh ( x i ) ≤ i ≤ I with mesh size h by x i := x min + ih, i ∈ I = { , . . . , I } , where h := x max − x min I + 1 . Hereafter, we denote by u a numerical approximation of v , the solution of (1), i.e. u ki ∼ v ( t k , x i ) . For each time step t k , the unkowns are the values u ki for i = 1 , . . . , I .Standard Dirichlet boundary conditions use the knowledge of the values at theboundary, v ( t, x min ) and v ( t, x max ) . Here, as a consequence of the size of the stencil forthe spatial BDF2 scheme below, we will assume that values at the two left- and right-most mesh points are given, that is, v ( t, x j ) for j ∈ {− , } as well as j ∈ { I +1 , I +2 } are known (corresponding to the values at the points ( x − , x , x I +1 , x I +2 ) ≡ ( x min − h, x min , x max , x max + h ) ). We then consider the following scheme, for k ≥ , i ∈ I , S ( τ,h ) ( t k , x i , u ki , [ u ] ki ) = (3) u ki − u k − i + u k − i τ + sup a ∈ Λ (cid:110) L a [ u k ]( t k , x i ) + r ( t k , x i , a ) u ki + (cid:96) ( t k , x i , a ) (cid:111) = 0 , where we denote as usual by [ u ] ki the numerical solution excluding at ( t k , x i ) , and L a [ u ]( t k , x i ) := − σ ( t k , x i , a ) D u i + b + ( t k , x i , a ) D , − u i − b − ( t k , x i , a ) D , + u i ,D u i := u i − − u i + u i +1 h , (the usual second order approximation of v xx ), b + := max( b, and b − := max( − b, denote the positive and negative part of b , respectively, and where a second order left-or right-sided BDF approximation is used for the first derivative in space: D , − u i := 3 u i − u i − + u i − h and D , + u i := − (cid:18) u i − u i +1 + u i +2 h (cid:19) . (4) In practice, this means that a sufficiently accurate approximation of these “boundary values”has to be available. Boundary approximations with modified schemes are commonly used and arenot the focus of this paper; it is seen in [15] that the use of a lower order scheme in the vicinity ofthe boundary does not affect the global provable convergence order.
O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Note in particular the implicit form of the scheme (3). The existence of a uniquesolution of this nonlinear implicit scheme will be addressed later on.We will also define the numerical Hamiltonian associated with the scheme: H [ u ]( t k , x i ) := sup a ∈ Λ (cid:110) L a [ u ]( t k , x i ) + r ( t k , x i , a ) u i + (cid:96) ( t k , x i , a ) (cid:111) . As discussed above, the scheme is completed by the following boundary conditions: u ki := v ( t k , x i ) , ∀ i ∈ {− , } ∪ { I + 1 , I + 2 } . Since (3) is a two-step scheme, for the first time step k = 1 , i ∈ I , we use a backwardEuler step, S ( τ,h ) ( t , x i , u i , [ u ] i ) = (5) u i − u i τ + sup a ∈ Λ (cid:110) L a [ u ]( t , x i ) + r ( t , x i , a ) u i + (cid:96) ( t , x i , a ) (cid:111) = 0 , and u i = v ( x i ) , i ∈ I (6)is given by the initial condition (2b). Remark
1. As the backward Euler step is only used once, it does not affect theoverall second order of the scheme.
Remark
2. Most of our results also apply to the scheme obtained by replacing theBDF approximation (4) of the drift term by a centred finite difference approximation: (cid:101) D , ± u i := u i +1 − u i − h . (7)However, numerical tests (see Section 7.1) show that the BDF upwind approximationas 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 evenfor possibly vanishing diffusion (Section 5.2). In the remainder of this paper, we provevarious stability and convergence results for the scheme (3). We state in this sectionthe first main well-posedness and stability results.Let u denote the solution of (3) and let v be the solution of (1). The errorassociated with the scheme is then defined by E ki := u ki − v ( t k , x i ) . For any function φ we will also use the notation φ ki := φ ( t k , x i ) as well as φ k :=( φ ki ) ≤ i ≤ I and [ φ ] ki := ( φ mj ) ( j,m ) (cid:54) =( i,k ) , and the error vector at time t k is defined by E k := ( E k , . . . , E kI ) T = u k − v k . The consistency error will be denoted by E k ( φ ) := ( E ki ( φ )) ≤ i ≤ I ∈ R I and isdefined in the classical way as follows, for any smooth enough function φ : E ki ( φ ) := S ( τ,h ) ( t k , x i , φ ki , [ φ ] ki ) − (cid:18) φ t + sup a ∈ Λ (cid:110) L a [ φ ]( t k , x i ) + r ( t k , x i , a ) φ + (cid:96) ( t k , x i , a ) (cid:111)(cid:19) . (8) DF2 SCHEMES FOR HJB EQUATIONS v of (1), we will simply define E ki ( v ) := S ( τ,h ) ( t k , x i , v ki , [ v ] ki ) . (9)Note that (9) is well-defined for any continuous function.In particular for the scheme (3) it is clear that we have second order consistencyin space and time, that is, |E ki ( φ ) | ≤ c ( φ ) τ + c ( φ ) h (10)for sufficiently regular data φ .Throughout the paper, A will denote the finite difference matrix associated to thesecond order derivative, i.e. A := 1 h − − − . . . . . . − . . . . . . . . . . . . − − . (11)Let (cid:104) x, y (cid:105) A := (cid:104) x, Ay (cid:105) . Then we consider the A -norm defined as follows: | x | A := (cid:104) x, Ax (cid:105) = (cid:88) ≤ i ≤ I +1 (cid:18) x i − x i − h (cid:19) (12)(with the convention in (12) that x = x I +1 = 0 ). Hence, √ h | x | A approximates the H semi-norm in Ω . Similarly, we will consider later the standard Euclidean normdefined by (cid:107) x (cid:107) := (cid:104) x, x (cid:105) , such that √ h (cid:107) x (cid:107) approximates the L -norm.Our first result concerns the solvability of the numerical scheme S ( τ,h ) (seen asan equation for u k , with [ u ] k given) and is the following. Assumption (A1). σ, b and r are bounded functions. Theorem Let (A1) and the following CFL condition hold: (cid:107) b (cid:107) ∞ τh < C. (13) Then, for τ small enough and C = 3 / (resp. C = 1 ) there exists a unique solutionof the scheme (3) for k ≥ (resp. k = 1 , for scheme (5) ). The scheme is hence well-defined even if σ vanishes. A uniform ellipticity condi-tion for σ will be needed for proving the H stability of the scheme. Assumption (A2).
There exists η > such that inf t ∈ [0 ,T ] inf x ∈ Ω inf a ∈ Λ σ ( t, x, a ) ≥ η. We provide a relaxation of the ellipticity condition for stability in the Euclidean normin Section 5.2.Our main stability result is the following.
O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Theorem Assume (A1), (A2), as well as the CFL condition (13) . Then thereexists a constant C ≥ (independent of τ and h ) and τ > such that, for any τ ≤ τ , max ≤ k ≤ N | E k | A ≤ C (cid:16) | E | A + | E | A + τ (cid:88) ≤ k ≤ N |E k ( v ) | A (cid:17) . (14)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 reg-ularity assumptions on the boundary data, we can deduce that the scheme (3) is A -norm bounded: max ≤ k ≤ N | u k | A ≤ C, (15)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 ifthe solution to (2) is classical and the supremum is attained for each x and t (andsimilarly for each i and k in (3)), we cannot make any assumptions on the regularityof 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 bythe 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 ≥ such that | σ ( t, x ) − σ ( t, y ) | ≤ L | x − y | ∀ x, y ∈ Ω , t ∈ [0 , T ] . Theorem Assume (A1), (A2), (A3), as well as the CFL condition (13) . Thenthere exists C ≥ (independent of τ and h ) and τ > such that, for any τ ≤ τ , max ≤ k ≤ N (cid:107) E k (cid:107) ≤ C (cid:16) (cid:107) E (cid:107) + (cid:107) E (cid:107) + τ (cid:88) ≤ k ≤ N (cid:107)E k ( v ) (cid:107) (cid:17) . (16)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 ≥ ) can be written in the following form: sup a ∈ Λ ( M ka X − q ka ) = 0 , where q ka ∈ R I and M ka ∈ R I × I with the following non-zero entries: ( M ka ) i,i := 32 + τ (cid:26) σ h + 3 b + h + 3 b − h + r (cid:27) (17) ( M ka ) i,i +1 := τ (cid:26) − σ h − b − h (cid:27) , ( M ka ) i,i − := τ (cid:26) − σ h − b + h (cid:27) (18) ( M ka ) i,i +2 := τ b − h ( M ka ) i,i − := τ b + h (19)with σ ≡ σ ( t k , x i , a ) , b ± ≡ b ± ( t k , x i , a ) and r ≡ r ( t k , x i , a ) . For k = 1 , the termsare different but the form (and analysis) is similar. The fact that ( M a ) i,i ± are non-negative 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: DF2 SCHEMES FOR HJB EQUATIONS Lemma 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: ( i ) for all a ∈ Λ , ( M a ) ii > ii ) (a form of diagonal dominance) sup a ∈ Λ max i ∈ I (cid:80) j>i | ( M a ) ij || ( M a ) ii | − (cid:80) j
8. For a fixed a ∈ Λ , we have max i ∈ I (cid:80) j>i | ( M a ) ij || ( M a ) ii | − (cid:80) j . Moreover, if Λ is compact and a → M a is continuous, then (20) is equivalent to inf a ∈ Λ min i ∈ I (cid:18) | ( M a ) ii | − (cid:88) j (cid:54) = i | ( M a ) ij | (cid:19) > . 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 µ := (cid:88) j>i | ( M a ) ij | ≤ τ (cid:18) σ i h + 5 b − i h (cid:19) (omitting the dependency on k and a in σ, b ± , r ) and µ := | ( M a ) ii | − (cid:88) j
32 + τ (cid:18) σ i h − b + i h + 3 b − i h + r (cid:19) . By the CFL condition (13), there exists (cid:15) > such that τ (cid:107) b (cid:107) ∞ h ≤ − (cid:15) . This implies − (cid:15) τ (cid:18) − b + i h + 3 b − i h (cid:19) ≥ (cid:15) τ b − i h and therefore µ ≥ (cid:18) τ σ i h + (cid:15) τ r (cid:19) + (cid:18) τ b − i h + (cid:15) (cid:19) . Then by using a + a c + c ≤ max (cid:16) a c , a c (cid:17) for numbers a i , c i ≥ , we obtain µ µ ≤ max (cid:18) τ σ i h τ σ i h + (cid:15) + τ r , τ b − i h τ b − i h + (cid:15) (cid:19) . Taking τ small enough such that for instance (cid:15) + τ r ≥ (cid:15) , and since b ( . ) and σ ( . ) arebounded functions (by (A1)), we obtain the bound sup a ∈ A max i ∈ I (cid:80) j>i | ( M a ) ij || ( M a ) ii | − (cid:80) j
4. Proof of Theorem 4 (stability in the A -norm) . The proof consists ofthree 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 (Lemma10); finally, we show the stability estimate from this scalar recursion (Lemma 11).
First, we have the following:
Lemma Let u be the solution of scheme (3) and v the solution of equation (2) .There exist coefficients ˜ σ ki , (˜ b ± ) ki , ˜ r ki , such that the error E k = u k − v k satisfies E ki − E k − i + E k − i τ −
12 (˜ σ ) ki D E ki + (˜ b + ) ki D , − E ki − (˜ b − ) ki D , + E ki + ˜ r ki E ki = −E ki (22) for any k ≥ and i ∈ I , where (˜ σ ) ki , (˜ b ± ) ki , ˜ r ki belong, respectively, to the convexhulls co ( σ ( t k , x i , Λ)) , co ( b ± ( t k , x i , Λ)) , co ( r ( t k , x i , Λ) .Proof. By definition of the consistency error (9), one has (for k ≥ , ≤ i ≤ I ) v ki − v k − i + v k − i τ + H [ v k ]( t k , x i ) = E ki . (23)The scheme simply reads u ki − u k − i + u k − i τ + H [ u k ]( t k , x i ) = 0 . (24)Subtracting (23) from (24), denoting also H [ u k ] ≡ ( H [ u k ]( t k , x i )) ≤ i ≤ I , the followingrecursion is obtained for the error in R I : E k − E k − + E k − τ + H [ u k ] − H [ v k ] = −E k . (25)For simplicity of presentation, we first consider the case when b and r vanish, i.e. b ( . ) ≡ and r ( . ) ≡ . In this case, H [ u k ] i = sup a ∈ Λ (cid:110) − σ ( t k , x i , a )( D u k ) i + (cid:96) ( t k , x i , a ) (cid:111) . (26)To simplify the presentation, we will assume that σ and (cid:96) are continuous functions of a so that the supremum is attained. For each given k, i , let then ¯ a ki ∈ Λ denote anoptimal control in (26).In the same way, let ¯ b ki denote an optimal control for H [ v k ] i . By using the optimalityof ¯ a ki , it holds H [ u k ] i − H [ v k ] i = − σ ( t k , x i , ¯ a ki )( D u k ) i + (cid:96) ( t k , x i , ¯ a ki ) − sup a ∈ Λ (cid:110) − σ ( t k , x i , a )( D v k ) i + (cid:96) ( t k , x i , a ) (cid:111) ≤ − σ ( t k , x i , ¯ a ki )( D u k ) i − (cid:16) − σ ( t k , x i , ¯ a ki )( D v k ) i (cid:17) = − σ ( t k , x i , ¯ a ki )( D E k ) i (27) The general case is obtained easily by considering sequences of (cid:15) -optimal controls and letting (cid:15) → , such that (30) below still holds for a suitably defined ˜ σ , ˜ b + , ˜ b − , ˜ r .DF2 SCHEMES FOR HJB EQUATIONS H [ u k ] i − H [ v k ] i ≥ − σ ( t k , x i , ¯ b ki )( D E k ) i . (28)Therefore, combining (27) and (28), H [ u k ] i − H [ v k ] i is a convex combination of − σ ( t k , x i , ¯ a ki )( D E k ) i and − σ ( t k , x i , ¯ b ki )( D E k ) i . In particular, we can write H [ u k ] i − H [ v k ] i = −
12 ˜ σ ( t k , x i )( D E k ) i , (29)where ˜ σ ( t k , x i ) is a convex combination of σ ( t k , x i , ¯ a ki ) and σ ( t k , x i , ¯ b ki ) .In the general case (i.e. b, r (cid:54)≡ ) one gets similarly H [ u k ] i − H [ v k ] i = −
12 (˜ σ ) ki D E ki + (˜ b + ) ki D , − E ki − (˜ b − ) ki D , + E ki + ˜ r ki E ki , (30)where, for φ = σ , b, r , ˜ φ ki := γ ki φ ( t k , x i , ¯ a ki ) + (1 − γ ki ) φ ( t k , x i , ¯ b ki ) for some γ ki ∈ [0 , . The same technique used above to deal with the non-linear operator applies also to Isaacs equations, i.e. equations of the following form: v t + sup a ∈ Λ inf b ∈ Λ (cid:110) − L ( a,b ) [ v ]( t, x ) + r ( t, x, a, b ) v + (cid:96) ( t, x, a, b ) (cid:111) = 0 , (31)where ( t, x ) ∈ [0 , T ] × R d , Λ , Λ ⊂ R m are compact sets and L ( a,b ) [ v ]( t, x ) = 12 σ ( t, x, a, b ) v xx + b ( t, x, a, b ) v x . To simplify the presentation, let us consider again b, r ≡ , and now also (cid:96) ≡ .By analogous definitions and reasoning to above, we get (25), where, for φ = u, v , H [ φ k ] i = sup a ∈ Λ inf b ∈ Λ (cid:110) − σ ( t, x, a, b )( D x φ k ) i (cid:111) . (32)Let (¯ a ki , ¯ b ki ) ∈ Λ × Λ denote an optimal control in (32). One has H [ u k ] i = sup a ∈ Λ (cid:110) − σ ( t, x, a, ¯ b ki )( D x u k ) i (cid:111) = inf b ∈ Λ (cid:110) − σ ( t, x, ¯ a ki , b )( D x v k ) i (cid:111) . Therefore H [ u k ] i − H [ v k ] i = sup a ∈ Λ (cid:110) − σ ( t, x, a, ¯ b ki )( D x u k ) i (cid:111) − sup a ∈ Λ inf b ∈ Λ (cid:110) − σ ( t k , x i , a, b )( D v k ) i (cid:111) ≥ sup a ∈ Λ (cid:110) − σ ( t, x, a, ¯ b ki )( D x u k ) i (cid:111) − sup a ∈ Λ (cid:110) − σ ( t k , x i , a, ¯ b ki )( D v k ) i (cid:111) ≥ inf a ∈ Λ (cid:110) − σ ( t, x, a, ¯ b ki )( D x E k ) i (cid:111) . (33) Or, if not attained, use an approximation argument. O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Analogously, one can prove H [ u k ] i − H [ v k ] i ≤ sup b ∈ Λ (cid:110) − σ ( t, x, ¯ a ki , b )( D x E k ) i (cid:111) (34)(here, we also use inf( a ) − inf( b ) ≤ sup( a − b ) and sup( a ) − sup( b ) ≥ inf( a − b ) ). Atthis point, it is sufficient to take for ˆ a ki ∈ Λ and ˆ b ki ∈ Λ optimal controls in (33)and (34), respectively, to be able to write H [ u k ] i − H [ v k ] i as a convex combination of − σ ( t, x, ˆ a ki , ¯ b ki )( D x E k ) i and − σ ( t, x, ¯ a ki , ˆ b ki )( D x E k ) i .From this, an equation exactly as in (22) can be derived, with a suitable convexcombination ( ˜ σ ) ki of diffusion coefficients, and similar for the drift and other terms. From (22), we can derive the following:
Lemma
Let assumptions (A1) and (A2) in Theorem 4 be satisfied. Then thereexists a constant C ≥ such that (cid:16) (3 − Cτ ) | E k | A − | E k − | A + | E k − | A (cid:17) + | E k − E k − | A − | E k − − E k − | A ≤ τ | E k | A |E k | A . (35) Proof.
For simplicity of presentation we will assume that b has constant positivesign. The terms coming from the negative part of b can be treated in a similar way.We remark that for E ∈ R I , − D E = AE, where A is the finite difference matrixdefined in (11). By (22), we get the following: E k − E k − + E k − τ + ∆ k AE k + F k BE k + R k E k = −E k , (36)where ∆ k := diag((˜ σ ) ki ) , F k = diag(˜ b ki ) , R k = diag(˜ r ki ) and B = 12 h − − . . . . . . . . . . . . . . . . . . . . . − . We form the scalar product of (36) with AE k . By using the identity (cid:104) a − b, a (cid:105) A = | a | A + | a − b | A − | b | A , one has: (cid:104) E k − E k − + E k − , E k (cid:105) A = 4 (cid:104) E k − E k − , E k (cid:105) A − (cid:104) E k − E k − , E k (cid:105) A = 12 (cid:0) | E k | A + 4 | E k − E k − | A − | E k − | A (cid:1) − (cid:0) | E k | A + | E k − E k − | A − | E k − | A (cid:1) ≥ (cid:0) | E k | A − | E k − | A + | E k − | A (cid:1) + | E k − E k − | A − | E k − − E k − | A , (37)where we have also used | a + b | ≤ | a | + 2 | b | . From ( σ ) ki ≥ η > for all k, i : (cid:104) ∆ k AE k , AE k (cid:105) ≥ η (cid:107) AE k (cid:107) , (38) DF2 SCHEMES FOR HJB EQUATIONS (cid:107) · (cid:107) denotes the canonical Euclidean norm in R I .In order to estimate the drift component, let us introduce the notation δE := ( E i − E i − ) ≤ i ≤ I (39)with the convention that E i = 0 for all indices i which are not in I . It holds: |(cid:104) F k BE k , AE k (cid:105)| = (cid:12)(cid:12)(cid:12)(cid:12) h (cid:104) F k (3 E ki − E ki − + E ki − ) i ∈ I , AE k (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) h (cid:104) F k (cid:0) δE k − δ E k (cid:1) , AE k (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) ≤ h (cid:110) (cid:107) F k δE k (cid:107) (cid:107) AE k (cid:107) + (cid:107) F k δ E k (cid:107) (cid:107) AE k (cid:107) (cid:111) . By using the boundedness of the drift term, and (cid:107) δE k (cid:107) , (cid:107) δ E k (cid:107) ≤ h | E k | A , |(cid:104) F k BE k , AE k (cid:105)| ≤ (cid:107) b (cid:107) ∞ h (cid:110) (cid:107) AE k (cid:107)(cid:107) δE k (cid:107) + (cid:107) AE k (cid:107)(cid:107) δ E k (cid:107) (cid:111) ≤ (cid:107) b (cid:107) ∞ (cid:107) AE k (cid:107) | E k | A . (40)For the last term, using the boundedness of r and the Cauchy-Schwarz inequality, |(cid:104) R k E k , AE k (cid:105)| ≤ (cid:107) r (cid:107) ∞ (cid:107) E k (cid:107)(cid:107) AE k (cid:107) . (41)Therefore, putting (38), (40) and (41) together, (cid:104) ∆ k AE k + F k BE k + R k E k , AE k (cid:105)≥ η (cid:107) AE k (cid:107) − (cid:107) b (cid:107) ∞ (cid:107) AE k (cid:107) | E k | A − (cid:107) r (cid:107) ∞ (cid:107) AE k (cid:107)(cid:107) E k (cid:107) . (42)Easy calculus shows that the minimal eigenvalue of A is λ min ( A ) = h sin ( πh ) ≥ .Hence (cid:104) X, AX (cid:105) ≥ (cid:104) X, X (cid:105) and therefore (cid:107) X (cid:107) ≤ | X | A . In the same way, we havealso | X | A ≤ (cid:107) AX (cid:107) . Hence it holds (cid:104) ∆ k AE k + F k BE k + R k E k , AE k (cid:105) ≥ η (cid:107) AE k (cid:107) − C (cid:107) AE k (cid:107)| E k | A (43)with C := 2 (cid:107) b (cid:107) ∞ + (cid:107) r (cid:107) ∞ . By using C (cid:107) AE k (cid:107)| E k | A ≤ η (cid:107) AE k (cid:107) + η C | E k | A , (cid:104) ∆ k AE k + F k BE k + R k E k , AE k (cid:105) ≥ − η C | E k | A . (44)Then, combining (37) and (44), we obtain the desired inequality with C := η C . In the following, it is assumed that | · | isany vectorial norm. We will use the result for the canonical Euclidean norm | · | ≡ (cid:107) · (cid:107) and the A -norm | · | ≡ | · | A .In order to prove the following Lemma 11, we will exploit properties of the matrix M τ := (3 − Cτ ) − − Cτ ) − . . . . . . . . . . . . . . . . . . −
40 (3 − Cτ ) , (45)in particular the fact that ( M τ ) − ≥ for τ small enough (which we prove).2 O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Lemma
Assume that there exists a constant C ≥ such that ∀ k = 2 , . . . , N : (cid:16) (3 − Cτ ) | E k | − | E k − | + | E k − | (cid:17) + | E k − E k − | − | E k − − E k − | ≤ τ | E k | |E k | . (46) Then there exists a constant C ≥ and τ > such that ∀ < τ ≤ τ , ∀ n ≤ N : max ≤ k ≤ n | E k | ≤ C (cid:16) | E | + | E | + τ (cid:88) ≤ j ≤ n |E j | (cid:17) . (47) Proof.
Let us denote x k := | E k | and y k := | E k − E k − | , so that (46) reads (cid:16) (3 − Cτ ) x k − x k − + x k − (cid:17) ≤ y k − − y k ) + 4 τ | E k | |E k | . (48)For a given τ > and given k , let M τ ∈ R ( k − × ( k − as defined in (45). Let z, w ∈ R k − be defined by z := ( x k , x k − , . . . , x ) T and w := (2( y j − − y j ) + 4 τ | E j | |E j | ) j = k,..., . By (48), we have M τ z ≤ w. (49)We notice that M τ = (3 − Cτ ) I − J + J with J := tridiag(0 , , . Hence, with λ = 2 + √ Cτ and λ = 2 − √ Cτ , the roots of λ − λ + (3 − Cτ ) = 0 for − Cτ ≥ , we can write M τ = ( λ I − J )( λ I − J ) = λ λ (cid:18) I − Jλ (cid:19) (cid:18) I − Jλ (cid:19) . Furthermore, since J k − = 0 , it holds M − τ = 1 λ λ (cid:18) I − Jλ (cid:19) − (cid:18) I − Jλ (cid:19) − = 1 λ λ (cid:88) ≤ q ≤ k − (cid:18) Jλ (cid:19) q (cid:88) ≤ q ≤ k − (cid:18) Jλ (cid:19) q = k − (cid:88) p =0 a p J p , where a p := p (cid:88) j =0 λ j +11 λ p − j +12 = 1 λ p +22 p (cid:88) j =0 (cid:18) λ λ (cid:19) j +1 . Therefore M − τ ≥ componentwise (for τ < /C ), and using (49) it holds z ≤ M − τ w . DF2 SCHEMES FOR HJB EQUATIONS τ > and a constant C ≥ (dependingonly on T ) such that ∀ < τ ≤ τ and ∀ p ≥ : ≤ a p ≤ C and a p − a p − ≥ . (50)We postpone the proof of (50) to the end. For the first component of z , we deduce x k ≤ k − (cid:88) j =0 a j w j +1 ≤ k − (cid:88) j =0 a j ( y k − j − − y k − j ) + 4 C τ k (cid:88) j =2 | E j | |E j | (51) = − a y k + 2 k − (cid:88) j =0 ( a j − a j +1 ) y k − j +1 + 2 a k − y + 4 C τ k (cid:88) j =2 | E j | |E j | , for all k ≥ , where, for (51), we have used the fact that a p ≤ C . Since y j ≥ , ∀ j ,by definition, a k − ≤ C , a = λ λ ≥ and a j − a j − ≥ , ∀ j , we obtain x k ≤ C y + 4 C τ k (cid:88) j =2 | E j | |E j | . (52)Recalling the definition of x k and y k , for any ≤ k ≤ n one has: | E k | ≤ C | E − E | + 4 C τ k (cid:88) j =2 | E j | |E j |≤ C ( | E | + | E | ) + 4 C τ (cid:16) max ≤ k ≤ n | E k | (cid:17) n (cid:88) j =2 |E j |≤ C ( | E | + | E | ) + 12 (cid:16) max ≤ k ≤ n | E k | (cid:17) + 8 C τ (cid:16) n (cid:88) j =2 |E j | (cid:17) (where we made use of ab ≤ a K + Kb for any a, b ≥ and K > ). Hence, we obtain (cid:16) max ≤ k ≤ n | E k | (cid:17) ≤ C (cid:16) | E | + | E | + τ n (cid:88) j =2 |E j | (cid:17) with C := max(8 C , C T ) (we used (cid:16) (cid:80) nj =2 |E j | (cid:17) ≤ n (cid:80) nj =2 |E j | and nτ ≤ T ).It remains to prove (50). From the definition of a p one has a p = 1 λ p +22 p (cid:88) j =0 (cid:18) λ λ (cid:19) j +1 ≤ λ p +22 (cid:18) − λ λ (cid:19) − for p = 0 , . . . , k − . Observing that λ λ ≤ , it follows that a p ≤ λ p +22 ≤ − √ Cτ ) n . O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Notice that √ Cτ ≤ Cτ , and also that e − x ≤ − x/ , ∀ x ∈ [0 , . Hence (2 − √ Cτ ) n ≥ (2 − (1 + Cτ )) n = (1 − Cτ ) n ≥ ( e − Cτ ) n = e − Ct n for Cτ ≤ , andtherefore a p ≤ e Ct n . The desired result follows with C := e CT and τ := C .Moreover, one has a p − a p − = 1 λ p +12 λ p (cid:88) j =0 (cid:18) λ λ (cid:19) j +1 − p − (cid:88) j =0 (cid:18) λ λ (cid:19) j +1 , which is nonnegative for τ small enough thanks to the fact that λ , λ ≥ and λ ≤ .
5. Stability in the Euclidean norm.
The fundamental stability result givenby Lemma 11 applies to any vectorial norm. In this section, we discuss some specialcases where (46) can be obtained for the Euclidean norm | · | = (cid:107) · (cid:107) .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 Lipschitzcontinuous in the following proof. We considerthe scalar product of (36) directly with E k (instead of AE k previously used), againin the situation where b ≥ to simplify the argument. We obtain: (cid:104) E k , E k − E k − + E k − (cid:105) + 2 τ (cid:104) E k , ∆ k AE k + F k BE k + R k E k (cid:105) = − τ (cid:104) E k , E k (cid:105) . (53)As in Section 4.3, we have (cid:104) E k , E k − E k − + E k − (cid:105) (54) ≥ (cid:0) (cid:107) E k (cid:107) − (cid:107) E k − (cid:107) + (cid:107) E k − (cid:107) (cid:1) + (cid:107) E k − E k − (cid:107) − (cid:107) E k − − E k − (cid:107) . We now focus on bounding the other terms on the left-hand side of (53).By using the Lipschitz continuity of σ one has (cid:104) E k , ∆ k AE k (cid:105) = (cid:88) i ∈ I ( σ ki ) h ( − E ki +1 + 2 E ki − E ki − ) E ki = (cid:88) i ∈ I ( σ ki − ) h ( E ki − − E ki ) + (cid:88) i ∈ I (cid:32) ( σ ki − ) h − ( σ ki ) h (cid:33) ( E ki − − E ki ) E ki ≥ η h (cid:88) i ∈ I ( E ki − − E ki ) − L h (cid:88) i ∈ I | E ki − − E ki || E ki | . Therefore, by the Cauchy-Schwarz inequality, one obtains (cid:104) E k , ∆ k AE k (cid:105) ≥ η h (cid:107) δE k (cid:107) − L h (cid:107) δE k (cid:107)(cid:107) E k (cid:107) , (55)where δE k is defined by (39). Moreover, for the first order term one has (cid:104) E k , F k BE k (cid:105) = (cid:88) i ∈ I b i h (3 E ki − E ki − + E ki − ) E ki ≥ − (cid:107) b (cid:107) ∞ h (cid:88) i ∈ I | E ki − E ki − || E ki | − (cid:107) b (cid:107) ∞ h (cid:88) i ∈ I | E ki − − E ki − || E ki |≥ − (cid:107) b (cid:107) ∞ h (cid:107) δE k (cid:107)(cid:107) E k (cid:107) , (56) DF2 SCHEMES FOR HJB EQUATIONS (cid:107) δ E k (cid:107) ≤ (cid:107) δE k (cid:107) . Putting togetherestimates (55) and (56), using the fact that (cid:104) E k , R k E k (cid:105) ≥ −(cid:107) r (cid:107) ∞ (cid:107) E k (cid:107) , we get (cid:104) E k , ∆ k AE k + F k BE k + R k E k (cid:105) ≥ η h (cid:107) δE k (cid:107) − C h (cid:107) δE k (cid:107)(cid:107) E k (cid:107) − (cid:107) r (cid:107) ∞ (cid:107) E k (cid:107) ≥ η h (cid:107) δE k (cid:107) − (cid:18) C η + (cid:107) r (cid:107) ∞ (cid:19) (cid:107) E k (cid:107) , where we have denoted C := L + 4 (cid:107) b (cid:107) ∞ and have used again the Cauchy-Schwarzinequality. Hence, together with (54), this gives (46) with | · | = (cid:107) · (cid:107) and the constant C := 4( C η + (cid:107) r (cid:107) ∞ ) . By using Lemma 11, this concludes the proof of Theorem 6. (cid:50) The next result con-cerns the case of a possibly degenerate diffusion term. It will require more restrictiveassumptions on the drift and diffusion terms, and we shall assume that there is nocontrol here. Indeed, in this case, one cannot count on the positive term coming fromthe non-degenerate diffusion which, in the proof of Theorem 6, is used to compensatethe negative correction terms coming from the drift term. This leads us to considerthe following assumptions:
Assumption (A4). r is bounded. The drift and diffusion coefficients are inde-pendent of the control, i.e. b ≡ b ( t, x ) and σ ≡ σ ( t, x ) , and there exist L , L ≥ suchthat, for all t, x, h : | b ( t, x ) − b ( t, y ) | ≤ L | x − y | , (57) σ ( t, x − h ) − σ ( t, x ) + σ ( t, x + h ) h ≥ − L . (58)(The last condition is equivalent to ( σ ) xx ≥ − L in the differentiable case.) Proposition
Let assumption (A4) be satisfied. Then (46) holds for |·| = (cid:107)·(cid:107) .Proof. We consider again the scalar recursion (53). For any vector E = ( E i ) ≤ i ≤ I (with E j = 0 for j ∈ {− , , I + 1 , I + 2 } ), it holds: E i (2 E i − E i − − E i +1 ) ≥ | E i | −
12 ( | E i | + | E i − | ) −
12 ( | E i | + | E i +1 | ) ≥
12 (2 | E i | − | E i − | − | E i +1 | ) . Hence, by the semi-concavity assumption (58) on σ , (cid:104) E k , ∆ k AE k (cid:105) = (cid:88) ≤ i ≤ I σ i h E ki (2 E ki − E ki − − E ki +1 ) ≥ (cid:88) ≤ i ≤ I σ i h ( −| E ki − | + 2 | E ki | − | E ki +1 | ) ≥ (cid:88) ≤ i ≤ I (cid:18) − σ i − + 2 σ i − σ i +1 h (cid:19) | E ki | . ≥ − L (cid:107) E k (cid:107) . (59)6 O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Now we focus on a lower bound for (cid:104) E k , F k BE k (cid:105) . Let y ki = | E ki − E ki − | . First, (3 E ki − E ki − + E ki − ) E ki = 12 (3 | E ki | − | E ki − | + | E ki − | )+ 12 (4 | E ki − E ki − | − | E ki − E ki − | ) ≥
12 (3 | E ki | − | E ki − | + | E ki − | ) + 12 (2 y ki − y ki − ) . We assume again b i ≥ for all i to simplify the presentation. The case where b i ≤ for some i is similar. Then, the following bound holds: (cid:104) E k , F k BE k (cid:105) = I (cid:88) i =1 b i h (3 E ki − E ki − + E ki − ) E ki = I +2 (cid:88) i =1 b i h (3 E ki − E ki − + E ki − ) E ki ≥ I +2 (cid:88) i =1 b i h (3 | E ki | − | E ki − | + | E ki − | ) + I +2 (cid:88) i =1 b i h ( y ki − y ki − ) ≥ I (cid:88) i =1 (cid:18) b i − b i +1 + b i +2 h (cid:19) | E ki | + I +1 (cid:88) i =1 (cid:18) b i − b i +1 h (cid:19) y ki (where we have used y k = y kI +2 = 0 and (cid:80) ≤ i ≤ I +2 b i ( E ki − ) = (cid:80) ≤ i ≤ I b i +2 ( E ki ) aswell as (cid:80) ≤ i ≤ I +2 b i ( E ki − ) = (cid:80) ≤ i ≤ I +1 b i +1 ( E ki ) = (cid:80) ≤ i ≤ I b i +1 ( E ki ) ). Then, bythe Lipschitz continuity of b ( . ) and the bound y ki ≤ E ki ) + 2( E ki − ) , we have (cid:104) E k , F k BE k (cid:105) ≥ − L I (cid:88) i =1 | E ki | − L I +1 (cid:88) i =1 y ki ≥ − L (cid:107) E k (cid:107) . (60)By combining the bounds (59) and (60), we obtain (cid:104) E k , ∆ k AE k (cid:105) + (cid:104) E k , F k BE k (cid:105) + (cid:104) E k , R k E k (cid:105) ≥ − ( L L + (cid:107) r (cid:107) ∞ ) (cid:107) E k (cid:107) . Therefore, inequality (46) is obtained with C := 4( L + 3 L + (cid:107) r (cid:107) ∞ ) , which leads tothe desired stability estimate. Under suitable assumptions, theresult of Theorem 6 can be extended to multi-dimensional equations. The nonlinearitycan be treated exactly as in Section 4.1 (or 4.2), so that we can focus on the linearcase v t −
12 tr[Σ( t, x ) D x v ] + b ( t, x ) D x v + r ( t, x ) v + (cid:96) ( t, x ) = 0 for a positive definite matrix Σ and a drift vector b . For simplicity, we furthermoreconsider the two-dimensional case d = 2 , with r, (cid:96) ≡ , and omit the dependence ofthe coefficients on the time variable, then with Σ( x, y ) := (cid:18) σ ( x, y ) ρσ σ ( x, y ) ρσ σ ( x, y ) σ ( x, y ) (cid:19) and b ( x, y ) := (cid:18) b ( x, y ) b ( x, y ) (cid:19) , where σ , σ ≥ and ρ ∈ [ − , is the correlation parameter, the equation reads v t − σ ( x, y ) v xx − ρσ σ ( x, y ) v xy − σ ( x, y ) v yy + b ( x, y ) v x + b ( x, y ) v y = 0 . DF2 SCHEMES FOR HJB EQUATIONS
Ω := ( x min , x max ) × ( y min , y max ) . We introducethe discretization in space defined by the steps h x , h y > and we denote by G ( h x ,h y ) the associated mesh. In what follows, given any function φ of ( x, y ) ∈ Ω , we willdenote φ ij = φ ( x i , y j ) for ( i, j ) ∈ I := I × I , where I = { , . . . , I } , I = { , . . . , I } .Assuming that ρ ≥ everywhere (the case when ρ ≤ is similar), we consider a7-point stencil for the second order derivatives (see [12, Section 5.1.4]): v xx ∼ v i − ,j − v ij + v i +1 ,j h x =: D xx v ij , v yy ∼ v i,j − − v ij + v i,j +1 h y =: D yy v ij v xy ∼ − v i,j − − v i,j +1 − v i − ,j − v i +1 ,j + v i − ,j − + v i +1 ,j +1 + 2 v ij h x h y =: D xy v ij and the BDF approximation of the first order derivatives D , − x u ij := 3 u ij − u i − ,j + u i − ,j h x and D , + x u ij := − (cid:18) u ij − u i +1 ,j + u i +2 ,j h x (cid:19) ,D , − y u ij := 3 u ij − u i,j − + u i,j − h y and D , + y u ij := − (cid:18) u ij − u i,j +1 + u i,j +2 h y (cid:19) . The scheme is therefore defined, for k ≥ , by u kij − u k − ij + u k − ij τ (61) − σ ( x i , y j ) D xx u kij − ρσ σ ( x i , y j ) D xy u kij − σ ( x i , y j ) D yy u kij + b +1 ( x i , y j ) D , − x u kij − b − ( x i , y j ) D , + x u kij + b +2 ( x i , y j ) D , − y u kij − b − ( x i , y j ) D , + y u kij . A straightforward calculation shows that the second order term also reads σ ( x i , y j ) D xx u ij + 2 ρσ σ ( x i , y j ) D xy u ij + σ ( x i , y j ) D yy u ij = α ij D xx u ij + β ij D yy u ij + γ ij ( u i − ,j − − u ij + u i +1 ,j +1 ) , (62)with α ij := σ ( x i , y j ) h x (cid:18) σ ( x i , y j ) h x − ρσ ( x i , y j ) h y (cid:19) ,β ij := σ ( x i , y j ) h y (cid:18) σ ( x i , y j ) h y − ρσ ( x i , y j ) h x (cid:19) , γ ij := ρ ( x i , y j ) σ ( x i , y j ) σ ( x i , y j ) h y h x . The scheme is completed with the following boundary conditions: u ki,j = v ( t k , x i , y j ) , ∀ i ∈ {− , } ∪ { I + 1 , I + 2 } , j ∈ I ,u ki,j = v ( t k , x i , y j ) , ∀ j ∈ {− , } ∪ { I + 1 , I + 2 } , i ∈ I . For simplicity, assume h x = h y =: h . We consider the following assumptions: Assumptions(A1’): (cid:107) b i (cid:107) ∞ < ∞ for i = 1 , ; (A2’): ∃ η > , ∀ ( x, y ) ∈ Ω , ∀ i (cid:54) = j : σ i ( x, y ) − ρ ( x, y ) σ i ( x, y ) σ j ( x, y ) ≥ η ; (A3’): ∀ i, j = 1 , , σ i σ j is Lipschitz continuous on Ω .We then have the following result. The proof is similar to the one of Theorem 6,using (62) with α ij , β ij , γ ij ≥ by assumption (A2’), and is therefore omitted.8 O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Proposition
Let assumptions (A1’),(A2’) and (A3’) be satisfied. Then thestability estimate (47) holds for | · | = (cid:107) · (cid:107) .Remark ( i ) If h x (cid:54) = h y and for instance h y = Ch x for some C ≥ , (A2’) hasto hold with σ replaced by σ /C as a result of the scaling properties of the scheme. ( ii ) Observe that assumption (A2’) is equivalent to requiring strong diagonaldominance of the covariance matrix. ( iii ) When the strong diagonal dominance of the matrix Σ is not guaranteed, onecan consider the generalized finite difference scheme in [7]. However, determining theprecise set of assumptions on the coefficients needed to apply the previous argumentsdoes not seem easy from the construction in [7].
6. Error estimates.
In this section, we give detailed error estimates for theimplicit BDF2 scheme (3). We consider the following rescaled norms on R I : | u | := (cid:32)(cid:88) i ∈ I u i h (cid:33) / = (cid:107) u (cid:107)√ h, | u | := (cid:32)(cid:88) i ∈ I (cid:18) u i − u i − h (cid:19) h (cid:33) / = | u | A √ h, corresponding to discrete approximations of L (Ω) - and H (Ω) 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 ) : | w | C ,α ( I ) := sup (cid:26) | w ( x ) − w ( y ) || x − y | α , x (cid:54) = y, x, y ∈ I (cid:27) . For a given open subset Ω ∗ T of (0 , T ) × Ω , we define C k,(cid:96) (Ω ∗ T ) as the set of functions v : Ω ∗ T → R which admit continuous derivatives ( ∂ i v∂t i ) ≤ i ≤ k and ( ∂ j v∂x j ) ≤ j ≤ (cid:96) on Ω ∗ T .We also denote by C k,(cid:96)b (Ω ∗ T ) the subset of functions with bounded derivatives on Ω ∗ T . Assumption (A5). v ∈ C , ((0 , T ) × Ω) and for some constant C ≥ : sup x ∈ Ω (cid:107) v t ( ., x ) (cid:107) C ,δ ([0 ,T ]) ≤ C, sup t ∈ (0 ,T ) (cid:107) v xx ( t, . ) (cid:107) C ,δ (¯Ω) ≤ C. (63) Remark
15. By results in [10] and [14], assumption (A5) is satisfied for sufficientlysmooth data and given a uniform ellipticity condition.We have the following error estimates:
Theorem
We assume (A1), (A2), (A3), and the CFL condition (13) . ( i ) If v ∈ C , b ((0 , T ) × Ω) , then max ≤ k ≤ N | v k − u k | ≤ Ch , where C is a constant which depends on the derivatives of v of order 3 and 4in t and x , respectively. ( ii ) If (A5) holds for some δ ∈ (0 , , then the numerical solution u of (3), (5)converges to v in the L -norm with max ≤ k ≤ N | v k − u k | ≤ Ch δ , for some constant C (possibly different from the one in (A5)). DF2 SCHEMES FOR HJB EQUATIONS Proof.
We first prove ( ii ) . By Taylor expansion, we can write for instance, forsome θ , θ ∈ [0 , , (cid:12)(cid:12)(cid:12)(cid:12) v t ( t, x ) − v ( t, x ) − v ( t − τ, x ) τ (cid:12)(cid:12)(cid:12)(cid:12) = | v t ( t, x ) − v t ( t − θ τ, x ) | ≤ Cτ δ and (cid:12)(cid:12)(cid:12)(cid:12) v t ( t, x ) − v ( t, x ) − v ( t − τ, x ) + v ( t − τ, x )2 τ (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12) v t ( t, x ) −
12 (3 v t ( t − θ τ, x ) − v t ( t − (1 + θ ) τ, x )) (cid:12)(cid:12)(cid:12)(cid:12) ≤ | v t ( t, x ) − v t ( t − θ τ, x ) | + 12 | v t ( t − θ τ, x ) − v t ( t − (1 + θ ) τ, x ) |≤ Cτ δ + 12 C (2 τ ) δ ≤ Cτ δ . Similarly, using the higher spatial regularity, there exists a constant C ≥ such that (cid:12)(cid:12)(cid:12)(cid:12) v x ( t, x ) − v ( t, x ) − v ( t, x − h ) + v ( t, x − h )2 h (cid:12)(cid:12)(cid:12)(cid:12) ≤ C Ch δ +1 , (cid:12)(cid:12)(cid:12)(cid:12) v xx ( t, x ) − v ( t, x + h ) − v ( t, x ) + v ( t, x − h ) h (cid:12)(cid:12)(cid:12)(cid:12) ≤ C Ch δ . The result ( ii ) now follows directly by inserting the obtained truncation error into thestability estimate of Theorem 6.For the proof of ( i ) (smooth case), expansion up to order 3 and 4 gives thetruncation error of higher order for k ≥ , and we use the fact that the error fromthe first backward Euler step is bounded by (cid:107) E (cid:107) ≤ Cτ ( τ + h ) ; in particular, we usethat ( E − E ) /τ + (∆ A + F B + R ) E = −E , with (cid:107)E (cid:107) ≤ C ( τ + h ) , E = 0 and the bound is otherwise similar and simpler than that for k ≥ .The previous arguments can also be used to derive error estimates for piecewisesmooth solutions. In this case, we will need to limit the number of non-regular pointsthat may appear in the exact solution (assumption (A6) ( i ) is similar to [5]). Assumption (A6).
There exists an integer p ≥ and functions ( x ∗ j ( t )) ≤ j ≤ p for t ∈ [0 , T ] , such that, with Ω ∗ T := (Ω × (0 , T )) \ (cid:83) ≤ j ≤ p { ( t, x ∗ j ( t )) , t ∈ (0 , T ) } , thefollowing holds: ( i ) v ∈ C , b (Ω ∗ T ) ; ( ii ) ∀ j , t → x ∗ j ( t ) is Lipschitz regular.We give the following straightforward preliminary result without proof: Lemma
Assume (A6) and the CFL condition (13) . Then for all t Card { j, x → v ( t, x ) not regular in [ x j − , x j +2 ] } ≤ p and Card { j, θ → v ( θ, x j ) not regular in [ t − τ, t ] } ≤ Cp. for some constant C ≥ independent of τ, h ("not regular" meaning not C in thefirst case and not C in the second one). Such a situation will be illustrated in the numerical example of Section 7.2.0
O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Theorem
We assume (A1), (A2), (A3) and the CFL condition (13) . Let(A5) and (A6) hold for some δ ∈ (0 , , then the numerical solution u of (3), (5)converges to v in the L -norm with max ≤ k ≤ N | v k − u k | ≤ Ch / δ , where C is a constant independent of h .Proof. Let I k be the (finite) set of indices i such that v is not regular in { t k } × ( x i − h, x i + 2 h ) ∪ ( t k − τ, t k ) × { x i } . Then |E k | = (cid:88) i ∈ I |E ki | h = (cid:88) i ∈ I k |E ki | h + (cid:88) i ∈ I \ I k |E ki | h ≤ C | I k | ( τ δ + h δ ) h + C ( τ + h ) . We then use the fact that | I k | ≤ C for some (different) constant C by Lemma 17and that ( τ + h ) = O ( h ) = O ( h δ ) , τ δ + h δ = O ( h δ ) by the CFL condition (13),in order to obtain the desired result. Remark
19. (i) Similar results can be derived for errors in the A -norm, howeverderivatives of one order higher are required due to the derivative in the definitionof the norm.(ii) The estimates in Theorem 16 are not always sharp, as symmetries and thesmoothing behaviour of the scheme can result in higher order convergence. Wediscuss such special cases for Examples 1 and 2 in Section 7, Remarks 21 and22, respectively.(iii) These error estimates can be compared with [5], where an error bound of order h / was obtained for diffusion problems with an obstacle term, under the mainassumption that v xx is a.e. bounded with a finite number of singularities (insteadof (A5)) . In the present context it seems natural to assume the Hölder regularityof u t and u xx coming from the ellipticity assumption (see Remark 15).
7. Numerical tests.
We now compare the performance of the BDF2 schemewith other second order finite difference schemes on two examples.
The first example is based on a deterministiccontrol problem ( σ ≡ ) and motivates the choice of the BDF2 approximation for thedrift term in (4), compared to the more classical centered scheme (7). We consider (cid:26) v t + | v x | = 0 , x ∈ ( − , , t ∈ (0 , T ) ,v (0 , x ) = v ( x ) , x ∈ ( − , , with v ( x ) = max(0 , − x ) and T = 0 . . The initial datum is shown in Figure 1(dashed line). The exact solution is v ( t, x ) = min( v ( x − t ) , v ( x + t )) . Remark
20. The Eikonal equation can be written as v t +max a ∈{− , } ( av x ) = 0 inHJB 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. DF2 SCHEMES FOR HJB EQUATIONS Fig. 1 . Test 1: Initial data (dashed line) and numerical solution at time T = 0 . computedfor I + 1 = 200 and N = 20 ( τ/h = 0 . ) using BDF in time and centred approximation of the drift(left), BDF in time and space (right). In Figure 1, we show the results obtained at the terminal time T = 0 . usingschemes (3)-(7) (left) and (3)-(4) (right) with τ /h = 0 . . We numerically observethat the centered approximation generates undesirable oscillations, whereas the BDF2scheme is stable.As stated in Theorem 3, in case of a degenerate diffusion, a CFL condition ofthe form τ ≤ Ch has to be satisfied for well-posedness of the BDF2 scheme. Table 1shows numerical convergence of order in both time and space, although the solutionis globally only Lipschitz. N I + 1 H norm L -norm L ∞ norm CPU (s)error order error order error order5 10 5.35E-01 - 1.25E-01 - 1.36E-01 - 0.09410 20 2.42E-01 1.14 4.51E-02 1.47 6.83E-02 0.99 0.09620 40 8.25E-02 1.55 1.55E-02 1.55 2.01E-02 1.77 0.12640 80 2.38E-02 1.80 4.32E-03 1.84 5.23E-03 1.94 0.14780 160 6.26E-03 1.92 1.11E-03 1.96 1.31E-03 2.00 0.194160 320 1.61E-03 1.96 2.79E-04 1.99 3.24E-04 2.01 0.335320 640 4.09E-04 1.98 7.10E-05 1.99 8.19E-05 2.00 0.759640 1280 1.03E-04 1.99 1.78E-05 2.00 2.05E-05 2.00 2.306 Table 1
Test 1. Error and convergence rate to the exact solution for the BDF2 scheme with τ/h = 0 . and initial data v ( x ) = max(0 , − x ) . Remark
21. The full convergence order here is due to the particular symmetry ofthe solution. To confirm this, we report in Table 2 the results obtained for the sameequation with initial data v (0 , x ) = − max(0 , − x ) (see also Figure 2). In this case, there is no such symmetry around the two singularpoints and as a result the full convergence order is lost: the scheme is globally onlyof order in the H norm and roughly . in the L and L ∞ norm.2 O. BOKANOWSKI, A. PICARELLI, C. REISINGER
Fig. 2 . Test 1: Initial data (dashed line) v ( x ) = − max(0 , − x ) and numerical solu-tion at time T = 0 . computed for I + 1 = 200 and N = 20 ( τ/h = 0 . ) using the BDF2scheme. The convergence rates for this exam-ple are reported in Table 2. N I + 1 H norm L norm L ∞ norm CPU (s)error order error order error order5 10 5.84E-01 - 1.62E-01 - 1.51E-01 - 0.00610 20 2.69E-01 1.12 5.23E-02 1.63 6.20E-02 1.28 0.00820 40 1.45E-01 0.89 1.86E-02 1.49 2.08E-02 1.58 0.01840 80 6.74E-02 1.10 5.95E-03 1.64 7.89E-03 1.40 0.03980 160 3.20E-02 1.08 1.81E-03 1.72 3.57E-03 1.15 0.093160 320 1.60E-02 1.00 5.44E-04 1.73 1.51E-03 1.24 0.233320 640 8.16E-03 0.97 1.65E-04 1.72 6.33E-04 1.25 0.695640 1280 4.20E-03 0.96 5.09E-05 1.70 2.64E-04 1.26 2.163 Table 2
Test 1. Error and convergence rate to the exact solution for the BDF2 scheme with τ/h = 0 . and initial data v ( x ) = − max(0 , − x ) . The secondtest we propose is a problem with controlled diffusion. We consider (cid:40) v t + sup σ ∈{ σ ,σ } (cid:16) − σ v xx (cid:17) = 0 , x ∈ ( − , , t ∈ (0 , T ) ,v (0 , x ) = sin( πx ) , x ∈ ( − , , with parameters σ = 0 . , σ = 0 . , T = 0 . .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 similaroptimal control problem. The BDF2 scheme, in contrast, has shown good performancefor that same problem in [6].Figure 3 (top row) shows the initial data and the value function at terminal timecomputed using the BDF2 scheme. The error and convergence rate in different normsare reported in Table 3. Here an accurate numerical solution computed by an implicitEuler 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 someinstability 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 , the CN scheme shows second order of convergence. Remark
22. In this example, due to the strict ellipticity, Assumption (A5) isguaranteed for some δ > (see Remark 15). Then Theorem 16 gives convergence DF2 SCHEMES FOR HJB EQUATIONS Fig. 3 . Test 2: Initial data (top, left), numerical solution at time T = 0 . (top, right) computedby the BDF2 scheme, second order derivative computed with CN scheme (bottom, left) and BDF2(bottom, right) for N = 256 and I + 1 = 5120 . N I + 1 H norm L norm L ∞ norm CPU (s)error order error order error order1 20 1.54E-01 - 5.11E-02 - 7.24E-02 - 0.1312 40 5.53E-02 1.48 1.88E-02 1.45 2.63E-02 1.46 0.1124 80 1.47E-02 1.91 5.17E-03 1.86 6.99E-03 1.91 0.1118 160 3.59E-03 2.04 1.27E-03 2.03 1.66E-03 2.08 0.12216 320 8.98E-04 2.00 3.14E-04 2.02 4.09E-04 2.02 0.14632 640 2.26E-04 1.99 7.84E-05 2.00 1.02E-04 2.00 0.18364 1280 5.65E-05 2.00 1.96E-05 2.00 2.56E-05 2.00 0.267128 2560 1.42E-05 2.00 4.90E-06 2.00 6.42E-06 2.00 0.598256 5120 1.21E-06 2.01 1.21E-06 2.01 1.59E-06 2.01 1.879 Table 3
Test 2. Error and convergence rate for the BDF2 scheme with high CFL number τ = 5 h . Areference solution computed by the implicit Euler scheme (6) with I + 1 = 20 × , N = 2 is used. 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, weconjecture that Assumption (A6) is satisfied, such that Theorem 18 would give thehigher order / δ . In the test, in fact the full order 2 is observed (see Table 3).
8. Conclusion.
We have proved the well-posedness and stability in L and H norms of a second order BDF scheme for HJB equations with enough regularity ofthe coefficients. The significance of the results is that this was achieved for a second4 O. BOKANOWSKI, A. PICARELLI, C. REISINGER
N I + 1 H norm L norm L ∞ norm CPU (s)error order error order error order1 20 4.11E-02 - 7.01E-03 - 9.44E-03 - 0.1492 40 7.82E-03 2.39 1.45E-03 2.27 2.29E-03 2.04 0.1134 80 1.97E-03 1.99 3.87E-04 1.91 5.62E-04 2.03 0.1118 160 5.16E-04 1.94 1.02E-04 1.92 1.45E-04 1.95 0.12816 320 1.09E-04 2.24 2.67E-05 1.94 3.77E-05 1.95 0.16632 640 2.96E-05 1.88 7.15E-06 1.90 9.87E-06 1.93 0.18864 1280 7.64E-06 1.96 2.03E-06 1.82 2.61E-06 1.92 0.310128 2560 9.50E-05 -3.64 1.98E-05 -3.29 3.49E-05 -3.74 0.992256 5120 7.18E-04 -2.92 8.40E-05 -2.08 1.62E-04 -2.22 4.251 Table 4
Test 2. Error and convergence rate for the CN scheme with high CFL number τ = 5 h . Areference solution computed by the implicit Euler scheme (6) with I + 1 = 20 × , N = 2 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 thenumerical solution in terms of the truncation error of the scheme. The latter dependson 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 thetheoretical results, and often better, due to symmetries of the solution or smoothingproperties of the equation and the scheme. This is in contrast to some alternativesecond order schemes, such as the central spatial difference in the case of a first orderequation, 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 anduniqueness of a solution to (21), we consider a fixed-point approach. The initialproblem (21) can be written as follows: sup a ∈ Λ ( L a X − ( q a − U a X )) = 0 , (64)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 i ≥ j , and U a the remaining upper triangular part, ( U a ) ij := ( M a ) ij i Math. Comput. , 74(260):1861–1893, 2007.[2] G. Barles and P.E. Souganidis. Convergence of approximation schemes for fully nonlinearsecond order equations. Asymptotic Anal. , 4:271–283, 1991.[3] T. Beale. Smoothing properties of implicit finite difference methods for a diffusion equation inmaximum norm. SIAM J. Numer. Anal. , 47(4):2476–2495, 2009.[4] J. Becker. A second order backward difference method with variable steps for a parabolicproblem. BIT Numer. Math. , 38(4):644–662, 1998.[5] O. Bokanowski and K. Debrabant. High order finite difference schemes for some nonlineardiffusion equations with an obstacle term. HAL preprint hal-01686742.[6] O. Bokanowski, A. Picarelli, and C. Reisinger. High-order filtered schemes for time-dependentsecond order HJB equations. ESAIM Math. Model. Numer. Anal. , 2017. Forthcoming.[7] J.F. Bonnans and H. Zidani. Consistency of generalized finite difference schemes for the stochas-tic HJB equation. SIAM J. Numer. Anal. , 41(3):1008–1021, 2003.[8] M.G. Crandall, H. Ishii, and P.L. Lions. User’s guide to viscosity solutions of second orderpartial differential equations. Bull. Amer. Math. Soc. , 27(1):1–67, 1992.[9] E. Emmrich. Stability and error of the variable two-step BDF for semilinear parabolic problems. J. Appl. Math. & Computing , 19(1-2):33–55, 2005.[10] L.C. Evans and S. Lenhart. The parabolic Bellman equation. Nonlinear Anal. , 5(7):765–773,1981.[11] S.K. Godunov. A difference method for numerical calculation of discontinuous solutions of theequations of hydrodynamics. Matematicheskii Sbornik , 89(3):271–306, 1959.[12] W. Hackbusch. Elliptic Differential Equations: Theory and Numerical Treatment , volume 18of Springer Series in Computational Mathematics . Springer, 2010.[13] A. Hill and E. Süli. Approximation of the global attractor for the incompressible Navier–Stokesequations. IMA J. Numer. Anal. , 20(4):633–667, 2000.[14] N.V. Krylov. Boundedly nonhomogeneous elliptic and parabolic equations. Izvestiya RossiiskoiAkademii Nauk. Seriya Matematicheskaya , 46(3):487–523, 1982.[15] A. Picarelli, C. Reisinger, and J. Rotaetxe. Error bounds for monotone schemes forparabolic Hamilton-Jacobi-Bellman equations in bounded domains. arXiv preprintarXiv:1710.11284 , 2017.[16] D.M. Pooley, P.A. Forsyth, and K.R. Vetzal. Numerical convergence properties of option pricingpdes with uncertain volatility. IMA J. Numer. Anal. , 23(2):241–267, 2003.[17] I. Smears and E. Süli. Discontinuous Galerkin finite element methods for time-dependentHamilton–Jacobi–Bellman equations with Cordes coefficients. Numer. Math. , 133(1):141–176, 2016.[18] J. Stoer and R. Bulirsch. Introduction to Numerical Analysis , volume 12 of Texts in Ap-plied Mathematics . Springer-Verlag, New York, second edition, 1993. Translated from theGerman by R. Bartels, W. Gautschi and C. Witzgall.[19] E. Süli and D.F. Mayers.