# Numerical viscosity solutions to Hamilton-Jacobi equations via a Carleman
estimate and the convexification method

Michael Klibanov Department of Mathematics and Statistics, University of North
Carolina at Charlotte, Charlotte, NC, 28223, USA<EMAIL_ADDRESS>Loc H.
Nguyen Department of Mathematics and Statistics, University of North Carolina
at Charlotte, Charlotte, NC, 28223, USA<EMAIL_ADDRESS>Hung V. Tran
Department of Mathematics, University of Wisconsin Madison, Madison, WI,
53706, USA<EMAIL_ADDRESS>(corresponding author).

###### Abstract

We propose a globally convergent numerical method, called the convexification,
to numerically compute the viscosity solution to first-order Hamilton-Jacobi
equations through the vanishing viscosity process where the viscosity
parameter is a fixed small number. By convexification, we mean that we employ
a suitable Carleman weight function to convexify the cost functional defined
directly from the form of the Hamilton-Jacobi equation under consideration.
The strict convexity of this functional is rigorously proved using a new
Carleman estimate. We also prove that the unique minimizer of the this
strictly convex functional can be reached by the gradient descent method.
Moreover, we show that the minimizer well approximates the viscosity solution
of the Hamilton-Jacobi equation as the noise contained in the boundary data
tends to zero. Some interesting numerical illustrations are presented.

Key words: numerical methods; convexification method; gradient descent method;
viscosity solutions; Hamilton-Jacobi equations; vanishing viscosity process;
boundary value problems.

AMS subject classification: 35D40, 35F30, 65N06.

##  1 Introduction

The aim of this paper is to compute viscosity solutions to a large class of
Hamilton-Jacobi equations possibly involving nonconvex Hamiltonians. The key
ingredient for us to reach this achievement is the use of a new Carleman
estimate and the convexification method. This method is only applicable when a
viscosity term is added to the Hamilton-Jacobi equation under consideration.
The idea of adding the viscosity term and passing to the limit to obtain
viscosity solutions is due to the seminal works [15, 14]. Let $R>0$ and
$\Omega=(-R,R)^{d}$ where $d\geq 1$ is the spatial dimension. Let
$F:\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$ and
$f:\overline{\Omega}\to\mathbb{R}$ be functions of the class $C^{2}$. In this
paper, we propose a globally convergent numerical method to solve the
Hamilton-Jacobi equation

$F({\bf x},u({\bf x}),\nabla u({\bf x}))=0\quad\mbox{for all }{\bf
x}\in\Omega$ (1.1)

with the Dirichlet boundary condition

$u({\bf x})=f({\bf x})\quad\mbox{for all }{\bf x}\in\partial\Omega.$ (1.2)

The smoothness condition imposing on $F$ and $f$ is for the simplicity. We
need it when analytically establishing the convergence of the proposed method.
However, this technical condition can be relaxed in our numerical study. In
this paper, we are interested in computing the viscosity solution to
(1.1)–(1.2). We only deal with the case that the Dirichlet boundary condition
holds in the classical sense in this paper. We refer the readers to [15, 14,
40, 5, 4, 56] and the references therein for the theory of viscosity solutions
to (1.1)–(1.2). It is worth mentioning that a number of different extremely
efficient and fast numerical approaches and techniques (many of which are of
high orders) have been developed for Hamilton-Jacobi equations. For finite
difference monotone and consistent schemes of first-order equations and
applications, see [16, 53, 6, 50, 45] for details and recent developments. If
$F=F({\bf x},s,{\bf p})$ is convex in ${\bf p}$ and satisfies some appropriate
conditions, it is possible to construct some semi-Lagrangian approximations by
the discretization of the Dynamical Programming Principle associated to the
problem, see [18, 19] and the references therein. For a non-exhaustive list of
results along these directions, see [46, 47, 58, 1, 2, 7, 51, 57, 23, 49, 11,
13, 10, 44, 39, 21]. Another approach to solve (1.1)–(1.2) is based on
optimization [22, 54, 38, 17]. However, due to the nonlinearity of the
function $F$ in (1.1), the least squares cost functional is nonconvex and
might have multiple local minima and ravines. Hence, the methods based on
optimization only provide reliable numerical solutions if good initial guesses
of the true solutions are given.

Unlike the mentioned optimization approach, we propose to use the
convexification method, which does not rely on the above assumptions. This
method is globally convergent in the sense that

  1. 1.

it delivers good approximation of the true solution without knowing any
advance knowledge of the true solution even when the given data is noisy;

  2. 2.

the claim in # 1 above is rigorously proved and numerically verified.

To the best of our knowledge, our method is new in the context of viscosity
solutions to Hamilton-Jacobi equations. It has two major advantages in solving
(1.1)–(1.2) numerically. Firstly, it works for some quite general $F({\bf
x},u,\nabla u)$, which might be nonconvex in $\nabla u$, and it does not
require a lot of specific structures on $F$. In particular, $F$ might be
dependent on $u$ and $\nabla u$ in a rather complicated way. Secondly, it is
quite stable and robust even with some noise level on the boundary data, which
occurs naturally in applications.

The main idea of the convexification method is to employ a suitable Carleman
weight function to convexify the mismatch functional derived from the given
boundary value problem. Several versions of the convexification method have
been developed since it was first introduced in [30] for a coefficient inverse
problem for a hyperbolic equation. We cite here [28, 26, 29, 27, 3, 25, 24,
52, 34] and references therein for some important works in this area and their
real-world applications in bio-medical imaging, non-destructed testing, travel
time tomography, identifying anti-personnel explosive devices buried under the
ground, etc. The crucial mathematical ingredients that guarantee the strict
convexity of this functional are the Carleman estimates. The original idea of
applying Carleman estimates to prove the uniqueness for a large class of
important nonlinear mathematical problems was first published in [8]. It was
discovered later in [30, 33], that the idea of [8] can be successfully
modified to develop globally convergent numerical methods for coefficient
inverse problems using the convexification.

In this paper, it is the first time we use a Carleman weight function to
numerically solve Hamilton-Jacobi equations. One of the strengths of the
convexification method is that it does not require the convexity of $F({\bf
x},s,{\bf p})$ with respect to ${\bf p}$. Still, it has a drawback. The theory
of the convexification method requires an additional information about $u_{z}$
on a part of $\partial\Omega$. In this paper, that part is

$\Gamma^{+}=\big{\\{}{\bf x}=(x_{1},x_{2},\dots,x_{d-1},z=R):|x_{i}|\leq
R,1\leq i\leq d-1\big{\\}}\subset\partial\Omega.$ (1.3)

For ${\bf x}\in\mathbb{R}^{d}$, write ${\bf x}=(x_{1},x_{2},\dots,x_{d-1},z)$.
More precisely, we impose the following condition.

######  Assumption 1.1.

Write $\nabla u=(u_{x_{1}},\dots,u_{x_{d-1}},u_{z})$. We assume that $u_{z}$
on $\Gamma^{+}$ is known.

In general, the additional knowledge of $u_{z}$ on $\Gamma^{+}$ in Assumption
1.1 makes the problem of computing solutions to Hamilton-Jacobi equations with
the Dirichlet data on $\partial\Omega$ and the Neumann data on $\Gamma^{+}$
over-determined. However, in many real-world circumstances, we are able to
compute $u_{z}$ on $\Gamma^{+}$ from the knowledge of $u$ on $\Gamma^{+}$
without further measurement. We provide here a classical example arising from
the traveling time problem. Denote by $c({\bf x})$, ${\bf
x}\in\overline{\Omega}$, the velocity of the light at the point ${\bf x}$.
Here, $c\in C(\overline{\Omega},(0,\infty))$ is a given function. Let $u({\bf
x})$ be the minimal time for the light to travel from $\partial\Omega$ to
${\bf x}\in\Omega$. This function is governed by the boundary value problem
for the eikonal equation

$\left\\{\begin{array}[]{rcll}c({\bf x})^{2}|\nabla u({\bf x})|^{2}&=&1&{\bf
x}\in\Omega,\\\ u({\bf x})&=&0&{\bf x}\in\partial\Omega.\end{array}\right.$
(1.4)

Since $u=0$ on $\partial\Omega$, in particular, $u=0$ on $\Gamma^{+}.$ Hence
$u_{x_{i}}=0$ on $\Gamma^{+}$ for all $1\leq i\leq d-1$. The function $u_{z}$
on $\Gamma^{+}$ is given by

$u_{z}({\bf x})=-\frac{1}{c({\bf x})}\quad\mbox{for all }{\bf
x}\in\Gamma^{+}.$

Above, the case $u_{z}({\bf x})=\frac{1}{c({\bf x})}$ is negligible since in
reality, $u\geq 0$ on $\overline{\Omega}$ and therefore, its partial
derivative with respect to $z$ on $\Gamma^{+}$ is non-positive.

######  Remark 1.1 (Reducing Assumption 1.1).

We have the following points.

  * 1.

It follows from the example above that in general, since $u$, and; therefore,
$u_{x_{1}},$ $\dots,$ $u_{x_{d-1}}$ on $\Gamma^{+}$ are known, to verify
Assumption 1.1, we simply solve the equation

$F({\bf x},u,u_{x_{1}},\dots,u_{x_{d-1}},u_{z})=0$ (1.5)

for $u_{z}$. So, a condition on $F$ such that Assumption 1.1 holds true is
that (1.5) is uniquely solvable for $u_{z}$. For example,

$\frac{\partial}{\partial z}F({\bf
x},u,u_{x_{1}},\dots,u_{x_{d-1}},u_{z})\not=0\quad\mbox{for all }{\bf
x}\in\Gamma^{+}.$

  * 2.

In many important Hamilton-Jacobi equations; for e.g., the eikonal equation
$F(x,u,\nabla u)=c({\bf x})^{2}|\nabla u|^{2}-1$ in (1.4) or $F({\bf
x},u,\nabla u)=|u_{x_{1}}|-|u_{z}|-g$ (the case $d=2$), for some function $g$,
equation (1.5) only provides $|u_{z}|$ rather than $u_{z}$ on $\Gamma^{+}$. In
this case, the sign of $u_{z}$ on $\Gamma^{+}$ is required.

  * 3.

We provide here an example in which the sign of $u_{z}$ is known. Let $d=3$
and $\mathbf{x}_{0}$ be a point in $\mathbb{R}^{3}\setminus\Omega$. In travel
time tomography, the Hamilton-Jacobi equation that describes the the travel
time of light traveling from the source $\mathbf{x}_{0}$ to a point
$\mathbf{x}\in\mathbb{R}^{3},$ has the form
$c^{2}(\mathbf{x})|\nabla_{\mathbf{x}}u(\mathbf{x},\mathbf{x}_{0})|^{2}=1$ for
all $\mathbf{x}\in\mathbb{R}^{3}\setminus\\{\mathbf{x}_{0}\\}$ with
$u(\mathbf{x}_{0},\mathbf{x}_{0})=0.$ Here $c(\mathbf{x})$ is the speed of the
light at $\mathbf{x}$. It was proved in [29, Lemma 4.1] that if
$c(\mathbf{x})$, $\mathbf{x}\in\Omega$ is an increasing function with respect
to $z$ and the source $x_{0}\in\left\\{z<-R\right\\},$ then the function
$u(x,x_{0})$ is strictly increasing in the $z-$direction for $x\in\Omega$
implying $u_{z}>0$ on $\Gamma^{+}.$ This result can be extended to all
dimensions by repeating the proof in [29, Lemma 4.1].

  * 4.

In the case when Assumption 1.1 cannot be verified or even when it might not
hold true, for e.g.,

$F({\bf x},u,\nabla u)=10u+\min\\{|\nabla u|,||\nabla u|-8|+6\\}-g$

or

$F({\bf x},u,\nabla u)=u+|\nabla u|-V\cdot\nabla u$

for some function $g$ and vector valued function $V$, the convexification
method still provides good numerical solutions, see Test 4 and Test 5 in
Section 5. However, the rigorous theorem that guarantees the efficiency of
this method is missing in this paper.

Let us give a brief description of the main results in the paper. We consider
the vanishing viscosity process (equations (2.1) and (4.1)) and aim at
computing $u^{\epsilon_{0}}$ for $\epsilon_{0}>0$ sufficiently small, which is
a good approximation of $u$, the viscosity solution to (1.1)–(1.2). The
convexification is developed to compute this $u^{\epsilon_{0}}$. Firstly, we
obtain a new Carleman inequality in Theorem 3.1: For $\beta>1$, $r>R+1$, and
$b>R+r$, we can find two numbers $\lambda_{0}=\lambda_{0}(\beta,r,R,b,d)>0$,
$C=C(r,R,b,d)>0$ such that for all $\lambda>\lambda_{0}$ and for all $u\in
C^{2}(\overline{\Omega})$ with $u=0$ on $\partial\Omega$ and $u_{z}=0$ on
$\Gamma^{+}$, we have

$\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta u|^{2}d{\bf x}\geq
C\lambda^{3}\beta^{2}(\beta-1)b^{-3\beta}(-R+r)^{2\beta}\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}d{\bf
x}\\\
+C\lambda(\beta-1)b^{-\beta}\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla
u|^{2}\,d{\bf x}.$

We then use this Carleman estimate to show in Theorem 4.1 that the functional

$J_{\lambda,\beta,\eta}(u)=\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\big{|}-\epsilon_{0}\Delta
u+F({\bf x},u,\nabla u)\big{|}^{2}\,d{\bf x}+\eta\|u\|^{2}_{H^{p}(\Omega)}$
(1.6)

