# Finite Element Approximation of Hamilton-Jacobi-Bellman equations with
nonlinear mixed boundary conditions

Bartosz Jaroszkowski and Max Jensen
Department of Mathematics, University of Sussex, Brighton, UK Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We show strong uniform convergence of monotone P1 finite element methods to
the viscosity solution of isotropic parabolic Hamilton-Jacobi-Bellman
equations with mixed boundary conditions on unstructured meshes and for
possibly degenerate diffusions. Boundary operators can generally be
discontinuous across face-boundaries and type changes. Robin-type boundary
conditions are discretised via a lower Dini derivative. In time the Bellman
equation is approximated through IMEX schemes. Existence and uniqueness of
numerical solutions follows through Howard’s algorithm. Finite element method;
Hamilton-Jacobi-Bellman equation; Mixed boundary conditions; Fully nonlinear
equation; Viscosity solution

##  1 Introduction

The value function of an optimal control problem is, under suitable
assumptions, the solution of a Hamilton-Jacobi-Bellman (HJB) equation. In this
work we consider the numerical solution of HJB equations with mixed boundary
conditions of the form:

$\displaystyle-\partial_{t}v+\sup_{\alpha\in A}(L^{\alpha}\;\;v-f^{\alpha})$
$\displaystyle=0$ $\displaystyle\qquad\text{in }[0,T)\times{\Omega},$ (1a)
$\displaystyle-\partial_{t}v+\sup_{\alpha\in
A}\bigl{(}L^{\alpha}_{\partial{\Omega}}v-g^{\alpha}\bigr{)}$ $\displaystyle=0$
$\displaystyle\qquad\text{on }[0,T)\times\partial{\Omega}_{t},$ (1b)
$\displaystyle\sup_{\alpha\in
A}\bigl{(}L^{\alpha}_{\partial{\Omega}}v-g^{\alpha}\bigr{)}$ $\displaystyle=0$
$\displaystyle\qquad\text{on }[0,T)\times\partial{\Omega}_{R},$ (1c)
$\displaystyle v-g\;\;\;$ $\displaystyle=0$ $\displaystyle\qquad\text{on
}[0,T)\times\partial{\Omega}_{D},$ (1d) $\displaystyle v-v_{T}\,$
$\displaystyle=0$ $\displaystyle\qquad\text{on
}\\{T\\}\times\overline{{\Omega}}.$ (1e)

Here $L^{\alpha}$ and $L^{\alpha}_{\partial{\Omega}}$ denote operators on the
domain $\Omega$ and its boundary, respectively. The sets
$\partial{\Omega}_{t}$, $\partial{\Omega}_{R}$ and $\partial{\Omega}_{D}$ form
a decomposition of $\partial\Omega$. While leaving further details of the
notation to the next section, it is already apparent how the basic fully
nonlinear structure of the PDE operator, meaning the left-hand side of (1a),
is mirrored in the Robin-type boundary conditions (1b) and (1c). But there is
a crucial, additional complication of the boundary operators
$L^{\alpha}_{\partial{\Omega}}$: They will in general depend on the full
gradient $\nabla v$ and not just on the tangential gradient
$\nabla_{\\!\partial{\Omega}}v$, meaning that $L^{\alpha}_{\partial{\Omega}}v$
cannot be evaluated with knowledge of $v|_{\partial{\Omega}}$ only.

Recalling the connection between optimal control and HJB equations, Bellman-
type equations as in (1b) naturally arise on sections $\partial{\Omega}_{t}$
of the boundary. Indeed, the boundary condition (1b) expresses the possibility
of controlling the particle or agent on the boundary through processes which
are implicitly described by the $L^{\alpha}_{\partial{\Omega}}$. In contrast,
Dirichlet conditions (1d) are appropriate for those sections
$\partial{\Omega}_{D}$ of $\partial{\Omega}$ where the possibility to control
may cease in exchange for the reward or cost of $g_{D}$.

Boundary conditions of type (1c) arise from the Skorokhod control problem,
which models particle reflection at the boundary [Lio85, KD01, Ser03].
Moreover, they have recently been used for the numerical solution of optimal
transport problems in the setting of Monge-Ampère equations, where the
transport boundary conditions are examined in Hamilton-Jacobi form [BFO12,
Kaw19].

For the authors the problem of primary interest is the Heston model of
financial interest rates with an uncertain market price of volatility risk
[JJ21]. The Heston equation is most naturally posed on an unbounded domain,
where already with a certain market price of volatility risk it appears with
mixed boundary terms corresponding to (1b), (1c) as well as (1d). All those
types of boundary conditions remain when introducing uncertainty and when
truncating the domain for the purposes of numerical approximation.

The aim of this work is to introduce a finite element method capable of
computing approximations to viscosity solutions for the aforementioned
problems. This paper extends results of [JS13] with the inclusion of mixed,
fully nonlinear boundary conditions. The presented method permits degenerate
diffusions. Boundary operators may exhibit discontinuities across face
boundaries and where the type of boundary condition changes.

A challenge for problems of this type is the discretisation of the first-order
directional derivatives in (1b) and (1c) which needs to be simultaneously
consistent and monotone. On the one hand establishing monotonicity with an
artificial diffusion approximating the Laplace-Beltrami operator of
$\partial{\Omega}$ would not be sufficient because of the normal component in
the directional derivatives of (1b) and (1c). On the other hand an artificial
diffusion approximating the Laplace operator of $\Omega$ would not vanish
under refinement due to different scaling of boundary and domain terms, thus
leading to an inconsistent method. Our formulation is based on the observation
that lower Dini directional derivatives exist for all functions in the P1
approximation space whenever the direction in question is in the tangent cone
of $\Omega$ at the position of interest.

A benefit of the finite element approach is that besides $L^{\infty}$ also
$L^{2}(H^{1})$ convergence can be established on unstructured meshes, as was
shown in [JS13, Jen17] for the Dirichlet problem. The $L^{2}(H^{1})$
convergence is for instance important for the above mentioned Heston model
[JJ21] as Delta hedging requires knowledge of partial derivatives of the value
function.

The numerical analysis of HJB equations with Neumann and Robin conditions
encompasses only few works. First results were provided by the finite
difference community; we refer to the text book [KD01]. More recently, the
transport boundary conditions of optimal transport were in [BFO12]
approximated with a filtered wide stencil scheme. In [AF12] a nonlinear
Neumann boundary operator is approximated by extending the boundary into a
strip of positive thickness, allowing the boundary conditions to be treated
like a PDE operator. Within the finite element setting one line of research
has developed around the approximation of Cordes solutions, in [Gal19] with a
mixed, non-conforming finite element method while in [Kaw19] with a
discontinuous Galerkin finite element method. Both [Gal19] and [Kaw19]
concentrate on the linear setting in non-divergence form. For a general review
of the approximation of fully nonlinear equations with other types of boundary
conditions we refer to [FGN13, NSZ17].

The structure of this article is as follows. In Section 2 we formulate the HJB
problem with mixed boundary conditions. In Section 3 we define the numerical
method. In Section 4 we prove monotonicity properties of the discretised
operators. In Section 5 we show the existence and uniqueness of numerical
solutions. In Sections 6 and 7 we establish consistency and stability,
respectively, leading us to our main result of convergence in Section 8.
Finally, we present numerical experiments in Section 9.

##  2 Mixed final time boundary value Bellman problem

In this section we introduce time-dependent Bellman equations with mixed
boundary conditions. We consider a polytopic domain
${\Omega}\subset\mathbb{R}^{d}$ with $d\geq 2$, i.e. a bounded, connected,
closed domain, whose interior is non-empty and whose boundary is formed of
flat faces. We allow $\Omega$ to be non-convex. Let $\mathcal{F}_{k}$ denote
set of open $k$-dimensional faces of $\Omega$ contained in $\partial{\Omega}$.

We consider Dirichlet and Robin boundary conditions on disjoint subsets of
boundary $\partial{\Omega}$. Additionally, the region of the Robin boundary
conditions breaks into two parts, one with and one without time derivative. We
denote those three disjoint regions as $\partial{\Omega}_{D}$,
$\partial{\Omega}_{R}$ and $\partial{\Omega}_{t}$, respectively. Therefore
$\partial{\Omega}_{D}\cap\partial{\Omega}_{R}\cap\partial{\Omega}_{t}=\emptyset$
and
$\partial{\Omega}_{D}\cup\partial{\Omega}_{R}\cup\partial{\Omega}_{t}=\partial{\Omega}$.

It is convenient to define the notion of a generalised face as the
intersection of an $\omega^{\prime}\in\mathcal{F}_{d-1}$ and a region linked
to a boundary condition:

$\mathcal{F}:=\bigl{\\{}\omega\subset\partial{\Omega}:\,\omega=\omega^{\prime}\cap\partial{\Omega}_{X}\text{
where
}\omega^{\prime}\in\mathcal{F}_{d-1},\partial{\Omega}_{X}\in\\{\partial{\Omega}_{D},\partial{\Omega}_{R},\partial{\Omega}_{t}\\}\bigr{\\}}$

We assume that the boundary conditions are continuous on each $\omega$;
however, discontinuities across (generalised) face boundaries may occur.

We introduce the normed space of piecewise continuous functions

$PC(\partial{\Omega},\mathbb{R}^{k}):=\\{g\in
L^{\infty}(\partial{\Omega},\mathbb{R}^{k}):\;g|_{\omega}\in
C(\omega,\mathbb{R}^{k})\quad\forall\,\omega\in\mathcal{F}\\},$

equipped with the $L^{\infty}(\partial{\Omega},\mathbb{R}^{k})$ norm. If
$k=1$, we simply write $PC(\partial{\Omega})$. We denote the standard inner
product of $L^{2}({\Omega})$ and $L^{2}({\Omega},\mathbb{R}^{d})$ by
$\langle\cdot,\cdot\rangle$.

Let $A$ be a compact metric space, $\alpha\in A$ and let $L^{\alpha}$ be a
linear operators of the following form:

$L^{\alpha}:\;\overline{{\Omega}}\times\mathbb{R}\times\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R},\;(x,q,p,s)\mapsto-a^{\alpha}(x)\,q-b^{\alpha}(x)\cdot
p+c^{\alpha}(x)\,s.$

The interpretation as differential operator follows with $q=\Delta w(x)$,
$p=\nabla w(x)$ and $s=w(x)$ for $w\in C^{2}(\overline{{\Omega}})$. The
mapping

$A\to C(\overline{{\Omega}})\times C(\overline{{\Omega}},\mathbb{R}^{d})\times
C(\overline{{\Omega}})\times C(\overline{{\Omega}}),\\\
\alpha\mapsto(a^{\alpha},b^{\alpha},c^{\alpha},f^{\alpha}),$

is assumed to be continuous such that the families of functions
$\\{a^{\alpha}\\}_{\alpha\in A}$, $\\{b^{\alpha}\\}_{\alpha\in A}$,
$\\{c^{\alpha}\\}_{\alpha\in A}$ and $\\{f^{\alpha}\\}_{\alpha\in A}$ are
equicontinuous. We require that $a^{\alpha}(x)\geq 0$ for all $\alpha\in A$ so
that all $L^{\alpha}$ are degenerate elliptic. Frequently we abbreviate
$L^{\alpha}(x,\Delta w(x),\nabla w(x),w(x))$ by $L^{\alpha}w(x)$ and $x\mapsto
L^{\alpha}w(x)$ by $L^{\alpha}w$.

For $\alpha\in A$ the Robin operators $L^{\alpha}_{\partial{\Omega}}$ are
defined as

$\displaystyle
L^{\alpha}_{\partial{\Omega}}:\partial{\Omega}\times\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R},\;(x,p,s)\mapsto-
b_{\partial{\Omega}}^{\alpha}(x)\cdot p+c_{\partial{\Omega}}^{\alpha}(x)s$ (2)

with $b_{\partial{\Omega}}^{\alpha}\in PC(\partial{\Omega},\mathbb{R}^{d})$,
$c_{\partial{\Omega}}^{\alpha}\in PC(\partial{\Omega})$. The abbreviations
$L^{\alpha}_{\partial{\Omega}}w(x)$ and $L^{\alpha}_{\partial{\Omega}}w$ are
used analogously to $L^{\alpha}$.

We can now pose Hamilton-Jacobi-Bellman (HJB) problem, whose numerical
solution is the subject of this paper:

$\displaystyle-\partial_{t}v+\sup_{\alpha\in A}(L^{\alpha}\;\;v-f^{\alpha})$
$\displaystyle=0$ $\displaystyle\qquad\text{in }[0,T)\times{\Omega},$ (3a)
$\displaystyle-\partial_{t}v+\sup_{\alpha\in
A}\bigl{(}L^{\alpha}_{\partial{\Omega}}v-g^{\alpha}\bigr{)}$ $\displaystyle=0$
$\displaystyle\qquad\text{on }[0,T)\times\partial{\Omega}_{t},$ (3b)
$\displaystyle\sup_{\alpha\in
A}\bigl{(}L^{\alpha}_{\partial{\Omega}}v-g^{\alpha}\bigr{)}$ $\displaystyle=0$
$\displaystyle\qquad\text{on }[0,T)\times\partial{\Omega}_{R},$ (3c)
$\displaystyle v-g\;\;\;$ $\displaystyle=0$ $\displaystyle\qquad\text{on
}[0,T)\times\partial{\Omega}_{D},$ (3d) $\displaystyle v-v_{T}\,$
$\displaystyle=0$ $\displaystyle\qquad\text{on
}\\{T\\}\times\overline{{\Omega}},$ (3e)

with $g\in C(\partial{\Omega})$, $g^{\alpha}\in PC(\partial{\Omega})$,
$v_{T}\in C(\overline{{\Omega}})$ and $T\in(0,\infty)$. The suprema are
applied pointwise. An interpretation of (3) in the context of optimal control
is given in Appendix A. Observe that the data terms $g^{\alpha}$ of the Robin
conditions are $\alpha$-dependent, while the corresponding Dirichlet data $g$
are not. We assume equicontinuity of mapping

$A\to PC(\partial{\Omega},\mathbb{R}^{d})\times PC(\partial{\Omega})\times
PC(\partial{\Omega}),\;\alpha\mapsto(b_{\partial{\Omega}}^{\alpha},c_{\partial{\Omega}}^{\alpha},g^{\alpha})$

in $\alpha\in A$. Additionally, we require the sign-conditions
$v_{T},\,c_{{\Omega}}^{\alpha},\,c_{\partial{\Omega}}^{\alpha},\,g,\,g^{\alpha}\geq
0$. It follows from continuity that

$\displaystyle\sup_{\alpha\in
A}\|\,(b_{\partial{\Omega}}^{\alpha},c_{\partial{\Omega}}^{\alpha},g^{\alpha})\,\|_{L^{\infty}(\partial{\Omega},\mathbb{R}^{d})\times
L^{\infty}(\partial{\Omega})\times L^{\infty}(\partial{\Omega})}<\infty.$ (4)

We require $v_{T}$ to satisfy the Dirichlet boundary conditions on
$\partial{\Omega}_{D}$.

It is useful to formulate the operator used in (3) more succinctly as

