# Finite Element Methods for Isotropic Isaacs Equations with Viscosity and
Strong Dirichlet Boundary Conditions††thanks: Bartosz Jaroszkowski
<EMAIL_ADDRESS>acknowledges the support of the EPSRC grant
1816514. Max Jensen<EMAIL_ADDRESS>acknowledges the support of the Dr
Perry James Browne Research Centre.

Bartosz Jaroszkowski  Max Jensen

(Department of Mathematics, University of Sussex, Brighton, UK)

###### Abstract

We study monotone P1 finite element methods on unstructured meshes for fully
non-linear, degenerately parabolic Isaacs equations with isotropic diffusions
arising from stochastic game theory and optimal control and show uniform
convergence to the viscosity solution. Elliptic projections are used to manage
singular behaviour at the boundary and to treat a violation of the consistency
conditions from the framework by Barles and Souganidis by the numerical
operators. Boundary conditions may be imposed in the viscosity or in the
strong sense, or in a combination thereof. The presented monotone numerical
method has well-posed finite dimensional systems, which can be solved
efficiently with Howard’s method.

##  1 Introduction

Hamilton–Jacobi–Isaacs equations, in short Isaacs equations, arise from
stochastic game theory and optimal control, in particular from stochastic two
player zero-sum games [Sou99]. The solution of the Isaacs equation corresponds
under appropriate assumptions to a value function of the game. The equations
are of the structure

$\displaystyle-\partial_{t}v+\adjustlimits{\inf}_{\beta}{\sup}_{\alpha}(L^{(\alpha,\beta)}v-f^{(\alpha,\beta)})=0,$
(1)

where $\inf\sup$ is taken over a family of second-order linear degenerately
elliptic differential operators $L^{\alpha,\beta}$.

Isaacs equations combine a number of features which are challenging from the
numerical point of view. They are fully nonlinear in non-divergence form and
the nonlinearity is in general non-convex. Their spatial differential
operators may only be degenerately elliptic and often do not exhibit
smoothness and structure properties to make linearisation a usable tool.
Furthermore, the theory of well-posedness relies, in the setting of viscosity
solutions, on comparison principles, which in turn make it natural to impose
monotonicity requirements on the numerical methods. Boundary conditions, such
as of Dirichlet type, can be posed in multiple non-equivalent forms and
comparison principles may rely on a particular form to ensure that they hold.

In this work we present a finite element method to solve time-dependent Isaacs
equations with isotropic, possibly degenerate diffusions. We show the
convergence to the viscosity solution on Lipschitz domains. Dirichlet boundary
conditions are imposed in the viscosity or in the strong sense, or in a
combination thereof. The presented analysis is a generalisation of [JS13],
where Bellman equations with homogeneous, strong Dirichlet boundary conditions
were approximated.

Compared to the literature for typical quasilinear differential equations, the
scientific discourse on the numerical approximation of second-order Isaacs
equations consists of a small number of papers. A two-scale finite element
method for Isaacs equations was proposed in [SZ19], where
Alexandrov–Bakelman–Pucci estimates are used to prove the convergence with
rates of numerical approximations to solutions of uniformly elliptic Isaacs
equations on convex domains. Such two-scale finite element methods are closely
related to semi-Lagrangian schemes [LN18], for which convergence rates for
Isaacs operators on unbounded domains were given in [DJ13], see also the
series of works [Jak04], [Jak06], [JK02], [JK05] as well as [MZB06] in this
context. Finite difference methods for uniformly parabolic Isaacs equations on
smooth bounded domains with Dirichlet boundary conditions were shown to
converge with rates in [Kry15] and [Tur15]. Deterministic two-player zero-sum
games lead to first-order Isaacs equations. We refer for the numerical
approximation of deterministic problems to [BCJ19], [BHT11], [Fal06], [FK14]
and the references therein. An important field of application is $H^{\infty}$
control, for which we highlight [BM98], [FRS10], [KKK20], [Sor99].

The main contributions of the paper are the following:

Convergence: We present a finite element method which is guaranteed to
converge to the viscosity solution of the Isaacs equation, including in the
degenerate case. Domains only need to be Lipschitz regular.

Boundary conditions: Dirichlet boundary conditions are treated in the
viscosity and strong sense, so that the numerical analysis can follow the
requirements arising from the construction of comparison principles.

Barrier functions: Barrier functions are used to ensure that upper and lower
semi-continuous envelopes of the numerical solutions locally satisfy the
boundary conditions in the strong sense. The construction of barrier functions
is characterized in the degenerate setting and on non-convex domains.

Stability: Boundary data is mapped onto the approximation space with elliptic
projections in order to avoid instabilities, which may arise for instance in
the presence of re-entrant corners of the domain.

This paper is organised as follows. In Section 2 we specify the class of
problems under consideration, introduce a notion of viscosity solution and
discuss the two types of Dirichlet boundary conditions. In Section 3 we define
the numerical scheme. In Section 4 we establish monotonicity properties; in
Section 5 we show the existence and uniqueness of the numerical solutions; in
Section 6 we ensure stability; in Section 7 we verify that the upper and lower
semi-continuous envelopes of numerical solutions are sub- and supersolutions,
respectively. In Section 8 we return to the different notions of the boundary
conditions. In particular, we introduce barrier functions which ensure that
boundary conditions can be satisfied in the strong sense at least on a part of
the boundary. With that we can state the convergence result in Section 9. In
Section 10 we show the construction of barrier functions in uniformly
parabolic as well as in degenerately parabolic cases. Finally, in Section 11
we conclude with numerical experiments verifying the convergence of the scheme
and we show an application to a two-player stochastic game.

##  2 Isotropic Isaacs equations

We consider Isaacs equations on bounded Lipschitz domains
${\Omega}\in\mathds{R}^{d}$ with $d\geq 2$. Let $A$ and $B$ be compact metric
spaces and $(\alpha,\beta)\in A\times B$. We introduce the linear operators

$L^{(\alpha,\beta)}:\;H^{2}({\Omega})\to
L^{2}({\Omega}),\;w\mapsto-a^{(\alpha,\beta)}\,\Delta
w-b^{(\alpha,\beta)}\cdot\nabla w+c^{(\alpha,\beta)}\,w$ (2)

and data terms $f^{(\alpha,\beta)},v_{T}\in C(\overline{{\Omega}})$,
$g\in\\{w\in W^{1,\infty}(\mathds{R}^{d}):\Delta w\in
L^{\infty}({\Omega})\\}$. The mapping

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

is assumed to be continuous such that the families of functions
$\\{a^{(\alpha,\beta)}\\}_{(\alpha,\beta)\in A\times B}$,
$\\{b^{(\alpha,\beta)}\\}_{(\alpha,\beta)\in A\times B}$,
$\\{c^{(\alpha,\beta)}\\}_{(\alpha,\beta)\in A\times B}$ and
$\\{f^{(\alpha,\beta)}\\}_{(\alpha,\beta)\in A\times B}$ are equicontinuous.
Moreover, $a^{(\alpha,\beta)}(x)\geq 0$ for any $(\alpha,\beta)\in A\times B$
and $x\in{\Omega}$. It follows that all $L^{\alpha,\beta}$ are degenerate
elliptic and that

$\sup_{(\alpha,\beta)\in A\times
B}\|\,(a^{(\alpha,\beta)},b^{(\alpha,\beta)},c^{(\alpha,\beta)},f^{(\alpha,\beta)})\,\|_{C(\overline{{\Omega}})\times
C(\overline{{\Omega}},\mathds{R}^{d})\times C(\overline{{\Omega}})\times
C(\overline{{\Omega}})}<\infty.$

For any sufficiently smooth $w$ we define the Hamiltonian

$Hw:=\adjustlimits{\inf}_{\beta\in B}{\sup}_{\alpha\in
A}(L^{(\alpha,\beta)}w-f^{(\alpha,\beta)}),$

assuming that the supremum and the infimum are applied pointwise. We wish to
study the Isaacs problems of the following form:

$\displaystyle-\partial_{t}v+Hv$ $\displaystyle=0$
$\displaystyle\qquad\text{in }(0,T)\times{\Omega},$ (3a) $\displaystyle v$
$\displaystyle=g$ $\displaystyle\qquad\text{on }(0,T)\times\partial{\Omega},$
(3b) $\displaystyle v$ $\displaystyle=v_{T}$ $\displaystyle\qquad\text{on
}\\{T\\}\times\overline{{\Omega}},$ (3c)

where $v_{T}\in L^{\infty}({\Omega})$ and
$v_{T}|_{\partial{\Omega}}=g|_{\partial{\Omega}}$. Even though only the values
of $g|_{\partial{\Omega}}$ appear in (3) and
$v_{T}|_{\partial{\Omega}}=g|_{\partial{\Omega}}$ it remains valuable to
distinguish $g$ and $v_{T}$: While $v_{T}$ is only assumed to be uniformly
bounded, we require that $g$ can be written as the trace of a function from
the more regular space $\\{w\in W^{1,\infty}(\mathds{R}^{d}):\Delta w\in
L^{\infty}({\Omega})\\}$.

Typically there does not exist a smooth solution allowing a classical
interpretation of (3). Our aim is to construct a finite element method which
approximates the viscosity solution of the Isaacs problem (3). We now
formalise what is meant by a viscosity solution throughout this chapter.
Firstly, let us reformulate (3) as $F=0$ where $F$ is the differential
operator