is strictly convex for $u\in H\cap\overline{B(M)}$. Here, $p>\lceil
d/2\rceil+2$ is such that $H^{p}\hookrightarrow C^{2}(\overline{\Omega})$, and

$H=\\{u\in H^{p}(\Omega):u|_{\partial\Omega}=f\mbox{ and
}u_{z}|_{\Gamma^{+}}=g\\},\quad{B(M)}=\\{u\in
H^{p}(\Omega):\|u\|_{H^{p}(\Omega)}<M\\}.$

Then, we use a gradient descent method (Theorem 4.14) to compute the minimizer
$u_{\rm min}$ of this functional. Assuming that $u_{\rm min}\in H\cap B(M/3)$,
we can start the gradient descent method at $u^{(0)}\in H\cap B(M/3)$ and
iterate

$u^{(k)}=u^{(k-1)}-\kappa DJ_{\lambda,\beta,\eta}(u^{(k-1)})\quad\text{ for
}k\in\mathbb{N}.$

Here, $\kappa\in(0,\kappa_{0})$ where $\kappa_{0}\in(0,1)$ depends only on
$\lambda,$ $\beta$, $R,$ $r,$ $b$, $d$, $M$ and $\epsilon_{0}$. We are able to
obtain

$\big{\|}u^{(k)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}\leq\theta^{k/2}\big{\|}u^{(0)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}\quad\text{ for }k\in\mathbb{N},$

for some $\theta\in(0,1)$ depending only on $\kappa,$ $\lambda,$ $\beta$, $R,$
$r,$ $b$, $d$, $M$, $F$ and $\epsilon_{0}$. Then, in Theorem 4.18, we show
that, even if there is a noise of size $\delta>0$, we still have a nice bound

$\|u_{\min}^{\delta}-u^{\epsilon_{0}}\|^{2}_{H^{1}(\Omega)}\leq
C(\eta\|u^{\epsilon_{0}}\|^{2}_{H^{p}(\Omega)}+\delta^{2}).$

Here, $u^{\delta}_{\min}$ is the minimizer of $J_{\lambda,\beta,\eta}(u)$ with
noisy data $u|_{\partial\Omega}=f^{\delta}$, $u_{z}|_{\Gamma^{+}}=g^{\delta}$.
Here, by saying that $\delta$ is the noise level, there exists an “error”
function $\mathcal{E}$ satisfying $\|\mathcal{E}\|_{H^{p}(\Omega)}\leq\delta$,
$\mathcal{E}|_{\partial\Omega}=f^{\delta}-f$,
$\mathcal{E}_{z}|_{\Gamma^{+}}=g^{\delta}-g$. Combining Theorem 4.14 and
Theorem 4.18, we have for each $k\geq 1$,

$\|u^{(k)}-u^{\epsilon_{0}}\|_{H^{1}(\Omega)}\leq
C(\sqrt{\eta}\|u^{\epsilon_{0}}\|_{H^{p}(\Omega)}+\delta)+\theta^{k/2}\|u^{(0)}-u^{\delta}_{\min}\|_{H^{p}(\Omega)}.$

This inequality shows the stability of our method with respect to noise. If
$\theta^{k}$ and $\eta$ are $O(\delta^{2})$ as $\delta$ tends to $0$, then the
convergence rate is Lipschitz. Finally, in Section 5, we implement the
convexification method based on the finite difference method and obtain
interesting numerical results in two dimensions.

We now address a bit further some state of the art numerical methods in
solving (1.1)–(1.2) in the literature. If $F$ is generically convex in $\nabla
u$, there have been extremely powerful approaches to compute the solutions
such as monotone numerical Hamiltonian based finite difference methods (see
[46, 50, 51, 57, 49] and the references therein). When $F$ is nonconvex in
$\nabla u$, the Lax–Friedrichs schemes ([47, 1, 45]) and the Lax–Friedrichs
sweeping algorithm ([23, 39]), in which numerical viscosity terms appear
naturally, are very efficient and accurate. Moreover, all the mentioned
methods have very quick running times with not too many iterations. Similar to
the Lax–Friedrichs schemes, the addition of a viscosity term is natural in our
approach as we deal with general $F$, which is possibly nonconvex in $\nabla
u$.

The paper is organized as follows. In Section 2, we give some preliminaries
about viscosity solutions to Hamilton-Jacobi equations, which are rather well-
known in the literature. We state and prove a Carleman estimate in Theorem 3.1
in Section 3. Section 4 is devoted to the theoretical results of the
convexification (Theorems 4.1–4.18), which is our main focus in this current
paper. Then, in Section 5, we implement the convexification method based on
the finite difference method and obtain interesting numerical results in two
dimensions.

##  2 Some preliminaries about viscosity solutions to Hamilton-Jacobi
equations

######  Definition 2.1 (Viscosity solutions of (1.1)–(1.2)).

Let $u\in C(\overline{\Omega})$.

  * (a)

We say that $u$ is a viscosity subsolution to (1.1)–(1.2) if for any test
function $\varphi\in C^{1}(\overline{\Omega})$ such that $u-\varphi$ has a
strict maximum at ${\bf x}_{0}\in\overline{\Omega}$, then

$F({\bf x}_{0},u({\bf x}_{0}),\nabla\varphi({\bf x}_{0}))\leq 0\quad\text{ if
}{\bf x}_{0}\in\Omega,$

or

$\min\left\\{F({\bf x}_{0},u({\bf x}_{0}),\nabla\varphi({\bf x}_{0})),u({\bf
x}_{0})-f({\bf x}_{0})\right\\}\leq 0\quad\text{ if }{\bf
x}_{0}\in\partial\Omega.$

  * (b)

We say that $u$ is a viscosity supersolution to (1.1)–(1.2) if for any test
function $\varphi\in C^{1}(\overline{\Omega})$ such that $u-\varphi$ has a
strict minimum at ${\bf x}_{0}\in\overline{\Omega}$, then

$F({\bf x}_{0},u({\bf x}_{0}),\nabla\varphi({\bf x}_{0}))\geq 0\quad\text{ if
}{\bf x}_{0}\in\Omega,$

or

$\max\left\\{F({\bf x}_{0},u({\bf x}_{0}),\nabla\varphi({\bf x}_{0})),u({\bf
x}_{0})-f({\bf x}_{0})\right\\}\geq 0\quad\text{ if }{\bf
x}_{0}\in\partial\Omega.$

  * (c)

We say that $u$ is a viscosity solution to (1.1)–(1.2) if it is both its
viscosity subsolution and its viscosity supersolution.

It is worth noting that the Dirichlet boundary condition might not hold in the
classical sense for viscosity solutions to (1.1)–(1.2) (see [56, Appendix E]
for example). In the following, we impose some compatibility conditions to
make sure that $u=f$ on $\partial\Omega$ classically.

We write $F=F({\bf x},s,{\bf p})$. Let $\nabla_{{\bf x}}F$ and $\nabla_{{\bf
p}}F$ the gradient vector of $F$ with respect to the first and the third
variables, respectively, and $\partial_{s}F=F_{s}$ the partial derivative of
$F$ with respect to the second variable. Here is one set of conditions that
often occurs in the literature.

######  Assumption 2.1.

1\. There exists $\alpha>0$ such that

$\alpha\leq F_{s}({\bf x},s,{\bf p})\leq\frac{1}{\alpha}\quad\text{ for all
}({\bf x},s,{\bf p})\in\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d}.$

2\. $F$ is coercive in ${\bf p}$ in the sense that

$\lim_{|{\bf p}|\to\infty}\inf_{({\bf
x},s)\in\overline{\Omega}\times\mathbb{R}}\left(F({\bf x},s,{\bf
p})^{2}+\partial_{s}F({\bf x},s,{\bf p})|{\bf p}|^{2}+\nabla_{{\bf x}}F({\bf
x},s,{\bf p})\cdot{\bf p}\right)=+\infty.$

3\. There exists $\phi\in C^{2}(\overline{\Omega})$ such that $\phi=f$ on
$\partial\Omega$, and

$F({\bf x},\phi({\bf x}),\nabla\phi({\bf x}))<0\quad\text{ in }\Omega.$

######  Theorem 2.1.

Suppose that Assumption 2.1 holds. Then, problem (1.1)–(1.2) has a unique
viscosity solution, denoted by $u$. Moreover, for $\epsilon>0,$ the following
problem

$\left\\{\begin{array}[]{rcll}-\epsilon\Delta u^{\epsilon}+F({\bf
x},u^{\epsilon},\nabla u^{\epsilon})&=&0&{\bf x}\in\Omega,\\\
u^{\epsilon}&=&f&{\bf x}\in\partial\Omega.\end{array}\right.$ (2.1)

has a unique solution in $u^{\epsilon}\in C^{2}(\overline{\Omega})$ and

$\lim_{\epsilon\to 0^{+}}\|u^{\epsilon}-u\|_{L^{\infty}(\Omega)}=0.$ (2.2)

Theorem 2.2 holds true under other appropriate sets of assumptions too. We
here just give one prototypical set of conditions, Assumption 2.1, for
demonstration. See [15, 14, 40, 56] for its proof. For other set of
appropriate conditions on $F$, see references listed in Section 5 in each
numerical test.

######  Remark 2.1.

In fact, under Assumption 2.1, we are able to quantify the rate of convergence
of $u^{\epsilon}$ to $u$ in $L^{\infty}(\Omega)$. There exists $C>0$
independent of $\epsilon\in(0,1)$ such that

$\|u^{\epsilon}-u\|_{L^{\infty}(\Omega)}\leq C\epsilon^{1/2}.$

See [16, 9, 56] for more details.

##  3 A Carleman estimate

We prove a Carleman estimate that plays an important role in our proof for the
convergence of the convexification method. In the proof of the Carleman
estimate, we will need the notation

$\Gamma^{-}=\big{\\{}{\bf x}=(x_{1},x_{2},\dots,x_{d-1},z=-R):|x_{i}|\leq
R,1\leq i\leq d-1\big{\\}}\subset\partial\Omega.$ (3.1)

######  Theorem 3.1 (Carleman estimate).

For $\beta>1$, $r>R+1$, and $b>R+r$, we can find a number
$\lambda_{0}=\lambda_{0}(\beta,r,R,b,d)$ such that for all
$\lambda>\lambda_{0}$ and for all $u\in C^{2}(\overline{\Omega})$ with $u=0$
on $\partial\Omega$ and $u_{z}=0$ on $\Gamma^{+}$, we have

$\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta u|^{2}d{\bf x}\geq
C\lambda^{3}\beta^{2}(\beta-1)b^{-3\beta}(-R+r)^{2\beta}\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}d{\bf
x}\\\
+C\lambda(\beta-1)b^{-\beta}\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla
u|^{2}\,d{\bf x}.$ (3.2)

where $C=C(r,R,b,d)$ is a constant. Here, $\lambda_{0}$ and $C$ depend only on
the listed parameters.

###### Proof.

In the proof below, $C_{1},$ $C_{2}$ and $C_{3}$ are constants depending only
on $r$, $R$, $b$, and $d$. We split the proof into several steps.

Step 1. Define the function

$v=e^{\lambda(\frac{z+r}{b})^{\beta}}u.$ (3.3)

For ${\bf x}=(x_{1},\dots,x_{d-1},z)$, $i=1,\dots,d-1$, we have

$\displaystyle u$ $\displaystyle=e^{-\lambda(\frac{z+r}{b})^{\beta}}v,$
$\displaystyle u_{x_{i}}$
$\displaystyle=e^{-\lambda(\frac{z+r}{b})^{\beta}}v_{x_{i}},$ $\displaystyle
u_{x_{i}x_{i}}$
$\displaystyle=e^{-\lambda(\frac{z+r}{b})^{\beta}}v_{x_{i}x_{i}},$
$\displaystyle u_{z}$
$\displaystyle=e^{-\lambda(\frac{z+r}{b})^{\beta}}\Big{(}-\lambda\beta
b^{-\beta}(z+r)^{\beta-1}v+v_{z}\Big{)},$

and

$u_{zz}=e^{-\lambda(\frac{z+r}{b})^{\beta}}\Big{(}\Big{(}\lambda^{2}\beta^{2}b^{-2\beta}(z+r)^{2\beta-2}-\lambda\beta(\beta-1)b^{-\beta}(z+r)^{\beta-2}\Big{)}v-2\lambda\beta
b^{-\beta}(z+r)^{\beta-1}v_{z}+v_{zz}\Big{)}.$

Therefore, using the inequality $(a-b+c)^{2}\geq-2ba-2bc$, we have

$\displaystyle e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta u|^{2}$
$\displaystyle=\ $
$\displaystyle\Big{[}\big{(}\lambda^{2}\beta^{2}b^{-2\beta}(z+r)^{2\beta-2}-\lambda\beta(\beta-1)b^{-\beta}(z+r)^{\beta-2}\big{)}v-2\lambda\beta
b^{-\beta}(z+r)^{\beta-1}v_{z}+\Delta v\Big{]}^{2}$ $\displaystyle\geq\ $
$\displaystyle-4\lambda^{2}\beta^{2}b^{-2\beta}(z+r)^{2\beta-3}\big{(}\lambda\beta
b^{-\beta}(z+r)^{\beta}-(\beta-1)\big{)}v_{z}v-4\lambda\beta
b^{-\beta}(z+r)^{\beta-1}v_{z}\Delta v.$ (3.4)

Dividing both sides of (3.4) by $2\lambda\beta b^{-\beta}(z+r)^{\beta-1}$ and
integrating the resulting equation, we have