$\displaystyle F(t,x,q,p,r,s)=$
$\displaystyle\left\\{\begin{array}[]{rllll}-r+\sup_{\alpha}\bigl{(}L^{\alpha}(x,q,p,s)-f^{\alpha}(x)\bigr{)}&\text{if
}(t,x)\in[0,T)\times\overline{{\Omega}},\\\
-r+\sup_{\alpha}\bigl{(}L^{\alpha}_{\partial{\Omega}}(x,p,s)-g^{\alpha}(x)\bigr{)}&\text{if
}(t,x)\in[0,T)\times\partial{\Omega}_{t},\\\
\phantom{-r+\;}\sup_{\alpha}\bigl{(}L^{\alpha}_{\partial{\Omega}}(x,p,s)-g^{\alpha}(x)\bigr{)}&\text{if
}(t,x)\in[0,T)\times\partial{\Omega}_{R},\\\ s-g(x)&\text{if
}(t,x)\in[0,T)\times\partial{\Omega}_{D},\\\ s-v_{T}(x)&\text{if
}(t,x)\in\\{T\\}\times\overline{{\Omega}}.\end{array}\right.$

We conclude the section with a definition of a viscosity solution similar to
the setting of [BS91] which will be used throughout the chapter. To this end,
let us consider a bounded function
$v:[0,T]\times\overline{{\Omega}}\to\mathbb{R}$ and its upper and lower semi-
continuous envelopes, defined respectively as

$v^{*}(t,x):=\limsup_{\begin{subarray}{c}(s,y)\to(t,x)\\\
(s,y)\in[0,T]\times\overline{{\Omega}}\end{subarray}}v(s,y)$

and

$v_{*}(t,x):=\liminf_{\begin{subarray}{c}(s,y)\to(t,x)\\\
(s,y)\in[0,T]\times\overline{{\Omega}}\end{subarray}}v(s,y).$

We analogously extend the definition of lower- and upper semicontinuous
envelopes to $F$.

######  Definition 1.

A bounded function $v$ is a viscosity supersolution (respectively,
subsolution) of (3) if, for any test function $\psi\in
C^{2}(\mathbb{R}\times\mathbb{R}^{d})$,

$F^{*}(t,x,\Delta\psi(t,x),\nabla\psi(t,x),\partial_{t}\psi(t,x),v_{*}(t,x))\geq
0,$

(respectively,

$\left.F_{*}(t,x,\Delta\psi(t,x),\nabla\psi(t,x),\partial_{t}\psi(t,x),v^{*}(t,x))\leq
0,\right)$

provided that $v_{*}-\psi$ attains a local minimum (respectively, $v^{*}-\psi$
attains a local maximum) at $(t,x)\in[0,T]\times\overline{{\Omega}}$. Finally,
we call $v:[0,T]\times\overline{{\Omega}}\to\mathbb{R}$ a viscosity solution
of (3) if it is simultaneously a viscosity sub- and supersolution of (3).

Here the value of

$F^{*}(t,x,\Delta\psi(t,x),\nabla\psi(t,x),\partial_{t}\psi(t,x),v_{*}(t,x))$

should only depend on $\psi$’s restriction to
$[0,T]\times\overline{{\Omega}}$. However, generally for lower-dimensional
faces $F\in\mathcal{F}_{d-2}$ one finds test functions $\psi,\phi\in
C^{2}(\mathbb{R}\times\mathbb{R}^{d})$ with
$\psi|_{[0,T]\times\overline{{\Omega}}}=\phi|_{[0,T]\times\overline{{\Omega}}}$
such that $\nabla\psi(t,x)\neq\nabla\phi(t,x)$ for $x\in F$.

We therefore demand that the coefficient $b_{\partial{\Omega}}^{\alpha}(x)$ of
(2) belongs to the tangent cone:

$\displaystyle b_{\partial{\Omega}}^{\alpha}(x)\in
K(x)\qquad\forall\,\alpha\in A,x\in\partial{\Omega}.$ (5)

Here, because of the polytopic nature of the domain, we define the tangent
cone $K(x)$ as

$\displaystyle
K(x):=\bigl{\\{}x^{\prime}\in\mathbb{R}^{d}\big{|}\;\exists\,\Lambda\in(0,\infty)\;\forall\,\lambda\in[0,\Lambda]:x+\lambda\,x^{\prime}\in\overline{{\Omega}}\bigr{\\}},$

i.e. $x^{\prime}$ is in the cone if there is a line segment from $x$ in the
direction of $x^{\prime}$ which is contained in $\overline{{\Omega}}$. Indeed,
for $b_{\partial{\Omega}}^{\alpha}(x)\in K(x)\setminus\\{0\\}$ we observe how

$\displaystyle-
b_{\partial{\Omega}}^{\alpha}(x)\cdot\nabla\psi(t,x)=\partial_{-b_{\partial{\Omega}}^{\alpha}(x)}\psi(t,x)=\lim_{\begin{subarray}{c}\lambda\to
0\\\\[0.85358pt]
\lambda>0\end{subarray}}\frac{\psi(t,x)-\psi(t,x+\lambda\,b_{\partial{\Omega}}^{\alpha}(x))}{\lambda}$
(6)

is expressed only referring to $\psi$ on $[0,T]\times\overline{{\Omega}}$ and
thus independently of $\psi$’s extension to $\mathbb{R}\times\mathbb{R}^{d}$.
The limit on the right-hand side of (6) is known as the lower Dini derivative
of $\psi$ in direction $-b_{\partial{\Omega}}^{\alpha}(x)$. On smooth sections
of the boundary and at outward pointing corners the requirement (5)
corresponds to an outflow condition, while at re-entrant corners (5) may
permit an inflow term. In that sense, (5) is less restrictive than oblique
boundary conditions such as [CIL92, (7.35)]. We remark, however, that
strengthened versions such as [CIL92, (7.35)] may be necessary to ensure the
existence of a comparison principle for the final time boundary value problem
of the specific application of interest.

##  3 Numerical scheme

For the discretisation of (3) we consider a sequence $V_{i},i\in\mathbb{N}$,
of piecewise linear, simplicial, shape-regular finite element spaces. Let
$\mathcal{T}_{i}$ be the mesh corresponding to the finite element space
$V_{i}$. The boundary mesh $\mathcal{B}_{i}$ consists of the
$(d-1)$-dimensional faces $F$ of elements $K\in\mathcal{T}_{i}$ with
$F\subset\partial{\Omega}$. We make the assumption that $\mathcal{B}_{i}$ is
subordinate to $\mathcal{F}$, i.e. every open $F\in\mathcal{B}_{i}$ is
contained entirely in exactly one generalised face $\omega\in\mathcal{F}$.

Let $V_{i}^{g}\subset V_{i}$ be the affine subspace of functions which
interpolate the Dirichlet boundary data on $\partial{\Omega}_{D}$ and
$V_{i}^{0}\subset V_{i}$ be the vector subspace of functions which interpolate
$0$ on $\partial{\Omega}_{D}$. The nodes of the finite element mesh are
denoted by ${y}_{i}^{\ell}$. Here the index $\ell$ ranges over the nodes in
the interior first, then the nodes on $\partial{\Omega}_{t}$, then
$\partial{\Omega}_{R}$ and finally $\partial{\Omega}_{D}$. Therefore,
${y}_{i}^{\ell}\in{\Omega}$ for $\ell\leq N_{i}^{{\Omega}}$ for some
$N_{i}^{{\Omega}}\in\mathbb{N}$ denoting number of interior nodes,
${y}_{i}^{\ell}\in{\Omega}\cup\partial{\Omega}_{t}$ for $\ell\leq N_{i}^{t}$
for some $N_{i}^{t}\in\mathbb{N}$ and, lastly,
${y}_{i}^{\ell}\in{\Omega}\cup\partial{\Omega}_{t}\cup\partial{\Omega}_{R}$
for $\ell\leq N_{i}:=\dim V_{i}^{g}$. These nodes
${y}_{i}^{\ell}\in{\Omega}\cup\partial{\Omega}_{t}\cup\partial{\Omega}_{R}$
are called non-Dirichlet nodes.

The associated hat functions ${\phi}_{i}^{\ell}\in V_{i}$ are chosen so that
${\phi}_{i}^{\ell}({y}_{i}^{\ell})=1$ while ${\phi}_{i}^{\ell}(y_{i}^{s})=0$
for $\ell\neq s$. Set
$\hat{\phi}_{i}^{\ell}:={\phi}_{i}^{\ell}/\|{\phi}_{i}^{\ell}\|_{L^{1}({\Omega})}$.
Therefore, the ${\phi}_{i}^{\ell}$ are normalised in the
$L^{\infty}(\overline{{\Omega}})$ norm whilst the $\hat{\phi}_{i}^{\ell}$ are
normalised in the $L^{1}(\overline{{\Omega}})$ norm.

The mesh size, i.e. the largest diameter of an element, is denoted $\Delta
x_{i}$. It is assumed that $\Delta x_{i}\to 0$ as $i\to\infty$. The uniform
time step size is denoted $h_{i}$ with the constraint that
$T/h_{i}\in\mathbb{N}$. It is assumed that $h_{i}\to 0$ as $i\to\infty$. Let
${s}_{i}^{k}$ be the $k$th time step at the refinement level $i$. Then the set
of time steps is $S_{i}:=\bigl{\\{}{s}_{i}^{k}:k=T/h_{i},\dots,0\bigr{\\}}$.

We introduce the operator $d_{i}$, which approximates the time derivative on
$\Omega$ and $\partial{\Omega}_{t}$ but which is $0$ on the remaining
boundary. More precisely, we let the $\ell$th entry of
$d_{i}w({s}_{i}^{k},\cdot)$ be

$(d_{i}w({s}_{i}^{k},\cdot))_{\ell}=\begin{cases}\frac{w({s}_{i}^{k+1},{y}_{i}^{\ell})-w({s}_{i}^{k},{y}_{i}^{\ell})}{h_{i}}\quad&\ell\leq
N_{i}^{t},\\\ \hskip 52.03227pt0&\textrm{otherwise.}\end{cases}$

Observe how $(d_{i}w({s}_{i}^{k},\cdot))_{\ell}=0$ for nodes
${y}_{i}^{\ell}\in\partial{\Omega}_{R}$ is consistent with the structure of
(3c).

For each discretisation of (3) we allow a splitting of $L^{\alpha}$ and
$L^{\alpha}_{\partial{\Omega}}$ into an explicit and an implicit part. For
each $\alpha$ and for each $i$, we introduce the explicit operator
${E}_{{\Omega},i}^{\alpha}$ and the implicit operator
${I}_{{\Omega},i}^{\alpha}$ such that

$\displaystyle{E}_{{\Omega},i}^{\alpha}$
$\displaystyle:\;C^{2}(\overline{{\Omega}})\to
C(\overline{{\Omega}}),\;w\mapsto-\bar{a}_{{\Omega},i}^{\alpha}\,\Delta
w\,-\bar{b}_{{\Omega},i}^{\alpha}\cdot\,\nabla
w+\bar{c}_{{\Omega},i}^{\alpha}\,w,$ $\displaystyle{I}_{{\Omega},i}^{\alpha}$
$\displaystyle:\;C^{2}(\overline{{\Omega}})\to
C(\overline{{\Omega}}),\;w\mapsto-\bar{\bar{a}}_{{\Omega},i}^{\alpha}\,\Delta
w\,-\bar{\bar{b}}_{{\Omega},i}^{\alpha}\cdot\,\nabla
w+\bar{\bar{c}}_{{\Omega},i}^{\alpha}\,w,$

where
$\bar{a}_{{\Omega},i}^{\alpha},\bar{\bar{a}}_{{\Omega},i}^{\alpha},\bar{c}_{{\Omega},i}^{\alpha},\bar{\bar{c}}_{{\Omega},i}^{\alpha}\in
C(\overline{{\Omega}})$ and
$\bar{b}_{{\Omega},i}^{\alpha},\bar{\bar{b}}_{{\Omega},i}^{\alpha}\in
C(\overline{{\Omega}},\mathbb{R}^{d})$. Assumption 1 below shows that the
explicit and implicit operators are chosen such that
${I}_{{\Omega},i}^{\alpha}+{E}_{{\Omega},i}^{\alpha}$ approximates
$L^{\alpha}$. Analogously, we introduce non-negative $f_{i}^{\alpha}\in
C(\overline{{\Omega}})$ which approximate $f^{\alpha}$.

The discretisations of ${E}_{{\Omega},i}^{\alpha}$ and
${I}_{{\Omega},i}^{\alpha}$ are the mappings from $V_{i}$ to
$\mathbb{R}^{N_{i}}$ which are given by

$\displaystyle({\mathsf{E}}_{{\Omega},i}^{\alpha}w)_{\ell}:=$
$\displaystyle\,\bar{a}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell})\langle\nabla
w,\nabla\hat{\phi}_{i}^{\ell}\rangle+\langle-\bar{b}_{{\Omega},i}^{\alpha}\cdot\nabla
w+\bar{c}_{{\Omega},i}^{\alpha}\,w,\hat{\phi}_{i}^{\ell}\rangle,$ (7a)
$\displaystyle({\mathsf{I}}_{{\Omega},i}^{\alpha}w)_{\ell}:=$
$\displaystyle\,\bar{\bar{a}}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell})\langle\nabla
w,\nabla\hat{\phi}_{i}^{\ell}\rangle+\langle-\bar{\bar{b}}_{{\Omega},i}^{\alpha}\cdot\nabla
w+\bar{\bar{c}}_{{\Omega},i}^{\alpha}w,\hat{\phi}_{i}^{\ell}\rangle,$ (7b)
$\displaystyle({\mathsf{F}}_{{\Omega},i}^{\alpha})_{\ell}:=$
$\displaystyle\,\langle{f}_{i}^{\alpha},\hat{\phi}_{i}^{\ell}\rangle,$ (7c)

where $\ell$ ranges over all internal nodes, i.e. $\ell\leq N^{\Omega}_{i}$.
For boundary nodes $\ell>N^{\Omega}_{i}$ we set
$({\mathsf{E}}_{{\Omega},i}^{\alpha}w)_{\ell}=({\mathsf{I}}_{{\Omega},i}^{\alpha}w)_{\ell}=({\mathsf{F}}_{{\Omega},i}^{\alpha})_{\ell}=0$.
Because of the scaling of the $\hat{\phi}_{i}^{\ell}$, integration-by-parts
gives for smooth $w$ and large $i$ that $\langle\nabla
w,\nabla\hat{\phi}_{i}^{\ell}\rangle\approx-\Delta w(x)$ if
$x\approx{y}_{i}^{\ell}$ away from the boundary.

Similarly, we define operators ${\mathsf{E}}^{\alpha}_{\partial{\Omega},i}$
and ${\mathsf{I}}^{\alpha}_{\partial{\Omega},i}$ on the boundary to discretise
$L^{\alpha}_{\partial{\Omega}}$ as the sum of an explicit and implicit part.
Starting point is the observation that the directional derivative
$\partial_{-b_{\partial{\Omega}}^{\alpha}}w$ is well-defined in the sense of
(6) for functions $w\in V_{i}$ even though $w$ is in general not
differentiable.

For interior nodes with index $0\leq\ell\leq N_{i}^{\Omega}$ we set
$({\mathsf{E}}^{\alpha}_{\partial{\Omega},i}w)_{\ell}=({\mathsf{I}}^{\alpha}_{\partial{\Omega},i}w)_{\ell}=({\mathsf{F}}^{\alpha}_{\partial{\Omega},i})_{\ell}=0$.
More interestingly, for $N_{i}^{\Omega}<\ell\leq N_{i}$ ranging over the nodes
of the Robin boundary conditions, we define the mappings from $V_{i}$ to
$\mathbb{R}^{N_{i}}$ by

$\displaystyle({\mathsf{E}}^{\alpha}_{\partial{\Omega},i}w)_{\ell}$
$\displaystyle:=\partial_{-\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})}w({y}_{i}^{\ell})+\bar{c}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})\,w({y}_{i}^{\ell}),$
(8a) $\displaystyle({\mathsf{I}}^{\alpha}_{\partial{\Omega},i}w)_{\ell}$
$\displaystyle:=\partial_{-\bar{\bar{b}}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})}w({y}_{i}^{\ell})+\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})\,w({y}_{i}^{\ell}),$
(8b) $\displaystyle({\mathsf{F}}^{\alpha}_{\partial{\Omega},i})_{\ell}$
$\displaystyle:={g}_{i}^{\alpha}({y}_{i}^{\ell}),$ (8c)