$\displaystyle F(t,x,q,p,r,s)=$
$\displaystyle\left\\{\begin{array}[]{rllll}\\!\\!\\!-r+{\displaystyle\inf_{\beta}\sup_{\alpha}}\bigl{(}-a^{(\alpha,\beta)}(x)\,q-b^{(\alpha,\beta)}(x)\\!\cdot\\!p+c^{(\alpha,\beta)}(x)\,s-f^{(\alpha,\beta)}(x)\bigr{)}\\!&:x\in{\Omega},t<T,\\\
s-g(x)\\!&:x\in\partial{\Omega},t<T,\\\\[5.69054pt]
s-v_{T}(x)\\!&:t=T.\end{array}\right.$

Given a bounded function $v:[0,T]\times\overline{{\Omega}}\to\mathds{R}$ we
define its upper and lower semi-continuous envelopes, 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 2.2 below imposes Dirichlet boundary conditions in the viscosity
sense on all of $\partial{\Omega}$ in the flavour of the original Barles-
Souganidis proof, see also [BS91, JS18]. For boundary value problems arising
from concrete applications, such viscosity boundary conditions may not be
strong enough to establish a comparison principle. We therefore introduce a
subset $\omega\subset\partial{\Omega}$ on which boundary conditions hold in a
stronger sense.

######  Definition 2.1.

Let $g:\partial{\Omega}\to\mathds{R}$ and
$v_{T}:\overline{{\Omega}}\to\mathds{R}$ be bounded functions. We say that a
bounded function $v:[0,T]\times\overline{{\Omega}}\to\mathds{R}$ satisfies the
Dirichlet boundary conditions $v(t,x)=g(x)$ at $x\in\partial{\Omega}$ strongly
if, for all $t\in[0,T]$,

$v^{*}(t,x)=g(x)\qquad\text{and}\qquad v_{*}(t,x)=g(x).$

Similarly, we say that $v$ satisfies the final time conditions at
$x\in\overline{{\Omega}}$ strongly if

$v^{*}(T,x)=v_{T}(x)\qquad\text{and}\qquad v_{*}(T,x)=v_{T}(x).$

Depending on the setting $\omega=\emptyset$, $\omega=\partial{\Omega}$ or
$\emptyset\subsetneqq\omega\subsetneqq\partial{\Omega}$ may be most
appropriate choices.

######  Definition 2.2.

A bounded function $v$ is a viscosity supersolution (respectively,
subsolution) of (3) strongly on $\omega\subset\partial{\Omega}$ if, for any
test function $\psi\in C^{2}(\mathds{R}\times\mathds{R}^{d})$,

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

(respectively,

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

provided that $v_{*}-\psi$ attains at $(t,x)\in[0,T)\times\overline{{\Omega}}$
a local minimum relative to $[0,T)\times\overline{{\Omega}}$ (respectively,
$v^{*}-\psi$ attains a local maximum) and additionally $v_{*}$ (respectively
$v^{*}$) satisfies the final time conditions and the boundary conditions on
$\omega$ in the strong sense. A function which is simultaneously a viscosity
super- and subsolution of (3) is called a viscosity solution.

The motivation for incorporating strong boundary conditions in the definition
of viscosity solutions can be illustrated with a simple example.

######  Example 2.1.

We investigate on ${\Omega}=(-1,1)$ the Bellman equation

$\displaystyle\sup_{\beta\in\mathds{R}}\beta\cdot
v^{\prime}(x)-\beta^{2}/4-(x^{2}-x^{4})^{2}=0$ (6)

with homogeneous Dirichlet boundary conditions. Identifying the structure of a
Fenchel conjugate in the $\beta$-dependent terms, we can rewrite the PDE as
$|v^{\prime}(x)|^{2}-(x^{2}-x^{4})^{2}=0$. Taking the square root, we may
simplify the expression to the eikonal equation $|v^{\prime}(x)|=x^{2}-x^{4}$
with a double-well potential on the right-hand side.

An elementary calculation shows that

$v(x)=\frac{|x|^{5}}{5}-\frac{|x|^{3}}{3}+\frac{2}{15}$

is twice continuously differentiable, solves (6) classically and satisfies the
boundary conditions strongly.

For $c\leq 0$ we define

$v_{c}:\;\overline{{\Omega}},\;x\mapsto\begin{cases}v(x)\quad&:x\in{\Omega},\\\
c\quad&:x\in\partial{\Omega}.\end{cases}$

We initially assume that $\omega$ of Definition 2.2 is equal to the empty set
and show that $v_{c}$ is in that case a viscosity solution for all $c\leq 0$.
In particular the viscosity solution is then not unique as required by the
construction in [BS91] and there cannot exists a corresponding comparison
principle.

We observe that $(v_{c})^{*}=v$ on $\overline{{\Omega}}$, whereas
$(v_{c})_{*}=v_{c}$ since $v_{c}$ is lower semi-continuous. Thus it is clear
that $v_{c}$ is a viscosity subsolution of the problem.

We now prove that $v_{c}$ is also a supersolution. We must show that for all
$\psi\in C^{2}$ such that $(v_{c})_{*}-\psi$ has a local minimum $x$ we have
as in [CIL92, equation (7.10)]

$\displaystyle\textstyle\sup_{\beta}\beta\cdot\psi^{\prime}(x)-\beta^{2}/4-(x^{2}-x^{4})^{2}\phantom{\\}}$
$\displaystyle\geq 0,$ $\displaystyle x\in{\Omega},$ (7a)
$\displaystyle\max\\{(v_{c})_{*},\textstyle\sup_{\beta}\beta\cdot\psi^{\prime}(x)-\beta^{2}/4-(x^{2}-x^{4})^{2}\\}$
$\displaystyle\geq 0,$ $\displaystyle x\in\partial{\Omega}.$ (7b)

Inequality (7a) holds because $v_{c}$ is a classical solution on ${\Omega}$.
Remembering that $x^{2}-x^{4}$ vanishes on $\partial{\Omega}$ and the
simplification to the eikonal form, the boundary case (7b) can be simplified
to

$\displaystyle\max\\{c,|\psi^{\prime}(x)|^{2}\\}\geq 0,\qquad
x\in\partial{\Omega},$

which holds trivially.

Instead, setting $\omega=\partial{\Omega}$ we enforce $(v_{c})_{*}=v_{c}=0$ on
the boundary, which filters out all the spurious solutions with $c<0$ and
uniqueness is regained.

Similar examples can also be constructed with second-order equations arising
from diffusive problems and in higher dimensions, e.g. see [JS18, Section 2].
While the loss of uniqueness of solutions directly interacts with the ability
to formulate comparison principles, the flexibility of viscosity solutions to
break continuity in the vicinity of the boundary can be helpful in some
circumstances, e.g. when examining singularly perturbed diffusion-advection
equations which exhibit a boundary layers for small diffusion coefficients
$\varepsilon>0$ and may fully detach from the Dirichlet data when
$\varepsilon\to 0$.

We conclude this section by relating our definition of strong boundary
conditions to that in [CIL92, Section 7.A]. There an upper semi-continuous
function $v:\overline{{\Omega}}\to\mathds{R}$ is a subsolution of the boundary
conditions $B(x,v(x),Du(x),D^{2}u(x))=0$ at $x\in\partial{\Omega}$ if

$\displaystyle B(x,v(x),p,X)\leq
0\quad\text{for}\quad(p,X)\in\overline{J}_{\overline{{\Omega}}}^{2,+}v(x).$
(8)

Similarly, a lower semi-continuous $v:\overline{{\Omega}}\to\mathds{R}$ is a
subsolution of the boundary conditions at $x\in\partial{\Omega}$ if

$\displaystyle B(x,v(x),p,X)\geq
0\quad\text{for}\quad(p,X)\in\overline{J}_{\overline{{\Omega}}}^{2,-}v(x).$
(9)

In our case $B(x,q,p,s)=s-g(x)$ and we refer for the definition of the closed
second-order jets $\overline{J}_{\overline{{\Omega}}}^{2,+}u(x)$,
$\overline{J}_{\overline{{\Omega}}}^{2,-}u(x)$ to [CIL92].

Comparing (8) and (9) to Definitions 2.1 and 2.2, we first of all make an
observation which is not directly related to the boundary conditions: we adopt
in this paper the notion of bounded, but not necessarily semi-continuous sub-
and supersolutions, following [BP88] and [BS91], see also [FS06, Chapter VII].
Hence, to bridge this difference, we can adapt the definition of [CIL92] to

$\displaystyle B(x,v^{*}(x),p,X)\leq
0\quad\text{for}\quad(p,X)\in\overline{J}_{\overline{{\Omega}}}^{2,+}v(x)$
(10)

and

$\displaystyle B(x,v_{*}(x),p,X)\geq
0\quad\text{for}\quad(p,X)\in\overline{J}_{\overline{{\Omega}}}^{2,-}v(x).$
(11)

Now if the closed second-order jets of $v$ are non-empty, then (10) implies
$v^{*}(x)\leq g(x)$ and (11) implies $v_{*}(x)\geq g(x)$. Since $v^{*}\geq
v_{*}$ by construction, we arrive at $v^{*}(x)=v_{*}(x)=g(x)$, mirroring
Definition 2.1.

##  3 The numerical method

We consider a sequence $V_{i},i\in\mathds{N}$, of piecewise linear shape-
regular finite element spaces. Let us denote the nodes of the finite element
mesh by ${y}_{i}^{\ell}$ where $\ell$ is the index over the set of nodes. The
associated hat functions are denoted ${\phi}_{i}^{\ell}$, i.e.
${\phi}_{i}^{\ell}\in V_{i}$ such that $\phi_{i}^{\ell}(y_{i}^{\ell})=1$ while
${\phi}_{i}^{\ell}(y_{i}^{s})=0$ for all $\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}$ norm
whilst the $\hat{\phi}_{i}^{\ell}$ are normalised in the $L^{1}$ norm.

Let $V_{i}^{0}\subset V_{i}$ be the subspace of functions which satisfy the
homogeneous Dirichlet conditions on $\partial{\Omega}$. Throughout the text we
let the index $\ell$ range over the boundary nodes first, in other words,
${y}_{i}^{\ell}\in{\Omega}$ for $\ell\leq N_{i}:=\dim V_{i}^{0}$.

In order to construct a finite element method which is both stable and
consistent, care needs to be taken when imposing non-homogeneous boundary
conditions. Due to the more delicate stability properties of Isaacs operators
compared to Bellman operators, on non-convex domains a nodal interpolant is
not suitable to map the boundary data onto the approximation space.

Given $w\in C(H^{1}({\Omega}))$, we consider therefore linear mappings
${P_{i}}$ which map $w$ into $V_{i}$ such that, for all
$\hat{\phi}_{i}^{\ell}\in V_{i}^{0}$,

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

######  Assumption 3.1.

There are linear mappings ${P_{i}}$ satisfying (12) and there is a constant
$C\geq 0$ such that for every $w\in C^{\infty}(\mathds{R}^{d})$ and
$i\in\mathds{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\rightarrow\infty}\left\|{P_{i}}w-w\right\|_{W^{1,\infty}({\Omega})}=0.$
(13)

It is shown in [DLSW12] that Assumption 3.1 holds if we chose ${P_{i}}w$ such
that it interpolates $w$ on the boundary, (3.1) is satisfied, ${\Omega}$ is a
bounded convex polyhedral domain in $\mathds{R}^{d}$, $d\in\\{2,3\\}$, and the
mesh satisfies a local quasi-uniformity condition.

To apply the result for non-convex domains ${\Omega}$ and general $w\in
C^{\infty}(\mathds{R}^{d})$, we consider a convex polyhedral domain $B$
containing ${\Omega}$ and assume there is a locally quasi-uniform mesh on $B$
which coincides with the original mesh on ${\Omega}$. Let $\eta$ be a smooth
cut-off function with relatively compact support in $B$ such that $\eta\equiv
1$ on ${\Omega}$. Then the classical elliptic projection on $B$, acting on
$\eta w:B\to\mathds{R}$, has the required properties. Given this construction
for ${P_{i}}$, it is natural to refer to it as an elliptic projection, see
also [JS13, Section 4]. This construction provides on non-convex domains an
approximation of the boundary data, which permits the proof of stability of
numerical solutions, see Theorem 6.1 below, as well as of consistency, see
Lemma 3.1. In conclusion, let $V_{i}^{g}\subset V_{i}$ be the subspace of
functions which attain ${P_{i}}g$ on the boundary.

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\mathds{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{\\}}$.

The time derivative is approximated by $d_{i}$ for which we let the $\ell$th
entry of $d_{i}w({s}_{i}^{k},\cdot)$ be

$(d_{i}w({s}_{i}^{k},\cdot))_{\ell}=\frac{w({s}_{i}^{k+1},{y}_{i}^{\ell})-w({s}_{i}^{k},{y}_{i}^{\ell})}{h_{i}}.$

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

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

For each $i$ we also seek a discretisation ${f}_{i}^{(\alpha,\beta)}$ of
$f^{(\alpha,\beta)}$. We now want to make the conceptual statements
$L^{(\alpha,\beta)}\approx{E}_{i}^{(\alpha,\beta)}+{I}_{i}^{(\alpha,\beta)}$
and $f^{(a,b)}\approx{f}_{i}^{(\alpha,\beta)}$ more precise.

######  Assumption 3.2.

The splitting is consistent in the sense that

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

The coefficients $\bar{c}_{i}^{(\alpha,\beta)}$ and
$\bar{\bar{c}}_{i}^{(\alpha,\beta)}$ are non-negative and that there exists
$\gamma\in\mathds{R}$ such that

$\displaystyle\|\bar{c}_{i}^{(\alpha,\beta)}\|_{L^{\infty}}+\|\bar{\bar{c}}_{i}^{(\alpha,\beta)}\|_{L^{\infty}}\leq\gamma,\qquad\forall\,i\in\mathds{N},\;\forall\,(\alpha,\beta)\in
A\times B.$ (14)

The family of mappings

$\displaystyle\begin{array}[]{ll}A\times B\to
L^{\infty},\;(\alpha,\beta)\mapsto(\bar{a}_{i}^{(\alpha,\beta)},\bar{b}_{i}^{(\alpha,\beta)},\bar{c}_{i}^{(\alpha,\beta)},\bar{\bar{a}}_{i}^{(\alpha,\beta)},\bar{\bar{b}}_{i}^{(\alpha,\beta)},\bar{\bar{c}}_{i}^{(\alpha,\beta)},f_{i}^{(\alpha,\beta)})\end{array}$

is continuous for each $i$.

The splitting into explicit and implicit part is used to define the internal
numerical operators ${\mathsf{E}}_{{\Omega},i}^{(\alpha,\beta)}$ and
${\mathsf{I}}_{{\Omega},i}^{(\alpha,\beta)}$ as mappings from
$H^{1}({\Omega})$ to $\mathds{R}^{\dim V_{i}}$:

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

where $\ell$ ranges over all internal nodes, i.e. $\ell\leq N_{i}$. For
boundary nodes $\ell>N_{i}$ we set
$({\mathsf{E}}_{{\Omega},i}^{(\alpha,\beta)}w)_{\ell}=({\mathsf{I}}_{{\Omega},i}^{(\alpha,\beta)}w)_{\ell}=({\mathsf{F}}_{{\Omega},i}^{(\alpha,\beta)})_{\ell}=0$.
Likewise we define the boundary numerical operators
${\mathsf{E}}^{(\alpha,\beta)}_{\partial{\Omega},i}$ and
${\mathsf{I}}^{(\alpha,\beta)}_{\partial{\Omega},i}$ from $H^{1}({\Omega})$ to
$\mathds{R}^{\dim V_{i}}$:

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

where $\ell$ ranges over all boundary nodes, i.e. $\ell>N_{i}$. For internal
nodes $\ell\leq N_{i}$ we set
$({\mathsf{E}}^{(\alpha,\beta)}_{\partial{\Omega},i}w)_{\ell}=({\mathsf{I}}^{(\alpha,\beta)}_{\partial{\Omega},i}w)_{\ell}=({\mathsf{F}}^{(\alpha,\beta)}_{\partial{\Omega},i})_{\ell}=0$.
Finally, we obtain the numerical operators ${\mathsf{E}}_{i}^{(\alpha,\beta)}$
and ${\mathsf{I}}_{i}^{(\alpha,\beta)}$ by adding the internal and the
boundary operators, i.e.
${\mathsf{E}}_{i}^{(\alpha,\beta)}:={\mathsf{E}}_{{\Omega},i}^{(\alpha,\beta)}+{\mathsf{E}}^{(\alpha,\beta)}_{\partial{\Omega},i}={\mathsf{E}}_{{\Omega},i}^{(\alpha,\beta)}$,
${\mathsf{I}}_{i}^{(\alpha,\beta)}:={\mathsf{I}}_{{\Omega},i}^{(\alpha,\beta)}+{\mathsf{I}}^{(\alpha,\beta)}_{\partial{\Omega},i}$.
Similarly,
${\mathsf{F}}_{i}^{(\alpha,\beta)}:={\mathsf{F}}_{{\Omega},i}^{(\alpha,\beta)}+{\mathsf{F}}^{(\alpha,\beta)}_{\partial{\Omega},i}$.

When restricted to the domain $V_{i}$, the numerical operators have matrix
representations with respect to the nodal bases
$\\{{\phi}_{i}^{\ell}\\}_{\ell}$, which we also denote by
${\mathsf{E}}_{i}^{(\alpha,\beta)}$ and ${\mathsf{I}}_{i}^{(\alpha,\beta)}$.

Throughout the paper, we make use of the partial ordering of $\mathds{R}^{n}$:
for $x,\,y\in\mathds{R}^{n}$, we write $x\geq y$ if and only if $x_{\ell}\geq
y_{\ell}$ for all $\ell\in\left\\{1,\dots,n\right\\}$. For collections
$\left\\{x^{\alpha}\right\\}_{\alpha}\subset\mathds{R}^{n}$ and
$\left\\{y^{\beta}\right\\}_{\beta}\subset\mathds{R}^{n}$, we define the
operators $\sup_{\alpha}$ and $\inf_{\beta}$ componentwise:

$\bigl{(}\sup_{\alpha}x^{\alpha}\bigr{)}_{\ell}=\sup_{\alpha}x^{\alpha}_{\ell}\quad\textrm{and}\quad\bigl{(}\inf_{\beta}y^{\beta}\bigr{)}_{\ell}=\inf_{\beta}y^{\beta}_{\ell}.$