$\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{2\lambda\beta b^{-\beta}(z+r)^{\beta-1}}\,d{\bf x}\geq I_{1}+I_{2},$
(3.5)

where

$I_{1}=-2\lambda\beta\int_{\Omega}b^{-\beta}(z+r)^{\beta-2}\big{(}\lambda\beta
b^{-\beta}(z+r)^{\beta}-(\beta-1)\big{)}v_{z}v\,d{\bf x}$ (3.6)

and

$I_{2}=-2\int_{\Omega}v_{z}\Delta v\,d{\bf x}.$ (3.7)

Step 2. In this step, we estimate $I_{2}$. Since $u=0$ on $\partial\Omega$,
$v=0$ on $\partial\Omega$. Since $u=u_{z}=0$ on $\Gamma^{+}$, $v_{z}=0$ on
$\Gamma^{+}$. Hence, $v_{z}=0$ on $\partial\Omega\setminus\Gamma^{-}$. Using
integration by parts, we have

$\displaystyle I_{2}$ $\displaystyle=-2\int_{\Omega}{\rm div}(v_{z}\nabla
v)\,d{\bf x}+2\int_{\Omega}\nabla v_{z}\cdot\nabla v\,d{\bf
x}=-2\int_{\Gamma^{+}}|v_{z}|^{2}\,d\sigma+2\int_{\Gamma^{-}}|v_{z}|^{2}\,d\sigma+\int_{\Omega}(|\nabla
v|^{2})_{z}\,d{\bf x}$
$\displaystyle=2\int_{\Gamma^{-}}|v_{z}|^{2}\,d\sigma+\int_{\Gamma^{+}}|\nabla
v|^{2}\,d\sigma-\int_{\Gamma^{-}}|\nabla v|^{2}\,d\sigma.$

Since $v_{x_{i}}=0$ on $\Gamma^{+}\cup\Gamma^{-}$, $i=1,\dots,d-1$, $|\nabla
v|=|v_{z}|$ on this set. Using the fact that $v_{z}=0$ on
$\partial\Omega\setminus\Gamma^{-}$ again, we have

$I_{2}=\int_{\Gamma^{+}}|v_{z}|^{2}\,d\sigma+\int_{\Gamma^{-}}|v_{z}|^{2}\,d\sigma\geq
0.$ (3.8)

Step 3. We estimate $I_{1}.$ It follows from (3.6) that

$\displaystyle I_{1}$
$\displaystyle=-\lambda\beta\int_{\Omega}b^{-\beta}(z+r)^{\beta-2}\big{(}\lambda\beta
b^{-\beta}(z+r)^{\beta}-(\beta-1)\big{)}(|v|^{2})_{z}\,d{\bf x}$
$\displaystyle=-\lambda\beta\int_{\Omega}\big{[}b^{-\beta}(z+r)^{\beta-2}\big{(}\lambda\beta
b^{-\beta}(z+r)^{\beta}-(\beta-1)\big{)}|v|^{2}\big{]}_{z}\,d{\bf x}$
$\displaystyle\hskip
170.71652pt+\lambda\beta\int_{\Omega}\big{[}b^{-\beta}(z+r)^{\beta-2}\big{(}\lambda\beta
b^{-\beta}(z+r)^{\beta}-(\beta-1)\big{)}\big{]}_{z}|v|^{2}\,d{\bf x}$
$\displaystyle=\lambda\beta\int_{\Omega}\big{[}2\lambda\beta(\beta-1)b^{-2\beta}(z+r)^{2\beta-3}-(\beta-1)(\beta-2)b^{-\beta}(z+r)^{\beta-3}\big{]}|v|^{2}\,d{\bf
x}.$

Now, for fixed $\beta>1$, we can find $\lambda_{0}\geq 1$, only depending on
$\beta$, $r$ and $R$, sufficiently large such that for all
$\lambda>\lambda_{0}$,

$I_{1}\geq
C_{1}\lambda^{2}\beta^{2}(\beta-1)b^{-2\beta}\int_{\Omega}(z+r)^{2\beta}|v|^{2}d{\bf
x}.$ (3.9)

Step 3. Combining (3.5), (3.8) and (3.9), we have

$\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{2\lambda\beta b^{-\beta}(z+r)^{\beta-1}}d{\bf x}\geq
C_{1}\lambda^{2}\beta^{2}(\beta-1)b^{-2\beta}\int_{\Omega}(z+r)^{2\beta}|v|^{2}d{\bf
x}.$ (3.10)

Recall that $v=e^{\lambda(\frac{z+r}{b})^{\beta}}u$. It follows from (3.10)
that

$\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{2\lambda\beta b^{-\beta}(z+r)^{\beta-1}}d{\bf x}\geq
C_{1}\lambda^{2}\beta^{2}(\beta-1)b^{-2\beta}\int_{\Omega}(z+r)^{2\beta}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}d{\bf
x}.$

or equivalently,

$\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta u|^{2}}{\beta
b^{-\beta}(z+r)^{\beta-1}}d{\bf x}\geq
C_{1}\lambda^{3}\beta^{2}(\beta-1)b^{-2\beta}\int_{\Omega}(z+r)^{2\beta}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}d{\bf
x}.$ (3.11)

Step 4. We estimate the term
$\displaystyle\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}u\Delta u\,d{\bf
x}$. Using integration by parts, since $u=0$ on $\partial\Omega,$ we have

$\displaystyle-\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}u\Delta
u\,d{\bf x}$ $\displaystyle=-\int_{\Omega}{\rm
div}\big{[}e^{2\lambda(\frac{z+r}{b})^{\beta}}u\nabla u\big{]}\,d{\bf
x}+\int_{\Omega}\nabla(e^{2\lambda(\frac{z+r}{b})^{\beta}}u)\cdot\nabla
u\,d{\bf x}$
$\displaystyle=\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\big{[}2\lambda\beta
b^{-\beta}(z+r)^{\beta-1}u+u_{z}\big{]}u_{z}\,d{\bf
x}+\sum_{i=1}^{d-1}\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u_{x_{i}}|^{2}\,d{\bf
x}$
$\displaystyle=\lambda\beta\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}b^{-\beta}(z+r)^{\beta-1}(|u|^{2})_{z}\,d{\bf
x}+\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla u|^{2}\,d{\bf x}$
$\displaystyle=\lambda\beta\int_{\Omega}\Big{[}e^{2\lambda(\frac{z+r}{b})^{\beta}}b^{-\beta}(z+r)^{\beta-1}|u|^{2}\Big{]}_{z}\,d{\bf
x}-\lambda\beta\int_{\Omega}\Big{[}e^{2\lambda(\frac{z+r}{b})^{\beta}}b^{-\beta}(z+r)^{\beta-1}\Big{]}_{z}|u|^{2}\,d{\bf
x}$ $\displaystyle\hskip
170.71652pt+\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}b^{-\beta}}|\nabla
u|^{2}\,d{\bf x}.$

Due to the fact that $u=0$ on $\partial\Omega$, the first integral in the last
row above vanishes. We have

$-\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}u\Delta u\,d{\bf
x}=-\lambda\beta\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\Big{(}\lambda\beta
b^{-2\beta}(z+r)^{2\beta-2}+(\beta-1)b^{-\beta}(z+r)^{-2+\beta}\Big{)}|u|^{2}\,d{\bf
x}\\\ +\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}b^{-\beta}}|\nabla
u|^{2}\,d{\bf x}.$

Hence, for $\lambda>\lambda_{0}$,

$-\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}u\Delta u\,d{\bf x}\geq-
C_{2}\lambda^{2}\beta^{2}b^{-2\beta}\int_{\Omega}(z+r)^{2\beta}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}\,d{\bf
x}+\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla u|^{2}\,d{\bf x}.$
(3.12)

Step 5. We complete the proof in this step. Using the inequality
$x^{2}+y^{2}\geq-2xy$, we have

$\lambda^{5/2}\beta
b^{-\beta}\int_{\Omega}(z+r)^{\beta-1}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}d{\bf
x}+\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{4\lambda^{1/2}\beta b^{-\beta}(z+r)^{\beta-1}}d{\bf x}\\\
\geq-\lambda\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}u\Delta u\,d{\bf
x}.$ (3.13)

Combining (3.12) and (3.13), we have

$\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{4\lambda^{1/2}\beta b^{-\beta}(z+r)^{\beta-1}}d{\bf
x}\geq-\lambda^{5/2}\beta
b^{-\beta}\int_{\Omega}(z+r)^{\beta-1}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}d{\bf
x}\\\
-C_{2}\lambda^{3}\beta^{2}b^{-2\beta}\int_{\Omega}(z+r)^{2\beta}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}\,d{\bf
x}+\lambda\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla
u|^{2}\,d{\bf x},$

which implies

$\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{4\lambda^{1/2}\beta b^{-\beta}(z+r)^{\beta-1}}d{\bf x}\geq-
C_{3}\lambda^{3}\beta^{2}b^{-2\beta}\int_{\Omega}(z+r)^{2\beta}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}\,d{\bf
x}\\\ +\lambda\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla
u|^{2}\,d{\bf x}.$ (3.14)

Multiply $\frac{C_{1}(\beta-1)}{2C_{3}}$ to both sides of (3.14) where $C_{1}$
is the constant in (3.11), we yield

$\frac{C_{1}(\beta-1)}{2C_{3}}\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{4\lambda^{1/2}\beta b^{-\beta}(z+r)^{\beta-1}}d{\bf
x}\geq-\frac{C_{1}\lambda^{3}\beta^{2}(\beta-1)b^{-2\beta}}{2}\int_{\Omega}(z+r)^{2\beta}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}\,d{\bf
x}\\\
+\frac{C_{1}\lambda(\beta-1)}{2C_{3}}\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla
u|^{2}\,d{\bf x}.$ (3.15)

Adding (3.11) and (3.15), we have

$\Big{(}1+\frac{C_{1}(\beta-1)}{2C_{3}\lambda^{1/2}}\Big{)}\int_{\Omega}\frac{e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta
u|^{2}}{\beta b^{-\beta}(z+r)^{\beta-1}}d{\bf
x}\geq\frac{C_{1}\lambda^{3}\beta^{2}(\beta-1)b^{-2\beta}}{2}\int_{\Omega}(z+r)^{2\beta}e^{2\lambda(\frac{z+r}{b})^{\beta}}|u|^{2}\,d{\bf
x}\\\
+\frac{C_{1}\lambda(\beta-1)}{2C_{3}}\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\nabla
u|^{2}\,d{\bf x}.$ (3.16)

The desired Carleman estimate (3.2) follows from (3.16).  $\square$

######  Remark 3.1.

In our previous publications, see e.g. [43, 36], the number $\beta$ must be
large. The reason for (3.2) holds true when $\beta>1$ is our trick of dividing
both sides of (3.4) by $2\lambda\beta b^{-\beta}(z+r)^{\beta-1}$ so that the
corresponding $I_{2}$, defined in (3.7), is nonnegative. This trick was used
in [36, Theorem 3.1]. However, since the Carleman weight in [36, Theorem 3.1]
is different from the one in (3.2), the parameter $\beta$ in [36, Theorem 3.1]
must be large. In the case when the principal differential operator in (3.2)
is replaced by the general elliptic operator, this trick is not applicable.

##  4 The convexification method to compute the viscosity solution to
Hamilton-Jacobi equations

In this section, we propose to use the convexification method to solve (1.1)
together with the boundary condition (1.2) supposing that Assumption 1.1 holds
true. That means we know the boundary value $u=f$ on $\partial\Omega$ and the
function $u_{z}|_{\Gamma^{+}}$ can be computed, say $u_{z}=g$ on $\Gamma^{+}.$
Due to Theorem 2.2, it is natural to try to approximate the solution $u$ by a
function $u^{\epsilon_{0}}$ that satisfies