where
$\bar{c}_{\partial{\Omega},i}^{\alpha},\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha},{g}_{i}^{\alpha}\in
PC(\partial{\Omega})$ and
$\bar{b}_{\partial{\Omega},i}^{\alpha},\bar{\bar{b}}_{\partial{\Omega},i}^{\alpha}\in
PC(\partial{\Omega},\mathbb{R}^{d})$. Here
$\partial_{-\bar{b}_{\partial{\Omega},i}^{\alpha}}$ and
$\partial_{-\bar{\bar{b}}_{\partial{\Omega},i}^{\alpha}}$ are understood as
lower Dini derivatives as in (6).

On the Dirichlet boundary the mappings
${\mathsf{E}}^{\alpha}_{\partial{\Omega},i}$,
${\mathsf{I}}^{\alpha}_{\partial{\Omega},i}$ and
${\mathsf{F}}^{\alpha}_{\partial{\Omega},i}$ implement nodal interpolation.
For $\ell>N_{i}$ we set

$\displaystyle({\mathsf{E}}^{\alpha}_{\partial{\Omega},i}w)_{\ell}$
$\displaystyle:=0,$ (9a)
$\displaystyle({\mathsf{I}}^{\alpha}_{\partial{\Omega},i}w)_{\ell}$
$\displaystyle:=w({y}_{i}^{\ell}),$ (9b)
$\displaystyle({\mathsf{F}}^{\alpha}_{\partial{\Omega},i})_{\ell}$
$\displaystyle:=g({y}_{i}^{\ell}).$ (9c)

We assume a fully implicit discretisation of the region
$\partial{\Omega}_{R}$; additionally, suppose that
$\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha}$ is chosen positive on
$\partial{\Omega}_{R}$, even if $c_{\partial{\Omega}}^{\alpha}=0$.

In summary, we require that the following assumption holds.

######  Assumption 1.

The coefficients satisfy

$\displaystyle\lim_{i\to\infty}\sup_{\alpha\in A}\Bigl{(}$
$\displaystyle\sup_{0\leq\ell\leq
N_{i}}\bigl{\|}a^{\alpha}-\bigl{(}\bar{a}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell})+\bar{\bar{a}}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell})\bigr{)}\bigr{\|}_{L^{\infty}({\rm
supp}\,\hat{\phi}_{i}^{\ell})}$
$\displaystyle+\bigl{\|}b^{\alpha}-\bigl{(}\bar{b}_{{\Omega},i}^{\alpha}+\bar{\bar{b}}_{{\Omega},i}^{\alpha}\bigr{)}\bigr{\|}_{L^{\infty}({\Omega},\mathbb{R}^{d})}+\bigl{\|}c^{\alpha}-\bigl{(}\bar{c}_{{\Omega},i}^{\alpha}+\bar{\bar{c}}_{{\Omega},i}^{\alpha}\bigr{)}\bigr{\|}_{L^{\infty}({\Omega})}$
$\displaystyle+\bigl{\|}f^{\alpha}-{f}_{i}^{\alpha}\bigr{\|}_{L^{\infty}({\Omega})}\Bigr{)}=0$

and

$\displaystyle\lim_{i\to\infty}\sup_{\alpha\in A}\Bigl{(}$
$\displaystyle\bigl{\|}b_{\partial{\Omega}}^{\alpha}-\bigl{(}\bar{b}_{\partial{\Omega},i}^{\alpha}+\bar{\bar{b}}_{\partial{\Omega},i}^{\alpha}\bigr{)}\bigr{\|}_{L^{\infty}(\partial{\Omega})}+\bigl{\|}c_{\partial{\Omega}}^{\alpha}-\bigl{(}\bar{c}_{\partial{\Omega},i}^{\alpha}+\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha}\bigr{)}\bigr{\|}_{L^{\infty}(\partial{\Omega})}$
$\displaystyle+\bigl{\|}f^{\alpha}-f_{i}^{\alpha}\bigr{\|}_{L^{\infty}(\partial{\Omega})}+\bigl{\|}g-g_{i}\bigr{\|}_{L^{\infty}(\partial{\Omega})}+\bigl{\|}g^{\alpha}-g_{i}^{\alpha}\bigr{\|}_{L^{\infty}(\partial{\Omega})}\Bigr{)}=0.$

We require that the family

$\displaystyle\begin{array}[]{ll}\\{(\bar{a}_{{\Omega},i}^{\alpha},\bar{b}_{{\Omega},i}^{\alpha},\bar{c}_{{\Omega},i}^{\alpha},\bar{b}_{\partial{\Omega},i}^{\alpha},\bar{c}_{\partial{\Omega},i}^{\alpha},\bar{\bar{a}}_{{\Omega},i}^{\alpha},\bar{\bar{b}}_{{\Omega},i}^{\alpha},\bar{\bar{c}}_{{\Omega},i}^{\alpha},\bar{\bar{b}}_{\partial{\Omega},i}^{\alpha},\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha},f_{i}^{\alpha},g_{i}^{\alpha})\\}_{\alpha\in
A}\end{array}$

is equicontinuous and depends continuously on $\alpha$. We impose
$\bar{a}_{{\Omega},i}^{\alpha}=\bar{c}_{{\Omega},i}^{\alpha}=\bar{c}_{\partial{\Omega},i}^{\alpha}=0\in\mathbb{R}$
and
$\bar{b}_{{\Omega},i}^{\alpha}=\bar{b}_{\partial{\Omega},i}^{\alpha}=0\in\mathbb{R}^{d}$
as well as $\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha}>0$ on the restriction
to $\partial{\Omega}_{R}$, $i\in\mathbb{N}$.

We define

${\mathsf{E}}_{i}^{\alpha}={\mathsf{E}}_{{\Omega},i}^{\alpha}+{\mathsf{E}}^{\alpha}_{\partial{\Omega},i},\quad{\mathsf{I}}_{i}^{\alpha}={\mathsf{I}}_{{\Omega},i}^{\alpha}+{\mathsf{I}}^{\alpha}_{\partial{\Omega},i},\quad{\mathsf{F}}_{i}^{\alpha}={\mathsf{F}}_{{\Omega},i}^{\alpha}+{\mathsf{F}}^{\alpha}_{\partial{\Omega},i}.$

We also use the notation ${\mathsf{I}}_{i}^{\alpha}$,
${\mathsf{E}}_{i}^{\alpha}$ and ${\mathsf{F}}_{i}^{\alpha}$ for the matrix
representations of exactly these ${\mathsf{I}}_{i}^{\alpha}$,
${\mathsf{E}}_{i}^{\alpha}$ and ${\mathsf{F}}_{i}^{\alpha}$ with respect to
the nodal basis $\\{{\phi}_{i}^{\ell}\\}_{\ell}$ for the trial functions.
Moreover, we assume that the supremum operator is applied componentwise, i.e.
$(\sup_{\alpha}v^{\alpha})_{\ell}=\sup_{\alpha}v_{\ell}^{\alpha}$ for
$v\in\mathbb{R}^{n}$. The expression $a\lesssim b$ means that there exists a
generic constant $C>0$, independent of $i$ and $\alpha$, such that $a\leq Cb$.
Relation $a\gtrsim b$ is defined analogously.

We can now state the numerical scheme used to approximate the solution of (3).
We initialise the scheme by the nodal interpolation of $v_{T}$ so that
$v_{i}(T,\cdot)\in V_{i}^{g}$. Then, in order to find the numerical solution
$v_{i}({s}_{i}^{k},\cdot)\in V_{i}^{g}$, we proceed inductively over the
remaining timesteps $k\in\left\\{T/h_{i}-1,\dots,1,0\right\\}$:

$\displaystyle-d_{i}v_{i}({s}_{i}^{k},\cdot)+\sup_{\alpha\in
A}\bigl{(}{\mathsf{E}}_{i}^{\alpha}v_{i}({s}_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{\alpha}v_{i}({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\bigr{)}=0.$
(10)

We also use an alternative formulation of the numerical scheme. The matrices
${\mathsf{E}}_{i}^{k,w}$, ${\mathsf{I}}_{i}^{k,w}$ and
${\mathsf{F}}_{i}^{k,w}$ are constructed row-wise out of the matrices
${\mathsf{E}}_{i}^{\alpha}$, ${\mathsf{I}}_{i}^{\alpha}$ and
${\mathsf{F}}_{i}^{\alpha}$. More precisely, given a node ${y}_{i}^{\ell}$,
timestep ${s}_{i}^{k}$ and a function $w({s}_{i}^{k},\cdot)\in
H^{1}(\overline{{\Omega}})$, let $\hat{\alpha}$ be a maximiser of

$\sup_{\alpha\in
A}\bigl{(}{\mathsf{E}}_{i}^{\alpha}w({s}_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{\alpha}w({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\bigr{)}_{\ell}=0.$

Note that choice of $\hat{\alpha}$ is not necessarily unique; the analysis is
valid for any choice of such $\hat{\alpha}$. We let the $\ell$th row of
${\mathsf{E}}_{i}^{k,w}$, ${\mathsf{I}}_{i}^{k,w}$ and
${\mathsf{F}}_{i}^{k,w}$ be equal to the $\ell$th row of
$\mathsf{E}_{i}^{\hat{\alpha}}$, $\mathsf{I}_{i}^{\hat{\alpha}}$ and
$\mathsf{F}_{i}^{\hat{\alpha}}$, respectively. In a non-ambiguous case we will
omit explicit mention of $k$ and simply write ${\mathsf{E}}_{i}^{w}$,
${\mathsf{I}}_{i}^{w}$ and ${\mathsf{F}}_{i}^{w}$. We can now reformulate (10)
using the newly constructed operators. We initialize the scheme with the
interpolant $v_{i}(T,\cdot)$. Then $v_{i}\in V_{i}^{g}$ for each
$k\in\\{T/h_{i}-1,\dots,1,0\\}$ and for $0\leq\ell\leq N_{i}^{t}$ solves

$\Bigl{(}(h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf
Id})\,v_{i}({s}_{i}^{k},\cdot)+(h_{i}{\mathsf{E}}_{i}^{k,v_{i}}-{\sf
Id})\,v_{i}({s}_{i}^{k+1},\cdot)-h_{i}{\mathsf{F}}_{i}^{k,v_{i}}\Bigr{)}_{\ell}=0\\\
$ (11a) and for each $k\in\\{T/h_{i}-1,\dots,1,0\\}$ and for $N_{i}^{t}<\ell$
solves
$\Bigl{(}h_{i}{\mathsf{I}}_{i}^{k,v_{i}}\,v_{i}({s}_{i}^{k},\cdot)-h_{i}{\mathsf{F}}_{i}^{k,v_{i}}\Bigr{)}_{\ell}=0,$
(11b)

recalling the implicit discretisation on
$\partial{\Omega}_{R}\cup\partial{\Omega}_{D}$, enforced through (9) and
Assumption 1.

For the sake of convenience let us also introduce the operators
$\hat{\mathsf{I}}_{i}^{k,v_{i}}$, $\hat{\mathsf{E}}_{i}^{k,v_{i}}$ and
$\hat{\mathsf{F}}_{i}$ which combine spatial and temporal terms. The $\ell$th
row of $\hat{\mathsf{I}}_{i}^{k,v_{i}}$, $\hat{\mathsf{E}}_{i}^{k,v_{i}}$ and
$\hat{\mathsf{F}}_{i}$ is equal to that of
$(h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf Id})$,
$(h_{i}{\mathsf{E}}_{i}^{k,v_{i}}-{\sf Id})$ and
$h_{i}{\mathsf{F}}_{i}^{k,v_{i}}$, respectively, if $\leq\ell\leq N_{i}^{t}$.
If $N_{i}^{t}<\ell$, the $\ell$th row is equal to
$(h_{i}{\mathsf{I}}_{i}^{k,v_{i}})$, a zero vector and
$h_{i}{\mathsf{F}}_{i}^{k,v_{i}}$, respectively. For a fixed control $\alpha$,
the operators ${\hat{\mathsf{I}}}_{i}^{\alpha}$,
${\hat{\mathsf{E}}}_{i}^{\alpha}$ and ${\hat{\mathsf{F}}}_{i}^{\alpha}$ are
constructed in an analogous manner. Then for all timesteps ${s}_{i}^{k}$ and
each node ${y}_{i}^{\ell}$ solution $v_{i}$ of (11) solves also:

$\hat{\mathsf{I}}_{i}^{k,v_{i}}\,v_{i}({s}_{i}^{k},{y}_{i}^{\ell})+\hat{\mathsf{E}}_{i}^{k,v_{i}}\,v_{i}({s}_{i}^{k+1},{y}_{i}^{\ell})-\hat{\mathsf{F}}_{i}^{k,v_{i}}=0.$
(12)

######  Remark 1.

To implement the lower Dini derivative
$\partial_{-\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})}w({y}_{i}^{\ell})$
in a computer code, we note that for $\lambda>0$ sufficiently small there is
an element $K\in\mathcal{T}_{i}$ whose closure contains both ${y}_{i}^{\ell}$
and
${y}_{i}^{\ell}+\lambda\,\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})$.
Indeed, choosing $\lambda$ such that
$\lambda\,\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})$ is smaller
than the smallest element edge diameter for all $N_{i}^{{\Omega}}<\ell\leq
N_{i}$ achieves this. We then have

$\partial_{-\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})}w({y}_{i}^{\ell})=\frac{w({y}_{i}^{\ell})-w({y}_{i}^{\ell}+\lambda\,\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell}))}{\lambda}.$

Importantly, because $w\in V_{i}$ is affine on $\overline{K}$, even without
taking a limit $\lambda\to 0$ as on the right-hand side of (6) the Dini
derivative is obtained exactly.

##  4 Monotonicity

In this section we consider monotonicity properties of the discrete
differential operators defined in the previous section. Monotonicity is
crucial for proving the existence of a unique numerical solution of (10) as
well as for establishing convergence to the viscosity solution.

######  Definition 2.

Let us consider $v\in V_{i}$ that has a local non-positive minimum at a node
${y}_{i}^{\ell}$. We say that an operator $F$ satisfies a Local Monotonicity
Property (LMP) if for any such $v$ it follows that $(Fv)_{\ell}\leq 0$.
Additionally, the operator $F$ satisfies the weak Discrete Maximum Principle
(wDMP) provided that, for any $v\in V_{i}$,

$\bigl{(}Fv\bigr{)}_{\ell}\geq
0\;\;\forall\,\ell\in\\{1,\dots,N_{i}\\}\;\;\implies\;\;\min_{{\Omega}\cup\partial{\Omega}_{R}\cup\partial{\Omega}_{t}}\\!\\!v\geq\min\\{\min_{\partial{\Omega}_{D}}v,0\\}.$
(13)

We now describe a method for choosing the artificial diffusion coefficients to
impose the LMP on the matrices ${\mathsf{I}}_{i}^{\alpha}$ and
${\mathsf{E}}_{i}^{\alpha}$. It is based on the assumption of strict acuteness
on the mesh. Consider an element $K\in\mathcal{T}_{i}$ with diameter $\Delta
x_{K}$. For a bounded function $g:\overline{{\Omega}}\to\mathbb{R}^{d}$ we
define $g$’s norm on the restriction to $K$ as