We can now state the numerical scheme for finding an approximate solution to
(3). We initialise the scheme with the elliptic projection
$v_{i}(T,\cdot)={P_{i}}v_{T}$. Then, in order to find the numerical solution
$v_{i}({s}_{i}^{k},\cdot)\in V^{g}_{i}$, we proceed inductively over the
remaining time steps $k\in\left\\{T/h_{i}-1,\dots,0\right\\}$:

$\displaystyle-d_{i}v_{i}({s}_{i}^{k},\cdot)+\adjustlimits{\inf}_{\beta\in
B}{\sup}_{\alpha\in
A}\bigl{(}{\mathsf{E}}_{i}^{(\alpha,\beta)}v_{i}({s}_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{(\alpha,\beta)}v_{i}({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{(\alpha,\beta)}\bigr{)}=0.$
(17)

If all ${\mathsf{I}}_{i}^{(\alpha,\beta)}$ vanish then (17) is an explicit
scheme, otherwise it is implicit. Note that the terms $d_{i}v_{i}$
representing the discretized time derivative vanish on the boundary nodes
because $g$ is time independent.

To show the existence and uniqueness of numerical solutions, it is useful to
provide an equivalent reformulation of the scheme. Let
$w:S_{i}\times\Omega\to\mathds{R}$ be a function that satisfies
$w({s}_{i}^{k},\cdot)\in H^{1}({\Omega})$ for all ${s}_{i}^{k}\in S_{i}$.
Given $\beta$ we denote by $\alpha_{i}^{k,\ell,\beta}(w)$ a maximiser of (18a)
below. Similarly $\beta_{i}^{k,\ell}(w)$ is a minimiser of (18b):

$\displaystyle\sup_{\alpha\in
A}\left({\mathsf{E}}_{i}^{(\alpha,\beta)}w({s}_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{(\alpha,\beta)}w({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{(\alpha,\beta)}\right)_{\ell},$
(18a) $\displaystyle\inf_{\beta\in B}\left({\sf
E}_{i}^{(\alpha_{i}^{k,\ell,\beta}(w),\beta)}w({s}_{i}^{k+1},\cdot)+{\sf
I}_{i}^{(\alpha_{i}^{k,\ell,\beta}(w),\beta)}w({s}_{i}^{k},\cdot)-{\sf
F}_{i}^{(\alpha_{i}^{k,\ell,\beta}(w),\beta)}\right)_{\ell}.$ (18b)

Finally, we write
$\alpha_{i}^{k,\ell}(w)=\alpha_{i}^{k,\ell,\beta_{i}^{k,\ell}(w)}(w)$, i.e.
drop the reference to $\beta$ to imply that $\beta$ is chosen optimally.

Let ${\mathsf{I}}_{i}^{k,w}$ and ${\mathsf{E}}_{i}^{k,w}$ be the matrices
whose $\ell$th row at $k$th time step is equal to that of

${\mathsf{I}}_{i}^{(\alpha_{i}^{k,\ell}(w),\beta_{i}^{k,\ell}(w))}\quad\mbox{and}\quad{\mathsf{E}}_{i}^{(\alpha_{i}^{k,\ell}(w),\beta_{i}^{k,\ell}(w))},$

respectively. Also let the $\ell$th entry of ${\mathsf{F}}_{i}^{k,w}$ be

$\left({\mathsf{F}}_{i}^{(\alpha_{i}^{k,\ell}(w),\beta_{i}^{k,\ell}(w))}\right)_{\ell}.$

Notice that control $(\alpha_{i}^{k,\ell}(w),\beta_{i}^{k,\ell}(w))$ is not
necessarily unique. The subsequent analysis is valid regardless of the choice
of the control.

We may reformulate the numerical scheme (17) with ${\mathsf{E}}_{i}^{k,w}$,
${\mathsf{I}}_{i}^{k,w}$ and ${\mathsf{F}}_{i}^{k,w}$ as follows. Initialise
the scheme with the elliptic projection so that $v_{i}(T,\cdot)=P_{i}v_{T}$.
Then for each $k\in\\{T/h_{i}-1,\dots,1,0\\}$, $v_{i}\in V_{i}^{g}$ solve

$(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}}=0.$ (19)

###  3.1 Consistency properties of elliptic projections

We conclude the section with the consistency properties of linear operators
for fixed $(\alpha,\beta)\in A\times B$. The result extends [JS13] even in the
Bellman setting by including non-homogenous boundary conditions. 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}):=(h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf
Id})\,w({s}_{i}^{k},\cdot)+(h_{i}{\mathsf{E}}_{i}^{k,w}-{\sf
Id})\,w({s}_{i}^{k+1},\cdot)-h_{i}{\mathsf{F}}_{i}^{k,w}$ (20)

for $w({s}_{i}^{k},\cdot)\in V^{g}_{i}$. Note that despite of the similar
notation, the operator ${\mathsf{F}}_{i}^{k,w}$ has a quite different meaning
in that it denotes the forcing term vector, as defined in the first part of
this section.

######  Lemma 3.1.

Let $\psi\in C^{2}(\mathds{R}\times\mathds{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))$
(21)

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)).$
(22)

###### Proof.

We only prove (21) since the result for (22) 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

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

Step 2: If ${y}_{i}^{\ell}\in{\Omega}$ then the proof is analogous to the
Hamilton–Jacobi–Bellman case described in [JS13], once the control $\alpha\in
A$ is replaced by the pair of controls $(\alpha,\beta)\in A\times B$. Both
Assumptions 3.1 and 3.2 are used here.

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

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

Step 4: 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}$ or
$\partial{\Omega}$, respectively. Then the conclusions of Steps 1 to 3 above
may be applied to the individual subsequences to arrive at (21). ∎

##  4 Monotonicity

Monotonicity properties of the numerical scheme play a crucial role in
ensuring convergence to the viscosity solution.

######  Definition 4.1.

An operator $F:V_{i}\to\mathds{R}^{N_{i}}$ is said to satisfy the local
monotonicity property (LMP) if for all $v\in V_{i}$ such that $v$ has a non-
positive local minimum at the internal node ${y}_{i}^{\ell}$, we have
$(Fv)_{\ell}\leq 0$. The operator $F$ satisfies the weak discrete maximum
principle (wDMP) provided that for any $v\in V_{i}$,

$\text{if \,}\bigl{(}Fv\bigr{)}_{\ell}\geq 0\text{ for all
}\ell\in\\{1,\dots,N_{i}\\},\quad\text{ then
}\min_{{\Omega}}v\geq\min\\{\min_{\partial{\Omega}}v,0\\}.$ (25)

Suppose an operator $F$ satisfies the LMP. If $v\in V_{i}$ has a negative
local minimum at an internal node ${y}_{i}^{\ell}$, then $((F+\epsilon{\sf
Id})v)_{\ell}<0$ for any positive $\epsilon$. In particular, the
contrapositive of (25) holds. Therefore $F+\epsilon{\sf Id}$ satisfies the
wDMP for any positive $\epsilon$.

######  Assumption 4.1.

For each $(\alpha,\beta)\in A\times B$, we assume that
${\mathsf{E}}_{i}^{(\alpha,\beta)}$, restricted to $V_{i}$, has non-positive
off-diagonal entries. We assume that $h_{i}$ is small enough so that
$h_{i}{\mathsf{E}}_{i}^{(\alpha,\beta)}-{\sf Id}$ is monotone for every
$(\alpha,\beta)$, i.e. so that all entries of all
$h_{i}{\mathsf{E}}_{i}^{(\alpha,\beta)}-{\sf Id}$ are non-positive. For each
$(\alpha,\beta)\in A\times B$, we suppose that
${\mathsf{I}}_{i}^{(\alpha,\beta)}$ satisfies the LMP.

We now turn to the matrices $h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf Id}$ and
$h_{i}{\mathsf{E}}_{i}^{k,w}-{\sf Id}$ which will later be used in the proof
of the well-posedness of the scheme (19).

######  Lemma 4.1.

Consider a $w:S_{i}\times\Omega\to\mathds{R}$ so that $w({s}_{i}^{k},\cdot)\in
H^{1}({\Omega})$ for all ${s}_{i}^{k}\in S_{i}$. Then:

  1. 1.

The matrices $h_{i}{\mathsf{E}}_{i}^{k,w}-{\sf Id}$ restricted to $V_{i}^{0}$
are monotone.

  2. 2.

The mapping $w\mapsto-(h_{i}{\mathsf{E}}_{i}^{k,w}-{\sf Id})\,w$ is positive:
if $w\geq 0$ then $-(h_{i}{\mathsf{E}}_{i}^{k,w}-{\sf Id})\,w\geq 0$.

  3. 3.

The matrices of $h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf Id}$ restricted to $V_{i}$
are strictly diagonally dominant M-matrices.

  4. 4.

For fixed $w$, the operators $v\mapsto{\mathsf{I}}_{i}^{k,w}v$ and
$v\mapsto(h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf Id})\,v$ satisfy, respectively, the
LMP and wDMP.

###### Proof.

The proof for all but 3. is analogous to that of [JS13, Lemma 2.3], once the
control $\alpha\in A$ is replaced by the pair of controls $(\alpha,\beta)\in
A\times B$, using Assumption 4.1.

For 3. the strict diagonal dominance of $h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf Id}$
for rows with indices $\ell\leq N_{i}$ is a direct consequence of its LMP by
the same argument as in the case of restriction to $V^{0}_{i}$ discussed in
[JS13]. For indices $\ell>N_{i}$ rows of matrices
$h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf Id}$ are equivalent to that of identity
matrix which is strictly diagonally dominant by definition. M-matrix property
follows from [BP94, Chapter 6, Theorem 2.3, $(M_{35})$] after noticing that
$h_{i}{\mathsf{I}}_{i}^{k,w}+{\sf Id}$ are Z-matrices. ∎

###  4.1 Monotonicity through artificial diffusion

We show now that, using the method of artificial diffusion one can ensure that
${\mathsf{E}}_{i}^{(\alpha,\beta)}$ and ${\mathsf{I}}_{i}^{(\alpha,\beta)}$
satisfy Assumption 4.1, provided the meshes are strictly acute.

Let $\mathcal{T}_{i}$ be the mesh corresponding to the finite element space
$V_{i}$. Given a function $g:\Omega\to\mathds{R}^{d}$ and an element $K$ of
$\mathcal{T}_{i}$, we denote

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

We say that the meshes $\mathcal{T}_{i}$ are strictly acute if there exists
${\vartheta}\in(0,\pi/2)$ such that for all $i\in\mathds{N}$:

$\nabla{\phi}_{i}^{\ell}\cdot\nabla\phi_{i}^{l}\bigl{|}_{K}\leq-\,\sin({\vartheta})\;|\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}.$ (26)

Consider a splitting of the form

$\displaystyle a^{(\alpha,\beta)}$
$\displaystyle=\tilde{a}_{i}^{(\alpha,\beta)}+\tilde{\tilde{a}}_{i}^{(\alpha,\beta)},$
$\displaystyle b^{(\alpha,\beta)}$
$\displaystyle=\bar{b}_{i}^{(\alpha,\beta)}+\bar{\bar{b}}_{i}^{(\alpha,\beta)},$
$\displaystyle c^{(\alpha,\beta)}$
$\displaystyle=\bar{c}_{i}^{(\alpha,\beta)}+\bar{\bar{c}}_{i}^{(\alpha,\beta)},$
$\displaystyle f^{(\alpha,\beta)}$ $\displaystyle={f}_{i}^{(\alpha,\beta)},$

where all terms are in $C(\overline{{\Omega}})$ or
$C(\overline{{\Omega}},\mathds{R}^{d})$. Notice the tilde in the notation of
the diffusion coefficients. Choose non-negative artificial diffusion
coefficients $\bar{\nu}_{i}^{(\alpha,\beta),\ell}$ and
$\bar{\bar{\nu}}_{i}^{(\alpha,\beta),\ell}$ such that for all $K$ that have
${y}_{i}^{\ell}$ as vertex:

$\displaystyle\bigl{(}|\bar{b}_{i}^{\alpha}|_{K}\,+\Delta
x_{K}\|\bar{c}_{i}^{\alpha}\|_{L^{\infty}(K)}\bigr{)}\leq\,$
$\displaystyle\bar{\nu}_{i}^{(\alpha,\beta),\ell}\,\sin({\vartheta})\,|\nabla\hat{\phi}_{i}^{\ell}|_{K}\,{\rm
vol}(K),$ (27a)
$\displaystyle\bigl{(}|\bar{\bar{b}}_{i}^{\alpha}|_{K}\,+\Delta
x_{K}\|\bar{\bar{c}}_{i}^{\alpha}\|_{L^{\infty}(K)}\bigr{)}\leq\,$
$\displaystyle\bar{\bar{\nu}}_{i}^{(\alpha,\beta),\ell}\,\sin({\vartheta})\,|\nabla\hat{\phi}_{i}^{\ell}|_{K}\,{\rm
vol}(K).$ (27b)

Select $\bar{a}_{i}^{(\alpha,\beta)}$ and $\bar{\bar{a}}_{i}^{(\alpha,\beta)}$
both in $C(\overline{{\Omega}})$ such that

$\displaystyle\bar{a}_{i}^{(\alpha,\beta)}({y}_{i}^{\ell})$
$\displaystyle\geq\max\bigl{\\{}\tilde{a}_{i}^{(\alpha,\beta)}({y}_{i}^{\ell}),\bar{\nu}_{i}^{(\alpha,\beta),\ell}\bigr{\\}},$
(28a) $\displaystyle\bar{\bar{a}}_{i}^{(\alpha,\beta)}({y}_{i}^{\ell})$
$\displaystyle\geq\max\bigl{\\{}\tilde{\tilde{a}}_{i}^{(\alpha,\beta)}({y}_{i}^{\ell}),\bar{\bar{\nu}}_{i}^{(\alpha,\beta),\ell}\bigr{\\}}.$
(28b)