$\left\\{\begin{array}[]{rcll}-\epsilon_{0}\Delta u^{\epsilon_{0}}+F({\bf
x},u^{\epsilon_{0}},\nabla u^{\epsilon_{0}})&=&0&{\bf x}\in\Omega,\\\
u^{\epsilon_{0}}&=&f({\bf x})&{\bf x}\in\partial\Omega,\\\
u^{\epsilon_{0}}_{z}&=&g({\bf x})&{\bf x}\in\Gamma^{+}.\end{array}\right.$
(4.1)

######  Remark 4.1.

1\. In general, (4.1) is over-determined. It might have no solution;
especially when the boundary data contains noise. However, the convexification
method can deliver a function that “most fits” (4.1) and show that this
function is an approximation of the true solution to (1.1)–(1.2) when the
given boundary data are noiseless. Again, due to Assumption 1.1, the Neumann
condition imposed in (4.1) makes sense.

2\. On the other hand, (4.1) is not over-determined in the sense that

$\left\\{\begin{array}[]{rcll}-\epsilon_{0}\Delta u^{\epsilon_{0}}+F({\bf
x},u^{\epsilon_{0}},\nabla u^{\epsilon_{0}})&=&0&{\bf x}\in\Omega,\\\
u^{\epsilon_{0}}&=&f({\bf x})&{\bf x}\in\partial\Omega,\\\ \end{array}\right.$
(4.2)

might not be uniquely solvable if we do not impose Assumption 2.1 or other
sets of appropriate conditions. For example, when $F({\bf
x},u^{\epsilon_{0}},\nabla u^{\epsilon_{0}})=-u^{\epsilon_{0}}$, $f=0$, and
$\frac{1}{\epsilon_{0}}$ is an eigenvalue of $-\Delta$, (4.2) has multiple
solutions. Therefore, imposing the additional Neumann boundary condition on
$\Gamma^{+}$, in the general case, is necessary.

Let $p>\lceil d/2\rceil+2$ such that $H^{p}(\Omega)\hookrightarrow
C^{2}(\overline{\Omega}).$ Define

$H=\\{u\in H^{p}(\Omega):u|_{\partial\Omega}=f\mbox{ and
}u_{z}|_{\Gamma^{+}}=g\\}.$

Clearly, $H$ is a closed subset of $H^{p}(\Omega).$ We will also need the
following set of test functions

$H_{0}=\\{u\in H^{p}(\Omega):u|_{\partial\Omega}=0\mbox{ and
}u_{z}|_{\Gamma^{+}}=0\\}.$

Let $M>0$ be chosen later, and set

${B(M)}=\\{u\in H^{p}(\Omega):\|u\|_{H^{p}(\Omega)}<M\\}.$

We assume that $H\cap\overline{B(M)}\not=\emptyset$. For each $\lambda>1$,
$\beta>1$ and $\eta>0$, introduce the functional

$J_{\lambda,\beta,\eta}(u)=\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\big{|}-\epsilon_{0}\Delta
u+F({\bf x},u,\nabla u)\big{|}^{2}\,d{\bf
x}+\eta\|u\|^{2}_{H^{p}(\Omega)}\quad u\in H\cap\overline{B(M)}$ (4.3)

where $\epsilon_{0}$ is a fixed small positive number.

The convexification theorem is to prove that for each $\beta>1$, there is a
number $\lambda_{0}>1$ such that for all $\lambda\geq\lambda_{0}$ and for all
$\eta>0$, the function $J_{\lambda,\beta,\eta}$ is strictly convex in
$H\cap\overline{B(M)}.$ The word “convexification” is suggested by the fact
that the Carleman weight function $e^{2\lambda(\frac{z+r}{b})^{\beta}}$
convexifies this functional. The convexification theorem is stated below.

######  Theorem 4.1 (The convexification theorem).

1\. For all $\lambda,\beta>1$, and $\eta>0$, for all $u\in H^{p}(\Omega)$,
$h\in H_{0}$, we have

$\lim_{H_{0}\in h\to
0}\frac{|J_{\lambda,\beta,\eta}(u+h)-J_{\lambda,\beta,\eta}(u)-DJ_{\lambda,\beta,\eta}(u)h|}{\|h\|_{H^{p}(\Omega)}}=0$
(4.4)

where

$DJ_{\lambda,\beta,\eta}(u)h=2\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\Big{[}-\epsilon_{0}\Delta
u+F({\bf x},u,\nabla u)\Big{]}\Big{[}-\epsilon_{0}\Delta h+\partial_{s}F({\bf
x},u,\nabla u)h\\\ +\nabla_{{\bf p}}F({\bf x},u,\nabla u)\cdot\nabla
h\Big{]}\,d{\bf x}+2\eta\langle u,h\rangle_{H^{p}(\Omega)}.$ (4.5)

2\. Let $M$ be an arbitrarily large number. For each $\beta>1$,
$\lambda>\lambda_{0}=\lambda_{0}(\epsilon_{0},M,b,d,r,F,\beta)>1$, $\eta>0$,
$u\in H\cap\overline{B(M)}$ and $v\in H\cap\overline{B(M)}$, we have

$J_{\lambda,\beta,\eta}(v)-J_{\lambda,\beta,\eta}(u)-DJ_{\lambda,\beta,\eta}(u)h\geq
C\|v-u\|^{2}_{H^{1}(\Omega)}+\eta\|v-u\|^{2}_{H^{p}(\Omega)},$ (4.6)

and

$\langle
DJ_{\lambda,\beta,\eta}(v)-DJ_{\lambda,\beta,\eta}(u),(v-u)\rangle_{H^{p}(\Omega)}\geq
C\|v-u\|^{2}_{H^{1}(\Omega)}+\eta\|v-u\|^{2}_{H^{p}(\Omega)}.$ (4.7)

Here, the constant $C$ depends only on $\lambda,$ $\beta$, $R,$ $r,$ $b$, $d$,
$M$, $F$ and $\epsilon_{0}$. As a result, the functional
$J_{\lambda,\beta,\eta}$ has a unique minimizer in $\overline{B(M)}$.

The key point for us to successfully establish the inequalities (4.6) and
(4.7) is the presence of the Carleman weight function
$e^{2\lambda(\frac{z+r}{b})^{\beta}}$ and the use of the Carleman estimate
(3.2). A direct consequence of the inequalities (4.6) and (4.7) is the strict
convexity of $J_{\lambda,\beta,\eta}$ in $H\cap\overline{B(M)}.$ It is worth
mentioning that $C$ in (4.6)–(4.7) depends on the viscosity coefficient
$\epsilon_{0}$, and we need
$\lambda>\lambda_{0}=\lambda_{0}(\epsilon_{0},M,b,d,r,F,\beta)>1$ so that
$J_{\lambda,\beta,\eta}$ is convex in $H\cap\overline{B(M)}.$ On the other
hand, $D^{2}J_{\lambda,\beta,\eta}\geq\eta Id$ in $H\cap\overline{B(M)}$,
which means that the uniform convexity of $J_{\lambda,\beta,\eta}$ only
depends on $\eta$, not $\epsilon_{0}$. The proof of Theorem 4.1 is similar to
that of the convexification theorem in [3], which is originally designed to
solve highly nonlinear and severely ill-posed inverse problems. Since this is
the first time the convexification method is employed in the area of numerical
methods for Hamilton-Jacobi equations, we present the proof here for the
reader’s convenience.

######  Remark 4.2.

For all $u\in H^{p}(\Omega)$, since $DJ_{\lambda,\beta,\eta}(u)$ is a bounded
linear map from $H_{0}$ into $\mathbb{R}$, by the Riesz theorem there exists
uniquely the function $J^{\prime}_{\lambda,\beta,\nu}(u)\in H_{0}$ such that

$\langle
J^{\prime}_{\lambda,\beta,\nu}(u),h\rangle_{H^{p}(\Omega)}=DJ_{\lambda,\beta,\eta}(u)h$

for all $h\in H_{0}.$

###### Proof of Theorem 4.1.

Since $F$ is in the class
$C^{2}(\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d})$, for all
$u,h\in H^{p}(\Omega)$ we can write

$F({\bf x},u+h,\nabla(u+h))=F({\bf x},u,\nabla u)+\partial_{s}F({\bf
x},u,\nabla u)h+\nabla_{{\bf p}}F({\bf x},u,\nabla u)\cdot\nabla
h+O(|h|^{2})+O(|\nabla h|^{2}).$ (4.8)

Here, $O(s)$ is the quantity satisfying $|O(s)|\leq C|s|$ where $C$ is a
constant that might depend on an upper bound of
$\|u\|_{C^{1}(\overline{\Omega})}$ and the function $F$. Using (4.3) and
(4.8), we obtain

$\displaystyle J_{\lambda,\beta,\eta}($ $\displaystyle
u+h)-J_{\lambda,\beta,\eta}(u)-2\eta\langle
u,h\rangle_{H^{p}(\Omega)}-\eta\|h\|^{2}_{H^{p}(\Omega)}$
$\displaystyle=\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\Big{(}|-\epsilon_{0}\Delta(u+h)+F({\bf
x},u+h,\nabla(u+h))|^{2}-|-\epsilon_{0}\Delta u+F({\bf x},u,\nabla
u)|^{2}\Big{)}\,d{\bf x}$
$\displaystyle=\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\big{(}-\epsilon_{0}\Delta
h+F({\bf x},u+h,\nabla(u+h))-F({\bf x},u,\nabla
u)\big{)}\big{(}-\epsilon_{0}\Delta h-2\epsilon_{0}\Delta u$
$\displaystyle\hskip 227.62204pt+F({\bf x},u+h,\nabla(u+h))+F({\bf x},u,\nabla
u)\big{)}\,d{\bf x}$
$\displaystyle=\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\Big{[}-\epsilon_{0}\Delta
h+\partial_{s}F({\bf x},u,\nabla u)h+\nabla_{{\bf p}}F({\bf x},u,\nabla
u)\cdot\nabla h+O(|h|^{2})+O(|\nabla h|^{2})\Big{]}$ $\displaystyle\hskip
28.45274pt\Big{[}2\big{(}-\epsilon_{0}\Delta u+F({\bf x},u,\nabla
u)\big{)}-\epsilon_{0}\Delta h+\partial_{s}F({\bf x},u,\nabla u)h+\nabla_{{\bf
p}}F({\bf x},u,\nabla u)\cdot\nabla h$ $\displaystyle\hskip
284.52756pt+O(|h|^{2})+O(|\nabla h|^{2})\Big{]}\,d{\bf x}.$

Hence,

$J_{\lambda,\beta,\eta}(u+h)-J_{\lambda,\beta,\eta}(u)-2\eta\langle
u,h\rangle_{H^{p}(\Omega)}-\eta\|h\|^{2}_{H^{p}(\Omega)}-2\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\Big{[}-\epsilon_{0}\Delta
u+F({\bf x},u,\nabla u)\Big{]}\\\ \Big{[}-\epsilon_{0}\Delta
h+\partial_{s}F({\bf x},u,\nabla u)h+\nabla_{{\bf p}}F({\bf x},u,\nabla
u)\cdot\nabla h\Big{]}\,d{\bf x}\\\
=2\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\Big{[}\epsilon_{0}^{2}|\Delta
h|^{2}+O(|h|^{2})+O(|\nabla h|^{2})\Big{]}\,d{\bf x}.$ (4.9)

Defining $DJ_{\lambda,\beta,\eta}$ as in (4.5) and using (4.9), we have

$\lim_{H_{0}\ni h\to
0}\frac{|J_{\lambda,\beta,\eta}(u+h)-J_{\lambda,\beta,\eta}(u)-DJ_{\lambda,\beta,\eta}(u)h|}{\|h\|_{H^{p}(\Omega)}}\\\
\leq\lim_{H_{0}\ni h\to
0}\frac{|J_{\lambda,\beta,\eta}(u+h)-J_{\lambda,\beta,\eta}(u)-DJ_{\lambda,\beta,\eta}(u)h|}{\|h\|_{H^{2}(\Omega)}}=0.$

We have proved part 1 of this theorem. We next prove part 2.

For any $u$ and $v$ in $H\cap\overline{B(M)}$, set $h=v-u\in H_{0}$. It
follows from (4.9) that

$\displaystyle J_{\lambda,\beta,\eta}(v)-J_{\lambda,\beta,\eta}(u)-$
$\displaystyle
DJ_{\lambda,\beta,\eta}(u)(v-u)=J_{\lambda,\beta,\eta}(u+h)-J_{\lambda,\beta,\eta}(u)-DJ_{\lambda,\beta,\eta}(u)h$
$\displaystyle\geq\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}\Big{[}\epsilon_{0}^{2}|\Delta
h|^{2}-O(|h|^{2})-O(|\nabla h|^{2})\Big{]}\,d{\bf
x}+\eta\|h\|_{H^{p}(\Omega)}^{2}.$ (4.10)

Note that $h|_{\partial\Omega}=0$ and $h_{z}|_{\Gamma^{+}}=0$. Applying
Theorem 3.1, for each $\beta>1$, we can find $\lambda_{0}$ depending on $R$,
$M$, $F$, $\beta,$ $b$ and $d$ such that for all $\lambda>\lambda_{0}$,

$\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}|\Delta h|^{2}\,d{\bf x}\geq
C\lambda\int_{\Omega}e^{2\lambda(\frac{z+r}{b})^{\beta}}(|\nabla
h|^{2}+\lambda^{2}|h|^{2})\,d{\bf x}.$ (4.11)

Here, the constant $C$ is allowed to depend on $\beta.$ We now choose
$\lambda$ such that $\epsilon_{0}^{2}\lambda$ is sufficiently large. Combining
(4.10) and (4.11), we have

$J_{\lambda,\beta,\eta}(v)-J_{\lambda,\beta,\eta}(u)-DJ_{\lambda,\beta,\eta}(u)(v-u)\geq
C\|h\|^{2}_{H^{1}(\Omega)}+\eta\|h\|^{2}_{H^{p}(\Omega)}.$

We have proved (4.6). Here, $C$ depends on $\lambda,$ $\beta,$ $\epsilon_{0}$,
$R$, $M$, $b$ and $d$. Interchanging the roles of $u$ and $v$ in (4.6) and
adding the resulting estimate to (4.6), we obtain (4.7).

We next prove that if (4.7) holds true, $J_{\lambda,\beta,\eta}$ has a unique
minimizer. The existence of the minimizer is obvious. It follows from the fact
that $H\cap\overline{B(M)}$ is convex in $H^{p}(\Omega)$ and the compact
embedding of $H^{p}(\Omega)$ to $H^{2}(\Omega).$ An alternative way to obtain
the existence of the minimizer is to argue similarly to the proofs of Theorem
2.1 in [3] or Theorem 4.1 in [32]. We now prove the uniqueness. Let $u_{1}$
and $u_{2}$ be two local minimizers of $J_{\lambda,\beta,\eta}$ on
$H\cap\overline{B(M)}$. It is clear that, see [3, Lemma 2],

$DJ_{\lambda,\beta,\eta}(u_{1})(u_{2}-u_{1})\geq 0\quad\mbox{and}\quad
DJ_{\lambda,\beta,\eta}(u_{2})(u_{1}-u_{2})\geq 0.$

Thus,