$|g|_{K}:=\Bigl{(}\sum_{j=1}^{d}\Bigl{\|}g_{j}\Bigr{\|}_{L^{\infty}(K)}^{2}\Bigr{)}^{\frac{1}{2}}.$

Then by strict acuteness of the meshes we mean that there exists a
$\theta\in(0,\frac{\pi}{2})$ such that the following holds:

$\nabla{\phi}_{i}^{\ell}\cdot\nabla\phi_{i}^{l}\bigl{|}_{K}\leq-\,\sin(\theta)\;|\nabla{\phi}_{i}^{\ell}|_{K}\;|\nabla\phi^{l}_{i}|_{K}\qquad\forall\ell,l\leq
N_{i},\;\ell\neq l,\;\forall K\in\mathcal{T}_{i}.$ (14)

We say that the family of meshes $\\{\mathcal{T}_{i}\\}_{i}$ is uniformly
strictly acute if $\theta$ does not depend on $i$. As discussed in [BE02], for
$d=2$ and $d=3$ the angle $\theta$ can be interpreted geometrically as
$\frac{\pi}{2}$ minus the largest angle between the pairs of
$(d-1)$-dimensional faces of the element $K$.

###  4.1 The LMP of ${\mathsf{E}}_{{\Omega},i}^{\alpha}$,
${\mathsf{E}}^{\alpha}_{\partial{\Omega},i}$,
${\mathsf{E}}_{{\Omega},i}^{\alpha}$ and
${\mathsf{I}}^{\alpha}_{\partial{\Omega},i}$

Let the functions
$\tilde{a}_{{\Omega},i}^{\alpha},\tilde{\tilde{a}}_{{\Omega},i}^{\alpha},\bar{c}_{{\Omega},i}^{\alpha},\bar{\bar{c}}_{{\Omega},i}^{\alpha}\in
C(\overline{{\Omega}})$ and
$\bar{b}_{{\Omega},i}^{\alpha},\bar{\bar{b}}_{{\Omega},i}^{\alpha}\in
C(\overline{{\Omega}},\mathbb{R}^{d})$ be given. These functions may be chosen
freely as long as Assumption 1 holds. Conceptually
$\tilde{a}_{{\Omega},i}^{\alpha}+\tilde{\tilde{a}}_{{\Omega},i}^{\alpha}\approx
a^{\alpha}$ is the splitting of the second-order coefficients into explicit
and implicit part without the addition of artificial diffusion. With the
addition of artificial diffusion, the coefficients $\bar{a}_{i}^{\alpha}$ and
$\bar{\bar{a}}_{i}^{\alpha}$ of Assumption 1 are obtained.

Indeed, as
$\bar{b}_{{\Omega},i}^{\alpha},\bar{\bar{b}}_{{\Omega},i}^{\alpha},\bar{c}_{{\Omega},i}^{\alpha},\bar{\bar{c}}_{{\Omega},i}^{\alpha}$
are bounded, we can select non-negative artificial diffusion coefficients
$\bar{\nu}_{{\Omega},i}^{\alpha,\ell}$ and
$\bar{\bar{\nu}}_{{\Omega},i}^{\alpha,\ell}$ so that we have for all interior
nodes ${y}_{i}^{\ell}$ and mesh elements $K$ with ${y}_{i}^{\ell}$ as vertex
that

$\displaystyle\begin{array}[]{c}|\bar{b}_{{\Omega},i}^{\alpha}|_{K}\,+\Delta
x_{K}\|\bar{c}_{{\Omega},i}^{\alpha}\|_{L^{\infty}(K)}\leq\bar{\nu}_{{\Omega},i}^{\alpha,\ell}\,\sin(\theta)\,|\nabla\hat{\phi}_{i}^{\ell}|_{K}\,{\rm
vol}(K),\\\ |\bar{\bar{b}}_{{\Omega},i}^{\alpha}|_{K}\,+\Delta
x_{K}\|\bar{\bar{c}}_{{\Omega},i}^{\alpha}\|_{L^{\infty}(K)}\leq\bar{\bar{\nu}}_{{\Omega},i}^{\alpha,\ell}\,\sin(\theta)\,|\nabla\hat{\phi}_{i}^{\ell}|_{K}\,{\rm
vol}(K).\end{array}$ (17)

Now, choosing
$\bar{a}_{{\Omega},i}^{\alpha},\bar{\bar{a}}_{{\Omega},i}^{\alpha}\in
C(\overline{{\Omega}})$ such that

$\displaystyle\bar{a}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell})\geq\max\bigl{\\{}\tilde{a}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell}),\bar{\nu}_{{\Omega},i}^{\alpha,\ell}\bigr{\\}},\qquad\bar{\bar{a}}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell})\geq\max\bigl{\\{}\tilde{\tilde{a}}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell}),\bar{\bar{\nu}}_{{\Omega},i}^{\alpha,\ell}\bigr{\\}},$
(18)

we obtain our splitting of $L^{\alpha}$ into implicit and explicit part.

######  Lemma 1.

Suppose that the mesh $\mathcal{T}_{i}$ is strictly acute and that (17) holds.
Then ${\mathsf{E}}_{{\Omega},i}^{\alpha}$ and
${\mathsf{I}}_{{\Omega},i}^{\alpha}$ satisfy the LMP for all $\alpha$.

###### Proof.

The argument from [JS13, Section 8] for the Dirichlet problem carries over
unchanged for local minima at interior nodes of $v$ from Definition 2. At
boundary nodes the LMP is trivially satisfied as
${\mathsf{E}}_{{\Omega},i}^{\alpha}$ and ${\mathsf{I}}_{{\Omega},i}^{\alpha}$
vanish there. ∎

We now turn to the monotonicity of the discrete boundary operators.

######  Lemma 2.

The operators ${\mathsf{E}}^{\alpha}_{\partial{\Omega},i}$ and
${\mathsf{I}}^{\alpha}_{\partial{\Omega},i}$ satisfy the LMP for all $\alpha$.

###### Proof.

Let $w\in V_{i}$ have a local non-positive minimum at a node
${y}_{i}^{\ell}\in\partial{\Omega}_{t}\cup\partial{\Omega}_{R}$. Then we find
for the lower Dini derivative
$\partial_{-\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})}w({y}_{i}^{\ell})\leq
0$. Also
$\bar{c}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})\,w({y}_{i}^{\ell})\leq
0$ because $c({y}_{i}^{\ell})\geq 0$. Hence
${\mathsf{E}}^{\alpha}_{\partial{\Omega},i}$ admits the LMP. The argument for
${\mathsf{I}}^{\alpha}_{\partial{\Omega},i}$ is analogous. ∎

###  4.2 Monotonicity properties of the ${\mathsf{E}}_{i}^{k,w}$,
${\hat{\mathsf{E}}}_{i}^{k,w}$, ${\mathsf{I}}_{i}^{k,w}$ and
${\hat{\mathsf{I}}}_{i}^{k,w}$

Having examined the basic building blocks of the numerical scheme in the
previous two subsections, we can now analyse the monotonicity properties of
the derived operators ${\hat{\mathsf{E}}}_{i}^{k,w}$ and
${\hat{\mathsf{I}}}_{i}^{k,w}$ as they appear in formulation (12) of the
scheme. We summarise the assumptions made so far in the selection of the
artificial diffusion coefficients.

######  Assumption 2.

Suppose that $\mathcal{T}_{i}$ is strictly acute and that (17) and (18) hold.

First we examine the explicit terms.

######  Lemma 3.

Consider a fixed $w:S_{i}\times\overline{{\Omega}}\to\mathbb{R}$ such that
$w({s}_{i}^{k},\cdot)\in V_{i}$ for all ${s}_{i}^{k}\in S_{i}$. Then the
operators $v\mapsto{\mathsf{E}}_{i}^{k,w}v$ satisfy the LMP and their matrix
has non-positive off-diagonal entries. For $h_{i}$ small enough,
${\hat{\mathsf{E}}}_{i}^{k,w}$ is monotone, i.e. all entries of the matrix
representation are non-positive.

###### Proof.

For any $i$ and $\alpha$,
${\mathsf{E}}_{i}^{\alpha}={\mathsf{E}}_{{\Omega},i}^{\alpha}+{\mathsf{E}}^{\alpha}_{\partial{\Omega},i}$
satisfies the LMP because its summands do. Let us consider a $v\in V_{i}$ that
has a local non-positive minimum at a node ${y}_{i}^{\ell}$. There is an
$\alpha\in A$ such that
$({\mathsf{E}}_{i}^{k,w}v)_{\ell}=({\mathsf{E}}_{i}^{\alpha}v)_{\ell}$. We
know $({\mathsf{E}}_{i}^{\alpha}v)_{\ell}\leq 0$ and therefore that
$v\mapsto{\mathsf{E}}_{i}^{k,w}v$ satisfies the LMP.

For $j\neq\ell$ the hat function $\phi_{i}^{j}$ attains a non-positive minimum
at ${y}_{i}^{\ell}$. Thus, by the LMP, we have that
$({\mathsf{E}}_{i}^{k,w}\phi_{i}^{j})_{\ell}\leq 0$. Hence all the off-
diagonal entries of ${\mathsf{E}}_{i}^{k,w}$ are non-positive.

Owing to Assumption 1, the discretization on $\partial{\Omega}_{R}$ is fully
implicit. Thus the rows of ${\mathsf{E}}_{i}^{k,w}$ belonging to the
discretization on $\partial{\Omega}_{R}$ contain only zeros. Similarly, the
rows linked to $\partial{\Omega}_{D}$ vanish, see (9). All other rows include
a term arising from the time derivative; their structure is
$(h_{i}{\mathsf{E}}_{i}^{k,v_{i}}-{\sf Id})$. Therefore, if $h_{i}$ is
sufficiently small then ${\hat{\mathsf{E}}}_{i}^{k,w}$ is monotone. ∎

Now we turn to the implicit terms.

######  Lemma 4.

Consider a fixed $w:S_{i}\times\overline{{\Omega}}\to\mathbb{R}$ such that
$w({s}_{i}^{k},\cdot)\in V_{i}$ for all ${s}_{i}^{k}\in S_{i}$. Then the
operators $v\mapsto{\mathsf{I}}_{i}^{k,w}v$ satisfy the LMP. Moreover, the
$v\mapsto{\hat{\mathsf{I}}}_{i}^{k,w}\,v$ satisfy the wDMP and their matrix
representation restricted to $V_{i}$ are strictly diagonally dominant
$M$-matrices.

###### Proof.

Analogously to the proof of Lemma 3, the $v\mapsto{\mathsf{I}}_{i}^{k,w}v$
satisfy the LMP and their off-diagonal entries are non-positive.

Before showing the wDMP we verify strict diagonal dominance. By construction,
$v\equiv-1$ attains a non-positive local minimum at each node. Since
${\mathsf{I}}_{i}^{k,w}$ satisfies the LMP property, we have

$\displaystyle
0\geq\left({\mathsf{I}}_{i}^{k,w}v\right)_{\ell}=-\left({\mathsf{I}}_{i}^{k,w}\right)_{\ell\ell}-\sum_{j\neq\ell}\left({\mathsf{I}}_{i}^{k,w}\right)_{\ell
j}.$ (19)

As the off-diagonal entries of ${\mathsf{I}}_{i}^{k,w}$ are non-positive, we
conclude the weak diagonal dominance of the ${\mathsf{I}}_{i}^{k,w}$:

$\left({\mathsf{I}}_{i}^{k,w}\right)_{\ell\ell}-\sum_{j\neq\ell}\left|\left({\mathsf{I}}_{i}^{k,w}\right)_{\ell
j}\right|\geq 0.$

The rows of ${\hat{\mathsf{I}}}_{i}^{k,w}$ which discretise on $\Omega$ and
$\partial{\Omega}_{t}$ are equal to the respective rows of the strictly
diagonally dominant matrix $h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf Id}$.

By Assumption 1 we have $\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha}>0$ on
$\partial{\Omega}_{R}$. Then

$\displaystyle
0>\left({\mathsf{I}}_{i}^{k,w}v\right)_{\ell}=\left({\hat{\mathsf{I}}}_{i}^{k,w}v\right)_{\ell}\quad\forall\,{y}_{i}^{\ell}\in\partial{\Omega}_{R}.$
(20)

Using the same argument as above, but noting the strict inequality of (20)
compared to (19), we conclude the strict diagonal dominance for rows linked to
$\partial{\Omega}_{R}$. On $\partial{\Omega}_{D}$ the rows resemble an
identity matrix, giving also strict diagonal dominance. It follows that the
${\hat{\mathsf{I}}}_{i}^{k,w}$ are invertible $M$-matrices because [BP94,
Chapter 6, Theorem 2.3, $(M_{35})$] applies as ${\hat{\mathsf{I}}}_{i}^{k,w}$
is a $Z$-matrix.

Finally, consider a $v\in V_{i}$ with
$\min_{{\Omega}\cup\partial{\Omega}_{t}\cup\partial{\Omega}_{t}}\\!\\!v<\min\\{\min_{\partial{\Omega}_{D}}v,0\\}.$
Let ${y}_{i}^{\ell}$ be a non-Dirichlet node, where the negative, global
minimum of $v$ is attained. Since ${\hat{\mathsf{I}}}_{i}^{k,w}$ is a strictly
diagonal dominant $M$-matrix it follows that
$({\hat{\mathsf{I}}}_{i}^{k,w}v)_{\ell}<0$. Hence
${\hat{\mathsf{I}}}_{i}^{k,w}$ admits the wDMP. ∎

###  4.3 Scaling of the artificial diffusion coefficients

In order to achieve convergence of the numerical scheme we expect the
artificial diffusion coefficients $\bar{\nu}_{{\Omega},i}^{\alpha,\ell}$,
$\bar{\bar{\nu}}_{{\Omega},i}^{\alpha,\ell}$ to vanish in the limit
$i\to\infty$.

Suppose that (14) holds uniformly for some $\theta$. In this subsection we
suppose $\bar{\nu}_{{\Omega},i}^{\alpha,\ell}$ are chosen quasi-optimally with
regard to (17), meaning

$\displaystyle\bar{\nu}_{{\Omega},i}^{\alpha,\ell}$
$\displaystyle\lesssim\sup\Bigl{\\{}\frac{|\bar{b}_{{\Omega},i}^{\alpha}|_{K}\,+\Delta
x_{K}\|\bar{c}_{{\Omega},i}^{\alpha}\|_{L^{\infty}(K)}}{\sin(\theta)\,|\nabla\hat{\phi}_{i}^{\ell}|_{K}\,{\rm
vol}(K)}\;\Big{|}\;K\subset\textrm{supp}\;{\phi}_{i}^{\ell}\Bigr{\\}}$ (21)
$\displaystyle\bar{\bar{\nu}}_{{\Omega},i}^{\alpha,\ell}$
$\displaystyle\lesssim\sup\Bigl{\\{}\frac{|\bar{\bar{b}}_{{\Omega},i}^{\alpha}|_{K}\,+\Delta
x_{K}\|\bar{\bar{c}}_{{\Omega},i}^{\alpha}\|_{L^{\infty}(K)}}{\sin(\theta)\,|\nabla\hat{\phi}_{i}^{\ell}|_{K}\,{\rm
vol}(K)}\;\Big{|}\;K\subset\textrm{supp}\;{\phi}_{i}^{\ell}\Bigr{\\}}.$ (22)

Generally, in implementations of the algorithm quasi-optimally is more easily
achieved than optimality. Because of shape-regularity of the domain one has
$|\nabla\hat{\phi}_{i}^{\ell}|_{K}\,{\rm vol}(K)\gtrsim\frac{1}{\Delta
x_{K}}$. We conclude that quasi-optimal artificial diffusion coefficients
satisfy