The splittings of
$a^{(\alpha,\beta)}=\bar{a}_{i}^{(\alpha,\beta)}+\bar{\bar{a}}_{i}^{(\alpha,\beta)}$,
$b^{(\alpha,\beta)}=\bar{b}_{i}^{(\alpha,\beta)}+\bar{\bar{b}}_{i}^{(\alpha,\beta)}$,
$c^{(\alpha,\beta)}=\bar{c}_{i}^{(\alpha,\beta)}+\bar{\bar{c}}_{i}^{(\alpha,\beta)}$
and $f^{(\alpha,\beta)}={f}_{i}^{(\alpha,\beta)}$ can always be chosen so that
also Assumption 3.2 and (28) hold. Moreover, it is shown in [JS13] that (28)
implies Assumption 4.1.

##  5 Wellposedness and a solution algorithm

Algorithm 1 Howard’s method

1:$\eta>0$, $k\in\\{0,\dots,T/h_{i}-1\\}$, $w\in V_{i}^{g}$,
$(\alpha_{0},\beta_{0})\in A\times B$

2:$u\leftarrow w$

3:$\beta\leftarrow\beta_{0}$

4:while $\|\Psi^{(\alpha,\beta)}(u,w)\|\geq\eta$ do

5: $\alpha\leftarrow\alpha_{0}$

6: while $\|\Psi^{(\alpha,\beta)}(u,w)\|\geq\eta$ do

7: $u\leftarrow\Phi^{(\alpha,\beta)}(w)$

8: $\alpha\leftarrow{\rm argmax}_{\alpha^{\prime}\in
A}\Psi^{(\alpha^{\prime},\beta)}(u,w)$

9: end while

10: $u\leftarrow\Phi^{(\alpha,\beta)}(w)$

11: $\beta\leftarrow{\rm argmin}_{\beta^{\prime}\in
B}\Psi^{(\alpha,\beta^{\prime})}(u,w)$

12:end while

13:return u

In this section we address the existence and uniqueness of numerical solutions
and propose a method to compute them.

######  Theorem 5.1.

There exists a unique numerical solution $v_{i}\colon S_{i}\to V_{i}^{g}$ that
solves (17).

###### Proof.

Fix $k$ and $v_{i}(s_{i}^{k+1},\cdot)$. Then [BMZ09, Theorem 5.2] shows the
existence and uniqueness of a solution $v_{i}(s_{i}^{k},\cdot)$ of (17),
noting the continuity of $(\alpha,\beta)\to{\mathsf{I}}_{i}^{(\alpha,\beta)}$
follows from Assumption 3.2 and the monotonicity from Assumption 4.1. ∎

The proof of Theorem 5.2 of [BMZ09], which we used here, relies on the exact
solution of the inner Bellman equation (18a), which is reached by a Howard’s
algorithm in the limit. Using these exact solutions of (18a) in an outer
Howard loop yields a sequence whose limit is the solution of (17).

To ensure the completion of the algorithm in finite time we require a
termination criterion for the inner and outer loop. In Algorithm 1 (Howard’s
method) the termination criterion is posed in terms of a tolerance $\eta$. The
statement of the algorithm also refers to

$\Psi^{(\alpha,\beta)}(u,w):=(h_{i}{\mathsf{I}}_{i}^{(\alpha,\beta)}+{\sf
Id})\,u+(h_{i}{\mathsf{E}}_{i}^{(\alpha,\beta)}-{\sf
Id})w-h_{i}{\mathsf{F}}_{i}^{(\alpha,\beta)},$

given $u,w\in V_{i}^{g}$ and $(\alpha,\beta)\in A\times B$. Given just $w$, we
let $\Phi^{(\alpha,\beta)}(w)$ be the solution $u$ of
$\Psi^{(\alpha,\beta)}(u,w)=0$.

Suppose that $\sum_{\ell}\eta_{\ell}<\infty$ and that $w_{\ell}$ is the
function returned by Algorithm 1 with $\eta=\eta_{\ell}$,
$w=v_{i}({s}_{i}^{k+1},\cdot)$ and a fixed $(\alpha_{0},\beta_{0})\in A\times
B$. Then Theorem 5.4 of [BMZ09] ensures that the $w_{\ell}$ exist and converge
to the unique numerical solution $v_{i}(s_{i}^{k},\cdot)$ as $\ell\to\infty$.

##  6 Stability

For Hamilton–Jacobi–Bellman equations one can bound the value function by the
solution of any linear evolution problem associated to a fixed control. In
contrast the value function of an Isaacs equation may lie in part above and in
part below the solution of such a linear problem. This difference between the
Bellman and Isaacs equations extends to the proof of stability of numerical
solutions. We begin by adapting [JS13, Lemma 3.2].

######  Lemma 6.1.

One has $\lVert(h_{i}{\mathsf{I}}_{i}^{w}+{\sf Id})^{-1}\rVert_{\infty}\leq 1$
and $\lVert h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id}\rVert_{\infty}\leq 1$ for all
$i\in\mathds{N}$ and $w\in V_{i}$, where the norms are the matrix
$\infty$-norms of linear mappings from $V_{i}$ into $V_{i}$ and $V_{i}^{0}$
into $V_{i}^{0}$, respectively.

###### Proof.

Define $v=\sum_{\ell=1}^{\dim V_{i}}{\phi}_{i}^{\ell}\equiv 1$, and
$v_{0}=\sum_{\ell=1}^{N_{i}}{\phi}_{i}^{\ell}\in V_{i}^{0}$. By Lemma 4.1,
$h_{i}{\mathsf{I}}_{i}^{w}+{\sf Id}$ is an invertible M-matrix on $V_{i}$.
Thus, $(h_{i}{\mathsf{I}}_{i}^{w}+{\sf Id})^{-1}\geq 0$ entrywise on $V_{i}$.
Thus

$\begin{split}\lVert(h_{i}{\mathsf{I}}_{i}^{w}+{\sf
Id})^{-1}\rVert_{\infty}&=\max_{1\leq\ell\leq\dim V_{i}}\sum_{j=1}^{\dim
V_{i}}(h_{i}{\mathsf{I}}_{i}^{w}+{\sf Id})^{-1}_{\ell
j}=\max_{1\leq\ell\leq\dim V_{i}}\left((h_{i}{\mathsf{I}}_{i}^{w}+{\sf
Id})^{-1}\mathbf{1}\right)_{\ell},\end{split}$ (29)

where $\mathbf{1}\in\mathds{R}^{\dim V_{i}}$ is the vector with all entries
equal to $1$.

Since $\nabla v\equiv 0$ (as $v\equiv 1$) we have for each $1\leq\ell\leq
N_{i}$ that

$((h_{i}{\mathsf{I}}_{i}^{w}+{\sf
Id})v)_{\ell}=1+h_{i}\langle\bar{\bar{c}}_{i}^{(\alpha_{i}^{\ell}(w),\beta^{\ell}(w))},\hat{\phi}_{i}^{\ell}\rangle\geq
1,$

where we have used non-negativity of $\bar{\bar{c}}_{i}^{\alpha}$ and defined
$(\alpha_{i}^{\ell}(w),\beta^{\ell}(w))$ analogously to
$(\alpha_{i}^{k,\ell}(w),\beta_{i}^{k,\ell}(w))$ in (18). Similarly, for each
boundary node index $N_{i}<\ell\leq\dim V_{i}$

$((h_{i}{\mathsf{I}}_{i}^{w}+{\sf Id})v)_{\ell}=1+h_{i}\geq 1.$

It follows that $(h_{i}{\mathsf{I}}_{i}^{w}+{\sf Id})v\geq\mathbf{1}$.
Combining positive invertibility of $(h_{i}{\mathsf{I}}_{i}^{w}+{\sf Id})$
with (29), we conclude that $\lVert(h_{i}{\mathsf{I}}_{i}^{w}+{\sf
Id})^{-1}\rVert_{\infty}\leq 1$, as required.

Turning our attention to explicit operators, one has

$\lVert h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id}\rVert_{\infty}=\max_{1\leq\ell\leq
N_{i}}\bigl{(}-(h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id})\,v_{0}\bigr{)}_{\ell}$

because all entries of the matrix $h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id}$ are
non-positive. For each $1\leq\ell\leq N_{i}$,

$\bigl{(}(h_{i}{\mathsf{E}}_{i}^{w}-{\sf
Id})\,v\bigr{)}_{\ell}=\bigl{(}(h_{i}{\mathsf{E}}_{i}^{w}-{\sf
Id})\,v_{0}\bigr{)}_{\ell}+\sum_{j=N_{i}+1}^{\dim
V_{i}}\left(h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id}\right)_{\ell
j}\leq\bigl{(}(h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id})\,v_{0}\bigr{)}_{\ell},$

so

$(-(h_{i}{\mathsf{E}}_{i}^{w}-{\sf
Id})\,v_{0})_{\ell}\leq(-(h_{i}{\mathsf{E}}_{i}^{w}-{\sf
Id})\,v)_{\ell}=1-h_{i}\langle\bar{c}_{i}^{(\alpha_{i}^{\ell}(w),\beta^{\ell}(w))},\hat{\phi}_{i}^{\ell}\rangle\leq
1$

because $\bar{c}_{i}^{(\alpha_{i}^{\ell}(w),\beta^{\ell}(w))}\geq 0$ where
$\bar{c}_{i}^{(\alpha_{i}^{\ell}(w),\beta^{\ell}(w))}$ is defined analogously
to $\bar{\bar{c}}_{i}^{(\alpha_{i}^{\ell}(w),\beta^{\ell}(w))}$ above.
Therefore, $-(h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id})v_{0}\leq\mathbf{1}$. So
$\lVert h_{i}{\mathsf{E}}_{i}^{w}-{\sf Id}\rVert_{\infty}\leq 1$. ∎

######  Theorem 6.1.

The numerical solutions $v_{i}$ are uniformly bounded in the $L^{\infty}$
norm:

$\displaystyle\left\|v_{i}\right\|_{L^{\infty}(S_{i}\times{\Omega})}\leq$
$\displaystyle\left\|{P_{i}}v_{T}\right\|_{L^{\infty}({\Omega})}+2\lVert
P_{i}g\rVert_{L^{\infty}({\Omega})}$
$\displaystyle+T\sup_{(\alpha,\beta)}\Bigl{[}\left(\lVert\bar{a}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}+\lVert\bar{\bar{a}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}\right)\lVert\Delta
g\rVert_{L^{\infty}({\Omega})}$
$\displaystyle+\left(\lVert\bar{b}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}+\lVert\bar{\bar{b}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}\right)\lVert{P_{i}}g\rVert_{W^{1,\infty}({\Omega})}$
$\displaystyle+\left(\lVert\bar{c}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}+\lVert\bar{\bar{c}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}\right)\lVert{P_{i}}g\rVert_{L^{\infty}({\Omega})}+\lVert
f_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}\Bigr{]}.$

###### Proof.

We split the numerical solution into two parts

$v_{i}({s}_{i}^{k},\cdot)=v_{i}^{0}({s}_{i}^{k},\cdot)+g_{i},\qquad{s}_{i}^{k}\in
S_{i},$

where $v_{i}^{0}({s}_{i}^{k},\cdot)\in V_{i}^{0}$ and $g_{i}:={P_{i}}g\in
V_{i}^{g}$. Then (19) becomes

$\displaystyle(h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf
Id})\,(v_{i}^{0}({s}_{i}^{k},\cdot)+g_{i})+(h_{i}{\mathsf{E}}_{i}^{k,v_{i}}-{\sf
Id})\,(v_{i}^{0}({s}_{i}^{k+1},\cdot)+g_{i})-h_{i}{\mathsf{F}}_{i}^{k,v_{i}}=0$

or equivalently

$\displaystyle
v_{i}^{0}({s}_{i}^{k},\cdot)=(h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf
Id})^{-1}\bigl{(}h_{i}{\mathsf{F}}_{i}^{k,v_{i}}-h_{i}({\mathsf{I}}_{i}^{k,v_{i}}+{\mathsf{E}}_{i}^{k,v_{i}})\,g_{i}-(h_{i}{\mathsf{E}}_{i}^{k,v_{i}}-{\sf
Id})\,v_{i}^{0}({s}_{i}^{k+1},\cdot)\bigr{)}.$

Using (12), we find

$\displaystyle\|({\mathsf{I}}_{i}^{k,v_{i}}+{\mathsf{E}}_{i}^{k,v_{i}})\,g_{i}\|_{L^{\infty}({\Omega})}\leq$
$\displaystyle\,\left(\lVert\bar{a}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}+\lVert\bar{\bar{a}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}\right)\lVert\Delta
g\rVert_{L^{\infty}({\Omega})}$
$\displaystyle+\left(\lVert\bar{b}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}+\lVert\bar{\bar{b}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}\right)\lVert
g_{i}\rVert_{W^{1,\infty}({\Omega})}$
$\displaystyle+\left(\lVert\bar{c}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}+\lVert\bar{\bar{c}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}\right)\lVert
g_{i}\rVert_{L^{\infty}({\Omega})}.$

Now an application of Lemma 6.1 and an inductive argument over the timesteps
complete the proof. ∎

##  7 Sub- and supersolutions of the PDE and the viscosity BCs

Recall the definition of the upper and lower semi-continuous envelopes $v^{*}$
and $v_{*}$ of a function $v$ and consider their numerical equivalent defined
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}).$
(30)

Owing to Theorem 5.1 and Theorem 6.1, $\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}$. The
proof of the next theorem closely follows that Bellman setting analysed in
[JS13]. We outline briefly the main steps to clarify the changes due to the
additional $\inf$ operation, control set $B$ and the more general boundary
conditions.

At this point we delay the proof that the envelopes satisfy boundary
conditions in the strong sense. We can express that by assuming momentarily
that $\omega=\emptyset$. The extension to $\omega\neq\emptyset$ follows in the
next section.

######  Lemma 7.1.

Let $\omega=\emptyset$. Then $\overline{v}$ is a viscosity subsolution of (3)
and $\underline{v}$ is a viscosity supersolution of (3).

###### Proof.