$(DJ_{\lambda,\beta,\eta}(u_{1})-DJ_{\lambda,\beta,\eta}(u_{2}))(u_{1}-u_{2})\leq
0$ (4.12)

Combining (4.7) for $u_{1}$ and $u_{2}$ and (4.12), we have

$C\|u_{1}-u_{2}\|^{2}_{H^{1}(\Omega)}+\eta\|u_{1}-u_{2}\|^{2}_{H^{p}(\Omega)}\leq
0.$

The proof is complete.  $\square$

######  Theorem 4.2 (The convergence of the gradient descent method).

Let $\lambda,\beta$ and $\eta$ as in part 2 of Theorem 4.1 and let
$J^{\prime}_{\lambda,\beta,\eta}$ be as in Remark 4.2. Let $u^{(0)}$ be any
function in $H\cap B(M/3)$. For $k\in\mathbb{N}$, define

$u^{(k)}=u^{(k-1)}-\kappa J_{\lambda,\beta,\eta}^{\prime}(u^{(k-1)})$ (4.13)

for all $\kappa\in(0,\kappa_{0})$ where $\kappa_{0}\in(0,1)$ is a number that
depends only on $\lambda,$ $\beta$, $R,$ $r,$ $b$, $d$, $M$ and
$\epsilon_{0}$. Then, if the minimizer $u_{\rm min}$ of
$J_{\lambda,\beta,\eta}$ is in $H\cap B(M/3)$ then there is a number
$\theta\in(0,1)$ depending only on $\kappa,$ $\lambda,$ $\beta$, $R,$ $r,$
$b$, $d$, $M$, $F$ and $\epsilon_{0}$ such that for all $k\geq 1$

$\big{\|}u^{(k)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}\leq\theta^{k/2}\big{\|}u^{(0)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}.$ (4.14)

Theorem 4.14 and the estimate (4.14) guarantee that the minimizer of
$J_{\lambda,\beta,\eta}$ can be found by the popular gradient descent method.
The success is due to the hypothesis that the desired minimizer is in the
interior of $B(M/3)$. We do not experience any difficulty due to not checking
this condition in the numerical study. However, to be more rigorous, in
general, if this condition cannot be verified, one can use the projected
gradient method as in [3, 25] to find the minimizer. In this paper, we choose
the gradient descent method because the implementation of the projection in
the projected gradient method is more complicated while there are many ready-
to-use packages; for e.g., the optimization toolbox of Matlab, for the
gradient descent method. In other words, Theorem 4.14 significantly reduces
our efforts in implementation. Although the proof of Theorem 4.14 is similar
to the proofs of [31, Theorem 6] in 1D case and [37, Theorem 2.2] in higher
dimensions, we briefly present the proof of Theorem 4.14 here for the
convenience of the reader.

###### Proof of Theorem 4.14.

Since $F$ is a function in the class
$C^{2}(\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d})$, it is obvious
that $\nabla_{{\bf x}}F$, $\partial_{s}F$, $\nabla_{{\bf p}}F$ are all
Lipschitz continuous in any bounded subdomain of
$\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d}$. As a result,
$DJ_{\lambda,\beta,\eta}$, see (4.5) for its definition, is Lipschitz
continuous on $H\cap\overline{B(M)}$. Hence, there is a positive number $L$
such that

$\|J_{\lambda,\beta,\eta}^{\prime}(u_{1})-J_{\lambda,\beta,\eta}^{\prime}(u_{2})\|_{H^{p}(\Omega)}\leq
L\|u_{1}-u_{2}\|_{H^{p}(\Omega)}$ (4.15)

for all $u_{1},u_{2}\in H\cap B(M).$ We claim that for any $k\geq 0$,

$\|u^{(k)}-u_{\rm min}\|_{H^{p}(\Omega)}<\frac{2M}{3}.$ (4.16)

This is true when $k=0$. Assume (4.16) is true for some $k$, we will prove
that (4.16) holds true for $k+1.$ Due to (4.16), $u^{(k)}\in H\cap B(M)$.
Since $u_{\rm min}$ is in $H\cap B(M/3)$,
$J_{\lambda,\beta,\eta}^{\prime}(u_{\rm min})=0$. Using (4.7) and (4.15), we
have

$\displaystyle\big{\|}u^{(k+1)}$ $\displaystyle-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}^{2}=\big{\|}u^{(k)}-u_{\rm
min}-\kappa[J_{\lambda,\beta,\eta}^{\prime}(u^{(k)})-J_{\lambda,\beta,\eta}^{\prime}(u_{\rm
min})]\big{\|}_{H^{p}(\Omega)}^{2}$ $\displaystyle=\big{\|}u^{(k)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}^{2}-2\kappa\big{\langle}J_{\lambda,\beta,\eta}^{\prime}(u^{(k)})-J_{\lambda,\beta,\eta}^{\prime}(u_{\rm
min}),u^{(k)}-u_{\rm min}\big{\rangle}_{H^{p}(\Omega)}$ $\displaystyle\hskip
199.16928pt+\kappa^{2}\big{\|}J_{\lambda,\beta,\eta}^{\prime}(u^{(k)})-J_{\lambda,\beta,\eta}^{\prime}(u_{\rm
min})\big{\|}_{H^{p}(\Omega)}^{2}$
$\displaystyle\leq(1-2\kappa\eta+\kappa^{2}L^{2})\big{\|}u^{(k)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}^{2}.$

Choosing $\kappa\in(0,\kappa_{0})$ where $\kappa_{0}=\eta L^{-2}$, we have
$\theta=1-2\kappa\eta+\kappa^{2}L^{2}\in(0,1)$ and

$\big{\|}u^{(k+1)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}\leq\theta^{1/2}\big{\|}u^{(k)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}<\big{\|}u^{(k)}-u_{\rm
min}\big{\|}_{H^{p}(\Omega)}<\frac{2M}{3}.$ (4.17)

We have proved (4.16). The estimate (4.14) follows from (4.17).  $\square$

Consider the case when the boundary data $f=u|_{\partial\Omega}$ contains
noise with the noise level $\delta>0$. Since the knowledge of
$g=u_{z}|_{\Gamma^{+}}$ can be computed from the knowledge of $u$ on this set,
see Assumption 1.1 and Remark 1.1, the Neumann data $u_{z}|_{\Gamma^{+}}$ also
contains noise. We assume that the noise level is still $\delta$. Denote by
$f^{\delta}$ and $g^{\delta}$ the noisy boundary data and denote the
corresponding noiseless data $f^{*}$ and $g^{*}$. Here, by saying that
$\delta$ is the noise level, there exists an “error” function $\mathcal{E}$
satisfying

$\left\\{\begin{array}[]{ll}\|\mathcal{E}\|_{H^{p}(\Omega)}\leq\delta,\\\
\mathcal{E}|_{\partial\Omega}=f^{\delta}-f^{*},\\\
\mathcal{E}_{z}|_{\Gamma^{+}}=g^{\delta}-g^{*}.\end{array}\right.$

The following theorem guarantees that the minimizer of
$J_{\lambda,\beta,\eta}$ subject to the boundary constraints defined by the
noisy data is an approximation of the true solution to (4.1).

######  Theorem 4.3.

Assume that the set

$H^{\delta}=\\{u\in H^{p}(\Omega):u|_{\partial\Omega}=f^{\delta}\mbox{ and
}u_{z}|_{\Gamma^{+}}=g^{\delta}\\}$

is nonempty and let $u^{\delta}_{\min}$ be the minimizer of
$J_{\lambda,\beta,\eta}$ in $H^{\delta}.$ Assume further that Problem (4.1),
in which the Dirichlet and Neumann data $f$ and $g$ are replaced by $f^{*}$
and $g^{*}$, respectively, has the unique solution $u^{\epsilon_{0}}\in
B(M-\delta)$. Then, for all $\lambda$ and $\beta$ such that (4.6) holds true,
we have

$\|u_{\min}^{\delta}-u^{\epsilon_{0}}\|^{2}_{H^{1}(\Omega)}\leq
C(\eta\|u^{\epsilon_{0}}\|^{2}_{H^{p}(\Omega)}+\delta^{2}).$ (4.18)

###### Proof.

Define $u=u^{\epsilon_{0}}+\mathcal{E}.$ We have $u\in H^{\delta}\cap B(M).$
Since $u^{\delta}_{\min}$ is the minimizer of $J_{\lambda,\beta,\eta}$ in
$H^{\delta}$, by [3, Lemma 2],

$\langle
J_{\lambda,\beta,\eta}^{\prime}(u^{\delta}_{\min}),u^{\delta}_{\min}-u^{\epsilon_{0}}-\mathcal{E}\rangle_{H^{p}(\Omega)}\leq
0.$

Applying (4.6) for $u$ and $u_{\min}$ gives

$\displaystyle J_{\lambda,\beta,\eta}(u^{\epsilon_{0}}+\mathcal{E})$
$\displaystyle\geq
J_{\lambda,\beta,\eta}(u^{\epsilon_{0}}+\mathcal{E})-J_{\lambda,\beta,\eta}(u^{\delta}_{\rm
min})-\langle J_{\lambda,\beta,\eta}^{\prime}(u^{\delta}_{\rm
min}),u^{\epsilon_{0}}+\mathcal{E}-u^{\delta}_{\rm
min}\rangle_{H^{p}(\Omega)}$ $\displaystyle\geq
C\|u^{\epsilon_{0}}+\mathcal{E}-u^{\delta}_{\min}\|^{2}_{H^{1}(\Omega)}+\eta\|u^{\epsilon_{0}}+\mathcal{E}-u^{\delta}_{\min}\|^{2}_{H^{p}(\Omega)}$
$\displaystyle\geq
C\|u^{\epsilon_{0}}+\mathcal{E}-u^{\delta}_{\min}\|^{2}_{H^{1}(\Omega)}.$
(4.19)

Applying (4.4), since $u^{\epsilon_{0}}$ solves (4.1), we have

$\displaystyle J_{\lambda,\beta,\eta}(u^{\epsilon_{0}}+\mathcal{E})$
$\displaystyle=J_{\lambda,\beta,\eta}(u^{\epsilon_{0}})+J_{\lambda,\beta,\eta}^{\prime}(u^{\epsilon_{0}})(\mathcal{E})+C\delta^{2}=\eta\|u^{\epsilon_{0}}\|^{2}_{H^{p}(\Omega)}+2\eta\langle
u^{\epsilon_{0}},\mathcal{E}\rangle_{H^{p}(\Omega)}+C\delta^{2}$
$\displaystyle\leq\eta\|u^{\epsilon_{0}}\|^{2}_{H^{p}(\Omega)}+\eta^{2}\|u^{\epsilon_{0}}\|_{H^{p}(\Omega)}^{2}+C\delta^{2}\leq
2\eta\|u^{\epsilon_{0}}\|^{2}_{H^{p}(\Omega)}+C\delta^{2}.$ (4.20)

Using the inequality $(a-b)^{2}\geq\frac{1}{2}a^{2}-b^{2}$ and combining
(4.19) and (4.20), we have

$\displaystyle\frac{1}{2}\|u^{\epsilon_{0}}-u^{\delta}_{\min}\|^{2}_{H^{1}(\Omega)}-\|\mathcal{E}\|^{2}_{H^{1}(\Omega)}\leq\|u^{\epsilon_{0}}+\mathcal{E}-u^{\delta}_{\min}\|^{2}_{H^{1}(\Omega)}\leq
C(\eta\|u^{\epsilon_{0}}\|^{2}_{H^{p}(\Omega)}+\delta^{2})$

The estimate (4.18) follows.  $\square$

Combining Theorem 4.14 and Theorem 4.18, we have for each $k\geq 1$,

$\|u^{(k)}-u^{\epsilon_{0}}\|_{H^{1}(\Omega)}\leq
C(\sqrt{\eta}\|u^{\epsilon_{0}}\|_{H^{p}(\Omega)}+\delta)+\theta^{k/2}\|u^{(0)}-u^{\delta}_{\min}\|_{H^{p}(\Omega)}$

for some constant $C>0$ and $\theta\in(0,1)$ where $u^{(k)}$ is the minimizing
sequence defined in Theorem 4.14. This inequality shows the stability of our
method with respect to noise. If $\theta^{k}$ and $\eta$ are $O(\delta^{2})$,
then the convergence rate is Lipschitz.

######  Remark 4.3.

In the case when the function $F$ is such that the comparison principle for
$-\epsilon_{0}\Delta u+F({\bf x},u,\nabla u)=0$ holds true; for e.g., $F$ is
strictly increasing with respect to its second variable (Assumption 2.1), one
can prove the stability of $u$ with respect to noise without imposing
Assumption 1.1 by the fact that $u-\delta$ and $u+\delta$ are a subsolution
and supersolution of $-\epsilon_{0}\Delta u+F({\bf x},u,\nabla u)=0$,
respectively. The main reason for us to successfully establish the stability
without assuming the comparison principle is due to the presence of the
Neumann data in Assumption 1.1.

In the next section, we present the numerical implementation and some
interesting numerical results.

##  5 Numerical study

We implement the convexification method based on the finite difference method.
For the simplicity, in this section, we only consider the case $d=2.$ On
$\overline{\Omega}=[-R,R]^{2}$, we arrange a uniform grid of points

$\mathcal{G}=\big{\\{}{\bf
x}_{ij}=(x_{i},z_{j}):x_{i}=-R+(i-1)\delta_{x},z_{j}=-R+(j-1)\delta_{z}),1\leq
i,j\leq N\big{\\}}$ (5.1)

