1. Introduction
In this work we deal with the numerical approximation of the following parabolic HamiltonJacobiBellman (HJB) equation
(1.1) 
In the system above, , is a nonempty smooth bounded open set and and are nonlinear functions having the form
(1.2)  
(1.3) 
where and are nonempty compact sets, , with , , , , with being an open set containing , , and .
If and , for some and , and , with
being the unit outward normal vector to
at , then (1.1) reduces to a standard linear parabolic equation with Neumann boundary conditions. In the general case, and after a simple change of the time variable in order to write (1.1) in backward form, the HJB equation (1.1) appears in the study of optimal control of diffusion processes with controlled reflection on the boundary (see e.g. [27] for the first order case, i.e. , and [26, 11] for the general case). Since the HJB equation (1.1) is possibly degenerate parabolic, one cannot expect the existence of classical solutions and we have to rely on the notion of viscosity solution (see e.g. [16]). Moreover, as it has been noticed in [25, 27], in general the boundary condition in (1.1) does not hold in the pointwise sense and we have to consider a suitable weak formulation of it. We refer the reader to [27, 4] and [16, 6, 7, 24, 12], respectively, for wellposedness results for HJB equations with oblique boundary condition in the first and second order cases.The study of the numerical approximation of solutions to HJB and, more generally, fully nonlinear second order Partial Differential Equations (PDEs), has made important progress over the last few decades. Most of the related literature consider the case where
, or where a Dirichlet boundary condition is imposed on the boundary . We refer the reader to [19, 20, 30] and the references therein for the state of the art on this topic. By contrast, the numerical approximation of solutions to (1.1) has been much less explored. Indeed, to the best of our knowledge only the methods in [31, 1] can be applied to approximate (1.1) in the particular first order case (). Moreover, in [31], where a finite difference scheme is proposed, the function defining the boundary condition has the particular form . On the other hand, both references consider Hamiltonians which are not necessarily convex with respect to . Let us also mention the reference [2], where, in the context of mean curvature motion with nonlinear Neumann boundary conditions, the authors propose a discretization that combines a SemiLagrangian (SL) scheme in the main part of the domain with a finite difference scheme near the boundary.The main purpose of this article is to provide a consistent, stable, monotone and convergent SL scheme to approximate the unique viscosity solution to (1.1). By the results in [6], the latter is wellposed in under the assumptions in Sect. 2 below. SemiLagrangian schemes to approximate the solution to (1.1) when (see e.g. [13, 17]) can be derived from the optimal control interpretation of (1.1) and a suitable discretization of the underlying controlled trajectories. These schemes enjoy the feature that they are explicit and stable under an inverse CourantFriedrichsLewy (CFL) condition and, consequentely, they allow large time steps. A second important feature is that they permit a simple treatement of the possibly degenerate second order term in . The scheme that we propose for preserves these two properties and seems to be the first convergent scheme to approximate (1.1) with the rather general asumptions in Sect. 2. In particular, our results cover the stochastic and degenerate case. Consequently, from the stochastic control point of view, our scheme allows to approximate the socalled value function of the optimal control of a controlled diffusion process with possibly oblique reflection on the boundary (see [11]). The main difficulty in devising such a scheme is to be able to obtain a consistency type property at points in the space grid which are near the boundary while maintaining the stability. This is achieved by considering a discretization of the underlying controlled diffusion which suitably emulates its reflection at the boundary in the continuous case. We refer the reader to [28] for a related construction of a semidiscrete in time approximation of a second order nondegenerate linear parabolic equation.
The remainder of this paper is structured as follows. In Sect. 2 we state our assumptions, we recall the notion of viscosity solution to (1.1) and the wellposedness result. In Sect. 3 we provide the SL scheme as well as its probabilistic interpretation (in the spirit of [28]). The latter will play an important role in Sect. 4, which is devoted to show a consistency type property and the stability of the SL scheme. By using the halfrelaxed limits technique introduced in [5], we show in Sect. 5 our main result, which is the convergence of solutions to the SL scheme towards the unique viscosity solution to (1.1). The convergence is uniform in and holds under the same asymptotic condition between the space and time steps than in the case . Next, in Sect. 6 we first illustrate the numerical convergence of the SL scheme in the case of a onedimensional linear equation with homogeneous Neumann boundary conditions. In this case the numerical results confirm that the boundary condition in (1.1) is not satisfied at every , but it is satisfied in the viscosity sense recalled in Sect. 2 below. In a second example, we consider a two dimensional degenerate second order nonlinear equation on a circular domain with nonhomogeneous Neumann and oblique boundary conditions. In the last example, we consider a twodimensional nondegenerate nonlinear equation on a nonsmooth domain. Due to the lack of regularity of , our convergence result does not apply. However, the SL scheme can be successfully applied, which suggests that our theoretical findings could hold for more general domains. This extension as well as the corresponding study in the stationary framework remain as interesting subjects of future research. Finally, we provide in the Appendix of this work some theoretical results concerning oblique projections and the regularity of the distance to , which play a key role in the definition of the scheme and in the proof of its main properties.
Acknowledgements
The first two authors would like to thank the Italian Ministry of Instruction, University and Research (MIUR) for supporting this research with funds coming from the PRIN Project (KKJPX entitled “Innovative numerical methods for evolutionary partial differential equations and applications”). Xavier Dupuis thanks the support by the EIPHI Graduate School (contract ANR17EURE0002). Elisa Calzola, Elisabetta Carlini and Francisco J. Silva were partially supported by KAUST through the subaward agreement OSRCRG..
2. Preliminaries
As mentioned in the introduction, it will be simpler to describe our approximation scheme when (1.1) is written in backward form. This can be done by a simple change of the time variable and a possible modification of the time dependency of . Let us set and . We consider the HJB equation
(HJB) 
For notational convenience, throughout this article, we will write for all and . Our standing assumptions for the data in (HJB) are the following.