We show that $\overline{v}$ is a subsolution. That $\underline{v}$ is a
supersolution follows analogously up to a minor asymmetry in the sign of
$\gamma\left|\mu_{i}\right|$ in (33) below. Suppose that $w\in
C^{\infty}(\mathds{R}\times\mathds{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)$.
Then following the argument in [JS13] there exists a sequence
$\left\\{{s}_{i}^{k},{y}_{i}^{\ell}\right\\}_{i}$ such 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$
(31)

and

$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})\;\Leftrightarrow\;{P_{i}}w({s}_{i}^{\kappa},{y}_{i}^{\lambda})+\mu_{i}\geq
v_{i}({s}_{i}^{\kappa},{y}_{i}^{\lambda}),$ (32)

where $\kappa\in\\{k,k+1\\}$, ${y}_{i}^{\lambda}\in{\rm
supp}\,\hat{\phi}_{i}^{\ell}$ and
$\mu_{i}=v_{i}({s}_{i}^{k},{y}_{i}^{\ell})-{P_{i}}w({s}_{i}^{k},{y}_{i}^{\ell})$.
Notice that $\mu_{i}\to 0$ as $i\to\infty$ because of (31). Similar to [JS13]
we conclude from the monotonicity of the scheme (17) that

$\displaystyle 0=\,$ $\displaystyle-
d_{i}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})+\adjustlimits{\inf}_{\beta\in
B}{\sup}_{\alpha\in
A}\left({\mathsf{E}}_{i}^{(\alpha,\beta)}v_{i}(s_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{(\alpha,\beta)}v_{i}({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{(\alpha,\beta)}\right)_{\ell}$
$\displaystyle\geq\,$ $\displaystyle-
d_{i}{P_{i}}w({s}_{i}^{k},{y}_{i}^{\ell})+\adjustlimits{\inf}_{\beta\in
B}{\sup}_{\alpha\in
A}\left({\mathsf{E}}_{i}^{(\alpha,\beta)}{P_{i}}w(s_{i}^{k+1},\cdot)+{\mathsf{I}}_{i}^{(\alpha,\beta)}{P_{i}}w({s}_{i}^{k},\cdot)-{\mathsf{F}}_{i}^{(\alpha,\beta)}\right)_{\ell}$
$\displaystyle-\gamma\left|\mu_{i}\right|$ $\displaystyle=\,$ $\displaystyle
F_{i}{P_{i}}w({s}_{i}^{k},{y}_{i}^{\ell})-\gamma\left|\mu_{i}\right|$ (33)

Recalling Lemma 3.1, and $\lim_{i}\mu_{i}=0$, we take the limit
$i\rightarrow\infty$ in inequality (33) to conclude that

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

Hence $\overline{v}$ is a viscosity subsolution. ∎

##  8 Boundary and final time conditions

A direct consequence of Lemma 7.1 is that envelopes of the numerical solution
satisfy the Dirichlet boundary conditions in the viscosity sense on the entire
$\partial{\Omega}$. For the consistency with the strong boundary conditions on
$\omega$ we assume the existence of two families of barrier functions
$(\zeta_{y,\epsilon})_{\epsilon>0}$ and $(\xi_{y,\epsilon})_{\epsilon>0}$
corresponding to the sets of super- and subsolutions respectively.

######  Assumption 8.1.

Let $\omega\subseteq\partial{\Omega}$. Let us assume the following:

  1. 1.

Family of upper barriers:
For all $y\in\omega$, $\varepsilon>0$ there exists a barrier function
$\zeta_{y,\epsilon}\in W^{1,\infty}([0,T]\times\overline{{\Omega}})$ with:

    1. (a)

There exists an $\hat{i}=\hat{i}(\varepsilon)\in\mathds{N}$ so that for all
$i\geq\hat{i}$, ${s}_{i}^{k}\in S_{i}$ there exists a minimiser of
$v_{i}-P_{i}\zeta_{y,\epsilon}$ over
$\\{{s}_{i}^{k}\\}\times\overline{{\Omega}}$ which lies on
$\\{{s}_{i}^{k}\\}\times\partial{\Omega}$.

    2. (b)

For $t\in[0,T]$ let $q_{t,y,\epsilon}$ be a minimiser of
$g-\zeta_{y,\epsilon}$ over $\\{t\\}\times\partial{\Omega}$. Then,

$\lim_{\epsilon\rightarrow 0}\
q_{t,y,\epsilon}=y\quad\text{and}\quad\lim_{\epsilon\rightarrow
0}\zeta_{y,\epsilon}(t,y)-\zeta_{y,\epsilon}\ (t,q_{t,y,\epsilon})\geq 0.$

  2. 2.

Family of lower barriers:
For all $y\in\omega$, $\varepsilon>0$ there exists a barrier function
$\xi_{y,\epsilon}\in W^{1,\infty}([0,T]\times\overline{{\Omega}})$ with:

    1. (a)

There exists an $\hat{i}=\hat{i}(\varepsilon)\in\mathds{N}$ so that for all
$i\geq\hat{i}$, ${s}_{i}^{k}\in S_{i}$ there exists a maximiser of
$v_{i}-P_{i}\xi_{y,\epsilon}$ over
$\\{{s}_{i}^{k}\\}\times\overline{{\Omega}}$ which lies on
$\\{{s}_{i}^{k}\\}\times\partial{\Omega}$.

    2. (b)

For $t\in[0,T]$ let $r_{t,y,\epsilon}$ be a maximiser of $g-\xi_{y,\epsilon}$
over $\\{t\\}\times\partial{\Omega}$. Then,

$\lim_{\epsilon\rightarrow 0}\
r_{t,y,\epsilon}=y\quad\text{and}\quad\lim_{\epsilon\rightarrow
0}\zeta_{y,\epsilon}(t,y)-\xi_{y,\epsilon}\ (t,r_{t,y,\epsilon})\leq 0.$

In order to understand the connection between convergence of the envelopes to
the boundary conditions and the barrier functions let us consider the
following example.

######  Example 8.1.

We study a $1$-dimensional test problem with the homogeneous boundary
conditions:

$\displaystyle-\partial_{t}u+\nabla u-1$ $\displaystyle=0,$
$\displaystyle\text{in }[0,T)\times(0,1),$ $\displaystyle u$
$\displaystyle=x,$ $\displaystyle\text{on }\\{T\\}\times(0,1).$

In order to find solution we consider a fully implicit numerical scheme with
artificial diffusion. In this spirit, the exact solution of
$-u_{t}-\upsilon\Delta u+\nabla u-1=0$, interpreting $\upsilon$ as the
artificial diffusion coefficient, is:

$v_{\upsilon}(t,x)=x-1+\frac{e^{\frac{x}{\upsilon}}-e^{\frac{1}{\upsilon}}}{1-e^{\frac{1}{\upsilon}}}.$

For a fixed $\upsilon$ we can choose $\Delta x$ small enough such that the
numerical solution $\tilde{v}_{\upsilon}$ with artificial diffusion $\upsilon$
and uniform mesh with mesh size $\Delta x$ satisfies
$\left\|\tilde{v}_{\upsilon}-v_{\upsilon}\right\|_{W^{1,\infty}}<\upsilon$. We
note that $v_{\upsilon}$ attains its maximum at

$x_{\mathrm{max}}=\upsilon\log\left[\upsilon\left(e^{\frac{1}{\upsilon}}-1\right)\right]$

with maximal value equal to

$v_{\upsilon}(t,x_{\mathrm{max}})=\frac{1}{e^{\frac{1}{\upsilon}}-1}+\upsilon\log\left[\upsilon\left(e^{\frac{1}{\upsilon}}-1\right)\right].$

In particular, it follows that $x_{\mathrm{max}}\to 1$ and
$v_{\upsilon}(t,x_{\mathrm{max}})\to 1$ as $\upsilon\to 0$. Since

$\left\|\tilde{v}_{\upsilon}-v_{\upsilon}\right\|_{W^{1,\infty}}<\upsilon,$

we conclude that for a decreasing artificial diffusion coefficient the
sequence of numerical solutions $\tilde{v}_{\upsilon}$ has the lower semi-
continuous envelope

$\underline{v}(t,x)=\begin{cases}x&\quad x\neq 1,\\\ 0&\quad x=1.\end{cases}$
(35)

However, since $\underline{v}$ has a discontinuity at $1$, there cannot exist
a Lipschitz continuous barrier function as described in Assumption 8.1 – as we
decrease $\upsilon$, $v_{\upsilon}$ and hence $\tilde{v}_{\upsilon}$ exhibit
an arbitrarily large gradient in the vicinity of the boundary. As a result,
for any barrier function $\xi$ there exists $i$ for which the maximum of
$\tilde{v}_{\upsilon}-P_{i}\xi_{y,\epsilon}$ does not lie on the boundary for
all $\epsilon>0$. Note that we can construct both upper and lower barrier
functions at $x=0$, e. g. $\zeta_{0}=x^{2}$ and $\xi_{0}=(x+1)^{2}$.

######  Lemma 8.1.

Given Assumption 8.1, we have $\overline{v}(t,x)=\underline{v}(t,x)=g(t,x)$
for all $(t,x)\in[0,T]\times\omega$.

###### Proof.

We focus on the case of $v^{*}(t,x)$ as the other case follows analogously.
Let $y\in\omega$ and consider a sequence
$({s}_{i}^{k},{y}_{i}^{\ell})\rightarrow(t,y)$ as $i\rightarrow\infty$. We
have that

$\displaystyle\limsup_{i\rightarrow\infty}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})=$
$\displaystyle\lim_{i\rightarrow\infty}P_{i}\xi_{y,\epsilon}({s}_{i}^{k},{y}_{i}^{\ell})+\limsup_{i\rightarrow\infty}[v_{i}({s}_{i}^{k},{y}_{i}^{\ell})-P_{i}\xi_{y,\epsilon}({s}_{i}^{k},{y}_{i}^{\ell})]$
$\displaystyle\leq$
$\displaystyle\xi_{y,\epsilon}(t,y)+\limsup_{i\rightarrow\infty}\sup_{z\in\partial{\Omega}}[v_{i}({s}_{i}^{k},z)-P_{i}\xi_{z,\epsilon}({s}_{i}^{k},z)]$
$\displaystyle=$
$\displaystyle\xi_{y,\epsilon}(t,y)+\limsup_{i\rightarrow\infty}\sup_{z\in\partial{\Omega}}[{P_{i}}g({s}_{i}^{k},z)-P_{i}\xi_{y,\epsilon}({s}_{i}^{k},z)]$
$\displaystyle=$
$\displaystyle\xi_{y,\epsilon}(t,y)+\limsup_{i\rightarrow\infty}\sup_{z\in\partial{\Omega}}\Big{[}{P_{i}}g(t,z)-P_{i}\xi_{y,\epsilon}(t,z)\Big{.}$
$\displaystyle\left.+\left({P_{i}}g({s}_{i}^{k},z)-{P_{i}}g(t,z)\right)-\left(P_{i}\xi_{y,\epsilon}({s}_{i}^{k},z)-P_{i}\xi_{y,\epsilon}(t,z)\right)\right]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}$
$\displaystyle\xi_{y,\epsilon}(t,y)+\limsup_{i\rightarrow\infty}\sup_{z\in\partial{\Omega}}[{P_{i}}g(t,z)-P_{i}\xi_{y,\epsilon}(t,z)]$
$\displaystyle=$
$\displaystyle\xi_{y,\epsilon}(t,y)+\limsup_{i\rightarrow\infty}\sup_{z\in\partial{\Omega}}[\left(g(t,z)-\xi_{y,\epsilon}(t,z)\right)$
$\displaystyle+\left({P_{i}}g(t,z)-g(t,z)\right)-\left(P_{i}\xi_{y,\epsilon}(t,z)-\xi_{y,\epsilon}(t,z)\right)]$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}$
$\displaystyle\xi_{y,\epsilon}(t,y)+\limsup_{i\rightarrow\infty}\sup_{z\in\partial{\Omega}}[g(t,z)-\xi_{y,\epsilon}(t,z)]$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}$
$\displaystyle\xi_{y,\epsilon}(t,y)+g(t,r_{t,y,\epsilon})-\xi_{y,\epsilon}(t,r_{t,y,\epsilon}),$
(36)

where we get (a) due to Lipschitz continuity of $g$ and $\xi_{y,\epsilon}$ in
time, (b) due to $L^{\infty}$ convergence of ${P_{i}}g(t,z)$ and
$P_{i}\xi_{y,\epsilon}(t,z)$ as $i\to\infty$ and (c) due to the definition of
$r_{t,y,\epsilon}$. By Assumption 8.1 we have that $\lim_{\epsilon\rightarrow
0}\xi_{y,\epsilon}(t,y)-\xi_{y,\epsilon}\ (t,r_{t,y,\epsilon})\leq 0$ for any
$t\in[0,T]$, so we can conclude that

$\limsup_{i\rightarrow\infty}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})\leq g(t,y),$
(37)

as $\epsilon\rightarrow 0$. The proof concludes by completing a similar
calculation for
$\liminf_{i\rightarrow\infty}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})$. This gives us

$g(t,y)\geq\limsup_{i\rightarrow\infty}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})\geq\liminf_{i\rightarrow\infty}v_{i}({s}_{i}^{k},{y}_{i}^{\ell})\geq
g(t,y),$

and the final result follows. ∎

Finally, also the final time conditions are attained strongly.

######  Lemma 8.2.

The sub and supersolutions $\overline{v}$ and $\underline{v}$ satisfy

$\overline{v}(T,\cdot)=\underline{v}(T,\cdot)=v_{T}\quad\text{ on
}\,\overline{{\Omega}}.$ (38)

###### Proof.

The proof is identical to the Hamilton–Jacobi–Bellman case described in [JS13]
once the control $\alpha\in A$ is replaced by the pair of controls
$(\alpha^{v_{i}(\ell)},\beta^{v_{i}(\ell)})\in A\times B$ similarly as in the
proof of Lemma 6.1, the interpolation operator $\mathcal{I}_{i}$ is replaced
with elliptic projection ${P_{i}}$ and function space $V_{i}^{0}$ is replaced
with $V_{i}^{g}$. ∎