where $\delta_{x}=\delta_{z}=h=2R/(N-1)$. In our computation, $R=1$ and
$N=50.$ The finite difference version of the objective function
$J_{\lambda,\beta,\eta}(u)$ is given by

$J_{\lambda,\beta,\eta}(u)=h^{2}\sum_{i,j=2}^{N-1}e^{2\lambda\big{|}\frac{z_{j}+r)}{b}\big{|}^{\beta}}\Big{|}-\epsilon_{0}\Delta^{h}u(x_{i},z_{j})+F((x_{i},z_{j}),u(x_{i},z_{j}),\partial_{x}^{h}u(x_{i},z_{j}),\partial_{z}^{h}u(x_{i},z_{j}))\Big{|}^{2}\\\
+\eta
h^{2}\Big{(}\sum_{i,j=1}^{N}|u(x_{i},z_{j})|^{2}+\sum_{i,j=2}^{N-1}\big{(}|\partial_{x}^{h}u(x_{i},z_{j})|^{2}+|\partial_{z}^{h}u(x_{i},z_{j})|^{2}+|\Delta^{h}u(x_{i},z_{j})|^{2}\big{)}\Big{)}.$
(5.2)

where

$\displaystyle\Delta^{h}u(x_{i},z_{j})$
$\displaystyle=\frac{u(x_{i},z_{j-1})+u(x_{i},z_{j+1})+u(x_{i-1},z_{j})+u(x_{i+1},z_{j})-4u(x_{i},z_{j})}{h^{2}}$
$\displaystyle\partial_{x}^{h}u(x_{i},z_{j})$
$\displaystyle=\frac{u(x_{i+1},z_{j})-u(x_{i-1},z_{j})}{2h}$
$\displaystyle\partial_{z}^{h}u(x_{i},z_{j})$
$\displaystyle=\frac{u(x_{i},z_{j+1})-u(x_{i},z_{j-1})}{2h}.$

######  Remark 5.1.

1\. In our computation $R=1$, $N=50$, $\epsilon_{0}=10^{-3}$, $\lambda=2,$
$\beta=8$, $b=10$, $r=1.2$ and $\eta=10^{-4}$. Although in the theoretical
part, the parameters $\lambda$ and $\beta$ should be large, these values are
already good for the numerical part. For the simplicity, we use this set of
parameters for all tests below.

2\. In theory $p>3$ when $d=2$. However, in numerical study, we can reduce the
norm in the regularization term to $p=2$ to simplify the implementation and to
improve the speed of computation. We do not experience any difficulty with
this small change.

In our implementation, instead of writing the computational code for the
gradient descent method, we use the optimization toolbox of Matlab, in which
the gradient descent method is coded. More precisely, we use the command
“fmincon” to minimize the functional $J_{\lambda,\beta,\eta}$ subject to the
boundary constraint in (4.1). The command “fmincon” requires an initial
solution, that is the function $u^{(0)}$ in Theorem 4.14. The function
$u^{(0)}$ is naturally assigned to be the zero function. This function
$u^{(0)}$ does not satisfy the Dirichlet boundary condition on
$\partial\Omega$ and Neumann condition on $\Gamma^{+}$. However, the command
“fmincon” corrects this error automatically.

We present here five (5) numerical tests in which the given boundary data are
noisy with noise level $\delta=0\%$, $5\%$, and $10\%$, respectively. In each
test, $u_{\rm true}$ and $u_{\rm comp}$ denote the true and computed viscosity
solutions, respectively. Given functions $f,g$, the noisy versions of $f,g$
are given by

$f^{\delta}=f(1+\delta\cdot\mbox{rand})\quad\text{ and }\quad
g^{\delta}=g(1+\delta\cdot\mbox{rand}),$

where rand is the function that generates uniformly distributed random numbers
in $[-1,1].$ The relative computed error is defined as

${\rm err}(\delta)=\frac{\|u_{\rm comp}-u_{\rm
true}\|_{L^{\infty}(\Omega)}}{\|u_{\rm
true}\|_{L^{\infty}(\Omega)}}\quad\text{ for }\delta>0.$

Test 1. In this test, we find the viscosity solution to (1.1)–(1.2) where

$F({\bf x},s,{\bf p})=\frac{1}{150}s+|{\bf
p}|+\frac{1}{150}(x^{2}+z^{2})-2\sqrt{x^{2}+z^{2}}\quad\mbox{for all }{\bf
x}\in\Omega,s\in\mathbb{R},{\bf p}\in\mathbb{R}^{2}$ (5.3)

and the boundary data are given by

$u({\bf x})=f({\bf x})=-(x^{2}+z^{2})\quad\mbox{for all }{\bf
x}=(x,z)\in\partial\Omega$ (5.4)

and

$u_{z}({\bf x})=g({\bf x})=-2z\quad\mbox{for all }{\bf x}=(x,z)\in\Gamma^{+}.$
(5.5)

The true solution is $u_{\rm true}(x,z)=-(x^{2}+z^{2})$, which is smooth in
$\Omega$. We are here in a standard setting, and the convergence of
$u^{\epsilon}$, the solution to (2.1), to the solution $u_{\rm true}$ to
(1.1)–(1.2) is guaranteed in [15, 14, 40, 5, 4, 56].

(a) The true solution $u_{\rm true}=-(x^{2}+z^{2})$.

(b) The solution $u_{\rm comp}$, computed from noiseless boundary data.

(c) The solution $u_{\rm comp}$, computed from 5% noisy boundary data.

(d) The solution $u_{\rm comp}$, computed from 10% noisy boundary data.

(e) The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=0\%$.

(f) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=5\%$.

(g) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=10\%$.

Figure 1:  Test 1. The true and computed viscosity solutions with $\delta$ is
$0\%$, 5%, and 10% noisy boundary Dirichlet data on $\partial\Omega$ and
Neumann data on $\Gamma^{+}$.

We show in Figure 1 the graphs of this true solution and that of the computed
ones from noiseless and noisy boundary data. It is evident that we
successfully obtained computed viscosity solutions to (1.1)–(1.2). By adding
high level noise into the boundary data, we have numerically shown that the
convexification is stable. It is evident from Table 1 that the relative
computed error is about the noise level, which clearly illustrates the
Lipschitz stability in Theorem 4.18. In this test, the true solution is
smooth. The function $F({\bf x},s,{\bf p})$ is strictly increasing with
respect to $s$.

Table 1:  Test 1. The performance of the convexification method. The computational time is the time for a Precisions Workstations T7810 with 24 cores to compute the solution $u_{\rm comp}$. In this table, the relative $L^{\infty}(\Omega)$ error is $\|u_{\rm comp}-u_{\rm true}\|_{L^{\infty}(\Omega)}/\|u_{\rm true}\|_{L^{\infty}(\Omega)}$. Noise level | computational time | number of iterations | relative $L^{\infty}(\Omega)$ error
---|---|---|---
0% | 23.47 minutes | 279 | $0.24\%$
5% | 24.40 minutes | 292 | 4.51%
10% | 27.35 minutes | 329 | 9.95%

Test 2. We now solve the eikonal equation of which the function $F$ in (1.1)
is not in the class $C^{1}$. In this test, the function $F$ is given by

$F({\bf x},s,{\bf p})=|{\bf p}|-\sqrt{2}\quad\mbox{for all }{\bf
x}\in\Omega,s\in\mathbb{R},{\bf p}\in\mathbb{R}^{2}$ (5.6)

and the boundary data are

$u({\bf x})=f({\bf x})=-(|x|+|z|)\quad\mbox{for all }{\bf
x}=(x,z)\in\partial\Omega$ (5.7)

and

$u_{z}({\bf x})=g({\bf x})=\left\\{\begin{array}[]{rl}1&z<0,\\\
-1&z>0,\end{array}\right.\quad\mbox{for all }{\bf x}=(x,z)\in\Gamma^{+}.$
(5.8)

We claim that the true solution to (1.1)–(1.2) is $u_{\rm true}({\bf
x})=-(|x|+|z|)$ for all ${\bf x}=(x,z)\in\Omega.$ Intuitively, this claim
holds as the graph of $u_{\rm true}$ only has corners from above and $F$ is
convex in ${\bf p}$. Let us provide a rigorous verification here. If $x\neq 0$
and $z\neq 0$, then $u_{\rm true}$ is differentiable at ${\bf x}=(x,z)$, and

$\nabla u_{\rm true}({\bf
x})=\left(-\frac{x}{|x|},-\frac{z}{|z|}\right)\quad\Rightarrow\quad|\nabla
u_{\rm true}({\bf x})|=\sqrt{2}.$

If $x=0$ or $z=0$, then $u_{\rm true}$ is not differentiable at ${\bf
x}=(x,z)$. We can only find smooth test functions that touch $u_{\rm true}$
from above at ${\bf x}$, and we cannot find smooth test functions that touch
$u_{\rm true}$ from below at ${\bf x}$. Let $\phi$ be a smooth test function
that touches $u_{\rm true}$ from above at ${\bf x}$. Without loss of
generality, we only need to consider the case $x=0$. If $z\neq 0$, then we
have that

$\phi_{x}({\bf x})\in[-1,1],\ \phi_{z}({\bf x})=-\frac{z}{|z|}.$

If $z=0$, then we have that

$\phi_{x}({\bf x})\in[-1,1],\ \phi_{z}({\bf x})\in[-1,1].$

In both cases,

$|\nabla\phi({\bf x})|\leq\sqrt{2}.$

Thus, the subsolution test holds for $u_{\rm true}$ at ${\bf x}$. We conclude
that $u_{\rm true}$ is a viscosity solution to (1.1)–(1.2). This true solution
and its computed versions $u_{\rm comp}$ from noisy boundary data are
displayed in Figure 2. The convergence of $u^{\epsilon}$, the solution to
(2.1), to the solution $u_{\rm true}$ to (1.1)–(1.2) is guaranteed in [20, 55]
with convergence rate $O(\epsilon^{1/2})$.

(a) The true solution $u_{\rm true}=-(|x|+|z|)$.

(b) The solution $u_{\rm comp}$, computed from noiseless boundary data.

(c) The solution $u_{\rm comp}$, computed from 5% noisy boundary data.

(d) The solution $u_{\rm comp}$, computed from 10% noisy boundary data.

(e) The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=0\%$.

(f) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=5\%$.

(g) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=10\%$.

Figure 2:  Test 2. The true and computed viscosity solutions with $\delta$ is
$0\%$, 5%, and 10% noisy boundary Dirichlet data on $\partial\Omega$ and
Neumann data on $\Gamma^{+}$.

This test is more challenging than Test 1. In this test, the functions $F$ and
$f$ are both not smooth. So, the smoothness condition for the theoretical part
does not satisfy. However, the convexification method still provides reliable
solutions even when the given boundary is noisy with the noise level
$\delta=10\%.$ This shows the robustness of the convexification method. The
relative errors are acceptable, see Table 2.

Table 2:  Test 2. The performance of the convexification method. The computational time is the time for a Precisions Workstations T7810 with 24 cores to compute the solution $u_{\rm comp}$. In this table, the relative $L^{\infty}(\Omega)$ error is $\|u_{\rm comp}-u_{\rm true}\|_{L^{\infty}(\Omega)}/\|u_{\rm true}\|_{L^{\infty}(\Omega)}$. Noise level | computational time | number of iterations | relative $L^{\infty}(\Omega)$ error
---|---|---|---
0% | 48 minutes | 581 | 3.75%
5% | 89 minutes | 1071 | 5.23%
10% | 64 minutes | 766 | 12.22 %

Test 3. We next consider a more complicated Hamilton-Jacobi equation. Unlike
the previous two tests, the function $F$ in this test is not convex and; more
interestingly, not coercive and not continuous. It is given by

$F({\bf x},s,{\bf p})=20s+|p_{1}|-|p_{2}|\\\
-\Big{(}20\big{(}-|x-0.5|+e^{\sin(\pi(x^{2}+z^{2}))}\big{)}+G(x,z)-\Big{|}2\pi
z\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))}\Big{|}\Big{)}$ (5.9)

for all ${\bf x}=(x,z)\in\Omega,s\in\mathbb{R},{\bf
p}=(p_{1},p_{2})\in\mathbb{R}^{2}$ where

$G(x,z)=\left\\{\begin{array}[]{ll}\big{|}1+2\pi
x\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))}\big{|}&x<0.5,\\\
\big{|}-1+2\pi
x\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))}\big{|}&x>0.5.\end{array}\right.$

Note that $G$ and $F$ are discontinuous at $x=0.5$ in general. The boundary
data are given by

$u({\bf x})=f({\bf x})=-|x-0.5|+e^{\sin(\pi(x^{2}+z^{2}))}\quad\mbox{for all
}{\bf x}=(x,z)\in\partial\Omega$ (5.10)

and

$u_{z}({\bf x})=g({\bf x})=2\pi
z\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))}\quad\mbox{for all }{\bf
x}=(x,z)\in\Gamma^{+}.$ (5.11)

The function $u_{\rm true}(x,z)=-|x-0.5|+e^{\sin(\pi(x^{2}+z^{2}))}$ is the
true viscosity solution for this test. Indeed, if $x\neq 0.5$, then $u_{\rm
true}$ is differentiable at ${\bf x}=(x,z)$, and

$\nabla u_{\rm true}({\bf x})=\left(-\frac{x-0.5}{|x-0.5|}+2\pi
x\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))},2\pi
z\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))}\right),$