$\displaystyle{\mathsf{O}}(\bar{\nu}_{{\Omega},i}^{\alpha,\ell})={\mathsf{O}}(\bar{\bar{\nu}}_{{\Omega},i}^{\alpha,\ell})=\Delta
x_{K}.$ (23)

We now turn our attention to the time step restrictions imposed through the
quasi-optimality (21). Recall that in order for the explicit operators to be
monotone we require all their entries in matrix representation to be non-
positive. This is satisfied trivially for nodes on
$\partial{\Omega}_{R}\cup\partial{\Omega}_{D}$ where we use a fully implicit
scheme. Therefore let us consider non-positivity of the diagonal terms of
$h_{i}{\mathsf{E}}_{i}^{\alpha}-{\sf Id}$ on the complement
${\Omega}\cup\partial{\Omega}_{t}$. For ${y}_{i}^{\ell}\in{\Omega}$ this
translates into the condition

$\displaystyle 1\geq\,$ $\displaystyle
h_{i}\,\bigl{(}\bar{a}_{{\Omega},i}^{\alpha}({y}_{i}^{\ell})\langle\nabla{\phi}_{i}^{\ell},\nabla\hat{\phi}_{i}^{\ell}\rangle+\langle-\bar{b}_{{\Omega},i}^{\alpha}\cdot\nabla{\phi}_{i}^{\ell}+\bar{c}_{{\Omega},i}^{\alpha}\,{\phi}_{i}^{\ell},\hat{\phi}_{i}^{\ell}\rangle\bigr{)}$

and for ${y}_{i}^{\ell}\in\partial{\Omega}_{t}$

$1\geq\,h_{i}\left(\partial_{-\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})}{\phi}_{i}^{\ell}({y}_{i}^{\ell})+\bar{c}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})\right).$

Because

$\displaystyle\langle\nabla{\phi}_{i}^{\ell},\nabla\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle={\mathsf{O}}\bigl{(}(\Delta x_{K})^{-2}\bigr{)},$
$\displaystyle\langle\nabla{\phi}_{i}^{\ell},\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle={\mathsf{O}}\bigl{(}(\Delta x_{K})^{-1}\bigr{)},$
$\displaystyle\langle{\phi}_{i}^{\ell},\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle={\mathsf{O}}\bigl{(}1\bigr{)},$
$\displaystyle\partial_{-\bar{b}_{\partial{\Omega},i}^{\alpha}({y}_{i}^{\ell})}{\phi}_{i}^{\ell}({y}_{i}^{\ell})$
$\displaystyle={\mathsf{O}}\bigl{(}(\Delta x_{K})^{-1}\bigr{)},$

we find $h_{i}={\mathsf{O}}((\Delta x_{K})^{2})$ if
$\|\bar{a}_{{\Omega},i}^{\alpha}\|_{\infty}>0$. Otherwise if
$\|\bar{b}_{{\Omega},i}^{\alpha}\|_{\infty}>0$ or
$\|\bar{b}_{\partial{\Omega},i}^{\alpha}\|_{\infty}>0$ we have
$h_{i}={\mathsf{O}}(\Delta x_{K})$ and if $\bar{a}_{{\Omega},i}^{\alpha}$,
$\bar{b}_{\partial{\Omega},i}^{\alpha}$ and $\bar{b}_{{\Omega},i}^{\alpha}$
vanish but not $\bar{c}_{{\Omega},i}^{\alpha}$ or
$\bar{c}_{\partial{\Omega},i}^{\alpha}$, then $h_{i}={\mathsf{O}}(1)$. If also
$\bar{c}_{{\Omega},i}^{\alpha}=\bar{c}_{\partial{\Omega},i}^{\alpha}=0$ then
there is no restriction on $h_{i}$, i.e. fully implicit discretisations are
monotone for any $h_{i}>0$.

##  5 Existence of numerical solutions

The discrete non-linear problem (12) can be solved by a version of Howard’s
algorithm discussed in [BMZ09]. We now present its formulation in our setting.

######  Algorithm 1.

Given are timestep $k\in\\{0,\dots,T/h_{i}-1\\}$, solution
$v_{i}(s_{i}^{k+1},\cdot)\in V_{i}$ at timestep $k+1$ and an (arbitrary)
choice of $\alpha\in A$. Find $w_{0}\in V_{i}$ such that

${\hat{\mathsf{I}}}_{i}^{\alpha}w_{0}={\hat{\mathsf{F}}}_{i}^{\alpha}-{\hat{\mathsf{E}}}_{i}^{\alpha}v_{i}(s_{i}^{k+1},\,\cdot).$

Inductively over $m\in\mathbb{N}$, compute $w_{m+1}\in V_{i}$ such that

$\hat{\mathsf{I}}^{w_{m}}_{i}w_{m+1}=\hat{\mathsf{F}}^{w_{m}}_{i}-\hat{\mathsf{E}}^{w_{m}}_{i}v_{i}(s_{i}^{k+1},\cdot).$
(24)

To show the convergence of the sequence $(w_{m})_{m}$ to the solution of (10)
we appeal to an auxiliary problem: for some fixed control $\alpha\in A$ we
consider the linear evolution problem associated to it. More precisely, we
define ${v}_{i}^{\alpha}\colon S_{i}\to V_{i}$ to be such that
${v}_{i}^{\alpha}(T,\cdot)=v_{i}(T,\cdot)$, the interpolant of $v_{T}$, and
for each $k\in\\{0,\dots,T/h_{i}-1\\}$

${\hat{\mathsf{I}}}_{i}^{\alpha}\,{v}_{i}^{\alpha}({s}_{i}^{k},\cdot)+{\hat{\mathsf{E}}}_{i}^{\alpha}\,{v}_{i}^{\alpha}({s}_{i}^{k+1},\cdot)-{\hat{\mathsf{F}}}_{i}^{\alpha}=0.$
(25)

Notice that ${v}_{i}^{\alpha}$ is well-defined due to the invertibility of
${\hat{\mathsf{I}}}_{i}^{\alpha}$.

######  Theorem 1.

There exists a unique numerical solution $v_{i}\colon S_{i}\to V_{i}$ which
solves (10) and (12). Algorithm 1, provided with the inputs $k$,
$v_{i}(s_{i}^{k+1},\cdot)$ and $\alpha$, generates a sequence $(w_{m})_{m}$
which converges superlinearly to $v_{i}({s}_{i}^{k},\cdot)$ as $m\to\infty$.
Moreover, $0\leq v_{i}\leq{v}_{i}^{\alpha}$ for all $\alpha\in A$.

###### Proof.

For a fixed timestep $k$, the superlinear convergence of Algorithm 1 to the
unique solution $v_{i}({s}_{i}^{k},\cdot)$ is shown in [BMZ09, Theorem 2.1]
under Assumptions (H1) and (H2) stated therein. Condition (H1) requires the
inverse positivity of the operators $\hat{\mathsf{I}}^{w_{m}}_{i}$, which
holds because according to Lemma 4 every $\hat{\mathsf{I}}^{w_{m}}_{i}$ is a
non-singular $M$-matrix. Condition (H2) requires that $\alpha\in
A\mapsto-{\hat{\mathsf{I}}}_{i}^{\alpha}$ and $\alpha\in
A\mapsto{\hat{\mathsf{E}}}_{i}^{\alpha}\,v_{i}({s}_{i}^{k+1},\cdot)-{\hat{\mathsf{F}}}_{i}^{\alpha}$
are continuous, which follows from Assumption 1. Induction over timesteps $k$
gives existence and uniqueness of the solution $v_{i}$.

We now show that $v_{i}\geq 0$ on $S_{i}\times\overline{{\Omega}}$ by
induction over $k$. Firstly, we notice that $v_{i}(T,\cdot)\geq 0$ because we
assumed that $v_{T}\geq 0$ on $\overline{{\Omega}}$ and the same holds for its
interpolant. Let us assume $v_{i}(s_{i}^{k+1},\cdot)\geq 0$ on
$\overline{{\Omega}}$ for some ${s}_{i}^{k+1}\in S_{i}$. Due to the LMP, all
entries of $\hat{\mathsf{E}}_{i}^{v_{i}}$ are non-positive and by assumption
all entries of $\hat{\mathsf{F}}_{i}^{v_{i}}$ are non-negative. Therefore,
using (12) we have that

$\displaystyle\hat{\mathsf{I}}_{i}^{v_{i}}v_{i}({s}_{i}^{k},\cdot)$
$\displaystyle=-\hat{\mathsf{E}}_{i}^{v_{i}}v_{i}({s}_{i}^{k+1},\cdot)+\hat{\mathsf{F}}_{i}^{v_{i}}\geq
0.$

We conclude that $v_{i}(s_{i}^{k},\cdot)\geq 0$ on $\overline{{\Omega}}$ due
to the inverse positivity of $\hat{\mathsf{I}}_{i}^{v_{i}}$.

We now prove the last statement, namely $v_{i}\leq v_{i}^{\alpha}$ for all
$\alpha\in A$, by induction over $k$. Consider any $\alpha\in A$. At time $T$
both $v_{i}$ and ${v}_{i}^{\alpha}$ interpolate $v_{T}$ and hence are equal.
Let us assume that for some $k\leq T/h_{i}-1$, $v_{i}(s_{i}^{k+1},\cdot)\leq
v_{i}^{\alpha}(s_{i}^{k+1},\cdot)$. From (10),

${\hat{\mathsf{I}}}_{i}^{\alpha}v_{i}(s_{i}^{k},\cdot)\leq{\hat{\mathsf{F}}}_{i}^{\alpha}-{\hat{\mathsf{E}}}_{i}^{\alpha}v_{i}(s_{i}^{k+1},\cdot).$

Now subtracting (25) from the above inequality, together with the monotonicity
of ${\hat{\mathsf{E}}}_{i}^{\alpha}$, gives

$\displaystyle{\hat{\mathsf{I}}}_{i}^{\alpha}\left(v_{i}(s_{i}^{k},\cdot)-v_{i}^{\alpha}(s_{i}^{k},\cdot)\right)$
$\displaystyle\leq{\hat{\mathsf{E}}}_{i}^{\alpha}\left(v_{i}^{\alpha}(s_{i}^{k+1},\cdot)-v_{i}(s_{i}^{k+1},\cdot)\right)\leq
0.$

Using the inverse positivity of ${\hat{\mathsf{I}}}_{i}^{\alpha}$ gives us
$v_{i}(s_{i}^{k},\cdot)-v_{i}^{\alpha}(s_{i}^{k},\cdot)\leq 0$ on
$\overline{{\Omega}}$, as required. ∎

##  6 Consistency

We will assume existence of an elliptic projection $P_{i}$, described in
[JS13], with the properties required in the following assumption.

######  Assumption 3.

There are linear mappings $P_{i}:C(H^{1}({\Omega}))\to V_{i}$ satisfying for
all interior hat functions $\hat{\phi}_{i}^{\ell}$, $\ell\leq N_{i}^{\Omega}$,

$\langle\nabla P_{i}w,\nabla\hat{\phi}_{i}^{\ell}\rangle=\langle\nabla
w,\nabla\hat{\phi}_{i}^{\ell}\rangle.$ (26)

There is a constant $C\geq 0$ such that for every $w\in
C^{\infty}(\mathbb{R}^{d})$ and $i\in\mathbb{N}$,

$\left\|P_{i}w\right\|_{W^{1,\infty}({\Omega})}\leq
C\left\|w\right\|_{W^{1,\infty}({\Omega})}\quad\text{and}\quad\lim_{i\to\infty}\left\|P_{i}w-w\right\|_{W^{1,\infty}({\Omega})}=0.$
(27)

To state consistency it is convenient to abbreviate the operator of the
numerical scheme as

$\displaystyle
F_{i}\,w({s}_{i}^{k},{y}_{i}^{\ell}):=\hat{\mathsf{I}}_{i}^{k,w}\,w({s}_{i}^{k},{y}_{i}^{\ell})+\hat{\mathsf{E}}_{i}^{k,w}\,w({s}_{i}^{k+1},{y}_{i}^{\ell})-\hat{\mathsf{F}}_{i}$
(28)

for $w({s}_{i}^{k},\cdot)\in V_{i}$. Note that while $F_{i}$ is the discrete
operator approximating the continuous operator $F$ defined in section 2, the
notationally similar $\hat{\mathsf{F}}_{i}^{\alpha}$ represents the
approximation of $f^{\alpha}$, $g^{\alpha}$ and $g$ as explained at the end of
section 3.

######  Theorem 2.

Let $\psi\in C^{2}(\mathbb{R}\times\mathbb{R}^{d})$, $s_{i}^{k(i)}\to
t\in[0,T)$ and $y_{i}^{\ell(i)}\to x\in\overline{{\Omega}}$ as $i\to\infty$.
Here $s_{i}^{k(i)}$ is a time step and $y_{i}^{\ell(i)}$ a node of the $i$-th
refinement. Then

$\displaystyle\limsup_{i\to\infty}F_{i}P_{i}\psi(s_{i}^{k(i)},y_{i}^{\ell(i)})\leq
F^{*}(t,x,\Delta\psi(t,x),\nabla\psi(t,x),\partial_{t}\psi(t,x),\psi(t,x))$
(29)

and

$\displaystyle\liminf_{i\to\infty}F_{i}P_{i}\psi(s_{i}^{k(i)},y_{i}^{\ell(i)})\geq
F_{*}(t,x,\Delta\psi(t,x),\nabla\psi(t,x),\partial_{t}\psi(t,x),\psi(t,x)).$
(30)

###### Proof.

We prove (29). The result for (30) follows analogously. For ease of notation,
the dependence of $k$ and $\ell$ on $i$ is made implicit.

Step 1: Standard finite difference bounds ensure that if
${y}_{i}^{\ell}\in{\Omega}\cup\partial{\Omega}_{t}$ then

$\lim_{i\to\infty}d_{i}P_{i}\psi({s}_{i}^{k},{y}_{i}^{\ell})=\partial_{t}\psi(t,x).$
(31)

Otherwise, if ${y}_{i}^{\ell}\in\partial{\Omega}_{R}\cup\partial{\Omega}_{D}$
then

$\lim_{i\to\infty}d_{i}P_{i}\psi({s}_{i}^{k},{y}_{i}^{\ell})=0.$ (32)

Step 2: It is shown in [JS13, Section 4] that if ${y}_{i}^{\ell}\in{\Omega}$
then