##  9 Uniform convergence

Summarising the previous two sections, we have the sub- and supersolution
properties of the complete final boundary value problem.

######  Theorem 9.1.

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

The last ingredient to establish the convergence of the numerical scheme is a
crucial property of the final boundary value problem: a comparison principle.
Especially when working with differential operators which degenerate on a part
of the domain or classes of singularly perturbed operators, imposing strong
Dirichlet conditions on a part of the boundary while viscosity boundary
conditions on the remainder is appropriate to guarantee comparison and well-
posedness. This is reflected in our analysis in the following assumption in
combination with the formulation of Definition 2.2.

######  Assumption 9.1.

Let $u$ be a viscosity subsolution and $w$ be a viscosity supersolution. Then
$u\leq w$ on ${\Omega}$.

######  Theorem 9.2.

One has $\underline{v}=\overline{v}=v$, where $v$ is the unique viscosity
solution of equation (3). Furthermore,

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

###### Proof.

The proof is essentially identical to the Hamilton–Jacobi–Bellman case
described in [JS13]. ∎

##  10 Construction of barrier functions

In order to justify Assumption 8.1 we now present settings for the
construction of barrier functions. First, we introduce a method for ensuring
the existence of barrier functions for the simpler case of uniformly parabolic
operators on a convex domain for fully implicit numerical schemes. After that
we show the extension to general IMEX schemes and allow non-convex domains as
well as degenerate operators.

We will focus on constructing lower barrier functions $\xi$, since the
argument for upper barrier functions $\zeta$ is symmetric and follows by
changing the direction of inequalities and exchanging $\sup$ and $\inf$
operators where required.

###  10.1 Barrier functions for uniformly parabolic equations on convex
domains

We assume the existence of a function $\xi\in\\{w\in
W^{1,\infty}([0,T]\times\overline{{\Omega}}):\Delta\xi\in L^{\infty}\\}$ which
solves

$\displaystyle-\lambda\Delta\xi$
$\displaystyle=M+\gamma_{1}\lVert\nabla\xi\rVert_{L^{\infty}(\mathds{R}^{d},{\Omega})}+\gamma_{2}\lVert\xi\rVert_{L^{\infty}({\Omega})}$
$\displaystyle\text{on }{\Omega},$ (40a) $\displaystyle\xi$ $\displaystyle=g$
$\displaystyle\text{on }\partial{\Omega},$ (40b)

where

$\gamma_{1}:=\sup_{\begin{subarray}{c}i\in\mathds{N}\\\ (\alpha,\beta)\in
A\times
B\end{subarray}}\lVert\bar{b}_{i}^{(\alpha,\beta)}+\bar{\bar{b}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})},\quad\gamma_{2}:=\sup_{\begin{subarray}{c}i\in\mathds{N}\\\
(\alpha,\beta)\in A\times
B\end{subarray}}\lVert\bar{c}_{i}^{(\alpha,\beta)}+\bar{\bar{c}}_{i}^{(\alpha,\beta)}\rVert_{L^{\infty}({\Omega})}$

and

$\displaystyle\lambda:=\inf_{\begin{subarray}{c}i\in\mathds{N}\\\
(\alpha,\beta)\in A\times
B\end{subarray}}\bar{a}_{i}^{(\alpha,\beta)}+\bar{\bar{a}}_{i}^{(\alpha,\beta)},\quad
M:=\max\Bigl{\\{}$ $\displaystyle\lambda\lVert
v_{T}\rVert_{L^{\infty}({\Omega})},\sup_{i,\alpha,\beta}{\mathsf{F}}_{i}^{(\alpha,\beta)}({y}_{i}^{\ell})+1\Bigr{\\}}.$

The Dirichlet boundary conditions of (40b) are imposed in the strong sense.
Because $g$ does not depend on time it is clear that $\xi$ is constant in $t$.
We shall therefore often write $\xi(x)$ instead of $\xi(t,x)$. The function
$\xi$ is used as barrier for all $y\in\omega$ and $\varepsilon>0$:
$\xi_{y,\varepsilon}:=\xi$.

We assume in this subsection that ${\Omega}$ is convex. Then, as outlined
below Assumption 3.1, we can construct ${P_{i}}$ so that ${P_{i}}\xi$
interpolates $\xi$ on the boundary. Furthermore, we require for the
construction of this section that $\Delta v_{T}\in L^{\infty}({\Omega})$ and
$v_{T}=g$ on $\partial{\Omega}$.

To show that a maximum of

$v_{i}({s}_{i}^{k},\cdot)-P_{i}\xi(x)$

over ${{s}_{i}^{k}}\times\overline{{\Omega}}$ is attained on the boundary for
all ${{s}_{i}^{k}}$ we use that the $h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf Id}$
are M-matrices. Hence it is enough to prove that

$\displaystyle\left((h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf Id})(P_{i}\xi-
v_{i}({s}_{i}^{k},\cdot))\right)_{\ell}\geq 0\qquad\forall\,\ell\leq N_{i}.$
(41)

####  10.1.1 Fully implicit numerical scheme

For the sake of clarity, we begin with a fully implicit numerical scheme of
the form

$(h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf
Id})\,v_{i}({s}_{i}^{k},\cdot)-v_{i}({s}_{i}^{k+1},\cdot)-h_{i}{\mathsf{F}}_{i}^{k,v_{i}}=0.$
(42)

Recall how $(\alpha_{i}^{k,\ell}(v_{i}),\beta_{i}^{k,\ell}(v_{i}))\in A\times
B$ is a pair of controls such that $\alpha_{i}^{k,\ell}(v_{i})$ is a maximiser
of (18a) and $\beta_{i}^{k,\ell}(v_{i})$ is a minimiser of (18b) for
$w=v_{i}$. Noting that
$-\langle\Delta\xi,\hat{\phi}_{i}^{\ell}\rangle=\langle\nabla
P_{i}\xi,\nabla\hat{\phi}_{i}^{\ell}\rangle$ is positive due to (40), we find
that

$\displaystyle\left((h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf
Id})(P_{i}\xi)\right)_{\ell}$ $\displaystyle=\;\;$
$\displaystyle(P_{i}\xi)_{\ell}-h_{i}\,\bar{\bar{a}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}({y}_{i}^{\ell})\,\langle\Delta\xi,\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle+h_{i}\langle-\bar{\bar{b}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}\cdot\nabla
P_{i}\xi+\bar{\bar{c}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}{P_{i}}\xi,\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle=\;\;$
$\displaystyle(P_{i}\xi)_{\ell}-h_{i}\,\bar{\bar{a}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}({y}_{i}^{\ell})\,\langle\Delta\xi,\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle+h_{i}\langle-\bar{\bar{b}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}\cdot\nabla\xi+\bar{\bar{c}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}\xi,\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle+h_{i}\langle-\bar{\bar{b}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}\cdot\nabla(P_{i}\xi-\xi)+\bar{\bar{c}}_{i}^{(\alpha^{k,\ell}(v_{i}),\beta^{k,\ell}(v_{i}))}({P_{i}}\xi-\xi),\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle\geq\;\;$
$\displaystyle(P_{i}\xi)_{\ell}+h_{i}\langle-\lambda\Delta\xi-\gamma_{1}\lVert\nabla\xi\rVert_{L^{\infty}({\Omega})}-\gamma_{2}\lVert\xi\rVert_{L^{\infty}({\Omega})},\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle-
h_{i}\left(\gamma_{1}\lVert\nabla{P_{i}}\xi-\nabla\xi\rVert_{L^{\infty}({\Omega})}+\gamma_{2}\lVert{P_{i}}\xi-\xi\rVert_{L^{\infty}({\Omega})}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{xidef}}}{{=}}$
$\displaystyle(P_{i}\xi)_{\ell}+h_{i}M-h_{i}\left(\gamma_{1}\lVert\nabla{P_{i}}\xi-\nabla\xi\rVert_{L^{\infty}({\Omega})}+\gamma_{2}\lVert{P_{i}}\xi-\xi\rVert_{L^{\infty}({\Omega})}\right).$

Furthermore, there is a $\hat{i}\in\mathds{N}$ such that for all
$i\geq\hat{i}$.

$\sup_{i}\gamma_{1}\lVert\nabla{P_{i}}\xi-\nabla\xi\rVert_{L^{\infty}({\Omega})}+\gamma_{2}\lVert{P_{i}}\xi-\xi\rVert_{L^{\infty}({\Omega})}\leq
1$

Assuming $i\geq\hat{i}$ for the remainder of the section, and using the
definition of $M$

$\left((h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf
Id})(P_{i}\xi)\right)_{\ell}\geq({P_{i}}\xi)_{\ell}+h_{i}\,\sup_{j,\alpha,\beta}{\sf
F}_{j}^{(\alpha,\beta)}(y_{j}^{\ell}).$

Recall that $\xi(x)=\xi({s}_{i}^{k},x)=\xi({s}_{i}^{k+1},x)$. With the
numerical scheme (42) we obtain

$\displaystyle\left((h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf Id})(P_{i}\xi-
v_{i}({s}_{i}^{k},\cdot))\right)_{\ell}$
$\displaystyle\geq(P_{i}\xi)_{\ell}-v_{i}({s}_{i}^{k+1},{y}_{i}^{\ell}).$ (43)

Because both $(P_{i}\xi)_{\ell}$ and $v_{i}({s}_{i}^{k},{y}_{i}^{\ell})$
interpolate $g$ on $\partial{\Omega}$, it follows $P_{i}\xi-
v_{i}({s}_{i}^{k},\cdot)\in V_{i}^{0}$. Owing to the M-matrix property of
$h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf Id}$,

$\displaystyle(P_{i}\xi)_{\ell}-v_{i}({s}_{i}^{k+1},{y}_{i}^{\ell})\geq 0$
(44)

thus implies condition (41). Hence, by induction, (41) holds for all $k$ as
soon as (44) is shown at the final time ${s}_{i}^{k+1}=T$.

From (40) we know that $-\lambda\Delta\xi\geq M$. Then, using
$v_{i}(T,\cdot)={P_{i}}v_{T}$ and $\lambda>0$,

$\displaystyle\langle\nabla\left({P_{i}}\xi^{T}-v_{i}(T,\cdot)\right),\nabla\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle=-\langle\Delta(\xi-v_{T}),\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle\geq\langle M/\lambda+\Delta
v_{T},\hat{\phi}_{i}^{\ell}\rangle\geq M/\lambda-\lVert\Delta
v_{T}\rVert_{L^{\infty}({\Omega})}.$

The definition of $M$ ensures that $\frac{M}{\lambda}\geq\lVert\Delta
v_{T}\rVert_{L^{\infty}({\Omega})}$. Using M-matrix property of the discrete
Laplacian formed with respect to the nodal basis
$\\{\hat{\phi}_{i}^{\ell}\\}_{\ell}$ we obtain (44) at the final time.

At this point we proved that the maximum of
$v_{i}({s}_{i}^{k},\cdot)-P_{i}\xi$ lies on the boundary as required. It
remains to show that for $t\in[0,T]$ we can select a maximiser
$r_{t,y,\epsilon}$ of $g-\xi_{y,\epsilon}$ over
$\\{t\\}\times\partial{\Omega}$ such that $\lim_{\epsilon\rightarrow 0}\
r_{t,y,\epsilon}=y$ and $\lim_{\epsilon\rightarrow
0}\xi_{y,\epsilon}(t,y)-\xi_{y,\epsilon}\ (t,r_{t,y,\epsilon})\leq 0$. Because
$\xi$ attains $g$ on whole the boundary, $r_{t,y,\epsilon}:=y$ is already such
a choice.

####  10.1.2 IMEX numerical schemes

We now extend the above argument to general numerical schemes as defined in
(19). Choosing $\gamma_{1},\gamma_{2},\lambda$ and $M$ as above, the
inequality (43) generalises in the IMEX setting to

$\displaystyle\left((h_{i}{\mathsf{I}}_{i}^{k,v_{i}}+{\sf Id})(P_{i}\xi-
v_{i}({s}_{i}^{k},\cdot))\right)_{\ell}$
$\displaystyle\geq-\left(\\!(h_{i}{\mathsf{E}}_{i}^{k,v_{i}}-{\sf
Id})(P_{i}\xi-v_{i}({s}_{i}^{k+1},\cdot))\\!\right)_{\ell}.$

According to Lemma 4.1 we have $-(h_{i}{\mathsf{E}}_{i}^{k,v_{i}}-{\sf
Id})\geq 0$. Therefore, if

$P_{i}\xi-v_{i}({s}_{i}^{k+1},\cdot))\geq 0$

implies (41). At this point the induction of the previous subsection can be
adapted to show that maxima are attained on the boundary. Setting
$r_{t,y,\epsilon}:=y$ completes the construction.

###  10.2 Barrier functions for degenerate equations and general domains

We now want to remove the two main assumptions of section 10.1, namely the
requirements that the differential operator is uniformly parabolic and that
the domain is convex.

We form $\omega$ of the $y$ for which it is possible to ensure the existence
of strict supersolutions in the following sense: for all $\varepsilon>0$ one
can find a $\xi_{y,\varepsilon}$ satisfying

$\displaystyle\sup_{\alpha}-a^{(\alpha,\beta)}\,\Delta\xi_{y,\varepsilon}-b^{(\alpha,\beta)}\\!\cdot\\!\nabla\xi_{y,\varepsilon}$
$\displaystyle+c^{(\alpha,\beta)}\,\xi_{y,\varepsilon}-f^{(\alpha,\beta)}\geq\varepsilon$
$\displaystyle\text{on }{\Omega}$ (45)

for all $\beta\in C({\Omega};A)$ as well as

$\sup_{\varepsilon>0}(\|\Delta\xi_{y,\varepsilon}\|+\left\|\xi_{y,\varepsilon}\right\|_{W^{1,\infty}({\Omega})})<\infty$

and