which gives that $F({\bf x},u_{\rm true}({\bf x}),\nabla u_{\rm true}({\bf
x}))=0$. If $x=0.5$, then $u_{\rm true}$ is not differentiable at ${\bf
x}=(x,z)$. We can only find smooth test functions that touch $u_{\rm true}$
from above at ${\bf x}$, and we cannot find smooth test functions that touch
$u_{\rm true}$ from below at ${\bf x}$. Let $\phi$ be a smooth test function
that touches $u_{\rm true}$ from above at ${\bf x}$. Then,

$\phi_{x}({\bf x})\in 2\pi
x\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))}+[-1,1],\ \phi_{z}({\bf
x})=2\pi z\cos(\pi(x^{2}+z^{2}))e^{\sin(\pi(x^{2}+z^{2}))},$

which yields that $F_{*}({\bf x},u_{\rm true}({\bf x}),\nabla\phi({\bf
x}))\leq 0$. Here, $F_{*}$ is the lower semicontinuous envelope of $F$.
Therefore, the subsolution test holds for $u_{\rm true}$ at ${\bf x}$. We
imply that $u_{\rm true}$ is a viscosity solution to (1.1)–(1.2).

(a) The true solution $u_{\rm true}(x,z)=-|x-0.5|+e^{\sin(\pi(x^{2}+z^{2}))}$.

(b) The solution $u_{\rm comp}$, computed from noiseless boundary data.

(c) The solution $u_{\rm comp}$, computed from 5% noisy boundary data.

(d) The solution $u_{\rm comp}$, computed from 10% noisy boundary data.

(e) The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=0\%$.

(f) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=5\%$.

(g) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=10\%$.

Figure 3:  Test 3. The true and computed viscosity solutions with $\delta$ is
$0\%$, 5%, and 10% noisy boundary Dirichlet data on $\partial\Omega$ and
Neumann data on $\Gamma^{+}$.

The numerical results are given in Figure 3. As mentioned, this test is
interesting since the function $F$, see (5.9), is nonconvex, noncoercive, and
discontinuous. Solving the Hamilton-Jacobi equation with this Hamiltonian is
challenging. Some existing methods might not be applicable. In contrast, the
numerical results in Figure 3 are out of expectation. The errors of
computation are small, see Table 3, although the solution has complicated
structure. This kind of nonconvex Hamiltonian occurs in the context of two-
player zero-sum differential games (see [4, 56]).

Table 3: Test 3. The performance of the convexification method. The computational time is the time for a Precisions Workstations T7810 with 24 cores to compute the solution $u_{\rm comp}$. In this table, the relative $L^{\infty}(\Omega)$ error is $\|u_{\rm comp}-u_{\rm true}\|_{L^{\infty}(\Omega)}/\|u_{\rm true}\|_{L^{\infty}(\Omega)}$. Noise level | computational time | number of iterations | relative $L^{\infty}(\Omega)$ error
---|---|---|---
0% | 10.72 minutes | 129 | 1.32%
5% | 11.12 minutes | 131 | 2.04%
10% | 8.43 minutes | 101 | 3.97%

######  Remark 5.2.

In Tests 1, 2 and 3 above, we are in the context that the knowledge of $u_{z}$
on $\Gamma^{+}$ can be computed from the knowledge of $u$ and the form of the
Hamilton-Jacobi equation, see Assumption 1.1 and Remark 1.1. However, if the
given Hamilton-Jacobi equation is rather complicated as in Tests 4 and 5
below, solving $u_{z}|_{\Gamma^{+}}$ from $u|_{\Gamma^{+}}$ is impossible. In
this case, we minimize $J_{\lambda,\beta,\eta}$ on the set $\\{u\in
H^{p}(\Omega):u({\bf x})=f({\bf x})\mbox{ for }{\bf x}\in\partial\Omega\\}$.
The convexification method still gives us out of expectation numerical results
in these two tests. However, the proofs of the convexification theorem and the
convergence of the numerical scheme for the problem with only Dirichlet
boundary condition (1.1)–(1.2) are extremely challenging, and they are out of
the scope of this paper.

Test 4. We next consider the G-equation, which arises from instantaneous flame
position. We solve (1.1)–(1.2) when

$F({\bf x},s,{\bf p})=s+|{\bf p}|-xp_{1}$ (5.12)

for ${\bf x}\in\Omega,s\in\mathbb{R},{\bf p}=(p_{1},p_{2})\in\mathbb{R}^{2}$
and the boundary data is given by

$f({\bf x})=-|x|-1\quad{\bf x}=(x,z)\in\partial\Omega.$ (5.13)

The function $u_{\rm true}(x,z)=-|x|-1$ is the true viscosity solution for
this test. The verification of this is similar to that of Tests 2 and 3, and
is hence omitted here. The numerical results are given in Figure 4. Relative
errors are $0.91\%$, $4.97\%$ and $9.99\%$ when $\delta$ is 0%, 5% and $10\%$
respectively (see Table 4). The G-equation is quite popular in the combustion
science literature. We refer the readers to [12, 59, 41] for some recent
important mathematical developments.

(a) The true solution $u_{\rm true}=-|x|-1$.

(b) The solution $u_{\rm comp}$, computed from noiseless boundary data.

(c) The solution $u_{\rm comp}$, computed from 5% noisy boundary data.

(d) The solution $u_{\rm comp}$, computed from 10% noisy boundary data.

(e) The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=0\%$.

(f) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=5\%$.

(g) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=10\%$.

Figure 4:  Test 4. The true and computed viscosity solutions with $\delta$ is $0\%$, 5% and 10% noisy boundary Dirichlet data on $\partial\Omega$. Table 4: Test 4. The performance of the convexification method. The computational time is the time for a Precisions Workstations T7810 with 24 cores to compute the solution $u_{\rm comp}$. The relative $L^{\infty}(\Omega)$ error is $\|u_{\rm comp}-u_{\rm true}\|_{L^{\infty}(\Omega)}/\|u_{\rm true}\|_{L^{\infty}(\Omega)}$. Noise level | computational time | number of iterations | relative $L^{\infty}(\Omega)$ error
---|---|---|---
0% | 22.62 minutes | 280 | 0.91%
5% | 24.99 minutes | 311 | 4.97%
10% | 26.50 minutes | 322 | 9.99%

Test 5. We finally consider a quite complicated form of the function $F$. For
${\bf x}=(x,z)\in\Omega$, let