$\lim_{i\to\infty}\left({\mathsf{E}}_{i}^{\alpha}P_{i}\psi({s}_{i}^{\ell},\cdot)+{\mathsf{I}}_{i}^{\alpha}P_{i}\psi({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\right)_{\ell}=L^{\alpha}\psi(t,x)-f^{\alpha}(x),$
(33)

where convergence to the limit is uniform over all $\alpha\in A$. We remark
that the orthogonality (26) is used in this step.

Step 3: Now suppose that ${y}_{i}^{\ell}\in\partial{\Omega}_{D}$. Then it
follows from (27) that

$\lim_{i\to\infty}F_{i}P_{i}\psi({s}_{i}^{k},{y}_{i}^{\ell})=\psi(t,x)-g(x).$
(34)

Step 4: Let ${y}_{i}^{\ell}\in\partial{\Omega}_{t}\cup\partial{\Omega}_{R}$.
Just like the continuous operators the corresponding first-order terms of the
discrete Robin operators employ the lower Dini derivative, giving consistency
directly. Thus with

$\displaystyle\lim_{i\to\infty}$
$\displaystyle\left|c_{\partial{\Omega}}^{\alpha}\,P_{i}\psi(t,\cdot)-\bar{\bar{c}}_{\partial{\Omega},i}^{\alpha}\,P_{i}\psi(s_{i}^{k},\cdot)-\bar{c}_{\partial{\Omega},i}^{\alpha}\,P_{i}\psi(s_{i}^{k+1},\cdot)\right|=0,$
(35)

using Assumptions 1 and 3, we conclude that if
${y}_{i}^{\ell}\in\partial{\Omega}_{t}\cup\partial{\Omega}_{R}$ then

$\lim_{i\to\infty}\left({\mathsf{E}}_{i}^{\alpha}P_{i}\psi({s}_{i}^{\ell},\cdot)+{\mathsf{I}}_{i}^{\alpha}P_{i}\psi({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\right)_{\ell}=L_{\partial{\Omega}}^{\alpha}\psi(t,x)-g^{\alpha}(x).$
(36)

Step 5: Consider the sequence $\\{({s}_{i}^{k},{y}_{i}^{\ell})\\}_{i}$ as
specified in the statement of the theorem, in particular with
${y}_{i}^{\ell}\in\overline{{\Omega}}$. We decompose
$\\{({s}_{i}^{k},{y}_{i}^{\ell})\\}_{i}$ into the subsequences of the
$({s}_{i}^{k},{y}_{i}^{\ell})$ where ${y}_{i}^{\ell}$ belongs to ${\Omega}$,
$\partial{\Omega}_{t}$, $\partial{\Omega}_{R}$ and $\partial{\Omega}_{D}$,
respectively. Then the conclusions of Steps 1 to 4 above may be applied to the
individual subsequences. ∎

##  7 Stability

In this section we present a lemma which ensures $L^{\infty}$ stability of the
numerical scheme (12). The stability statement goes back to the boundedness of
a supersolution of the continuous linear problem for a fixed $\alpha$. Let

$\displaystyle F^{\alpha}(t,x,q,p,r,s)=$
$\displaystyle\left\\{\begin{array}[]{rllll}-r+L^{\alpha}(x,q,p,s)-f^{\alpha}(x)&\text{if
}(t,x)\in[0,T)\times\overline{{\Omega}},\\\
-r+L^{\alpha}_{\partial{\Omega}}(x,p,s)-g^{\alpha}(x)&\text{if
}(t,x)\in[0,T)\times\partial{\Omega}_{t},\\\
\phantom{-r+\;}L^{\alpha}_{\partial{\Omega}}(x,p,s)-g^{\alpha}(x)&\text{if
}(t,x)\in[0,T)\times\partial{\Omega}_{R},\\\ s-g(x)&\text{if
}(t,x)\in[0,T)\times\partial{\Omega}_{D},\\\ s-v_{T}(x)&\text{if
}(t,x)\in\\{T\\}\times\overline{{\Omega}}.\end{array}\right.$

######  Assumption 4.

There exists an $\alpha\in A$ and a $w(t,x)\in
C^{2}(\mathbb{R}\times\mathbb{R}^{d})$ which is a strict supersolution of the
associated linear problem. More precisely, there is an $\varepsilon>0$ such
that

$F^{\alpha}(t,x,\Delta w(t,x),\nabla
w(t,x),\partial_{t}w(t,x),w(t,x))\geq\varepsilon$

on $[0,T]\times\overline{{\Omega}}$.

The assumption is essentially fulfilled if the linear equation

$\displaystyle
F^{\alpha}(t,x,\Delta\psi(t,x),\nabla\psi(t,x),\partial_{t}\psi(t,x),\psi(t,x))=2\varepsilon$
(37)

is a well-posed problem in a suitable sense. For example, $\psi$ may be a weak
solution of (37) as in [OR73, Chapter 1] which admits a bounded extension to
$\mathbb{R}\times\mathbb{R}^{d}$. In such cases one may pass to a strict
supersolution in $C^{2}(\mathbb{R}\times\mathbb{R}^{d})$ through
mollification.

######  Lemma 5.

Let $\alpha$ be as in Assumption 4, $s_{i}^{k(i)}\to t\in[0,T)$ and
$y_{i}^{\ell(i)}\to x\in\overline{{\Omega}}$ as $i\to\infty$. Here
$s_{i}^{k(i)}$ is a time step and $y_{i}^{\ell(i)}$ a node of the $i$-th
refinement. Then

$\displaystyle\liminf_{i\to\infty}\left[{\hat{\mathsf{I}}}_{i}^{\alpha}P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})+{\hat{\mathsf{E}}}_{i}^{\alpha}P_{i}w({s}_{i}^{k+1},{y}_{i}^{\ell})-({\hat{\mathsf{F}}}_{i}^{\alpha})_{\ell}\right]\geq\varepsilon.$

###### Proof.

We use Theorem 2 for a singleton control set $A=\\{\alpha\\}$ to cover the
linear case. The result now follows because

$\displaystyle(F^{\alpha})_{*}(x,\Delta w(t,x),\nabla
w(t,x),\partial_{t}w(t,x),w(t,x))\geq\varepsilon,$

owing to Assumption 4. ∎

######  Theorem 3.

The numerical solutions $v_{i}$ are uniformly bounded in the $L^{\infty}$
norm. More precisely, there exists a finite constant $C>0$ such that

$\displaystyle\lVert
v_{i}\rVert_{L^{\infty}(S_{i}\times\overline{{\Omega}})}\leq
C\qquad\forall\,i\in\mathbb{N}.$ (38)

###### Proof.

Recall the solution $v_{i}^{\alpha}$ of the linear problem (25). We define

$\displaystyle{w}_{i}^{k}$
$\displaystyle:=P_{i}w^{\alpha}({s}_{i}^{k},\cdot),$
$\displaystyle{\tilde{v}}_{i}^{k}$
$\displaystyle:={w}_{i}^{k}-v_{i}^{\alpha}({s}_{i}^{k},\cdot).$

It is convenient to set

$\displaystyle{e}_{i}^{k}:={\hat{\mathsf{I}}}_{i}^{\alpha}{\tilde{v}}_{i}^{k}+{\hat{\mathsf{E}}}_{i}^{\alpha}{\tilde{v}}_{i}^{k+1}={\hat{\mathsf{I}}}_{i}^{\alpha}{w}_{i}^{k}+{\hat{\mathsf{E}}}_{i}^{\alpha}{w}_{i}^{k+1}-{\hat{\mathsf{F}}}_{i}^{\alpha}.$
(39)

Because of Assumptions 3 and 4, the $w_{i}$ are uniformly bounded in the
$L^{\infty}$ norm. Moreover, $0\leq v_{i}\leq v_{i}^{\alpha}$ due to Theorem
1. Thus the statement of the theorem is proved once we demonstrate that the
$v_{i}^{\alpha}$ are bounded from above independently of $i$. This is
equivalent to showing a lower bound for the ${\tilde{v}}_{i}^{k}$.

It follows from Lemma 5 that $e^{k}\geq 0$ for $i$ larger than some constant
$M$. It is trivial that there exists a constant $C$ such that the inequality
of (38) holds for all $i\leq M$. We therefore may assume w.l.o.g. that
$e^{k}\geq 0$ throughout. By possibly modifying $w$ through the addition of a
positive constant we can assume that

$w_{i}^{T/h_{i}}=P_{i}w(T,\cdot)\geq\|v_{T}(T,\cdot)\|_{L^{\infty}({\Omega})}$

while maintaining the supersolution property of $w$ because
$c_{{\Omega}}^{\alpha},c_{\partial{\Omega}}^{\alpha}\geq 0$. This implies
$\tilde{v}_{i}^{T/h_{i}}\geq 0$, i.e. the non-negativity at the final time.

Now suppose ${\tilde{v}}_{i}^{k+1}\geq 0$. Then

${\tilde{v}}_{i}^{k}=({\hat{\mathsf{I}}}_{i}^{\alpha})^{-1}\bigl{(}{e}_{i}^{k}-{\hat{\mathsf{E}}}_{i}^{\alpha}{\tilde{v}}_{i}^{k+1}\bigr{)}\geq
0$

because also $({\hat{\mathsf{I}}}_{i}^{\alpha})^{-1}\geq 0$ and
$-{\hat{\mathsf{E}}}_{i}^{\alpha}{\tilde{v}}_{i}^{k+1}\geq 0$. Now induction
in $k$ completes the proof. ∎

##  8 Uniform Convergence

Analogously to the envelopes of functions introduced in Section 2 we define
envelopes of the numerical solutions as follows

$\overline{v}(t,x)=\limsup_{i\to\infty}\sup_{({s}_{i}^{k},{y}_{i}^{\ell})\to(t,x)}v_{i}({s}_{i}^{k},{y}_{i}^{\ell}),\qquad\underline{v}(t,x)=\liminf_{i\to\infty}\inf_{({s}_{i}^{k},{y}_{i}^{\ell})\to(t,x)}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})$

where limits are taken over all sequences of nodes in
$[0,T]\times\overline{{\Omega}}$ which converge to
$(t,x)\in[0,T]\times\overline{{\Omega}}$. Owing to Theorem 3, $\overline{v}$
and $\underline{v}$ attain finite values. By construction, $\overline{v}$ is
upper and $\underline{v}$ lower semi-continuous and
$\underline{v}\leq\overline{v}$.

######  Theorem 4.

The function $\overline{v}$ is a viscosity subsolution and $\underline{v}$ is
a viscosity supersolution.

###### Proof.

Step 1 ($\overline{v}$ is a subsolution). To show that $\overline{v}$ is a
viscosity subsolution, suppose that $w\in
C^{\infty}(\mathbb{R}\times\mathbb{R}^{d})$ is a test function such that
$\overline{v}-w$ has a strict local maximum at
$(s,y)\in(0,T)\times\overline{{\Omega}}$, with $\overline{v}(s,y)=w(s,y)$.
Note that $(s,y)$ may be on the boundary. Consider a closed neighbourhood
$B:=\bigl{\\{}(t,x)\in(0,T)\times\overline{{\Omega}}\;:\;|t-s|+|x-y|\leq\delta\bigr{\\}}$
with $\delta>0$ such that

$\overline{v}(s,y)-w(s,y)>\overline{v}(t,x)-w(t,x)\quad\forall(t,x)\in
B\setminus(s,y).$

Choose $i$ sufficiently large for $B$ to contain nodes. As in [JS13] we choose
a sequence of nodes $\\{(s_{i(j)}^{k},y_{i(j)}^{\ell})\\}_{j}$ which maximise
$v_{i}(s_{i(j)}^{\kappa},y_{i(j)}^{\lambda})-P_{i}w(s_{i(j)}^{\kappa},y_{i(j)}^{\lambda})$
among all nodes $(s_{i(j)}^{\kappa},y_{i(j)}^{\lambda})\in B$ and converge to
$(s,y)$. It follows that

$\displaystyle
v_{i}({s}_{i}^{k},{y}_{i}^{\ell})-P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})\rightarrow\overline{v}(s,y)-w(s,y)=0.$
(40)

Moreover, because of $({s}_{i}^{k},{y}_{i}^{\ell})\rightarrow(s,y)$, the
neighbours of the $({s}_{i}^{k},{y}_{i}^{\ell})$ eventually also belong to
$B$: for $i$ sufficiently large, we have
$({s}_{i}^{\kappa},{y}_{i}^{\lambda})\in B$ if
$\kappa\in\left\\{k,k+1\right\\}$ and
${y}_{i}^{\lambda}\in\operatorname{supp}\hat{\phi}_{i}^{\ell}$, in which case

$\displaystyle
v_{i}({s}_{i}^{\kappa},{y}_{i}^{\lambda})-P_{i}w({s}_{i}^{\kappa},{y}_{i}^{\lambda})\leq
v_{i}({s}_{i}^{k},{y}_{i}^{\ell})-P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})$
$\displaystyle\Leftrightarrow\;P_{i}w({s}_{i}^{\kappa},{y}_{i}^{\lambda})+\mu_{i}\geq
v_{i}({s}_{i}^{\kappa},{y}_{i}^{\lambda}),$

with
$\mu_{i}=v_{i}({s}_{i}^{k},{y}_{i}^{\ell})-P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})$,
and $\mu_{i}\to 0$ as $i\to\infty$ because of (40).

Recall that the matrices ${\mathsf{E}}_{i}^{\alpha}$ have non-zero off
diagonal entries $\left({\mathsf{E}}_{i}^{\alpha}\right)_{\ell\lambda}$ only
if ${y}_{i}^{\lambda}\in\operatorname{supp}\hat{\phi}_{i}^{\ell}$ and that
$v_{i}(s_{i}^{k+1},\cdot)\leq P_{i}w(s_{i}^{k+1},\cdot)+\mu_{i}$ on
$\operatorname{supp}\hat{\phi}_{i}^{\ell}$. Therefore, monotonicity of
$h_{i}{\mathsf{E}}_{i}^{\alpha}-{\sf Id}$ for all $\alpha\in A$ implies that

$\left((h_{i}{\mathsf{E}}_{i}^{\alpha}-{\sf
Id})\left[P_{i}w(s^{k+1}_{i},\cdot)+\mu_{i}\right]\right)_{\ell}\leq\left((h_{i}{\mathsf{E}}_{i}^{\alpha}-{\sf
Id})v_{i}(s_{i}^{k+1},\cdot)\right)_{\ell}.$

Applying the LMP and linearity of ${\mathsf{I}}_{i}^{\alpha}$ to
$P_{i}w({s}_{i}^{k},\cdot)+\mu_{i}-v_{i}({s}_{i}^{k},\cdot)$, which has a non-
positive local minimum at ${y}_{i}^{\ell}$, yields

$\left((h_{i}{\mathsf{I}}_{i}^{\alpha}+{\sf
Id})\left[P_{i}w({s}_{i}^{k},\cdot)+\mu_{i}\right]\right)_{\ell}\leq\left((h_{i}{\mathsf{I}}_{i}^{\alpha}+{\sf
Id})v_{i}({s}_{i}^{k},\cdot)\right)_{\ell}.$

From the definition of the scheme, with
$\gamma:=\sup_{\alpha,i}\|\bar{c}_{i}^{\alpha}+\bar{\bar{c}}_{i}^{\alpha}\|_{\infty}$,

$\displaystyle 0=\,$ $\displaystyle-
d_{i}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})+\sup_{\alpha\in
A}\left({\mathsf{E}}_{i}^{\alpha}v_{i}(s_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{\alpha}v_{i}({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\right)_{\ell}$
$\displaystyle\geq\,$ $\displaystyle-
d_{i}\\!\left(P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})+\mu_{i}\right)+\sup_{\alpha\in
A}\\!\left({\mathsf{E}}_{i}^{\alpha}\\!\left(P_{i}w(s_{i}^{k+1},\cdot)+\mu_{i}\right)+{\mathsf{I}}_{i}^{\alpha}\\!\left(P_{i}w({s}_{i}^{k},\cdot)+\mu_{i}\right)-{\mathsf{F}}_{i}^{\alpha}\right)_{\ell}$
$\displaystyle=\,$ $\displaystyle-
d_{i}P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})+\sup_{\alpha\in
A}\left[\left({\mathsf{E}}_{i}^{\alpha}P_{i}w(s_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{\alpha}P_{i}w({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\right)_{\ell}+\mu_{i}\langle\bar{c}_{i}^{\alpha}+\bar{\bar{c}}_{i}^{\alpha},\hat{\phi}_{i}^{\ell}\rangle\right]$
$\displaystyle\geq\,$ $\displaystyle-
d_{i}P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})+\sup_{\alpha\in
A}\left({\mathsf{E}}_{i}^{\alpha}P_{i}w(s_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{\alpha}P_{i}w({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\right)_{\ell}-\gamma\left|\mu_{i}\right|$
$\displaystyle=\,$ $\displaystyle
F_{i}P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})-\gamma\left|\mu_{i}\right|.$ (41)

For a fixed $i$, evaluating $F_{i}P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})$ may
involve a boundary operator even if $(s,y)$ is internal and vice versa may
involve the PDE operator even if $(s,y)$ belongs to the boundary. Referring to
the semi-continuous envelope $F_{*}$, it now follows from (41),
$\lim_{i}\mu_{i}=0$ and Theorem 2 that