() is a nonemtpy, bounded domain with boundary of class .

The functions , , , and are continuous. Moreover, for every , the functions and are Lipschitz continuous, with Lipschitz constants independent of .

The function is of class . We also assume that
where, for every , we recall that denotes the unit outward normal vector to at .
We now recall the notion of viscosity solution to (see [6]). We need first to introduce some notation. Given a bounded function , its upper semicontinuous (resp. lower semicontinuous) envelope is defined by
(2.1) 
Definition 2.1.
[Viscosity solution]
(i) An upper semicontinuous function is a viscosity subsolution to (HJB) if for any and such that has a local maximum at , we have
(2.2) 
if ,
(2.3) 
if and,
(2.4) 
if .
(ii) A lower semicontinuous function is a viscosity supersolution to (HJB) if for any and such that has a local minimum at , we have
(2.5) 
if ,
(2.6) 
if and,
(2.7) 
if .
(iii) A bounded function is a viscosity solution to (HJB) if and , defined in (2.1), are, respectively, sub and supersolutions to (HJB).
Remark 2.1.
Theorem 2.1.
Assume (H1)(H3). Then there exists a unique viscosity solution to (HJB).
Remark 2.2.
(i) [Comparison principle and uniqueness] The existence of at most one solution to (HJB) follows from the following comparison principle (see [6, Theorem II.1] and also [11, Proposition 3.4]). If is a bounded viscosity subsolution to (HJB) and is a bounded viscosity supersolution to (HJB), then
(ii) [Existence] Once a comparison principle has been shown, the existence of a solution to (HJB) follows usually from the existence of sub and supersolutions to (HJB) and Perron’s method. In Sect. 5, we construct sub and supersolutions to (HJB) as suitable limits of solutions to the approximation scheme that we present in the next section. Together with the comparison principle, this yields an alternative existence proof of solutions to (HJB).
An different and interesting technique to show the existence of a solution to (HJB) is to consider a suitable stochastic optimal control problem, with controlled reflection of the state trajectory at the boundary , and to show that the associated value function is a viscosity solution to (HJB). This strategy has been followed in [11].
(iii) [Continuity] The continuity of the unique viscosity solution to (HJB) follows directly from the comparison principle and the continuity properties required in the definition of sub and supersolutions to (HJB). Notice that, as usual for parabolic problems with Neumann type boundary conditions, we do not require any compatibility condition between and the operator at the boundary .
3. The fully discrete scheme
We introduce in this section a fully discrete SL scheme that approximates the unique viscosity solution to . Throughout this section, we assume that (H1)(H3) are fulfilled.
3.1. Discretization of the space domain
Let us fix and consider a polyhedral domain such that
(3.1) 
for some . A construction of such a domain can be found in [8, Section 3] for or , which explain the dimension constraint in (H1). However, the results in the remainder of this article can be extended to , provided that a numerical domain satisfying (3.1) exists. Let be a triangulation of consisting of simplicial finite elements with vertices in (for some ). We assume that is the mesh size, i.e. the maximum of the diameters of , all the vertices on belong to , at most one face of each element , with vertices on , intersects , and satisfies the following regularity condition: there exists , independent of , such that each is contained in a ball of radius and contains a ball of radius . As in [18], we introduce an auxiliary exact triangulation of with vertices in . The boundary elements of are allowed to be curved and we have
Denoting by the projection on , the projection is defined by
if  
and the element has the same vertices than . 
Set and denote by the linear finite element basis function on . More precisely, for each , is a continuous function, affine on each , , , for all , with , and for all . For any
its linear interpolation
on the mesh is defined by(3.2) 
Lemma 3.1.
Let and denote by its restriction to . Then there exists a constant , independent of , such that
(3.3) 
Proof.
Let and let and be two elements having the same vertices and such that . By the triangular inequality
Using that is Lipschitz, we deduce from (3.1) the existence of , independent of and , such that
. In addition, by standard error estimates for
interpolation (see for instance [15]) and (3.2), there exists , independent of and , such that . Relation (3.3) follows from these two estimates. ∎3.2. A semiLagrangian scheme
Let , set , and . We define the time grid . Given , , and , we define the discrete characteristics
(3.4) 
Let and let be a fixed constant. For any we set
By Proposition 7.1 in the Appendix, there exist and two functions and , uniquely determined, such that
(3.5) 
Therefore, there exists such that for all , , , , and , the reflected characteristic
(3.6) 
is welldefined. In Figure 1 we illustrate how the reflected characteristic is computed from the projection of onto parallel to . Let us also set
(3.7)  
(3.8) 
Notice that if , then (3.5), (3.6), and (3.7) imply that
(3.9) 
For and , let us define by
(3.10) 
and set
(3.11) 
In the remainder of this work, we will consider the following fully discrete SL scheme to approximate the solution to .
(HJB) 
3.3. Probabilistic interpretation of the scheme
The fullydiscrete SL to approximate the solution to (HJB) in the unbounded case, i.e. , has a natural interpretation in terms of a discrete time, finite state, Markov control process (see e.g. [13, Section 3]). We show below that a similar interpretation holds for (HJB). The latter will play an important role in the stability analysis of (HJB) presented in the next section. Given and , , let us define the controlled transition law
(3.12) 
We say that is a policy if for all we have . The set of policies is denoted by . Let us fix and, for notational convenience, set . Associated to and
, there exists a probability measure
on (the powerset of) and a Markov chain
, with state space , such that(3.13) 
for . Now, consider a family of
valued independent random variables, which are also independent of
, and with common distribution given bywhere denotes the th canonical vector of . By a slight abuse of notation (see (3.4)), for , , and , let us set
(3.14) 
For , , , and , define the random variable
(3.15) 
For all and , let us define
where, for notational convenience, we have denoted, respectively, by and the first and the last coordinates of . Notice that, by construction and (3.10), we have that
Moreover, setting
for all , the dynamic programming principle (see e.g. [23, Theorem 12.1.5]) implies that satisfies (HJB). Since the latter has a unique solution, we deduce that for all and .
4. Properties of the fully discrete scheme
In this section, we establish some basic properties of (HJB).
Proposition 4.1.
The following hold:
(i) (Monotonicity) For all with , we have
(ii) (Commutation by constant) For any and ,
We show in Proposition 4.2 below a consistency result for (HJB). For this purpose, let us set
and for all , , , , , and , define
(4.3) 
Proposition 4.2 (Consistency).
Let and denote by its restriction to . Then the following hold:

For all , , , and , we have
where the set of constants is bounded, independently of .

For all and , we have
Comments
There are no comments yet.