$G(x,z)=\left\\{\begin{array}[]{ll}1+2\pi\cos(2\pi(x+z))&x>0,\\\
-1+2\pi\cos(2\pi(x+z))&x<0.\end{array}\right.$

We solve (1.1)–(1.2) when

$F({\bf x},s,{\bf p})=15s+\min\left\\{|{\bf p}|,||{\bf
p}|-10|+6\right\\}-\Big{[}15(-|x|+\sin(2\pi(x+z)))\\\
+\min\left\\{|\sqrt{|G(x,z)|^{2}+4\pi^{2}\cos(2\pi(x+z))}|,\big{|}\sqrt{|G(x,z)|^{2}+4\pi^{2}\cos(2\pi(x+z))}-10\big{|}+6\right\\}\Big{]}$
(5.14)

for ${\bf x}\in\Omega,s\in\mathbb{R},{\bf p}\in\mathbb{R}^{2}$, and the
Dirichlet boundary data is given by

$f({\bf x})=-|x|+\sin(2\pi(x+z))\quad\text{ for all }{\bf
x}=(x,z)\in\partial\Omega.$ (5.15)

It it worth mentioning that $F,G$ are not continuous in general at $x=0$. The
function $u_{\rm true}(x,z)=-|x|+\sin(2\pi(x+z))$ is the true viscosity
solution for this test. Its graph and the graphs of the computed solutions are
displayed in Figure 5. Let us give a rigorous verification here. If $x\neq 0$,
then $u_{\rm true}$ is differentiable at $(x,z)$, and

$\nabla u_{\rm
true}(x,z)=\left(-\frac{x}{|x|}+2\pi\cos(2\pi(x+z)),2\pi\cos(2\pi(x+z))\right),$

which implies that $F({\bf x},u({\bf x}),\nabla u_{\rm true}({\bf x}))=0$. If
$x=0$, then $u_{\rm true}$ is not differentiable at ${\bf x}=(x,z)$. We can
only find smooth test functions that touch $u_{\rm true}$ from above at ${\bf
x}$, and we cannot find smooth test functions that touch $u_{\rm true}$ from
below at ${\bf x}$. Let $\phi$ be a smooth test function that touches $u_{\rm
true}$ from above at ${\bf x}$. Then,

$\phi_{x}({\bf x})\in 2\pi\cos(2\pi z)+[-1,1],\ \phi_{z}({\bf
x})=2\pi\cos(2\pi z),$

which yields that

$|\nabla\phi({\bf x})|^{2}\leq
1+4\pi+4\pi^{2}\quad\Rightarrow\quad|\nabla\phi({\bf x})|\leq 1+2\pi<8.$

Therefore,

$\displaystyle F_{*}({\bf x},u_{\rm true}({\bf x}),\nabla\phi({\bf x}))$
$\displaystyle=10\sin(2\pi z)+|\nabla\phi({\bf x})|-10\sin(2\pi
z)-\sqrt{|G|^{*}(0,z)+4\pi^{2}\cos(2\pi z)}$ $\displaystyle=|\nabla\phi({\bf
x})|-\sqrt{1+2\pi|\cos(2\pi z)|+4\pi^{2}}\leq 0.$

Here, $F_{*}$ is the lower semicontinuous envelope of $F$, and $|G|^{*}$ is
the upper semicontinuous envelope of $|G|$. Thus, the subsolution test holds
for $u_{\rm true}$ at ${\bf x}$. We imply that $u_{\rm true}$ is a viscosity
solution to (1.1)–(1.2).

It is evidently clear that the convexification method delivers satisfactory
numerical results. The relative errors are provided in Table 5. This kind of
Hamiltonian was considered in [48] in the context of the periodic
homogenization theory.

(a) The true solution $u_{\rm true}=-|x|+\sin(2\pi(x+z))$.

(b) The solution $u_{\rm comp}$, computed from noiseless boundary data.

(c) The solution $u_{\rm comp}$, computed from 5% noisy boundary data.

(d) The solution $u_{\rm comp}$, computed from 10% noisy boundary data.

(e) The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=0\%$.

(f) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=5\%$.

(g) [The relative error $\frac{|u_{\rm comp}-u_{\rm true}|}{\|u_{\rm
true}\|_{L^{\infty}}}$, $\delta=10\%$.

Figure 5:  Test 5. The true and computed viscosity solutions with $\delta$ is $0\%$, 5%, and 10% noisy boundary Dirichlet data on $\partial\Omega$. The relative errors with these noise levels are given in Table 5. Table 5: Test 5. The performance of the convexification method. The computational time is the time for a Precisions Workstations T7810 with 24 cores to compute the solution $u_{\rm comp}$. In this table, the relative $L^{\infty}(\Omega)$ error is $\|u_{\rm comp}-u_{\rm true}\|_{L^{\infty}(\Omega)}/\|u_{\rm true}\|_{L^{\infty}(\Omega)}$. Noise level | computational time | number of iterations | relative $L^{\infty}(\Omega)$ error
---|---|---|---
0% | 11.83 minutes | 144 | 1.78%
5% | 10.06 minutes | 121 | 4.97%
10% | 11.54 minutes | 143 | 9.77%

######  Remark 5.3.

It follows from tests 3 and 5 that the convexification method is effective in
the interesting case when the function $F$ is nonconvex in ${\bf p}.$ The non-
convexity is illustrated in Figure 6. This numerically confirms the strength
of our method in solving Hamilton-Jacobi equations. We refer the readers to
[23, 39] for some other examples dealing with non-convex, discontinuous
Hamiltonians via the Lax-Friedrichs sweeping method.

(a) ${\bf p}\mapsto|p_{1}|-|p_{2}|$.

(b) ${\bf p}\mapsto\min\left\\{|{\bf p}|,||{\bf p}|-10|+6\right\\}$

Figure 6:  The dependence of the function $F({\bf x},s,{\bf p})$ on the third
variable ${\bf p}$ in Test 3 (a) and Test 5 (b). It is evident that in those
tests, the function $F$ is nonconvex in ${\bf p}.$

######  Remark 5.4.

The relative errors in computation are displayed in the last rows of Figures
1–5. It is evident from those figures that the errors in computation occur at
$\partial\Omega$ where the noise is added and at the places where the true
solution $u_{\rm true}$ is not differentiable. This again reflects the
strength of the convexification method.

It has been shown both analytically and numerically that the convexification
method is robust in solving a large class of Hamilton-Jacobi equations. The
strengths of the convexification method involves the facts

  1. 1.

that it does not require any special structure of the Hamiltonian; especially,
the convexity condition of the Hamiltonian with respect to the third variable
is relaxed;

  2. 2.

that it yields the satisfactory numerical solutions even when the given
boundary is noisy.

However, the convexification method has a drawback. It is time consuming in
comparison to the well-known methods for nonconvex Hamiltonians; for e.g., the
Lax–Friedrichs schemes ([47, 1, 45]) and the Lax–Friedrichs sweeping algorithm
([23, 39]). In our tests, it takes from 10 minutes to 90 minutes, depending on
the forms of the given Hamiltonians, for a Workstations T7810 with 24 cores to
compute the solutions (see Tables 1–5). The slow computational time is
acceptable in the sense that we consider the computational program as a “proof
of concept” to numerically confirm the analysis of the convexification method.
The convexification method is the first generalization of the numerical method
based on Carleman estimates to solve Hamilton-Jacobi equations. We expect to
improve the computational time in the next generation. The next generation
will be developed based on the fixed point iterative scheme similar to the
ones in [35, 36, 42] for quasi-linear elliptic and hyperbolic equations. The
rate of convergence in those papers is $O(\theta^{n})$ as $n\to\infty$ for
some $\theta\in(0,1).$ Hence, the success in reducing computational time is
very promising.

##  6 Concluding remarks

In this paper, we introduce a new method to solve highly nonlinear Hamilton-
Jacobi equations in a rectangular domain. This method is called the
convexification. The key idea of the convexification is to involve a Carleman
weight function into a cost functional defined directly from the equation
under consideration. Using a Carleman estimate, we established some important
theoretical results. The first theorem guarantees that this cost functional is
strictly convex and has a unique solution. Then, we proved in the second
theorem that the gradient descent method with sufficiently small step size
delivers a sequence converging to the unique minimizer. Then, in the third
theorem, we prove that the minimizer above converges to the solution in the
vanishing viscosity process, a good approximation of the viscosity solution to
the Hamilton-Jacobi equation, as the noise contained in the boundary data
tends to $0$. The rate of the convergence is Lipschitz. All theoretical
results are valid in the framework that we know the value of the solution on
the boundary of the domain and its normal derivative in a part of this
boundary. We have pointed out that this framework is acceptable in some real-
world applications. We have shown some interesting numerical tests in 2D.
These tests confirm the convergence of the convexification method even when
the Hamiltonian is not convex or discontinuous. Moreover, these numerical
results are out of expectation even when the solved equations are not in the
framework above.

As of now, the convexification method is more time consuming in comparison to
the well-known methods as addressed by the end of Section 5. We expect to
improve the computational time in the next generation. Besides, we also intend
to study Hamilton-Jacobi equations with other types of boundary conditions in
the near future.

## Acknowledgement

The works of MVK and LHN were supported by US Army Research Laboratory and US
Army Research Office grant W911NF-19-1-0044. The work of HT is supported in
part by NSF grant DMS-1664424 and NSF CAREER grant DMS-1843320.

## References

  * [1] R. Abgrall.  Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes.  Comm. Pure Appl. Math., 49(12):1339–1373, 1996.
  * [2] R. Abgrall.  Numerical discretization of boundary conditions for first order Hamilton-Jacobi equations.  SIAM J. Numer. Anal., 41(6):2233–2261, 2003.
  * [3] A. B. Bakushinskii, M. V. Klibanov, and N. A. Koshev.  Carleman weight functions for a globally convergent numerical method for ill-posed Cauchy problems for some quasilinear PDEs.  Nonlinear Anal. Real World Appl., 34:201–224, 2017.
  * [4] M. Bardi and I. Capuzzo-Dolcetta.  Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations.  Systems & Control: Foundations & Applications. Birkhäuser Boston, Inc., Boston, MA, 1997.  With appendices by Maurizio Falcone and Pierpaolo Soravia.
  * [5] G. Barles.  Solutions de viscosité des équations de Hamilton-Jacobi, volume 17 of Mathématiques & Applications (Berlin) [Mathematics & Applications].  Springer-Verlag, Paris, 1994.
  * [6] G. Barles and P. E. Souganidis.  Convergence of approximation schemes for fully nonlinear second order equations.  Asymptotic Anal., 4(3):271–283, 1991.
  * [7] S. Bryson and D. Levy.  High-order central WENO schemes for multidimensional Hamilton-Jacobi equations.  SIAM J. Numer. Anal., 41(4):1339–1369, 2003.
  * [8] A. L. Bukhgeim and Michael V. Klibanov.  Global uniqueness of a class of multidimensional inverse problems.  Soviet Math. Doklady, 17:244–247, 1981.
  * [9] F. Cagnetti, D. Gomes, and H. V. Tran.  Adjoint methods for obstacle problems and weakly coupled systems of PDE.  ESAIM Control Optim. Calc. Var., 19(3):754–779, 2013.
  * [10] F. Cagnetti, D. Gomes, and H. V. Tran.  Convergence of a semi-discretization scheme for the Hamilton-Jacobi equation: a new approach with the adjoint method.  Appl. Numer. Math., 73:2–15, 2013.
  * [11] F. Camilli, I. Capuzzo Dolcetta, and D. A. Gomes.  Error estimates for the approximation of the effective Hamiltonian.  Appl. Math. Optim., 57(1):30–57, 2008.
  * [12] P. Cardaliaguet, J. Nolen, and P. E. Souganidis.  Homogenization and enhancement for the $G$-equation.  Arch. Ration. Mech. Anal., 199(2):527–561, 2011.
  * [13] Bernardo Cockburn, Ivan Merev, and Jianliang Qian.  Local a posteriori error estimates for time-dependent Hamilton-Jacobi equations.  Math. Comp., 82(281):187–212, 2013.
  * [14] M. G. Crandall, L. C. Evans, and P.-L. Lions.  Some properties of viscosity solutions of Hamilton-Jacobi equations.  Trans. Amer. Math. Soc., 282(2):487–502, 1984.
  * [15] M. G. Crandall and P.-L. Lions.  Viscosity solutions of Hamilton-Jacobi equations.  Trans. Amer. Math. Soc., 277(1):1–42, 1983.
  * [16] M. G. Crandall and P.-L. Lions.  Two approximations of solutions of Hamilton-Jacobi equations.  Math. Comp., 43(167):1–19, 1984.
  * [17] P. Daniel and J.-D. Durou.  From deterministic to stochastic methods for shape from shading.  In In Proc. 4th Asian Conf. on Comp. Vis, pages 187–192, 2000.
  * [18] M. Falcone and R. Ferretti.  Semi-Lagrangian schemes for Hamilton-Jacobi equations, discrete representation formulae and Godunov methods.  J. Comput. Phys., 175(2):559–575, 2002.
  * [19] M. Falcone and R. Ferretti.  Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2014.
  * [20] W. H. Fleming and P. E. Souganidis.  Asymptotic series and the method of vanishing viscosity.  Indiana Univ. Math. J., 35(2):425–447, 1986.
  * [21] D. Gallistl, T. Sprekeler, and E. Süli.  Mixed finite element approximation of periodic Hamilton–Jacobi–Bellman problems with application to numerical homogenization, 2020.
  * [22] B. K. P. Horn and M. J. Brooks.  The variational approach to shape from shading.  Computer Vision, Graphics, and Image Processing, 33(2):174–208, 1986.
  * [23] C. Y. Kao, S. Osher, and J. Qian.  Lax-Friedrichs sweeping scheme for static Hamilton-Jacobi equations.  J. Comput. Phys., 196(1):367–391, 2004.
  * [24] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astratov.  Convexification and experimental data for a 3D inverse scattering problem with the moving point source.  Inverse Problems, 36:085007, 2020.
  * [25] V. A. Khoa, M. V. Klibanov, and L. H. Nguyen.  Convexification for a 3D inverse scattering problem with the moving point source.  SIAM J. Imaging Sci., 13(2):871–904, 2020.
  * [26] M. V. Klibanov.  Global convexity in a three-dimensional inverse acoustic problem.  SIAM J. Math. Anal., 28:1371–1388, 1997.
  * [27] M. V. Klibanov.  Global convexity in diffusion tomography.  Nonlinear World, 4:247–265, 1997.
  * [28] M. V. Klibanov.  Carleman weight functions for solving ill-posed Cauchy problems for quasilinear PDEs.  Inverse Problems, 31:125007, 2015.
  * [29] M. V. Klibanov.  Travel time tomography with formally determined incomplete data in 3D.  Inverse Problems and Imaging, 13:1367–1393, 2019.
  * [30] M. V. Klibanov and O. V. Ioussoupova.  Uniform strict convexity of a cost functional for three-dimensional inverse scattering problem.  SIAM J. Math. Anal., 26:147–179, 1995.
  * [31] M. V. Klibanov, V. A. Khoa, A. V. Smirnov, L. H. Nguyen, G. W. Bidney, L. Nguyen, A. Sullivan, and V. N. Astratov.  Nonlinear SAR imaging via a globally convergent inverse algorithm.  to appear on J. Applied and Industrial Mathematics, see also Arxiv.2103.10431, 2021.
  * [32] M. V. Klibanov, T. T. Le, L. H. Nguyen, A. Sullivan, and L. Nguyen.  Convexification-based globally convergent numerical method for a 1D coefficient inverse problem with experimental data.  to appear on Inverse Problems and Imaging, see also Arxiv:2104.11392, 2021.
  * [33] M. V. Klibanov and J. Li.  Inverse Problems and Carleman Estimates: Global Uniqueness, Global Convergence and Experimental Data.  to appear on De Gruyter, 2021.
  * [34] M. V. Klibanov, J. Li, and W. Zhang.  Convexification for an inverse parabolic problem.  Inverse Problems, 36:085008, 2020.
  * [35] T. T. Le, M. V. Klibanov, L. H. Nguyen, A. Sullivan, and L. Nguyen.  Carleman contraction mapping for a 1D inverse scattering problem with experimental time-dependent data.  preprint arXiv:2109.11098, 2021.
  * [36] T. T. Le and L. H. Nguyen.  A convergent numerical method to recover the initial condition of nonlinear parabolic equations from lateral Cauchy data.  Journal of Inverse and Ill-posed Problems, DOI: https://doi.org/10.1515/jiip-2020-0028, 2020.
  * [37] T. T. Le and L. H. Nguyen.  The gradient descent method for the convexification to solve boundary value problems of quasi-linear PDEs and a coefficient inverse problem.  preprint, Arxiv:2103.04159, 2021.
  * [38] Y. G. Leclerc and A. F. Bobick.  The direct computation of height from shading.  In Proceedings. 1991 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, pages 552–558, 1991.
  * [39] W. Li and J. Qian.  Newton-type Gauss-Seidel Lax-Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization.  Comput. Math. Appl., 79(4):1222–1239, 2020.
  * [40] P.-L. Lions.  Generalized solutions of Hamilton-Jacobi equations, volume 69 of Research Notes in Mathematics.  Pitman (Advanced Publishing Program), Boston, Mass.-London, 1982.
  * [41] Y.-Y. Liu, J. Xin, and Y. Yu.  Asymptotics for turbulent flame speeds of the viscous G-equation enhanced by cellular and shear flows.  Arch. Ration. Mech. Anal., 202(2):461–492, 2011.
  * [42] L. H. Nguyen and M. V. Klibanov.  Carleman estimates and the contraction principle for an inverse source problem for nonlinear hyperbolic equations.  preprint arXiv:2108.03500, 2021.
  * [43] L. H. Nguyen, Q. Li, and M. V. Klibanov.  A convergent numerical method for a multi-frequency inverse source problem in inhomogenous media.  Inverse Problems and Imaging, 13:1067–1094, 2019.
  * [44] A. M. Oberman and T. Salvador.  Filtered schemes for Hamilton-Jacobi equations: a simple construction of convergent accurate difference schemes.  J. Comput. Phys., 284:367–388, 2015.
  * [45] S. Osher and R. Fedkiw.  Level set methods and dynamic implicit surfaces, volume 153 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2003.
  * [46] S. Osher and J. A. Sethian.  Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations.  J. Comput. Phys., 79(1):12–49, 1988.
  * [47] S. Osher and C.-W. Shu.  High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations.  SIAM J. Numer. Anal., 28(4):907–922, 1991.
  * [48] J. Qian, H. V. Tran, and Y. Yu.  Min-max formulas and other properties of certain classes of nonconvex effective Hamiltonians.  Math. Ann., 372(1-2):91–123, 2018.
  * [49] J. Qian, Y.-T. Zhang, and H.-K. Zhao.  A fast sweeping method for static convex Hamilton-Jacobi equations.  J. Sci. Comput., 31(1-2):237–271, 2007.
  * [50] J. A. Sethian.  Level set methods and fast marching methods, volume 3 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, second edition, 1999.  Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science.
  * [51] J. A. Sethian and A. Vladimirsky.  Ordered upwind methods for static Hamilton-Jacobi equations: theory and algorithms.  SIAM J. Numer. Anal., 41(1):325–363, 2003.
  * [52] A. V. Smirnov, M. V. Klibanov, and L. H. Nguyen.  Convexification for a 1D hyperbolic coefficient inverse problem with single measurement data.  Inverse Probl. Imaging, 14(5):913–938, 2020.
  * [53] P. E. Souganidis.  Approximation schemes for viscosity solutions of Hamilton-Jacobi equations.  J. Differential Equations, 59(1):1–43, 1985.
  * [54] R. Szeliski.  Fast shape from shading.  In O. Faugeras, editor, Computer Vision — ECCV 90, pages 359–368, Berlin, Heidelberg, 1990. Springer Berlin Heidelberg.
  * [55] H. V. Tran.  Adjoint methods for static Hamilton-Jacobi equations.  Calc. Var. Partial Differential Equations, 41(3-4):301–319, 2011\.
  * [56] H.V. Tran.  Hamilton–Jacobi Equations: Theory and Applications.  Graduate Studies in Mathematics. American Mathematical Society, 2021.
  * [57] Y.-H. R. Tsai, L.-T. Cheng, S. Osher, and H.-K. Zhao.  Fast sweeping algorithms for a class of Hamilton-Jacobi equations.  SIAM J. Numer. Anal., 41(2):673–694, 2003.
  * [58] J. N. Tsitsiklis.  Efficient algorithms for globally optimal trajectories.  IEEE Trans. Automat. Control, 40(9):1528–1538, 1995.
  * [59] J. Xin and Y. Yu.  Periodic homogenization of the inviscid $G$-equation for incompressible flows.  Commun. Math. Sci., 8(4):1067–1078, 2010.