$\displaystyle v_{T}-\xi_{y,\varepsilon}$ $\displaystyle\leq-2\varepsilon$
$\displaystyle\text{on }\overline{{\Omega}}\setminus
B_{y}(\delta(\varepsilon)),$ (46a) $\displaystyle v_{T}-\xi_{y,\varepsilon}$
$\displaystyle\leq-\phantom{2}\varepsilon$ $\displaystyle\text{on
}\overline{{\Omega}}\cap B_{y}(\delta(\varepsilon)),$ (46b) $\displaystyle
v_{T}(y)-\xi_{y,\varepsilon}(y)$ $\displaystyle>-2\varepsilon,$ (46c)

where $B_{y}(\delta(\varepsilon))$ is the ball centred at $y$ with radius
$\delta(\varepsilon)>0$, which in turn is a positive parameter depending on
$\varepsilon$. Observe that $v_{T}|_{\partial{\Omega}}=g|_{\partial{\Omega}}$
and therefore (46) also provides control on the boundary, due to $g$ and
$\xi_{y,\varepsilon}$ being time independent. It ensures that the maximum of
$y-\xi_{y,\varepsilon}$ is attained in the vicinity of $y$ and is non-
positive. We can view (45) as generalisation of (40a) because $\xi$ of (40a)
is a strict supersolution of $L^{(\alpha,\beta)}w-f^{(\alpha,\beta)}=0$ for
all $\alpha\in A$, $\beta\in B$.

In light of Assumptions 3.1 and 3.2, we may choose $\hat{i}$ such that

$\displaystyle
C_{1}(\hat{i})\,(\|\Delta\xi_{y,\varepsilon}\|+\left\|\xi_{y,\varepsilon}\right\|_{W^{1,\infty}({\Omega})})+C_{2}\sup_{i\geq\hat{i}}\left\|{P_{i}}\xi_{y,\varepsilon}-\xi_{y,\varepsilon}\right\|_{W^{1,\infty}({\Omega})}\leq\varepsilon.$
(47)

where

$\displaystyle C_{1}(\hat{i}):=\sup_{\begin{subarray}{c}(\alpha,\beta)\in
A\times B\\\ i\geq\hat{i}\end{subarray}}$
$\displaystyle\bigl{(}\sup_{\ell}\bigl{\|}a^{(\alpha,\beta)}-\bigl{(}\bar{a}_{i}^{(\alpha,\beta)}({y}_{i}^{\ell})+\bar{\bar{a}}_{i}^{(\alpha,\beta)}({y}_{i}^{\ell})\bigr{)}\bigr{\|}_{L^{\infty}({\rm
supp}\,\hat{\phi}_{i}^{\ell})}$
$\displaystyle+\bigl{\|}b^{(\alpha,\beta)}-\bigl{(}\bar{b}_{i}^{(\alpha,\beta)}+\bar{\bar{b}}_{i}^{(\alpha,\beta)}\bigr{)}\bigr{\|}_{L^{\infty}({\Omega},\mathds{R}^{d})}$
$\displaystyle+\bigl{\|}c^{(\alpha,\beta)}-\bigl{(}\bar{c}_{i}^{(\alpha,\beta)}+\bar{\bar{c}}_{i}^{(\alpha,\beta)}\bigr{)}\bigr{\|}_{L^{\infty}({\Omega})}+\bigl{\|}f^{(\alpha,\beta)}-{f}_{i}^{(\alpha,\beta)}\bigr{\|}_{L^{\infty}({\Omega})}\bigr{)}$

and

$\displaystyle
C_{2}:=\max\bigl{\\{}1,\sup_{\begin{subarray}{c}(\alpha,\beta)\\\
i\in\mathds{N}\end{subarray}}\|\bar{a}_{i}^{(\alpha,\beta)}+\bar{\bar{a}}_{i}^{(\alpha,\beta)}\|_{\infty}+\|\bar{b}_{i}^{(\alpha,\beta)}+\bar{\bar{b}}_{i}^{(\alpha,\beta)}\|_{\infty}+\|\bar{c}_{i}^{(\alpha,\beta)}+\bar{\bar{c}}_{i}^{(\alpha,\beta)}\|_{\infty}\bigr{\\}}.$

Recall the definition of $(\alpha_{i}^{k,\ell}(w),\beta_{i}^{k,\ell}(w))\in
A\times B$ in (18). For the sake of readability, we write $\hat{\beta}$ in
place of $\beta_{i}^{k,\ell}(P_{i}\xi_{y,\varepsilon}-v_{i})$ and
$\hat{\alpha}$ in place of
$\alpha_{i}^{k,\ell,\hat{\beta}}(P_{i}\xi_{y,\varepsilon}-v_{i})$ in this
section, i.e. $\hat{\alpha}$ and $\hat{\beta}$ are optimal choices when
evaluating the numerical operator at $P_{i}\xi_{y,\varepsilon}-v_{i}$.

Then the numerical method, applied to $P_{i}\xi^{k+1}-v_{i}$ with the control
of the infimum frozen at $P_{i}\xi_{y,\varepsilon}-v_{i}$, returns at the node
${y}_{i}^{\ell}$

$\displaystyle\bigl{(}(h_{i}{\mathsf{I}}_{i}^{(\hat{\alpha},\hat{\beta})}+{\sf
Id})(P_{i}\xi_{y,\varepsilon}-v_{i}({s}_{i}^{k},\cdot))+(h_{i}{\mathsf{E}}_{i}^{(\hat{\alpha},\hat{\beta})}-{\sf
Id})(P_{i}\xi_{y,\varepsilon}-v_{i}({s}_{i}^{k+1},\cdot))\bigr{)}_{\ell}$
$\displaystyle=\,$ $\displaystyle
h_{i}\bigl{(}-(\bar{a}_{i}^{(\hat{\alpha},\hat{\beta})}({y}_{i}^{\ell})+\bar{\bar{a}}_{i}^{(\hat{\alpha},\hat{\beta})}({y}_{i}^{\ell}))\langle\Delta\xi_{y,\varepsilon},\hat{\phi}_{i}^{\ell}\rangle-\langle(\bar{b}_{i}^{(\hat{\alpha},\hat{\beta})}+\bar{\bar{b}}_{i}^{(\hat{\alpha},\hat{\beta})})\cdot\nabla{P_{i}}\xi_{y,\varepsilon},\hat{\phi}_{i}^{\ell}\rangle$
$\displaystyle\phantom{+\bigl{(}}+\langle(\bar{c}_{i}^{(\hat{\alpha},\hat{\beta})}+\bar{\bar{c}}_{i}^{(\hat{\alpha},\hat{\beta})}){P_{i}}\xi_{y,\varepsilon}-{f}_{i}^{(\hat{\alpha},\hat{\beta})},\hat{\phi}_{i}^{\ell}\rangle\bigr{)}$
$\displaystyle\geq\,$
$\displaystyle\;h_{i}\,\langle-{a}^{(\hat{\alpha},\hat{\beta})}\Delta\xi_{y,\varepsilon}-{b}^{(\hat{\alpha},\hat{\beta})}\cdot\nabla\xi_{y,\varepsilon}+{c}^{(\hat{\alpha},\hat{\beta})}\xi_{y,\varepsilon}-{f}^{(\hat{\alpha},\hat{\beta})},\hat{\phi}_{i}^{\ell}\rangle-
h_{i}\,\varepsilon$ $\displaystyle\geq\,$ $\displaystyle 0.$

We conclude with Lemma 4.1 that

$\displaystyle P_{i}\xi_{y,\varepsilon}-v_{i}({s}_{i}^{k+1},\cdot)\geq 0$ (48)

implies, for $\ell\leq N_{i}$,

$\displaystyle\bigl{(}(h_{i}{\mathsf{I}}_{i}^{(\hat{\alpha},\hat{\beta})}+{\sf
Id})(P_{i}\xi_{y,\varepsilon}-v_{i}({s}_{i}^{k},\cdot))\bigr{)}_{\ell}$
$\displaystyle\geq\bigl{(}-(h_{i}{\mathsf{E}}_{i}^{(\hat{\alpha},\hat{\beta})}-{\sf
Id})(P_{i}\xi_{y,\varepsilon}-v_{i}({s}_{i}^{k+1},\cdot))\bigr{)}_{\ell}\geq
0.$

Because unlike in the case of subsection 10.1 the functions
$P_{i}\xi_{y,\varepsilon}-v_{i}({s}_{i}^{k},\cdot)$ do not vanish on the
boundary we need to extend this result to the case when $N_{i}<\ell\leq\dim
V_{i}$. Owing to (47) and (46b) we know that (48) also holds on the boundary
for all time steps ${s}_{i}^{k+1}$. Furthermore, $\eqref{xiyeBCb}$ implies
that (48) is satisfied on all of $\overline{{\Omega}}$ at the final time
${s}_{i}^{k+1}=T$. In summary we have the induction step that if (48) on
$\overline{{\Omega}}$ then
$P_{i}\xi_{y,\varepsilon}-v_{i}({s}_{i}^{k},\cdot)\geq 0$ on
$\overline{{\Omega}}$ and the induction base to guarantee that the maximum of

$v_{i}({s}_{i}^{k},\cdot)-P_{i}\xi(x)$

over ${{s}_{i}^{k}}\times\overline{{\Omega}}$ is attained on the boundary for
all ${{s}_{i}^{k}}$.

It remains to show that for $t\in[0,T]$ we can select a maximiser
$r_{t,y,\epsilon}$ of $g-\xi_{y,\epsilon}$ over
$\\{t\\}\times\partial{\Omega}$ such that $\lim_{\epsilon\rightarrow 0}\
r_{t,y,\epsilon}=y$ and $\lim_{\epsilon\rightarrow
0}\xi_{y,\epsilon}(t,y)-\xi_{y,\epsilon}\ (t,r_{t,y,\epsilon})\leq 0$. The
former holds because of (46a) and (10.2), while the latter follows from (46b).

Similarly we assume the existence of strict subsolutions in the following
sense: For all $y\in\omega$ and $\varepsilon>0$ one can find a
$\zeta_{y,\varepsilon}$ satisfying

$\inf_{\beta}-a^{(\alpha,\beta)}\,\Delta\zeta_{y,\varepsilon}-b^{(\alpha,\beta)}\\!\cdot\\!\nabla\zeta_{y,\varepsilon}+c^{(\alpha,\beta)}\,\zeta_{y,\varepsilon}-f^{(\alpha,\beta)}\leq-\varepsilon\qquad\text{on
}{\Omega}$

for all $\alpha\in C({\Omega};A)$ as well as

$\sup_{\varepsilon>0}(\|\Delta\zeta_{y,\varepsilon}\|+\left\|\zeta_{y,\varepsilon}\right\|_{W^{1,\infty}({\Omega})})<\infty$

and

$\displaystyle v_{T}-\zeta_{y,\varepsilon}$ $\displaystyle\geq 2\varepsilon$
$\displaystyle\text{on }\overline{{\Omega}}\setminus
B_{y}(\delta(\varepsilon)),$ $\displaystyle v_{T}-\zeta_{y,\varepsilon}$
$\displaystyle\geq\phantom{2}\varepsilon$ $\displaystyle\text{on
}\overline{{\Omega}}\cap B_{y}(\delta(\varepsilon)),$ $\displaystyle
v_{T}(y)-\zeta_{y,\varepsilon}(y)$ $\displaystyle<2\varepsilon.$

##  11 Numerical experiments

In this section we present two numerical experiments showing the viability of
the presented method. In the first experiment, we analyse convergence rates of
a fully nonlinear second order Isaacs problem with a known solution and we
confirm at least linear convergence in $L^{2}$, $L^{\infty}$ and $H^{1}$
norms. In the second experiment we calculate the value function of a
stochastic tag-chase game with asymmetric velocities and vanishing diffusion
on a non-convex domain. The presented finite element method was implemented in
Python with FEniCS [LMW+12] and is available from the public repository
[Jar21] under the GNU Lesser General Public License.

###  11.1 Isaacs problem with exact solution

Let the spatial domain ${\Omega}$ in $\mathds{R}^{2}$ be the equilateral
triangle with vertices $(\pm\sqrt{3},\ \frac{1}{2})$ and $(0,\\-1)$. We will
study the following Isaacs problem:

$-v_{t}-\inf_{\beta\in[\frac{1}{4},\
\frac{1}{2}]}\Bigl{\\{}\beta\sqrt{\frac{x^{2}+y^{2}}{T-t+1}}\Delta
v+\frac{1}{2}\frac{1}{\sqrt{T-t+1}}\left|\nabla
v\right|\Bigr{\\}}=-\frac{1}{2}\frac{\sqrt{x^{2}+y^{2}}}{(T-t+1)^{3/2}}.$ (49)

Equation (49) is indeed of Isaacs type because the Euclidean norm of the
gradient may be written alternatively as

$\left|\nabla
v\right|=\sup_{\\{\alpha\in\mathds{R}^{2}:\left|\alpha\right|=1\\}}\\{\alpha\cdot\nabla
v\\}.$

One can now verify through direct calculation that the function

$\displaystyle
v(x,y,t)=\exp\left(-\sqrt{\frac{x^{2}+y^{2}}{T-t+1}}\right)+\sqrt{\frac{x^{2}+y^{2}}{T-t+1}}$
(50)

solves (49) exactly. The final and boundary data is also given by the right-
hand side of (50).

We split the scheme so that the advection term is treated explicitly. Note
that in order to ensure the monotonicity of the scheme we may also assign part
of the diffusion to the explicit operator. In order to improve the rate of
convergence this was done locally, i.e. on a node-wise basis. In case the
naturally occurring diffusion at the node is not sufficient, we introduce
artificial diffusion. Overall, this approach leads to artificial diffusion
near the origin where the differential operator is degenerately elliptic,
while for the majority of nodes of the mesh it is zero. We choose the largest
timestep guaranteeing the monotonicity of the method which leads to
$O(h_{i})=O(\Delta x_{i})$. Any amount of the diffusion left after ensuring
the monotonicity of the explicit operators is treated implicitly.

The convergence of the scheme is reflected on Figure 1 which plots errors at
$t=0$ for different mesh sizes. The rates of convergence are as follows:

$\Delta x$ | $L^{\infty}$ | Rate | $L^{2}$ | Rate | $H^{1}$ | Rate
---|---|---|---|---|---|---
0.4330 |  1.364e-02 | 0.66 |  7.062e-03 | 0.96 |  6.792e-02 | 0.91
0.2165 |  1.040e-02 | 0.91 |  3.695e-03 | 0.78 |  3.737e-02 | 0.98
0.1083 |  5.715e-03 | 1.01 |  2.361e-03 | 0.93 |  1.899e-02 | 1.01
0.0541 |  2.838e-03 | 1.07 |  1.266e-03 | 1.02 |  9.391e-03 | 1.03
0.0271 |  1.330e-03 | 1.10 |  6.234e-04 | 1.06 |  4.581e-03 | 1.03
0.0135 |  6.034e-04 | 1.11 |  2.941e-04 | 1.07 |  2.223e-03 | 1.03
0.0068 | 2.708e-04 |  | 1.374e-04 |  | 1.082e-03 |

$10^{1}$$10^{2}$$10^{-4}$$10^{-3}$$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

###  11.2 Tag-Chase Game with random noise

Imagine two players moving on the $R^{2}$ plane. One player is the pursuer (we
will denote them by P) and tries to catch the other one, who is the evader (we
will denote them by E). Both of them are allowed to choose their direction
freely: $\alpha,\beta\in[-\pi,\pi]$ and $\alpha$ is the direction chosen by
the evader, while $\beta$ is the direction chosen by the pursuer. The pursuer
moves with speed $v^{\beta}_{P}$ while the evader moves with speed
$v^{\alpha}_{E}$. In particular, the respective speeds are functions of the
angles. Additionally, some choices of a direction for the pursuer and the
evader may be subject to a random noise behaving like a standard Brownian
motion. Having specified the setting we are able to formulate the dynamics
explicitly as

$\displaystyle dx_{1}^{P}$
$\displaystyle=(v^{\beta}_{P})_{1}\,\sin(\beta)dt+\sigma^{\beta}_{P}(x_{1}^{P})\,dW_{(1,P)},$
$\displaystyle dx_{2}^{P}$
$\displaystyle=(v^{\beta}_{P})_{2}\,\cos(\beta)dt+\sigma^{\beta}_{P}(x_{2}^{P})\,dW_{(2,P)},$
$\displaystyle dx_{1}^{E}$
$\displaystyle=(v^{\alpha}_{E})_{1}\,\sin(\alpha)dt+\sigma^{\alpha}_{E}(x_{1}^{E})\,dW_{(1,E)},$
$\displaystyle dx_{2}^{E}$
$\displaystyle=(v^{\alpha}_{E})_{2}\,\cos(\alpha)dt+\sigma^{\alpha}_{E}(x_{2}^{E})\,dW_{(2,E)}.$

Additionally, we assume that

$\displaystyle W_{P}$ $\displaystyle=(W_{(1,P)},W_{(2,P)}),$ $\displaystyle
W_{E}$ $\displaystyle=(W_{(1,E)},W_{(2,E)}),$

are two $R^{2}$-valued, mutually independent standard Wiener processes. We
reduce this $4$-dimensional problem to a $2$-dimensional one by allowing the
origin of the coordinate system to move along with the pursuer. In this case
our dynamics are

$\displaystyle d\hat{x}_{1}$
$\displaystyle=((v^{\beta}_{P})_{1}\,\sin(\beta)-(v^{\alpha}_{E})_{1}\,\sin(\alpha))\,dt+\sigma^{(\alpha,\beta)}(\hat{x}_{1})\,d\hat{W}_{1},$
$\displaystyle d{\hat{x}_{2}}$
$\displaystyle=((v^{\beta}_{P})_{2}\,\cos(\beta)-(v^{\alpha}_{E})_{2}\,\cos(\alpha))\,dt+\sigma^{(\alpha,\beta)}(\hat{x}_{2})\,d\hat{W}_{2},$

where
$\sigma^{(\alpha,\beta)}(\hat{x}_{i})=\sqrt{\sigma^{\beta}_{P}(x^{P}_{i}))^{2}+(\sigma^{\alpha}_{E}(x^{E}_{i}))^{2}}$
and $\hat{W}$ is an $R^{2}$-valued standard Wiener process satisfying