$\displaystyle 0$
$\displaystyle\geq\liminf_{i\to\infty}F_{i}P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})$
$\displaystyle\geq F_{*}(t,y,\Delta w(s,y),\nabla
w(s,y),\partial_{t}w(s,y),\overline{v}(s,y)).$

Therefore $\overline{v}$ is a viscosity subsolution.

Step 2 ($\underline{v}$ is a supersolution). Arguments similar to those above
show that $\underline{v}$ is a viscosity supersolution, where the principal
change to the proof is that one considers $w\in
C^{\infty}(\mathbb{R}\times\mathbb{R}^{d})$ such that $\underline{v}-w$ has a
strict local minimum at some $(s,y)\in(0,T)\times{\Omega}$ with
$\underline{v}(s,y)=w(s,y)$. With analogous notation, the last line in (41)
corresponds to

$0\leq-d_{i}P_{i}w({s}_{i}^{k},{y}_{i}^{\ell})+\sup_{\alpha\in
A}\left({\mathsf{E}}_{i}^{\alpha}P_{i}w(s_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{\alpha}P_{i}w({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{\alpha}\right)_{\ell}+\gamma\left|\mu_{i}\right|,$

i.e. there is a slight asymmetry in the argument due to the last sign in (41).
Nevertheless, it is then deduced that

$0\leq F^{*}(t,x,\Delta w(t,x),\nabla
w(t,x),\partial_{t}w(t,x),\underline{v}(t,x)).$

Thus $\underline{v}$ is a viscosity supersolution. ∎

The above proof is an adaptation of the Barles-Souganidis argument [BS91] to
the finite element setting, in line with that in [JS13] but differing in the
treatment of the boundary conditions.

######  Assumption 5.

Let $\overline{v}$ be a lower semi-continuous supersolution and
$\underline{v}$ be an upper semi-continuous subsolution. Then
$\underline{v}\leq\overline{v}$.

######  Theorem 5.

One has $\underline{v}=\overline{v}=v$, where $v$ is the unique viscosity
solution with $v(T,\cdot)=v_{T}$. Furthermore

$\displaystyle\lim_{i\to\infty}\|v_{i}-v\|_{L^{\infty}((0,T)\times{\Omega})}=0.$
(42)

###### Proof.

Follows as in the proof of Theorem 6.2 in [JS13]. ∎

##  9 Numerical experiments

$10^{1}$$10^{2}$$10^{-2}$$10^{-1}$$\textrm{Inverse of mesh size}\ 1/\Delta
x$Error$L^{\infty}$$L^{2}$$H^{1}$ Figure 1: Approximation error of Experiment
1

The first experiment investigates rates of convergence for a known smooth
solution. The remaining experiments examine the approximation of solutions
with singularities near type changes of boundary conditions as well as the
solution behaviour in the vicinity of nonlinear boundary conditions. The code
is available from the public repository [Jar21] under the GNU Lesser General
Public License.

Experiment 1 (Rates for smooth known solution): We consider a final time
boundary value problem on the square domain ${\Omega}=[-1,1]^{2}$ with Robin
conditions on the right face $\partial{\Omega}_{t}=\\{1\\}\times(-1,1)$ and
Dirichlet conditions on the remaining faces
$\partial{\Omega}_{D}=\partial{\Omega}\setminus\partial{\Omega}_{t}$. We have
the control set $A=[0,1]$ and the final time $T=1$ for the system

$\displaystyle-\partial_{t}v+\sup_{\alpha\in A}\left(-(\alpha+|x|^{2}/2)\Delta
v+xv_{x}-f^{\alpha}\right)$ $\displaystyle=0$ $\displaystyle\quad\text{in
}[0,T)\times{\Omega},$ (43) $\displaystyle-\partial_{t}v+\sup_{\alpha\in
A}\bigl{(}\alpha v_{x}-g^{\alpha}\bigr{)}$ $\displaystyle=0$
$\displaystyle\quad\text{on }[0,T)\times\partial{\Omega}_{t},$ $\displaystyle
v$ $\displaystyle=0$ $\displaystyle\quad\text{on
}[0,T)\times\partial{\Omega}_{D},$ $\displaystyle v-(1-x^{2})(1-y^{2})$
$\displaystyle=0$ $\displaystyle\quad\text{on
}\\{T\\}\times\overline{{\Omega}}.$

We choose $g^{\alpha}$ and $f^{\alpha}$ such that

$v(x,y,t):=t(1-x^{2})(1-y^{2})+(1-t)\sin{(\pi x)}\cos{\left(\frac{\pi
y}{2}\right)}$

is the exact solution of (43).

The artificial diffusion coefficients are selected quasi-optimally, cf.
Section 4.3. The time dependent Robin boundary condition is treated fully
explicitly. The time step size is chosen to ensure monotonicity while
permitting a large time step, leading to $O(h_{i})=O(\Delta x_{i})$. The
$L^{2}$, $H^{1}$ and $L^{\infty}$ errors at time $t=0$, presented also in
Figure 1, obey in essence the same rates as those observed previously [JS13]
with Dirichlet conditions and $O(h_{i})=O(\Delta x_{i})$ scaling:

$\Delta x$ | $L^{2}$ | Rate | $L^{\infty}$ | Rate | $H^{1}$ | Rate
---|---|---|---|---|---|---
0.1165 |  1.186e-1 | 0.98 |  1.645e-1 | 0.92 |  5.089e-1 | 0.93
0.0583 |  6.044e-2 | 0.98 |  8.960e-2 | 0.95 |  2.743e-1 | 0.94
0.0291 |  3.072e-2 | 0.99 |  4.706e-2 | 0.97 |  1.457e-1 | 0.95
0.0146 |  1.558e-2 | 0.99 |  2.426e-2 | 0.98 |  7.696e-2 | 0.95
0.0073 |  7.894e-3 | 1.04 |  1.234e-2 | 1.03 |  4.058e-2 | 0.96
0.0036 |  3.806e-3 |  |  6.006e-3 |  |  2.120e-2 |

Figure 2: Value function of the Skorokhod problem for a coarse mesh size
$\Delta x\approx 0.12$

Experiment 2 (Skorokhod problem): The second numerical experiment is set on a
non-convex, less regular domain, which is depicted in Figure 2. The stochastic
controlled process is subject to a terminal cost of $10$ everywhere apart from
$\overline{\omega_{1}}$ where it is 0. There is no running cost. On
$\omega_{2}$ the particle is transported through a Skorokhod reflection
independently of the angle of incidence in the direction of the inner normal
vector. Ultimately, a particle can avoid penalisation only by reaching
$\overline{\omega_{1}}$ before the terminal time $T=1$. On $\Omega$ the
particle may only choose between an upwards drift and drift to the right:

$\displaystyle-\partial_{t}v+\sup_{\alpha}\left(-a^{\alpha}\Delta
v-b^{\alpha}\cdot\nabla v\right)$ $\displaystyle=0$
$\displaystyle\textstyle\text{in }[0,T)\times{\Omega},$ (44a) $\displaystyle-
b_{\partial{\Omega}}\cdot\nabla v$ $\displaystyle=0$
$\displaystyle\textstyle\text{on
}[0,T)\times\omega_{2}\cup\left\\{\left(\frac{1}{4},0\right)\right\\},$ (44b)
$\displaystyle v$ $\displaystyle=0$ $\displaystyle\textstyle\text{on
}[0,T)\times\overline{\omega_{1}},$ (44c) $\displaystyle v$ $\displaystyle=10$
$\displaystyle\textstyle\text{on
}[0,T)\times\omega_{0}\cup\overline{\omega_{3}}\cup\left\\{\left(\frac{1}{2},\frac{1}{4}\right)\right\\},$
(44d) $\displaystyle v$ $\displaystyle=v_{T}$ $\displaystyle\textstyle\text{on
}\\{T\\}\times\overline{{\Omega}},$ (44e)

where $a^{\alpha}=0.1(1-x_{2})\alpha$ and
$b^{\alpha}=(-2\alpha,2(\alpha-1))^{T}$ for $\alpha\in\\{0,1\\}$. Moreover,
$b_{\partial{\Omega}}=(1,-1)^{T}$ and

$\displaystyle
v_{T}(x)=\begin{cases}10\quad&x\in\overline{{\Omega}}\setminus\omega_{1},\\\
0\quad&x\in\omega_{1}.\end{cases}$ (45)

Hence when drifting to the right the particle is exposed to Brownian noise,
while the equation is degenerate when the upward drift is selected. The
numerical operators are given by

$\displaystyle({\mathsf{E}}_{{\Omega},i}^{\alpha}v)_{\ell}$
$\displaystyle:=\,\bar{\nu}_{{\Omega},i}^{\alpha,\ell}\langle\nabla
v,\nabla\hat{\phi}_{i}^{\ell}\rangle+\langle-b^{\alpha}\cdot\nabla
v,\hat{\phi}_{i}^{\ell}\rangle$ (46a)
$\displaystyle({\mathsf{I}}_{{\Omega},i}^{\alpha}v)_{\ell}$
$\displaystyle:=\,\max\left(a^{\alpha}-\bar{\nu}_{{\Omega},i}^{\alpha,\ell},0\right)\langle\nabla
v,\nabla\hat{\phi}_{i}^{\ell}\rangle.$ (46b)
$\displaystyle({\mathsf{E}}^{\alpha}_{\partial{\Omega},i}v)_{\ell}$
$\displaystyle:=0,$ (46c)
$\displaystyle({\mathsf{I}}^{\alpha}_{\partial{\Omega},i}v)_{\ell}$
$\displaystyle:=\frac{v(t,x_{1},x_{2})-v(t,x_{1}-\lambda,x_{2}+\lambda)}{\lambda}.$
(46d)

Figure 2 shows the approximation $v_{i}$ at time $t=0$ on a coarse mesh.
Notice how the introduction of a Skorokhod type boundary gives the particle
starting in the vicinity of $\omega_{2}$ a high probability of reaching the
penalty-free exit zone $\overline{\omega_{1}}$. The node at
$(\frac{1}{2},\frac{1}{4})$ already belongs to the Dirichlet boundary. We
observe that the penalty of $10$ in the boundary segment between $\omega_{1}$
and $\omega_{2}$ leads to a layer-like behaviour of the solution of only one
element thickness. The related numerical experiments on finer meshes depicted
in Figure 3 also exhibit this aspect of the numerical solution.

(a) Skorokhod problem

(b) Nonlinear boundary condition

Figure 3: Linear and nonlinear boundary conditions on $\omega_{2}$ with
$\Delta x=0.0035$.

Experiment 3 (Internal barrier and nonlinear boundary conditions): The third
experiment is an adaptation of the previous one to examine the effect of
nonlinear boundary conditions. We break the adaptation into two parts.

Part (a) (Internal barrier): The experiment is identical to the previous one
with exception of the drift terms on $\Omega$:

$b^{\alpha}(x)=\begin{cases}\bigl{(}-2\alpha,2(\alpha-1)\bigr{)}^{T}\quad&:|x-\frac{3}{8}|>\frac{1}{20},\\\
\phantom{-\alpha}\bigl{(}0,2(\alpha-1)\bigr{)}^{T}\quad&:\textrm{otherwise.}\end{cases}$

This means that the strip of all $x$ with $|x-\frac{3}{8}|\leq\frac{1}{20}$
acts as barrier in $\Omega$: Within this strip there is no process with drift
to the right. The only way a particle can cross the strip from the left to the
right in order to avoid penalisation is by adopting $\alpha=1$. In this case
the particle might cross the barrier by means of diffusion; however, there is
no drift term to aid the crossing.

The construction results in a value function which at first sight resembles a
piecewise constant function. It is close to 10 left of the barrier as the
particle is unlikely to reach the penalty-free exit zone $\omega_{1}$. It is
mostly close to 0 right of the barrier; however, reaches 10 at Dirichlet
boundary conditions on $\omega_{0}$ as already indicated in Experiment 2. At
the barrier there is an internal layer arising from the possibility of
crossing owing to diffusion.

This description of the value function is matched entirely by the numerical
computations with the scheme of this paper: see Figure 3(a), where the
internal layer as well as the boundary conditions right of the barrier are
well resolved.

Part (b) (Nonlinear boundary condition): Now the boundary condition on
$\omega_{2}$ is replaced by a nonlinear operator which corresponds to the
choice between the previously used Skorokhod reflection and instantaneous
transport along the boundary towards the right. In other words, (44b) is
replaced by

$\sup\\{-b^{0}_{\partial{\Omega}}\cdot\nabla
v,-b^{1}_{\partial{\Omega}}\cdot\nabla v\\}=0\quad\textstyle\text{on
}[0,T)\times\omega_{2}\cup\left\\{\left(\frac{1}{4},0\right)\right\\}.$

with $b^{0}_{\partial{\Omega}}=(1,-1)$ and $b^{1}_{\partial{\Omega}}=(-1,-1)$.

This modification has a striking impact on the behaviour of the controlled
system. Now an optimally controlled particle may drift onto the boundary
segment $\omega_{2}$ left of the barrier to be then transported on the
boundary past the barrier. Right of the barrier the control is changed to
$b^{0}_{\partial{\Omega}}$ in order to reflect the particle into $\Omega$ to
avoid the penalty at the point $(\frac{1}{2},\frac{1}{4})$ at the end of
$\omega_{2}$.

An approximation of the resulting value function at time $0$ is shown in
Figure 3(b), which illustrates how this time regions left of the barrier have
a value function close to zero since particles located there can now reach the
penalty-free exit zone $\omega_{1}$. Indeed, when computing the solutions for
earlier times $t<0$, one observes further growth of the blue region as more
time is available to arrive at $\omega_{1}$ before termination.

Experiment 4 (Reflection vs. termination): We consider a final time boundary
value problem on the same domain, but now with a nonlinear boundary condition
corresponding to a choice between a Skorokhod reflection and termination of
the process in exchange for an oscillatory cost $g^{\alpha}$:

$\displaystyle-\partial_{t}v+\sup_{\alpha}\left(-a^{\alpha}\Delta
v-b^{\alpha}\cdot\nabla v\right)$ $\displaystyle=0$
$\displaystyle\quad\text{in }[0,T)\times{\Omega},$
$\displaystyle\sup_{\alpha}\bigl{(}-b_{\partial{\Omega}}^{\alpha}\cdot\nabla
v+c^{\alpha}_{\partial{\Omega}}v-g^{\alpha}\bigr{)}$ $\displaystyle=0$
$\displaystyle\quad\textstyle\text{on
}[0,T)\times\omega_{3}\cup\left\\{\left(\frac{1}{4},1\right)\right\\},$
$\displaystyle v$ $\displaystyle=0$ $\displaystyle\quad\text{on
}[0,T)\times\overline{\omega_{1}},$ $\displaystyle v$ $\displaystyle=10$
$\displaystyle\textstyle\quad\text{on
}[0,T)\times\omega_{0}\cup\overline{\omega_{2}}\cup\left\\{\left(\frac{1}{2},\frac{3}{4}\right)\right\\},$
$\displaystyle v$ $\displaystyle=v_{T}$ $\displaystyle\quad\text{on
}\\{T\\}\times\overline{{\Omega}},$

where $\alpha\in\\{0,1\\}$, $v_{T}$ as in (45) and

$\displaystyle a^{\alpha}$
$\displaystyle=0.2(1-x_{2})(1-\alpha)+0.2(1-x_{1})\alpha,$ $\displaystyle
b^{\alpha}$ $\displaystyle=(-2\alpha,2(\alpha-1))^{T},\quad
b^{\alpha}_{\partial{\Omega}}=(\alpha,\alpha)^{T},$ $\displaystyle c^{\alpha}$
$\displaystyle=(1-\alpha),\quad
g^{\alpha}=-(10\cos(160x_{1}/\pi+4)(1-\alpha).$

Note that compared to the two previous examples the Robin type boundary has
moved to $\omega_{3}$ while $\omega_{2}$ is part of the Dirichlet region with
value $10$. Both operators $L^{\alpha}$ now have regions of degeneracy, when
either $x_{1}$ or $x_{2}$ is near $1$.

The PDE operator on $\Omega$ is discretised according to (46a)–(46b), while
the Robin operators are approximated implicitly, consistent with Assumption 1.
The behaviour of the value function $v$ in the vicinity of $\omega_{3}$ is
depicted in Figure 4. One observes how troughs of $g^{0}$ are attained by the
value function, while near peaks of $g^{0}$ the numerical scheme switches to
the reflection principle. Overall the experiment demonstrates how the
framework of the paper not only allows us to approximate nonlinear Robin
conditions, but also incorporates a nonlinear switching between Robin
conditions on the one hand and Dirichlet conditions on the other hand.

Figure 4: Value function with nonlinear boundary condition on $\omega_{3}$.

## Acknowledgements

Bartosz Jaroszkowski gratefully acknowledges the support of the EPSRC grant
1816514. Max Jensen gratefully acknowledges the support of the Dr Perry James
Browne Research Centre.

## References

  * [AF12] Yves Achdou and Maurizio Falcone.  A semi-Lagrangian scheme for mean curvature motion with nonlinear Neumann conditions.  Interfaces Free Bound., 14(4):455–485, 2012.
  * [BE02] Erik Burman and Alexandre Ern.  Nonlinear diffusion and discrete maximum principle for stabilized Galerkin approximations of the convection–diffusion-reaction equation.  Computer Methods in Applied Mechanics and Engineering, 191(35):3833–3855, 2002.
  * [BFO12] Jean-David Benamou, Brittany D Froese, and Adam M Oberman.  A viscosity solution approach to the Monge-Ampere formulation of the Optimal Transportation Problem.  arXiv, August 2012.
  * [BMZ09] Olivier Bokanowski, Stefania Maroso, and Hasnaa Zidani.  Some convergence results for Howard’s algorithm.  SIAM Journal on Numerical Analysis, 47(4):3001–3026, 2009.
  * [BP94] Abraham Berman and Robert J. Plemmons.  Nonnegative matrices in the mathematical sciences, volume 9 of Classics in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1994.  Revised reprint of the 1979 original.
  * [BS91] Guy Barles and Panagiotis E. Souganidis.  Convergence of approximation schemes for fully nonlinear second order equations.  Asymptotic Analysis, 4(3):271–283, 1991.
  * [CIL92] Michael G. Crandall, Hitoshi Ishii, and Pierre-Louis Lions.  User’s guide to viscosity solutions of second order partial differential equations.  Bulletin of the American Mathematical Society, 27(1):1–67, 1992\.
  * [FGN13] Xiaobing Feng, Roland Glowinski, and Michael Neilan.  Recent developments in numerical methods for fully nonlinear second order partial differential equations.  SIAM Rev., 55(2):205–267, 2013.
  * [Gal19] Dietmar Gallistl.  Numerical approximation of planar oblique derivative problems in nondivergence form.  Math. Comp., 88(317):1091–1119, 2019.
  * [Jar21] Bartosz Jaroszkowski.  FEISol (2021).  https://github.com/BartoszJaroszkowski/FEISol, 2021.
  * [Jen17] Max Jensen.  $L^{2}(H^{1}_{\gamma})$ finite element convergence for degenerate isotropic Hamilton–Jacobi–Bellman equations.  IMA Journal of Numerical Analysis, 37(3):1300–1316, 2017.
  * [JJ21] Bartosz Jaroszkowski and Max Jensen.  Valuation of European options under an uncertain market price of volatility risk.  submitted, preprint on arXiv, 2021.
  * [JS13] Max Jensen and Iain Smears.  On the convergence of finite element methods for Hamilton-Jacobi-Bellman equations.  SIAM Journal on Numerical Analysis, 51(1):137–162, 2013.
  * [Kaw19] Ellya L. Kawecki.  A discontinuous Galerkin finite element method for uniformly elliptic two dimensional oblique boundary-value problems.  SIAM J. Numer. Anal., 57(2):751–778, 2019.
  * [KD01] Harold J. Kushner and Paul G. Dupuis.  Numerical methods for stochastic control problems in continuous time.  In Applications of Mathematics, volume 24. Springer-Verlag, New York, 2 edition, 2001.
  * [Lio85] P L Lions.  Neumann type boundary conditions for Hamilton-Jacobi equations.  Duke Mathematical Journal, 52(4):793–820, 1985.
  * [NSZ17] Michael Neilan, Abner J. Salgado, and Wujun Zhang.  Numerical analysis of strongly nonlinear pdes.  Acta numer., 26:137–303, 2017.
  * [OR73] O. A. Oleĭnik and E. V. Radkevič.  Second order equations with nonnegative characteristic form.  Plenum Press, New York, London, 1973.
  * [Ser03] Oana-Silvia Serea.  On reflecting boundary problem for optimal control.  SIAM Journal on Control and Optimization, 42(2):559–575, 2003.

##  Appendix A Appendix: Interpretation of mixed boundary conditions

We briefly sketch in a simplified setting how mixed boundary conditions can
arise from an underlying optimal control problem. For a start we assume here
that the solution $v$ of (3) is smooth. In optimal control formulations the
coefficients $c^{\alpha}$ and $c_{\partial{\Omega}}^{\alpha}$ are typically
either $0$ or they all coincide with some constant in order to model a
discounting of cost. We shall assume the former.

We consider a particle, or agent, which occupies the state
$\mathbf{x}(t)\in\overline{{\Omega}}$ at time $t\in[0,T]$. Its movements are
described by the following rules:

  1. 1.

Suppose $\mathbf{x}(t)\in\Omega$, $t<T$ and the control $\alpha\in A$ is
selected. Then the particle’s immanent movement is described by the SDE

$\displaystyle\mathrm{d}\mathbf{x}=b^{\alpha}(\mathbf{x})\,\mathrm{d}t+\sqrt{2\,a^{\alpha}}\,\mathrm{d}W,$
(47)

where $W$ is a $d$-dimensional Brownian motion. While the particle follows
(47) it is subject to the cost $f^{\alpha}\mathrm{d}t$.

  2. 2.

Suppose $\mathbf{x}(t)\in\partial{\Omega}_{D}^{\circ}$, $t<T$ and the control
$\alpha\in A$ is selected. Here $\partial{\Omega}_{D}^{\circ}$ refers to the
interior of $\partial{\Omega}_{D}$ relative to $\partial{\Omega}$. Then, in
the context of viscosity boundary conditions, the particle may either follow
(47) at the running cost $f^{\alpha}\mathrm{d}t$ or terminate its movement at
a cost of $g(\mathbf{x}(t))$. If instead pointwise Dirichlet conditions were
imposed in Definition 1 then the boundary conditions would correspond to the
guarantee that the particle terminates its movement.

  3. 3.

Suppose $\mathbf{x}(t)\in\partial{\Omega}_{R}^{\circ}$, $t<T$ and the control
$\alpha\in A$ is selected. Then the Skorokhod reflection principle may apply.
Indeed, as described in [KD01, Sections 1.4, 3.1.3, …], upon reaching the
boundary the particle may instantaneously be transported a distance
$b_{\partial{\Omega}}^{\alpha}\,\delta$ where $\delta>0$ is small.
Alternatively, because of the context of viscosity boundary conditions, the
particle may continue to follow (47). To remain within the scope of [KD01] we
assume $g^{\alpha}=0$ on $\partial{\Omega}_{R}$.

  4. 4.

Suppose $\mathbf{x}(t)\in\partial{\Omega}_{t}^{\circ}$, $t<T$ and the control
$\alpha\in A$ is selected. Then the particle may either move according to

$\displaystyle\mathrm{d}\mathbf{x}=b_{\partial{\Omega}}^{\alpha}(\mathbf{x})\mathrm{d}t$
(48)

with running cost $g^{\alpha}\mathrm{d}t$ or according to (47) with running
cost $f^{\alpha}\mathrm{d}t$.

  5. 5.

Suppose $t=T$ and the particle’s movement has not yet terminated at a
Dirichlet boundary then the final time cost $v_{T}(\mathbf{x}(T))$ incurs.

  6. 6.

Suppose that
$\mathbf{x}(t)\in\partial(\partial{\Omega}_{D})\cup\partial(\partial{\Omega}_{R})\cup\partial(\partial{\Omega}_{t})$,
$t<T$. Here the outer $\partial$ of $\partial(\partial{\Omega}_{X})$ refers to
the boundary of $\partial{\Omega}_{X}$ relative to $\partial{\Omega}$ where
$X\in\\{D,R,t\\}$. Then the particle’s behaviour may be selected from multiple
of the above scenarios. E.g. if
$\mathbf{x}(t)\in\partial(\partial{\Omega}_{D})\cup\partial(\partial{\Omega}_{R})$
then the particle movement may terminate, the particle may be reflected or it
may be transported according to (47).

All these possible scenarios occur when representing uncertain market price of
volatility risk in a Heston model found in [JJ21].

We now link the above description of the particle by means of the value
function to the HJB final time boundary value problem (3). Let $:[0,T]\to A$
represent a choice of controls for each time $s\in[0,T]$. Similarly, let
$\xi_{\Omega},\xi_{\partial{\Omega}_{D}},\xi_{\partial{\Omega}_{R}},\xi_{\partial{\Omega}_{t}}:[0,T]\times\overline{{\Omega}}\to\\{0,1\\}$
be indicator functions such that
$\operatorname{supp}\xi_{\partial{\Omega}_{X}}\subset\overline{\partial{\Omega}_{X}}$
for $X\in\\{D,R,t\\}$, where $\overline{\partial{\Omega}_{X}}$ is the closure
of $\partial{\Omega}_{X}$ relative to $\partial{\Omega}$. Furthermore,

$\xi_{\Omega}+\xi_{\partial{\Omega}_{D}}+\xi_{\partial{\Omega}_{R}}+\xi_{\partial{\Omega}_{t}}\equiv
1.$

Where $\xi_{\Omega}=1$ the particle path $\mathbf{x}$ obeys (47), where
$\xi_{\partial{\Omega}_{D}}=1$ the particle terminates, where
$\xi_{\partial{\Omega}_{R}}=1$ the particle is reflected and where
$\xi_{\partial{\Omega}_{t}}=1$ the particle follows (48). Since the particle
terminates where $\xi_{\partial{\Omega}_{D}}=1$ we requite
$\xi_{\partial{\Omega}_{D}}(s_{1},x)=1\Rightarrow\xi_{\partial{\Omega}_{D}}(s_{2},x)=1$
for $s_{1}\leq s_{2}$. The value function $v$ at $(t,x)$ is the smallest cost
realised among all possible choices for and
${\boldsymbol{\xi}}=(\xi_{\Omega},\xi_{\partial{\Omega}_{D}},\xi_{\partial{\Omega}_{R}},\xi_{\partial{\Omega}_{t}})$:

$\displaystyle v(t,x)=\inf_{{\boldsymbol{\xi}}}{\bf E}_{xt}$
$\displaystyle\Bigl{(}\int_{t}^{\tau}\xi_{\Omega}(\mathbf{x}(s))\,f^{(s)}(\mathbf{x}(s))+\xi_{\partial{\Omega}_{t}}(\mathbf{x}(s))\,g^{(s)}(\mathbf{x}(s))\mathrm{d}s$
$\displaystyle\,+\xi_{\partial{\Omega}_{D}}(\mathbf{x}(\tau))\,g(\mathbf{x}(\tau))+\xi_{\Omega}(\mathbf{x}(\tau))\,v_{T}(\mathbf{x}(\tau))\Bigr{)},$

where $\tau$ is the exit time from $[0,T)\times\Omega$ of $\mathbf{x}$ and
${\bf E}_{xt}$ the expectation conditional to $\mathbf{x}(t)=x$.

We fix some ${\boldsymbol{\xi}}$ which are not necessarily optimal. Suppose
that $\mathbf{x}(s)\in\operatorname{supp}\xi_{\partial{\Omega}_{t}}$ for a
short duration $[t,t+\varepsilon)\ni s$. Then

$\displaystyle g^{(t)}(\mathbf{x}(t))$ $\displaystyle=\lim_{h\to
0}\frac{1}{h}\int_{t}^{t+h}g^{(s)}(\mathbf{x}(s))\mathrm{d}s\geq-\lim_{h\to
0}\frac{v(t+h,\mathbf{x}(t+h))-v(t,x)}{h}$
$\displaystyle=-\partial_{t}v(t,\mathbf{x})-\nabla
v(t,\mathbf{x})\cdot\dot{\mathbf{x}}=-\partial_{t}v(t,\mathbf{x})-\nabla
v(t,\mathbf{x})\cdot b_{\partial{\Omega}}^{\alpha}(\mathbf{x}).$

Here the first equality follows from continuity, the inequality from the
dynamic programming principle, the second equality from the chain rule and the
third from (48).

Now, suppose that $\mathbf{x}(s)\subset\operatorname{supp}\xi_{\Omega}$ for a
short duration $[t,t+\varepsilon)\ni s$. Then by a similar argument, detailed
in [KD01], one finds

$\displaystyle f^{(t)}(\mathbf{x}(t))\geq-\partial_{t}v(t,\mathbf{x})-\nabla
v(t,\mathbf{x})\cdot b^{\alpha}(\mathbf{x})-a^{\alpha}\nabla v(t,\mathbf{x}).$

On $\operatorname{supp}\xi_{\partial{\Omega}_{D}}$ we find
$g(\mathbf{x}(t))\geq v(t,\mathbf{x}(t))$ as the minimal cost cannot be more
than the cost of termination.

Suppose now that the particle is located at
$x\in\operatorname{supp}\xi_{\partial{\Omega}_{R}}$. It cannot be more
beneficial for the particle to be at
$x+b_{\partial{\Omega}}^{\alpha}\,\lambda$ as it will immanently be
transported there. Thus $v(t,x)\leq
v(t,x+b_{\partial{\Omega}}^{\alpha}\,\lambda)$ and therefore, with $\lambda\to
0$,

$-b_{\partial{\Omega}}^{\alpha}(x)\cdot\nabla v(t,x)\leq 0.$

When and where-ever the choice of ${\boldsymbol{\xi}}$ is optimal, the
respective above inequality turns into an equality. With the compactness of
$A$ and the continuous dependence of the coefficients on $\alpha$, such
optimal controls exist. Therefore, taking suprema over $A$, one obtains that
the value function solves (3), at least conceptually, with boundary conditions
in the viscosity sense. We refer here to the viscosity sense because the use
of semi-continuous envelopes in Definition 1 is interpreted as permitting (47)
as transport law on all of the closure $\overline{{\Omega}}$ and as offering
at the interfaces between boundary regions multiple boundary operators for the
choice of the optimal strategy, like indicated in scenario 6 of the above
list. We note that the choice between (47) and the various boundary operators
will in general be subject to some delicate restrictions, arising from the
sub- and superjets. At the boundary these jets are increased in size compared
to their counterparts in the domain interior [CIL92, Remark 2.7].