$\sigma^{(\alpha,\beta)}(\hat{x}_{i})\hat{W}_{i}(t)=\sigma^{\beta}_{P}(x_{i}^{P})W_{i,P}(t)-\sigma^{\alpha}_{E}(x_{i}^{E})W_{i,E}(t),\quad
t\geq 0,\;i=1,2.$

The pursuer catches the evader (and thus wins the game) if they manage to
reduce their distance from the evader to some value $r\in R$. The evader wins
the game when they manage to increase the distance to pursuer to some given
$R>r$ or if they manage to avoid the capture before some time $T$. Note that
in this case the spatial domain of the problem becomes
$\overline{{\Omega}}:=\overline{B}(0,R)\setminus B(0,r)$.

Mathematically, the evader receives the pay-out of $1$ whenever they win the
game and receive $0$ otherwise. We write for the expected pay-out to the
evader $\mathcal{J}(\vec{\alpha},\vec{\beta},x,t)$, where
$\vec{\alpha}:[0,T]\to A$, $\vec{\beta}:[0,T]\to B$ are functions of time so
that $\vec{\alpha}(t)$, $\vec{\beta}(t)$ are the controls chosen by the evader
and pursuer at time $t$, respectively.

The value function

$v(x,t):=\adjustlimits{\inf}_{\beta}{\sup}_{\alpha}\mathcal{J}(\alpha,\beta,x,t)$

then solves the second-order Isaacs equation

$\displaystyle-\partial_{t}v+\inf_{\beta}\sup_{\alpha}\big{(}-a^{(\alpha,\beta)}\Delta
v-b^{(\alpha,\beta)}\cdot\nabla v\big{)}=0\qquad$ $\displaystyle\text{on
}[0,T)\times{\Omega},$ $\displaystyle v=1\qquad$ $\displaystyle\text{on
}[0,T)\times\\{x:\|x\|=R\\},$ $\displaystyle v=0\qquad$ $\displaystyle\text{on
}[0,T)\times\\{x:\|x\|=r\\},$ $\displaystyle v=1\qquad$ $\displaystyle\text{on
}\\{T\\}\times\overline{{\Omega}},$

where

$a^{(\alpha,\beta)}:=\frac{(\sigma^{(\alpha,\beta)})^{2}}{2},\qquad
b^{\alpha,\beta}:=\begin{pmatrix}(v_{P}^{\beta})_{1}\,\sin(\beta)-(v_{E}^{\alpha})_{1}\,\sin(\alpha)\\\
(v_{P}^{\beta})_{2}\,\cos(\beta)-(v_{E}^{\alpha})_{2}\,\cos(\alpha)\end{pmatrix}.$

We shall assume that the pursuer is faster than the evader when moving in the
horizontal direction. Moreover, we assume that the diffusion in the upper part
of the domain scales with the vertical position $x_{2}$, while it is constant
where $x_{2}\leq 0.1$. Specifically, this corresponds to the following choice
of the diffusion and advection coefficients:

$\displaystyle a^{(\alpha,\beta)}=\max\\{x_{2},0.1\\},\qquad
b^{(\alpha,\beta)}$
$\displaystyle=\begin{pmatrix}4\sin(\beta)-0.5\sin(\alpha)\\\
\cos(\beta)-\cos(\alpha)\end{pmatrix}.$

The numerical approximation of the value function is displayed in Figure 2.
Note the asymmetric nature of the graph, due to the different speeds of the
pursuer and the evader in the horizontal direction as well as a larger effect
of the stochastic component of the equation in the upper part of the domain.

Figure 2: The value function at $t=0$ for asymmetric velocities and the
diffusion having a lower value in the lower part of ${\Omega}$.

## Declarations

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

Conflicts of interest: The authors declare that they have no competing
interest.

Code availability: At https://doi.org/10.5281/zenodo.4598310.

## References

  * [BCJ19] Imran H Biswas, Indranil Chowdhury, and Espen Robstad Jakobsen.  On the rate of convergence for monotone numerical schemes for nonlocal Isaacs equations.  SIAM Journal on Numerical Analysis, 57(2):799–827, 2019.
  * [BHT11] Nikolai D Botkin, Karl-Heinz Hoffmann, and Varvara L Turova.  Stable numerical schemes for solving Hamilton–Jacobi–Bellman–Isaacs equations.  SIAM Journal on Scientific Computing, 33(2):992–1007, 2011.
  * [BM98] Randal W. Beard and Timothy W. McLain.  Successive Galerkin approximation algorithms for nonlinear optimal and robust control.  International Journal of Control, 71(5):717–743, 1998.
  * [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.
  * [BP88] Guy Barles and Benoît Perthame.  Exit Time Problems in Optimal Control and Vanishing Viscosity Method.  SIAM Journal on cControl and Optimization, 26(5):1133–1148, September 1988.
  * [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\.
  * [DJ13] Kristian Debrabant and Espen R. Jakobsen.  Semi-Lagrangian schemes for linear and fully non-linear diffusion equations.  Mathematics of Computation, 82(283):1433–1462, 2013.
  * [DLSW12] Alan Demlow, Dmitriy Leykekhman, Alfred H. Schatz, and Lars B. Wahlbin.  Best approximation property in the $W^{1}_{\infty}$ norm for finite element methods on graded meshes.  Mathematics of Computation, 81(278):743–764, 2012.
  * [Fal06] Maurizio Falcone.  Numerical methods for differential games based on partial differential equations.  International Game Theory Review, 8(2):231–272, 2006.
  * [FK14] Maurizio Falcone and Dante Kalise.  A high-order semi-Lagrangian/finite volume scheme for Hamilton–Jacobi–Isaacs equations.  In Christian Pötzsche, Clemens Heuberger, Barbara Kaltenbacher, and Franz Rendl, editors, System Modeling and Optimization. CSMO 2013. IFIP Advances in Information and Communication Technology, volume 443, pages 105–117. Springer, Berlin, Heidelberg, 2014.
  * [FRS10] Henrique C Ferreira, Paulo H Rocha, and Roberto M Sales.  On the convergence of successive Galerkin approximation for nonlinear output feedback $H_{\infty}$ control.  Nonlinear Dynamics. An International Journal of Nonlinear Dynamics and Chaos in Engineering Systems, 60(4):651–660, 2010.
  * [FS06] Wendell H Fleming and Halil Mete Soner.  Controlled Markov Processes and Viscosity Solutions, volume 25 of Stochastic Modelling and Applied Probability.  Springer, New York, second edition, 2006.
  * [Jak04] Espen Robstad Jakobsen.  On error bounds for approximation schemes for non-convex degenerate elliptic equations.  BIT Numerical Mathematics, 44(2):269–285, 2004.
  * [Jak06] Espen Robstad Jakobsen.  On error bounds for monotone approximation schemes for multi-dimensional Isaacs equations.  Asymptotic Analysis, 49(3-4):249–273, 2006.
  * [Jar21] Bartosz Jaroszkowski.  FEISol https://doi.org/10.5281/zenodo.4598310, 2021.
  * [JK02] Espen Robstad Jakobsen and Kenneth H Karlsen.  Continuous dependence estimates for viscosity solutions of fully nonlinear degenerate parabolic equations.  Journal of Differential Equations, 183(2):497–525, August 2002\.
  * [JK05] Espen Robstad Jakobsen and Kenneth H Karlsen.  Continuous dependence estimates for viscosity solutions of integro-PDEs.  Journal of Differential Equations, 212(2):278–318, May 2005.
  * [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.
  * [JS18] Max Jensen and Iain Smears.  On the notion of boundary conditions in comparison principles for viscosity solutions.  In Hamilton–Jacobi–Bellman equations, pages 143–154. De Gruyter, Berlin, 2018.
  * [KKK20] Dante Kalise, Sudeep Kundu, and Karl Kunisch.  Robust feedback control of nonlinear PDEs by numerical approximation of high-dimensional Hamilton–Jacobi–Isaacs equations.  SIAM Journal on Applied Dynamical Systems, 19(2):1496–1524, 2020\.
  * [Kry15] Nicolai V Krylov.  To the theory of viscosity solutions for uniformly parabolic Isaacs equations.  Methods and Applications of Analysis, 22(3):259–280, 2015.
  * [LMW+12] Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al.  Automated Solution of Differential Equations by the Finite Element Method.  Springer, 2012.
  * [LN18] Wenbo Li and Ricardo H Nochetto.  Optimal pointwise error estimates for two-scale methods for the Monge-Ampére equation.  SIAM Journal on Numerical Analysis, 56(3):1915–1941, January 2018\.
  * [MZB06] Stefania Maroso, Housnaa Zidani, and J Frederic Bonnans.  Error estimates for stochastic differential games: the adverse stopping case.  IMA Journal of Numerical Analysis, 26(1):188–212, 2006.
  * [Sor99] Pierpaolo Soravia.  Equivalence between nonlinear $H_{\infty}$ control problems and existence of viscosity solutions of Hamilton–Jacobi–Isaacs equations.  Applied Mathematics and Optimization, 39:17–32, 1999.
  * [Sou99] Panagiotis E. Souganidis.  Two-player, zero-sum differential games and viscosity solutions.  In Martino Bardi, T.E.S. Raghavan, and Thiruvenkatachari Parthasarathy, editors, Stochastic and Differential Games. Annals of the International Society of Dynamic Games, volume 4, chapter 2, pages 69–104. Birkhäuser, Boston, MA, 1999.
  * [SZ19] Abner J. Salgado and Wujun Zhang.  Finite element approximation of the Isaacs equation.  ESAIM Mathematical Modelling and Numerical Analysis, 53(2):351–374, 2019.
  * [Tur15] Olga Turanova.  Error estimates for approximations of nonlinear uniformly parabolic equations.  Nonlinear Differential Equations and Applications NoDEA, 22(3):345–389, June 2015.