# Inviscid Burgers as a degenerate elliptic problem

Uditnarayan Kouskiya Amit Acharya Department of Civil & Environmental
Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, email:
<EMAIL_ADDRESS>of Civil & Environmental Engineering, and Center
for Nonlinear Analysis, Carnegie Mellon University, Pittsburgh, PA 15213,
email<EMAIL_ADDRESS>
###### Abstract

We demonstrate the feasibility of a scheme to obtain approximate weak
solutions to the (inviscid) Burgers equation in conservation and Hamilton-
Jacobi form, treated as degenerate elliptic problems. We show different
variants recover non-unique weak solutions as appropriate, and also specific
constructive approaches to recover the corresponding entropy solutions.

##  1 Introduction

In this paper we continue, following [1], the assessment of a recently
introduced duality based approach to solving differential equations involving
evolution in time. A nonlinear ODE system was considered in [1] (along with
some linear PDE); here we consider the (inviscid) Burgers equation in
conservation and Hamilton-Jacobi form.

Our essential idea (see App. A) is to treat the primal PDE under consideration
as constraints and invoke a more-or-less arbitrarily designable strictly
convex, auxiliary potential to be optimized. Then, a dual variational
principle for the Lagrange multiplier (dual) fields can be designed involving
a dual-to-primal (DtP) mapping (i.e. an adapted change of variables), with the
special property that its Euler-Lagrange equations are exactly the primal PDE
system, interpreted as equations for the dual fields using the DtP mapping,
and this, even though the primal system may not have the required symmetries
necessary to be the Euler-Lagrange equations of any objective functional of
the primal fields alone. We use a simple Galerkin Finite element
discretization of the dual Euler-Lagrange system (which is the primal system
of interest using a change of variables) to look for a computational
approximation of weak solutions to (inviscid) Burgers equation. The system is
degenerate-elliptic in space-time domains.

To our knowledge, our approach for generating approximate solutions to a
nonlinear hyperbolic problem (albeit scalar here, but seamlessly generalizable
to systems, formally at least) is new. Brenier [2, Sec. 4] solves the inviscid
Burgers problem without approximation through a maximization of a functional
related to ours, but without the use of ‘base states,’ (see Sec. 2 and App. B)
which we find crucial in making our ideas work; in fact, we have to employ a
sequence of dual functionals parametrized by evolving base states, the latter
self-consistently prescribed by the scheme. As very nicely explained by a
‘safe mountain climbing’ analogy, Brenier’s exact scheme [2, Sec. 4] is
different from our approach, and computational results are shown in [3,
pp.101-105]; he also established connections to an Optimal Transport based
method of attack, which can be utilized for computational approximations via
the celebrated Benamou-Brenier formulation of optimal transport theory, as
mentioned in [2, Sec. 4.1, p. 597].

An outline of the paper is as follows: Sec. 2 and Sec. 3 comprise the
development of the weak formulation for the dual inviscid Burgers equation in
conservation and Hamilton-Jacobi forms, respectively. Sec. 4 describes the
algorithm for the results computed in the paper. In Sec. 5, five selected
problems are solved to illustrate and evaluate the features of the formulation
and algorithm developed in Sections 2 and 3, as they relate to generating
approximate solutions to the (inviscid) Burgers equation. The paper also
contains five appendices supporting various sections of the main narrative. A
word on notation: we always use the Einstein summation convention for indices
unless otherwise mentioned and except when the indices are $x,t$.

##  2 Dual formulation of the Burgers Equation

Burgers equation [4, 5] is a (convection-dominated) convection-diffusion
partial differential equation given by

$\partial_{t}u+u\,\partial_{x}u=\nu\partial_{xx}u\qquad\text{in }\Omega,$ (1)

where $\Omega=(0,L)\times(0,T)\in\mathbb{R}\times\mathbb{R}^{+}$ is a domain
in (1-d)space-time, $u:\Omega\to\mathbb{R}$ and $\nu\in\mathbb{R}^{+}$ are the
flow velocity and viscosity, respectively, in the context of fluid mechanics.

Figure 1: Domain of interest

A schematic of the domain is shown in Fig. 1. Henceforth, $L=1$. The
boundaries $\partial\Omega_{l}$, $\partial\Omega_{r}$, $\partial\Omega_{0}$
and $\partial\Omega_{T}$ will be referred to as the left, right, bottom and
top boundaries, respectively.

We will be mainly interested in the inviscid Burgers equation, obtained by
setting $\nu=0$ in (1):

$\partial_{t}u+u\,\partial_{x}u=0\qquad\text{in }\Omega,$

and, as is standard, explore weak solutions of this quasilinear equation
expressed in conservation form:

$\partial_{t}u+\partial_{x}\left(\frac{1}{2}u^{2}\right)=0\qquad\text{in
}\Omega,$ (2)

with initial and the boundary conditions given by

$u(x,0)=u_{0}(x);\qquad u(0,t)=u_{l}(t).$ (3)

Following the ideas in Appendix A, we form the product of (2) with a dual
field $\lambda:\Omega\to\mathbb{R}$ (that acts like a Lagrange multiplier) and
integrate by parts in the space-time domain to define a pre-dual functional

$\displaystyle\widehat{S}_{H}\,[u,\lambda]$
$\displaystyle=\int_{\Omega}\,\left(-u\,\partial_{t}\lambda-\frac{u^{2}\,\partial_{x}\lambda}{2}+H(u,x,t)\right)\,dt\,dx-\int_{0}^{L}\,u_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,u_{l}(t)\,\lambda(0,t)\,dt$
$\displaystyle=\int_{\Omega}\,\mathcal{L}_{H}(u,\mathcal{D},x,t)\,dt\,dx-\int_{0}^{L}\,u_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,u_{l}(t)\,\lambda(0,t)\,dt,$

where $\mathcal{D}:=(\partial_{t}\lambda,\partial_{x}\lambda)$ and $H(u,x,t)$
is a ‘free’ auxiliary potential. The design of the functional involves
inclusion of only those space-time boundary terms that contain information
available from the primal problem [6, 1]. Additionally, we impose the
following Dirichlet space-time boundary conditions on $\lambda$:

$\displaystyle\lambda(L,t)$ $\displaystyle=\lambda_{r}(t);$ (4a)
$\displaystyle\lambda(x,T)$ $\displaystyle=\lambda_{T}(x),$ (4b)

where $\lambda_{r}(t)$ and $\lambda_{T}(x)$ can be specified arbitrarily. The
choice of the auxiliary potential $H$ needs to be such that it renders the
following equation ‘solvable’ for $u$ (Appendix A):

$\frac{\partial\mathcal{L}_{H}}{\partial
u}:-\partial_{t}\lambda-u\,\partial_{x}\lambda+\frac{\partial}{\partial
u}H(u,x,t)=0.$ (5)

We employ the following shifted quadratic form for $H(u,x,t)$

$H(u,x,t)=\frac{\beta_{u}}{2}\big{(}u-\bar{u}(x,t)\big{)}^{2},$ (6)

where $\beta_{u}\gg 1$. The function $\bar{u}:=\bar{u}(x,t)$, is a ‘free’
function to choose that is referred to as a base state (see [6, Sec. 5]-[1,
Sec. 4]). The choice of this base state is often crucial for the success of
the scheme (see Appendix B for a simple example) and their algorithmic use is
shown in [1, Sec. 5.3] and in Section 5 of this work, among others. Solving
for $u$ in (5), one obtains the following dual-to-primal mapping (DtP):

$\displaystyle u^{(H)}(\mathcal{D},x,t)$
$\displaystyle=\bar{u}+\frac{\bar{u}\,\partial_{x}\lambda+\partial_{t}\lambda}{\beta_{u}-\partial_{x}\lambda}$
(7) $\displaystyle\hat{u}(x,t)$
$\displaystyle:=u^{(H)}\big{(}\mathcal{D}(x,t),x,t\big{)}.$

The dual functional is now defined as

$\displaystyle
S_{H}[\lambda]:=\widehat{S}_{H}\left[u^{(H)}(\mathcal{D},x,t),\lambda\right]$
$\displaystyle=\int_{\Omega}\,\mathcal{L}_{H}\left(u^{(H)}(\mathcal{D},x,t),\mathcal{D},x,t\right)\,dt\,dx$
(8)
$\displaystyle-\int_{0}^{L}\,u_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,u_{l}(t)\,\lambda(0,t)\,dt.$

The above functional can be explicitly written in terms of $\lambda$ as

$S_{H}[\lambda]=\int_{\Omega}\,\mathbb{K}\big{|}_{\mathcal{D}}\,\left(\partial_{t}\lambda+\bar{u}\,\partial_{x}\lambda\right)^{2}\,dt\,dx-\int_{\Omega}\,\left(\bar{u}\,\partial_{t}\lambda-\bar{u}^{2}\,\frac{\partial_{x}\lambda}{2}\right)\,dt\,dx\\\
-\int_{0}^{L}\,u_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,u_{l}(t)\,\lambda(0,t)\,dt.$
(9)

where

$\mathbb{K}=\frac{-1}{2\beta_{u}\left(1-\frac{\partial_{x}\lambda}{\beta_{u}}\right)}.$

Hence, for $\beta_{u}\gg\left|\partial_{x}\lambda\right|$, the nonlinear part
of the bulk integrand is non-positive. The first variation of the above
functional in a direction $\delta\lambda$ yields

$\delta^{(1)}S[\lambda;\delta\lambda]=\bigintsss_{\Omega}\,\left(\frac{\partial\mathcal{L}_{H}}{\partial\hat{u}}\frac{\partial\hat{u}}{\partial\mathcal{D}_{k}}d\mathcal{D}_{k}+\frac{\partial\mathcal{L}_{H}}{\partial\mathcal{D}_{k}}d\mathcal{D}_{k}\right)\,dt\,dx-\int_{0}^{L}\,u_{0}(x)\,\delta\lambda(x,0)\,dx\\\
-\int_{0}^{T}\,u_{l}(t)\,\delta\lambda(0,t)\,dt.$

The first term in the domain integral vanishes due to (5). Since
$\mathcal{L}_{H}$ is necessarily affine in $\mathcal{D}$ one has

$\delta^{(1)}S[\lambda;\delta\lambda]=\int_{\Omega}\,\left(\partial_{t}\hat{u}+\frac{\partial}{\partial
x}\left(\frac{\hat{u}^{2}}{2}\right)\right)\delta\lambda\,dt\,dx-\int_{0}^{L}\,\hat{u}(x,T)\,\delta\lambda(x,T)\,dx-\int_{0}^{T}\,\hat{u}(L,t)\,\delta\lambda(L,t)\,dt\\\
\int_{0}^{L}\,(\hat{u}(x,0)-u_{0}(x))\,\delta\lambda(x,0)\,dx+\int_{0}^{T}\,(\hat{u}(0,t)-u_{l}(t))\,\delta\lambda(0,t)\,dt,$

and using $\delta\lambda$ consistent with (4), the E-L equations and side-
conditions are given by

$\displaystyle\partial_{t}\hat{u}+\frac{\partial}{\partial
x}\left(\frac{\hat{u}^{2}}{2}\right)=0\qquad\text{in }\Omega;$ (10a)
$\displaystyle\hat{u}(x,0)=u_{0}(x);\qquad\hat{u}(0,t)=u_{l}(t).$ (10b)

which are simply (2)-(3) with the replacement $u\to\hat{u}$.

Thus, solutions of the problem defined by (2), (3) can be generated by solving
(10) along with (4), where the DtP mapping (7) is utilized to bridge the
primal variable $u$ with the dual variable $\lambda$. In short, our scheme may
be interpreted as designing an adapated change of variables for solving
Burgers equation.

In [7], the degenerate ellipticity of the dual formulation of a significant
class of equations from continuum mechanics is analyzed. We apply those ideas
next to the dual formulation of Burgers equation given by (10) written as

$\partial_{t}\big{(}\mathcal{F}_{1}(\hat{u})\big{)}+\partial_{x}\big{(}\mathcal{F}_{2}(\hat{u})\big{)}=0,$

some of whose analytical properties are governed by the terms

$A_{ij}=\frac{\partial\,\mathcal{F}_{i}}{\partial(\nabla
D)_{j}}=\left.\frac{\partial\,\mathcal{F}_{i}}{\partial
u}\right|_{\hat{u}}\,\frac{\partial\,u}{\partial(\nabla D)_{j}},$

where $\nabla D:=(\partial_{t}\lambda,\partial_{x}\lambda)$, and

$\begin{gathered}A_{11}=\frac{1}{\beta_{u}-\partial_{x}\lambda};\quad
A_{12}=\frac{\hat{u}}{\beta_{u}-\partial_{x}\lambda};\quad
A_{21}=\frac{\hat{u}}{\beta_{u}-\partial_{x}\lambda};\quad
A_{22}=\frac{\hat{u}^{2}}{\beta_{u}-\partial_{x}\lambda}.\end{gathered}$

Consequently, for $\mathcal{D}=(0,0)$

$c_{i}\,A_{ij}\,c_{j}=\frac{(c_{1}+\bar{u}\,c_{2})^{2}}{\beta_{u}}\geq
0\qquad\forall(c_{1},c_{2})\in\mathbb{R}^{2}$

which establishes that $A_{ij}|_{\mathcal{D}=0}$ is positive semi-definite.
Let $\mathcal{N}$ be neighborhood around $\mathcal{D}=(0,0)$ and let
$\beta_{u}\gg|\partial_{x}\lambda|$. Then, the following expression holds:

$c_{i}\,A_{ij}\,c_{j}=\frac{(c_{1}+\hat{u}\,c_{2})^{2}}{\beta_{u}-\partial_{x}\lambda}=\frac{(c_{1}+\hat{u}\,c_{2})^{2}}{\beta_{u}\left(1-\frac{\partial_{x}\lambda}{\beta_{u}}\right)}\geq
0\qquad\forall(c_{1},c_{2})\in\mathbb{R}^{2}.$

This relation establishes the positive semi-definiteness of $A$ within this
specific neighborhood. Consequently, the equation (10) is locally degenerate
elliptic and for any surface with normal $(n_{1},n_{2})$ in the space-time
domain ($n_{1}$ and $n_{2}$ represents the projection of normal in the time
and space directions, respectively), ellipticity fails when

$(n_{1},n_{2})=\kappa(-\hat{u},1)\quad\text{for }\kappa\in\mathbb{R},$

where $\hat{u}$ is given by (7).

###  2.1 Weak formulation of the dual Burgers equation

A weak formulation of (10) with the side conditions (4) can be generated in
the usual way through integration by parts or by considering the first
variation of the functional (8) or (9) and setting it equal to zero. For any
field $\delta\lambda$ (within an appropriate class) satisfying the conditions
stated below, we intend to find the dual field $\lambda$ which satisfies the
following equations:

$\begin{gathered}R[\lambda;\delta\lambda]:=\int_{\Omega}\,\left(-\hat{u}\,\partial_{t}\delta\lambda-\frac{\hat{u}^{2}\,\partial_{x}\delta\lambda}{2}\right)\,dt\,dx-\int_{0}^{L}\,u_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,u_{l}(t)\,\lambda(0,t)\,dt=0;\\\
\delta\lambda(x,T)=0;\qquad\delta\lambda(L,t)=0;\\\
\lambda(x,T)=\lambda_{T}(T);\qquad\lambda(L,t)=\lambda_{r}(t),\end{gathered}$
(11)

where $\lambda_{T}(\cdot)$ and $\lambda_{r}(\cdot)$ are arbitrarily specified
functions satisfying $\lambda_{T}(L)=\lambda_{r}(T)$. The dual scheme formally
guarantees that the solution to (11) implies the solution to the set of
equations (10). We subsequently make use of the above weak form to numerically
compute an approximate solution for the dual field and utilize the DtP mapping
(7) to obtain the corresponding field for the primal problem i.e. the inviscid
Burgers equation.

##  3 Dual formulation of the (inviscid) Burgers Equation in Hamilton-Jacobi
form

Here we consider the Burgers equation (2) in Hamilton-Jacobi form:

$\partial_{t}Y+\frac{(\partial_{x}Y)^{2}}{2}=\nu\,\partial_{xx}Y\qquad\text{in
}\Omega.$ (12)

Differentiating (12) w.r.t. $x$, and defining $\partial_{x}Y=:u$ gives (1).
Throughout the following text, we will refer to this form of the Burgers
equation as Burgers-HJ and when $\nu=0$, we will refer to it as the inviscid
Burgers-HJ equation. As already stated, our primary interest is in the
inviscid case $\nu=0$. We will consider the $\nu\neq 0$ case to shed light on
the inviscid case, the motivation for which will be provided in the relevant
examples.

It is clear that (12) requires $Y$ to be specified at one point of the domain,
above and beyond information available from the corresponding Burgers
equation. If the function $u$ satisfies Burgers equation and $\partial_{x}Y=u$
then $\partial_{t}Y(x,t)+\frac{1}{2}u^{2}(x,t)-\nu\partial_{x}u(x,t)=f(x*,t)$
for an arbitrarily fixed $x^{*}$, where
$f(x^{*},t)=\partial_{t}Y(x^{*},t)+\frac{1}{2}u^{2}(x^{*},t)-\nu\partial_{x}u(x*,t)$.
In all problems considered in the text and Appendix, we will always work with
$f(x^{*},t)=0$ with knowledge of the Burgers solution of interest known at the
point $x^{*}$ (mostly a boundary point) from which $Y(x^{*},t)$ is determined,
with an initial condition $Y(x^{*},0)$ arbitrarily specified.

The dual formulation of (12) for $\nu=0$ is obtained by writing it in first-
order form [8, Sec. 6.2]:

$\displaystyle\partial_{t}Y=-\frac{u^{2}}{2}\quad\text{in }\Omega;$ (13a)
$\displaystyle\partial_{x}Y=u\qquad\text{in }\Omega.$ (13b)

The initial and the boundary conditions can be given as

$Y(x,0)=Y_{0}(x);\qquad Y(0,t)=Y_{l}(t).$ (14)

Corresponding to the equations (13a) and (13b) we introduce the dual fields
$\lambda$ and $\gamma$, respectively. Similar to the setup explained in Sec.
2, we employ a shifted quadratic form for the auxiliary potential $H$ in both
$Y$ and $u$:

$H(Y,u,x,t)=\frac{\beta_{Y}}{2}(Y-\bar{Y})^{2}+\frac{\beta_{u}}{2}(u-\bar{u})^{2},$
(15)

where $\bar{Y}$ and $\bar{u}$ correspond to the base states for $Y$ and $u$,
respectively.

With this choice of $H$, the DtP mapping for $Y$ and $u$, given by $Y^{(H)}$
and $u^{(H)}$, can be expressed in terms of the dual fields and their
derivatives as

$\displaystyle
Y^{(H)}(\mathcal{D},x,t)=\bar{Y}+\frac{\partial_{t}\lambda+\partial_{x}\gamma}{\beta_{Y}};$
(16a) $\displaystyle
u^{(H)}(\mathcal{D},x,t)=\bar{u}+\frac{\gamma-\lambda\bar{u}}{\beta_{u}+\lambda},$
(16b)

where $\mathcal{D}:=(\lambda,\gamma,\partial_{t}\lambda,\partial_{x}\gamma)$.
Additionally, the following boundary conditions are imposed:

$\displaystyle\gamma(L,t)$ $\displaystyle=\gamma_{r}(t);$ (17a)
$\displaystyle\lambda(x,T)$ $\displaystyle=\lambda_{T}(x),$ (17b)

where $\gamma_{r}(t)$ and $\lambda_{T}(x)$ can be specified arbitrarily. The
dual functional obtained using the above procedure is

$S_{H}[\lambda,\gamma]=\bigintsss_{\Omega}\,\Biggl{(}-\hat{Y}\,\partial_{t}\lambda\,+\frac{\hat{u}\,^{2}\lambda}{2}-\hat{Y}\,\partial_{x}\gamma-\hat{u}\,\gamma+H\left(\hat{Y},\hat{u},x,t\right)\Bigg{)}\,dt\,dx\\\
-\int_{0}^{L}\,Y_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,Y_{l}(t)\,\gamma(0,t)\,dt,$
(18)

where

$\hat{Y}(x,t):=Y^{(H)}\big{(}\mathcal{D}(x,t),x,t\big{)};\qquad\hat{u}(x,t):=u^{(H)}\big{(}\mathcal{D}(x,t),x,t\big{)}.$

The functional above can be written explicitly in terms of the dual variables
as

$S_{H}[\lambda,\gamma]=\int_{\Omega}\,\left(-\frac{\mathbb{K}_{1}}{\beta_{Y}}\,(\partial_{t}\lambda\,+\,\partial_{x}\gamma)^{2}\,-\frac{\mathbb{K}_{2}}{\beta_{u}}(\bar{u}\lambda-\gamma)^{2}\right)\,dt\,dx\\\
+\int_{\Omega}\,\left(\bar{Y}(\partial_{t}\lambda+\partial_{x}\gamma)+\bar{u}\left(\frac{\lambda\bar{u}}{2}-\gamma\right)\right)\,dt\,dx\\\
-\int_{0}^{L}\,Y_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,Y_{l}(t)\,\gamma(0,t)\,dt,$
(19)

where

$\mathbb{K}_{1}=\frac{1}{2};\qquad\mathbb{K}_{2}=\frac{1}{2\left(1+\frac{\lambda}{\beta_{u}}\right)},$

and for $\beta_{u}\gg|\lambda|$, $\mathbb{K}_{2}>0$. Hence, the integrand
containing nonlinear terms in the dual fields is negative semi-definite for
this range of $|\lambda|$. The Euler-Lagrange equations of the dual functional
(18), when extracted using the conditions (17), recover the primal equations
(13) and (14), expressed in terms of their dual counterparts:

$\begin{gathered}\partial_{t}\hat{Y}=-\frac{\bigl{(}\hat{u}\bigl{)}^{2}}{2}\quad\text{in
}\Omega;\\\ \partial_{x}\hat{Y}=\hat{u}\quad\text{in }\Omega;\end{gathered}$
(20a) $\hat{Y}(x,0)=Y_{0}(x);\qquad\hat{Y}(0,t)=Y_{l}(t).$ (20b)

Thus, a solution to the the above set of equations in dual variables implies a
solution the primal problem using the DtP map (16).

To examine the ellipticity of the dual equations (20b), we rewrite them as

$\begin{gathered}\partial_{t}\big{(}\mathcal{F}_{11}(\hat{Y},\hat{u})\big{)}+\partial_{x}\big{(}\mathcal{F}_{12}(\hat{Y},\hat{u})\big{)}+\mathcal{G}_{1}\big{(}\hat{Y},\hat{u}\big{)}=0;\\\
\partial_{t}\big{(}\mathcal{F}_{21}(\hat{Y},\hat{u})\big{)}+\partial_{x}\big{(}\mathcal{F}_{22}(\hat{Y},\hat{u})\big{)}+\mathcal{G}_{2}\big{(}\hat{Y},\hat{u}\big{)}=0,\end{gathered}$

where $\mathcal{F}_{12}=\mathcal{F}_{21}=0$. The behavior of the above set of
equations is now governed by matrix

$\begin{gathered}\mathbb{A}_{ijkl}=\frac{\partial\mathcal{F}_{ij}}{\partial(\nabla
D)_{kl}};\qquad\nabla
D=\begin{bmatrix}\partial_{t}\lambda&\partial_{x}\lambda\\\
\partial_{t}\gamma&\partial_{x}\gamma\\\ \end{bmatrix},\end{gathered}$

with each of the indices $i,j,k,l$ ranging from $1$ to $2$. Also,
$\mathbb{A}_{ijkl}=0$ except

$\mathbb{A}_{1111}=\mathbb{A}_{1122}=\mathbb{A}_{2211}=\mathbb{A}_{2222}=\frac{1}{\beta_{Y}}$

Equivalently, we represent $\mathbb{A}$ as a $2\times 2$ matrix $A$:

$A=\begin{bmatrix}\mathbb{A}_{1111}&\mathbb{A}_{1112}&\mathbb{A}_{1121}&\mathbb{A}_{1122}\\\
\mathbb{A}_{1211}&\mathbb{A}_{1212}&\mathbb{A}_{1221}&\mathbb{A}_{1222}\\\
\mathbb{A}_{2111}&\mathbb{A}_{2112}&\mathbb{A}_{2121}&\mathbb{A}_{2122}\\\
\mathbb{A}_{2211}&\mathbb{A}_{2212}&\mathbb{A}_{2221}&\mathbb{A}_{2222}\\\
\end{bmatrix}$ (21)

which then satisfies

$c_{i}A_{ij}c_{j}=\frac{(c_{1}+c_{4})^{2}}{\beta_{Y}}\geq 0\qquad\forall
c\in\mathbb{R}^{4},$

establishing that $\mathbb{A}_{ijkl}$ is positive semi-definite on the space
of $2\times 2$ matrices and

$C_{ij}\mathbb{A}_{ijkl}C_{kl}=0\quad\mbox{for}\quad\left\\{C\in\mathbb{R}^{2\times
2}\big{|}\ C_{11}+C_{22}=0\right\\}.$

Hence the equation set (20b) is degenerate elliptic. Furthermore, for any
rank-one matrix $C$ of the form $a\otimes n$ satisfying $C:\mathbb{A}C=0$, for
some $a\in\mathbb{R}^{2}$ and some $n\in\mathbb{R}^{2},|n|=1$, the latter
representing a unit normal in the space-time domain, the degenerate
ellipticity (positive semi-definiteness of $\mathbb{A}$) implies a loss of
ellipticity along the direction $n$. Clearly, for any arbitrarily fixed unit
normal $n$, there always exists a vector $a$ such that
$a_{1}n_{1}=-a_{2}n_{2}$. As a result, ellipticity, i.e.
$\det(\mathbb{A}_{ijkl}n_{j}n_{l})>0$ fails along every unit normal in the
domain (which is also apparent from the direct calculation
$\mathbb{A}_{ijkl}n_{j}n_{l}=\frac{1}{\beta_{Y}}n_{i}n_{j}$ which is a rank-
one matrix).

###  3.1 Weak formulation for the dual Hamilton-Jacobi form of inviscid
Burgers

We construct a weak form for the system of equations (20b) with the side
condition (17) in the usual way through integration by parts. As in the case
of Burgers equation, the same weak form can also be obtained from the first
variation of the dual functional (18) or (19). For any field $\delta\lambda$
satisfying the conditions stated below, we intend to find the dual field
$\lambda$ which satisfies the following equations:

$\begin{gathered}\begin{aligned}
R\,[\lambda,\gamma;\delta\lambda,\delta\gamma]:=\int_{\Omega}\,\bigg{(}-\hat{Y}\,\partial_{t}\delta\lambda\,+&\frac{\hat{u}\,^{2}\,\delta\lambda}{2}-\hat{Y}\,\partial_{x}\delta\gamma-\hat{u}\,\delta\gamma\bigg{)}\,dt\,dx\\\
&-\int_{0}^{L}\,Y_{0}(x)\,\big{(}\delta\lambda(x,0)\big{)}\,dx-\int_{0}^{T}\,Y_{l}(t)\,\big{(}\delta\gamma(0,t)\big{)}\,dt=0;\end{aligned}\\\
\delta\lambda(x,T)=0;\qquad\delta\gamma(L,t)=0;\\\
\lambda(x,T)=\lambda_{T}(T);\qquad\gamma(L,t)=\gamma_{r}(t),\end{gathered}$
(22)

where $\lambda_{T}(\cdot)$ and $\gamma_{r}(\cdot)$ are arbitrarily specified
functions.

##  4 Algorithm

As introduced in [1], we compute approximate solutions to the problems of
interest by solving them in a time-concatenated series of space-time
subdomains whose (closed) union forms (the closure of) the entire domain
$\Omega$ of interest. We solve a distinct, dual (space-time) bvp in each of
these stages, step-by-step, marching in time, where the initial condition and
base state in any stage depends on the output produced from the previous
stage, details of which will be explained in the description of example
problems to follow. Each of these subdomains, $\Omega^{(s)}$, is referred to
as a stage, indexed by $s=1,2,3,\ldots,N$,
$\overline{\Omega}=\cup_{s=1}^{n}\overline{\Omega^{(s)}}$.

A two-point Gauss quadrature scheme, in each of the $x$ and $t$ directions,
has been utilized to approximate all the domain and boundary integrals
appearing in the work. We will refer to the collection of nodes corresponding
to a specific time (falling on the nodes of discretization) as a nodal
timeline. Similarly, we will refer to the collection of all Gauss points
corresponding to a particular time as a Gauss timeline. The setup for a sample
mesh of $2\times 2$ has been shown in Fig. 2(a).

Due to the nature of the arbitrarily chosen (without loss of generality) $0$
b.c. utilized on the dual fields ((4) and (17)) at the final-time boundary of
a stage, strong gradients/boundary-layers can arise in the dual solution that
can require excessive refinement to resolve. Since the time-like extent of any
stage is an arbitrary choice subject only to computational expediency and
accuracy, in every stage, solutions obtained for a specified layer of elements
in the time-like direction near the final time boundary of the stage are
ignored. We refer to this operation as ‘truncation’. The following stage
initiates from a suitable time near the end of the current stage on the
retained mesh (either at a Gauss timeline or at a nodal timeline, see Fig.
2(a)), referred to as the ‘cutoff’ line. A visual representation illustrating
this concept is presented in the Fig. 2(c).

We employ the Newton Raphson method (NR) within each stage to approximate the
solutions. A linear span of globally continuous, piecewise smooth finite
element shape functions corresponding to an FE mesh for $\Omega^{(s)}$ is used
to achieve this discretization. These shape functions are represented by
$N^{(\cdot)}$, where $(\cdot)$ denotes the index of any node on the space-time
mesh. Let $D$ and $U$ represent the set of dual and primal variables for the
problem under consideration. So, for the Burgers equation,

$D:=(D_{1})=(\lambda);\quad U:=(U_{1})=(u),$

whereas for the Burgers-HJ,

$D:=(D_{1},D_{2})=(\lambda,\gamma);\quad U:=(U_{1},U_{2})=(Y,u).$

The discretized dual fields and test functions are expressed as

$D_{i}(x,t)=D^{A}_{i}\,N^{A}(x,t);\qquad\delta D_{i}(x,t)=\delta
D^{A}_{i}\,N^{A}(x,t);\qquad dD_{i}(x,t)=dD^{A}_{i}\,N^{A}(x,t),$ (23)

where $D^{A}_{i}$ denotes the finite element nodal degrees of freedom. The
discretized version of the appropriate weak form (11) or (22) generates a
discrete residual $R^{A}_{i}(\mathcal{D})$ (derived in Sec. 5.1 for the
Burgers equation and Sec. 5.2 for the Burgers-HJ) given by

$R[D;\delta D]=\delta D^{A}_{i}\,R^{A}_{i}(\mathcal{D})$

and its variation in the direction $dD$, the Jacobian

$J[D;\delta D,dD]=\delta D^{A}_{i}\,J^{AB}_{ij}(\mathcal{D})\,dD^{B}_{j}.$

Any stage begins with an initial guess of $D_{i}^{A(0)}=0$ followed by solving
the following matrix equation for the $k^{th}$ correction:

$\displaystyle-R^{A}_{i}\left(\mathcal{D}^{(k-1)}\right)$
$\displaystyle=J^{AB}_{ij}\left(\mathcal{D}^{(k-1)}\right)\,dD^{B}_{j}$ (24)
$\displaystyle D^{B(k)}_{j}$ $\displaystyle=D^{B(k-1)}_{j}+dD^{B}_{j},$

where $\mathcal{D}=(\partial_{t}\lambda,\partial_{x}\lambda)$ for the Burgers
equation and
$\mathcal{D}=(\lambda,\gamma,\partial_{t}\lambda,\partial_{x}\gamma)$ for
Burgers-HJ. The convergence criteria for NR is set as follows:

$\underset{A,i}{\max}\,\left|R^{A}_{i}\right|<tol,$

where $tol$ represents a user-defined threshold tolerance.

The following notation is used in the algorithm:

$(\cdot)^{(s)(k)}$ | value of $(\cdot)$ at $k^{th}$ NR iterate for stage $s$
---|---
$N_{x}$ | number of elements in $x$ direction
$N_{t}^{(s)}$ | number of elements in $t$ direction per stage
$t_{i}^{(s)}$ | time at the start of stage $s$
$T^{(s)}$ | time over which the stage $s$ is solved
$t_{f}^{(s)}$ | time at the cutoff line for stage $s$
$\Omega^{(s)}$ |  $(0,L)\times\left(t_{i}^{(s)},t_{i}^{(s)}+T^{(s)}\right)$; domain for stage $s$
$\bar{U}_{k}^{(s)}$ | base state for $U_{k}$ in stage $s$
$\beta_{k}$ | coefficient(s) defining $H$ in (6) or (15)
$U_{0}^{(s)}$ | primal initial condition for stage $s$
$\Omega^{(s)\times}$ | discarded domain for the stage $s$
$N_{c}\times N_{x}$ | number of discarded elements per stage

Algorithm

Initialization:

  1. 1.

Set $s=1$, $t_{f}^{(0)}=0$ and $U_{0}^{(1)}(x)$ from (3) or (14) depending on
the problem under consideration.

  2. 2.

Choose the values for $T^{(s)}$, $N_{x}$, $N_{t}^{(s)}$, $\beta_{k}$, $N_{c}$
and $tol$.

$\boldsymbol{s_{th}}$ stage:

  1. 1.

Set $t_{i}^{(s)}=t_{f}^{(s-1)}$ and generate $U_{0}^{(s)}$ ($s>1$, see Sec. 5
for details). Over the domain $\Omega^{(s)}$, set appropriate
$\bar{U}_{k}^{(s)}$.

  2. 2.

Initiate NR: Set $D_{k}^{A(s)(0)}=0$.

  3. For $j\geq 0$:

    1. i

Evaluate $R_{i}^{A(s)(j)}$.

    2. ii

Set $d^{(s)(j)}=\underset{A,i}{\max}\,\left|R^{A(s)(j)}_{i}\right|.$

    3. if $d^{(s)(j)}<tol$ then go to step 3 (exit loop).

    4. else continue.

    5. iii

Evaluate $J_{ij}^{AB(s)(j)}$.

    6. iv

Evaluate $D_{k}^{A(s)(j+1)}$ using (24).

    7. v

do  $j=j+1$ and go to step 2i

  4. 3.

$D_{k}^{(s)(j)}$ serves as the dual solution for the current stage, while the
corresponding $U_{k}$ serves as the solution for the primal problem. Discard
the results obtained in the region $\Omega^{(s)\times}$.

  5. 4.

Set an appropriate $t_{f}^{(s)}$ near the final time obtained on
$\Omega^{(s)}\backslash\Omega^{(s)\times}$.

  6. 5.

Set $s=s+1$ and repeat steps 1-4 until $t_{f}^{(s)}\geq T$

Table 1: Algorithm to solve Burgers equation or Burgers-HJ. The subscript
index $k$ takes the value of 1 for Burgers and ranges over $(1,2)$ for
Burgers-HJ.

(a) FEM element geometry

(b) Time marching

(c) $\Omega^{(s)}$

Figure 2: (a) A sample $2\times 2$ mesh with element numbers at center. Red
dashed lines, passing through the nodal points (red dots), represent nodal
timelines. $g_{1}$ and $g_{2}$ represent the first and the second Gauss
timelines passing through the Gauss points of the elements $1$ and $2$. (b)
represents the concatenated domains after truncation. Additionally,
$t_{f}^{(s)}=t_{i}^{(s+1)}$.

The algorithm has been presented in table 1. Additional details for the
algorithm are as follows:

  * •

Corresponding to the discretized dual fields $D$, the primal field $U$ can be
evaluated using the appropriate DtP mapping (7) or (16) on the Gauss points
within any element, which is further used in the evaluation of $R_{i}^{A}$ and
$J_{ij}^{AB}$.

  * •

For each stage, we discard the results over the final $T_{c}$ duration of time
in that stage (see Fig. 2(c)). Denoting by $T_{e}$ the edge-length of elements
in the time-like direction (see Fig. 2(a)), the truncation operation is
carried out over $N_{c}:=\frac{T_{c}}{T_{e}}$ space-like layers of elements.

  * •

For the Burgers equation, $U_{0}(x):=u_{0}(x)$, whereas for the Burgers-HJ,
$U_{0}(x):=Y_{0}(x)$. For the first stage, initial condition is used from the
problem definition (3)or (14), and for all subsequent stages, the initial
condition is generated from the output of the previous stage (see Sec. 4.1 and
Sec. 4.2). The primal boundary conditions obviously need no adjustment in the
course of dual time evolution. These initial and boundary conditions enter
into the dual formulation as natural boundary conditions in each stage.

  * •

To evaluate the quality of the results produced via the algorithm, the
corresponding exact entropy solution to each example is presented in Appendix
E.

The figures presented in this work are produced as follows: For each Gauss
timeline, we compute the average of the primal data obtained from the two
Gauss points within each element through which the Gauss timeline traverses.
This average is then assigned to the $x$-coordinate of the element center and
$t$-coordinate of the Gauss timeline under consideration.

We demonstrate and discuss the results of the computation of five selected
examples that are used to evaluate and understand our dual scheme applied to
the Burgers equation and its corresponding Hamilton-Jacobi form. Sec. 5.1 and
Sec. 5.2 deal with the inviscid case of the Burgers equation and Burgers-HJ,
respectively, and we often omit the adjective ‘inviscid.’ Specific details
related to each of these equations are presented below.

###  4.1 Algorithmic details for the Burgers equation

In order to evaluate solutions to the dual form of the Burgers equation (10),
we utilize the residual (11) and generate the corresponding Jacobian by
considering its variation in a direction $d\lambda$:

$J[\delta\lambda;d\lambda]=\bigintsss_{\Omega}\,\left(-\partial_{t}\delta\lambda\,\left(\frac{\partial\hat{u}}{\partial\mathcal{D}_{i}}d\mathcal{D}_{i}\right)-\hat{u}\,\partial_{x}\delta\lambda\,\left(\frac{\partial\hat{u}}{\partial\mathcal{D}_{i}}d\mathcal{D}_{i}\right)\right)\,dt\,dx,$

where $\mathcal{D}:=(\partial_{t}\lambda,\partial_{x}\lambda)$.

We employ the approximate dual field (23) in the residual (11) to generate a
discrete residual $R^{A}$ at each node $A$ given by

$R^{A}=\bigintsss_{\Omega}\,\left(-\hat{u}\,\partial_{t}N^{A}-\frac{\hat{u}^{2}\,\partial_{x}N^{A}}{2}\right)\,dt\,dx-\int_{0}^{L}u_{0}(x)\,N^{A}(x,0)\,dx-\int_{0}^{T}\,u_{l}(t)\,N^{A}(0,t)\,dt.$

However, $R^{A}=0$ is not imposed for all the nodes corresponding to the right
and top boundaries as a consequence of the Dirichlet b.cs (4). The discrete
Jacobian matrix associated with the pair of degrees of freedom $(A,B)$ (the
coefficient of the term $(\delta\lambda^{A}d\lambda^{B})$) is given by

$\begin{gathered}J^{AB}=\int_{\Omega}\,\Big{(}-\partial_{t}N^{A}-\hat{u}\,\partial_{x}N^{A}\Big{)}\,\left(\frac{\partial\hat{u}}{\partial\mathcal{D}_{i}}\partial_{i}N^{B}\right)\,dt\,dx;\\\
\frac{\partial\hat{u}}{\partial(\partial_{t}\lambda)}=\frac{1}{\beta_{u}-\partial_{x}\lambda};\qquad\frac{\partial\hat{u}}{\partial(\partial_{x}\lambda)}=\frac{\hat{u}}{\beta_{u}-\partial_{x}\lambda},\end{gathered}$

where $\partial_{1}(\cdot)=\partial_{t}(\cdot)$ and
$\partial_{2}(\cdot)=\partial_{x}(\cdot)$.

Based on the algorithm presented in Table. 1, after truncation, the cutoff
line for each stage (from which the next stage commences) is the $g_{2}$ Gauss
timeline of the last layer of retained elements with $t_{f}^{(s)}$ the
corresponding physical time. The data for $\hat{u}$ obtained at this time
(required only at Gauss points) is denoted by
$\hat{u}^{(s)}\left(x,t_{f}^{(s)}\right)$. The initial condition and the base
state for the next stage is set as

$\displaystyle u_{0}^{(s+1)}(x)$
$\displaystyle=\hat{u}^{(s)}\left(x,t_{f}^{(s)}\right);$
$\displaystyle\bar{u}^{(s+1)}(x,t)$
$\displaystyle=\mathcal{S}\left[\hat{u}^{(s)}\left(x,t_{f}^{(s)}\right)\right],$

where $\mathcal{S}$ denotes a smoothing operator, and the specific details can
be found in Appendix D. In essence, the primal fields generated in each stage,
relying on the dual fields and their derivatives, may exhibit oscillations,
particularly in the vicinity of any shock; the smoothing operation mitigates
these high-wave number oscillations in the base state for the subsequent stage
that remains constant throughout the subsequent cycle. We note that such
smoothing is not used to define the initial condition for the next stage.

###  4.2 Algorithmic details for the Burgers-HJ

To assess solutions to (20b), we employ the residual (22) and generate a
corresponding Jacobian by considering its variation in the direction
$dD:=(d\lambda,d\gamma)$, expressed as

$J[\delta\lambda,\delta\gamma;d\lambda,d\gamma]=\bigintsss_{\Omega}\,\left(\bigl{(}-\partial_{t}(\delta\lambda)-\partial_{x}(\delta\gamma)\bigl{)}\,\left(\frac{\partial\hat{Y}}{\partial\mathcal{D}_{i}}d\mathcal{D}_{i}\right)+\bigl{(}\hat{u}\,\delta\lambda-\delta\gamma\bigl{)}\,\left(\frac{\partial\hat{u}}{\partial\mathcal{D}_{i}}d\mathcal{D}_{i}\right)\right)\,dt\,dx,$
(25)

where $\mathcal{D}:=(\partial_{t}\lambda,\partial_{x}\gamma,\lambda,\gamma)$.
The discrete version for the residual (22), computed using the approximate
dual fields $D$ as defined in (23), for each node $A$ is given by

$R_{1}^{A}=\bigintsss_{\Omega}\,\Biggl{(}-\hat{Y}\,\partial_{t}N^{A}\,+\frac{\hat{u}\,^{2}N^{A}}{2}\Biggl{)}\,dt\,dx-\int_{0}^{L}\,Y_{0}(x)\,(N^{A}(x,0))\,dx;$
$R_{2}^{A}=\int_{\Omega}\,\bigl{(}-\hat{Y}\,\partial_{x}N^{A}+\hat{u}\,N^{A}\bigl{)}\,dt\,dx-\int_{0}^{T}\,Y_{l}(t)\,(N^{A}(0,t))\,dt.$

$R_{1}^{A}=0$ for all the nodes corresponding to the top boundary and
$R_{2}^{A}=0$ for all the nodes corresponding to the right boundary are not
imposed as a consequence of the Dirichlet b.cs (17). The discrete version of
the Jacobian corresponding to the degree of freedom pair $(A,i),(B,j)$ is
given as follows:

$\begin{gathered}J_{11}^{AB}=\bigintsss_{\Omega}\,\left(-\frac{\partial\hat{Y}}{\partial(\partial_{t}\lambda)}\,\partial_{t}N^{A}\,\partial_{t}N^{B}+\hat{u}\,\frac{\partial\hat{u}}{\partial\lambda}\,N^{A}\,N^{B}\right)\,dt\,dx;\\\
J_{12}^{AB}=\bigintsss_{\Omega}\,\left(-\frac{\partial\hat{Y}}{\partial(\partial_{x}\gamma)}\,\partial_{t}N^{A}\,\partial_{x}N^{B}+\hat{u}\,\frac{\partial\hat{u}}{\partial\gamma}\,N^{A}\,N^{B}\right)\,dt\,dx;\\\
J_{21}^{AB}=\bigintsss_{\Omega}\,\left(-\frac{\partial\hat{Y}}{\partial(\partial_{t}\lambda)}\,\partial_{x}N^{A}\,\partial_{t}N^{B}-\frac{\partial\hat{u}}{\partial\lambda}\,N^{A}\,N^{B}\right)\,dt\,dx;\\\
J_{22}^{AB}=\bigintsss_{\Omega}\,\left(-\frac{\partial\hat{Y}}{\partial(\partial_{x}\gamma)}\,\partial_{x}N^{A}\,\partial_{x}N^{B}-\frac{\partial\hat{u}}{\partial\gamma}\,N^{A}\,N^{B}\right)\,dt\,dx;\\\
\frac{\partial\hat{Y}}{\partial(\partial_{t}\lambda)}=\frac{1}{\beta_{Y}};\quad\frac{\partial\hat{Y}}{\partial(\partial_{x}\gamma)}=\frac{1}{\beta_{Y}};\\\
\frac{\partial\hat{u}}{\partial\lambda}=-\frac{\hat{u}}{\beta_{u}+\lambda};\quad\frac{\partial\hat{u}}{\partial\gamma}=\frac{1}{\beta_{u}+\lambda}.\end{gathered}$

In contrast to the examples associated with the Burgers equation, the
oscillations observed in the solution (for $\hat{Y}$ here) obtained on the
Gauss points are significantly more pronounced. To address this issue, in each
stage, the truncation is followed by an $L^{2}$ projection of $\hat{Y}$
(Appendix. C) along the $n_{2}$ nodal timeline of the last layer of elements
obtained on the retained mesh. Such an operation provides us with a continuous
$C^{0}$ approximant along the considered nodal timeline. This timeline also
serves as the cutoff line for the current stage. The initial condition and the
base states for the following stage are set as follows:

$\begin{gathered}Y_{0}^{(s+1)}(x)=L^{2}\left[Y^{(s)}\left(x,t_{f}^{(s)}\right)\right];\\\
\bar{Y}^{(s+1)}(x,t)=Y^{(s)}\left(x,t_{f}^{(s)}\right);\quad\bar{u}^{(s+1)}=\frac{d}{dx}Y^{(s)}\left(x,t_{f}^{(s)}\right),\end{gathered}$

where $L^{2}[\cdot]$ represents the $L^{2}$ projection operator. Unlike the
scheme for the Burgers equation, no smoothing for base states is used for the
inviscid Burgers-HJ algorithm.

##  5 Results

To evaluate the dual scheme, we demonstrate five examples for the two types of
equations, each presenting distinct levels of complexity. For each of the
examples introduced in Sec. 5.1, a related example is generated in Sec. 5.2.
This associated example is supplied with the initial and boundary conditions
in a manner that the exact entropy solution produced by the Burgers-HJ
equation should have its derivative equal to the exact entropy solution
obtained for the corresponding problem from the Burgers equation.

###  5.1 Burgers Equation

For each of the examples presented in this section, the following parameters
were set for each stage in the algorithm: $T^{(s)}=5\times 10^{-3}$,
$N_{x}=100$, $N_{t}^{(s)}=100$, $\beta_{k}=10^{6}$, $N_{c}=5$ and
$tol=10^{-16}$. This amounts to a fairly fine discretization of the problem,
but our main goal here is not to probe the computational efficiency of the
proposed elementary scheme, but to simply use it to evaluate and understand
the capabilities of our dual variational formulation for nonlinear PDE in this
specific context.

####  5.1.1 Expansion Fan

We consider the Riemann problem corresponding to an expansion fan solution
with initial condition given by:

$u_{0}(x)=\begin{cases}U_{L}=0&\mbox{for }x<0.5\\\ U_{R}=1&\mbox{for
}x>0.5.\end{cases}$ (26)

The problem has non-unique weak solutions (consisting of the continuous
rarefaction wave and an entire family of discontinuous solutions with shocks),
and a unique ‘entropy’ solution is obtained by imposing an entropy condition
[9, 10] (also see App. E). This entropy solution is a rarefaction wave, which
for $t>0$ in an infinite domain is given by (38).

We use the primal boundary condition $u_{l}(t)=0$ and apply the Dirichlet b.cs
$\lambda_{T}(x)=\lambda_{r}(t)=0$. The results for this setup are shown in
Fig. 3.

Comparing Fig. 3(a) against Fig. 3(b), it is surprising to see that the dual
scheme automatically picks up only the entropy solution without enforcing any
further conditions, and seems to be incapable of recovering the other possible
weak solutions for the equation.

(a) $u(x,t)$

(b) Exact $u(x,t)$ based on (38)

(c) Line plots for $u(x,t)$

Figure 3: DtP mapping generated primal field $u$ for the expansion fan. Fig.
(c) shows minor overshoots as the fan opens up which may be attributed to the
$C^{0}$ approximation of the dual fields.

####  5.1.2 Shock

We next consider the Riemann problem (26) with

$U_{L}=1;\qquad U_{R}=0.$

The exact entropy solution to this problem on an infinite domain is given by
(40).

We use the primal boundary condition $u_{l}(t)=1$ and apply the Dirichlet b.cs
$\lambda_{T}(x)=\lambda_{r}(t)=0$. The results for this setup has been shown
in Fig. 4. The over(under)shoot seen around the shock profile at any time in
Fig. 4(b) arise due to the $C^{0}$ FE interpolation of the dual fields.
Nevertheless, the shock profile effectively captures the accurate height and
speed, aligning closely with the exact entropy solution.

(a) $u(x,t)$

(b) Line plots for $u(x,t)$

Figure 4: DtP mapping generated primal field $u$ for the shock problem . In
Figure (a), the black asterisks represent the shock trajectory based on the
exact entropy solution, which is being superimposed on the original plot.

####  5.1.3 Double Shock

Here we assess the dual formulation in the setting of two traveling shocks
interacting with each other. Consider the following initial condition

$u_{0}(x)=\begin{cases}\begin{aligned} 1&\quad\text{for }x<0.25\\\
0.5&\quad\text{for }0.25<x<0.5\\\ 0&\quad\text{for
}x>0.5.\end{aligned}\end{cases}$

The solution to this problem on an infinite domain is given by (42).

We use the primal boundary condition $u_{l}(t)=1$ and apply the Dirichlet b.cs
$\lambda_{T}(x)=\lambda_{r}(t)=0$. The results for this setup are shown in
Fig. 5.

Evident form Fig. 5(b), the instabilities near the two shocks, caused by the
continuous interpolation, do not accumulate as the shocks merge. Instead, they
vanish, leaving no indication of a double shock presence. Based on Fig. 5(a),
it is also evident that the dual problem correctly captures the accurate shock
speeds for each individual shock.

(a) $u(x,t)$

(b) Line plots for $u(x,t)$

Figure 5: DtP mapping generated primal field $u$ for a double shock. In Figure
(a), the black asterisks represent the shock trajectories for the two shocks
based on the exact entropy solution, which is being superimposed on the
original plot. The two shocks merge at $(x,t)=(0.625,0.5)$

####  5.1.4 Half N-wave

(a) $u(x,t)$

(b) Line plots for $u(x,t)$

Figure 6: DtP mapping generated primal field $u$ for the half N-wave problem.
In Fig.(a), the black asterisks represent the shock trajectory based on the
exact entropy solution, which is being superimposed on the original plot. In
Figure (b), the adjacent black lines, aligned with the original plot lines at
various time points, depict the position and height of the shock based on the
exact entropy solution for those respective times.

The objective of this problem is to evaluate the dual formulation when a shock
wave decreases in size and experiences a change in speed during its
propagation Fig. E.4. We call this problem as a half N-wave, adapting standard
terminology in the literature (even though it looks like a ‘reflected’ half
$N$). Consider the following initial condition:

$\displaystyle u_{0}(x)=\begin{cases}0&\text{for }x<x_{0}\\\\[2.0pt]
\dfrac{h_{0}}{l_{0}}(x-x_{0})&\text{for }x_{0}\leq x<x_{0}+l_{0}\\\\[4.0pt]
0&\text{for }x>x_{0}+l_{0},\end{cases}$

where $x_{0}=0.25$, $l_{0}=0.25$ and $h_{0}=2$. The solution to this problem
on an infinite domain is given by (44).

We use the primal boundary condition $u_{l}(t)=0$ and apply the Dirichlet b.cs
$\lambda_{T}(x)=\lambda_{r}(t)=0$. The results for this setup are shown in
Fig. 6.

The geometry in this example is more intricate compared to previous examples,
featuring both a fan at $x=0.25$ and a traveling shock. As explained in
Appendix. E.4, conservation of the area under the right-angled triangle
dictates the speed and the height of the shock, which vary non-linearly in
time. The dual formulation successfully captures this phenomenon.

####  5.1.5 N-wave

The objective of this problem is to evaluate the dual formulation in the
context of two expansion fans in opposite directions converging to form a
standing shock whose magnitude dissipates over time (Fig. 19(a)). Consider the
following initial condition:

$\displaystyle u_{0}(x)=\begin{cases}0&\text{for }x<0.25\\\
-4h_{0}(x-0.5)&\text{for }0.25\leq x<0.75\\\ 0&\text{for }x>0.75;\end{cases}$
$h_{0}=2.$

The solution to this problem on an infinite domain is given by (48).

We use the primal boundary condition $u_{l}(t)=0$ and apply the Dirichlet b.cs
$\lambda_{T}(x)=\lambda_{r}(t)=0$.

Based on the discussion found in E.5 related to Fig. 19(b), each of the two
tips at $x=0.25$ and $x=0.75$ travel in opposite directions while maintaining
a constant speed and height. The information originating from the range
$x\in(0.25,0.75)$ at time $t=0$ ultimately converges to form a standing shock
at $x=0.5$ and $t_{m}=0.125$. Subsequently, the standing shock gradually self-
annihilates. The results based on the dual scheme for this setup, shown in
Fig. 7, manages to capture this whole sequence of events quite well.

(a) $u(x,t)$

(b) Exact $u(x,t)$ based on (48)

(c) Line plots for $u(x,t)$

Figure 7: DtP mapping generated primal field $u$ for the N-wave problem. In
Fig. 7(c), the dashed lines use colors matching those of the line plots at
specific times to represent the true height at those times according to the
exact entropy solution. The black dashed line corresponds to the height of the
tip before the shock formation occurs.

###  5.2 Inviscid Burgers-HJ

For the examples presented in this section, the following parameters were set
for each stage in the algorithm: $\beta_{Y}=\beta_{u}=10^{6}$, $N_{c}=5$ and
$tol=10^{-16}$. The parameters defining the mesh were
$N_{x}=50,N_{t}^{(s)}=10$ with $T^{(s)}=5\times 10^{-3}$. This is again a fine
discretization, and our concern here is simply to evaluate the dual
formulation of the problem.

####  5.2.1 Expansion Fan

Consider the following initial condition:

$Y_{0}(x)=\begin{cases}0&\mbox{for }x<0.5\\\ x&\mbox{for }x>0.5.\end{cases}$

The exact entropy solution for this setup on an infinite domain is given by
(39). To evaluate the dual solutions numerically on $\Omega$, we impose
$Y_{l}(t)=0$, naturally. The results obtained for this setup are show in Fig.
8, which shows an unexpected dip in the $Y$ profiles at $x=0.5$.

(a) Exact time plots for $Y(x,t)$

(b) Line plots for $Y(x,t)$ obtained using the dual scheme

Figure 8: The two plots have been drawn on same scale. Fig.(b) shows DtP
mapping generated primal field $Y$ for the Fan problem. The dots in various
colors at $x=1$ correspond to the exact value of $Y$ derived from the exact
entropy solution (39), corresponding to the time indicated by the color of
each respective dot.

This discrepancy indicates that the dual scheme does not capture the entropy
solution for this particular setup. Upon closer inspection one can notice that
the dip actually represents a standing shock which can be confirmed by
approximating the spatial gradient ($u$) to the immediate left and to the
right of $x=0.5$, and using the Rankine-Hugoniot jump condition such that

$c\,(x=0.5,t)=\frac{\partial_{x}Y(x=0.5^{+},t)+\partial_{x}Y(x=0.5^{-},t)}{2}\approx\frac{1+(-1)}{2}=0.$

Thus the evolution profiles indicates that a fan breaks up into a locally
standing shock which then propagates information in both directions. We note
that within an expansion fan, a shock is acceptable as a weak solution.

####  5.2.2 Shock

Consider the initial condition and the boundary condition given by

$\begin{gathered}Y_{0}(x)=\begin{cases}x&\mbox{for }x<0.5\\\ 0.5&\mbox{for
}x>0.5;\end{cases}\\\ Y_{l}(t)=-\frac{t}{2},\end{gathered}$

respectively. For these conditions, we expect to recover the solution given by
(41), the derivative of which recovers the shock presented in 5.1.2. The $Y$
data obtained at different times from the simulation of the dual problem has
been shown in Fig. 9.

Referring to Fig. 9(b), it is apparent that kink in the $Y$ profiles moves at
an approximate speed of 0.5. This kink signifies the presence of a shock in
the $u$ profile, where the slope of the $Y$ (equivalent to $u$) profile
changes from 1 to 0. As depicted in Fig. 9(b), minor oscillation-type
instabilities are still noticeable in the $Y$ profiles as time advances.

(a) $Y(x,t)$

(b) Line plots for $Y(x,t)$

Figure 9: DtP mapping generated primal field $Y$ for the Shock problem. In
Figure(a), the black dots represent the trajectory of the kink based on the
exact entropy solution (41). In Figure (b), the dots plotted at $Y=0.5$ with
colors matching those of the line plots represent the position of the kink
based on the exact $Y$ profile (41) at that specific time.

####  5.2.3 Double Shock

Consider the initial condition and the boundary condition given by

$\begin{gathered}Y_{0}(x)=\begin{aligned} \begin{cases}x&\quad\text{for
}x<0.25\\\ x/2+0.125&\quad\text{for }0.25<x<0.5\\\ 0.375&\quad\text{for
}x>0.5;\end{cases}\end{aligned}\\\ Y_{l}(t)=-\frac{t}{2},\end{gathered}$

respectively. For these conditions, we expect to recover the solution given by
(43), the derivative of which recovers the double shock presented in 5.1.3.
The results for this setup are shown in Fig. 10.

Similar to the shock example, this problem has two kinks emnating at $t=0$
corresponding to the two shocks, traveling at different speeds, merging at
$t=0.5$. Apparent from Fig. 10(a), the characteristics originating between
$x=0.25$ and $x=0.5$ move at a slower speed compared to the approaching kink
from the left, eventually merging with it. The dual formulation captures this
phenomenon, although (very minor) instabilities are again visible as evident
in Fig. 10(b)

(a) $Y(x,t)$

(b) Line plots for $Y(x,t)$

Figure 10: DtP mapping generated primal field $Y$ for the double shock. In
Fig.(a), the black dots represent the trajectory of the two kinks in the
space-time domain based on the exact expressions (43). The solid black lines
represent contour lines corresponding to the characteristics of $u$ derived
through dual scheme in the vicinity of the double shock. In Figure (b), the
dots plotted with colors matching those of the line plots represent the
position and height of the kinks based on the exact $Y$ profile (43) at that
specific time.

####  5.2.4 Half N-wave

Consider the following equation set:

$\displaystyle Y_{0}(x)=\begin{cases}0&\text{for }x<x_{0}\\\\[4.0pt]
\dfrac{h_{0}}{l_{0}}\left(\dfrac{x^{2}+x_{0}^{2}}{2}-x_{0}\,x\right)&\text{for
}x_{0}\leq x<x_{0}+l(t)\\\\[10.0pt] \dfrac{h_{0}\,l_{0}}{2}&\text{for
}x>x_{0}+l(t).\end{cases}$

The exact expression for $Y$ for this setup is given by (45), the spatial
derivative of which represents the Half N-wave presented in Sec. 5.1.4. The
results obtained via our dual scheme is shown in Fig. 11.

(a) Exact time plots for $Y(x,t)$

(b) Line plots for $Y(x,t)$ obtained using the dual scheme

Figure 11: The two plots have been drawn on same scale. Fig.(b) represents the
DtP mapping generated $Y$ profiles. The black circles along the kinks
represent the trajectory of one of the multiple kinks whose spatial
derivatives exhibit a shock-like behavior.

Comparing against the profiles for $Y$ based on the exact entropy solution
illustrated in 11(a), we immediately notice the standing shock behavior (as in
the case of the fan example) here at $x=0.25$. Additionally, we notice several
small shock-like structures nucleating one after other at different points in
space as the curved section of the profile between $x=0.25$ and $x=0.5$
advances. A ‘shock’ speed analysis similar to the fan example, but now on one
such structure characterized by a kink for which the spatial derivative ($u$)
approximately goes from 1 to 0 (shown by the circles in Fig. 11(b)), reveals
that the distance travelled by this kink in a given time approximates well the
distance suggested by the ‘shock speed’ of the kink within that time,
establishing that these structures indeed behave like shocks.

####  5.2.5 N-wave

As the last example, the following initial and boundary conditions are
considered:

$\begin{gathered}\begin{aligned} Y_{0}(x)=\begin{cases}0&\text{for }x<0.25\\\
-2h_{0}(x^{2}-x)-3h_{0}/8&\text{for }0.25\leq x<0.75\\\ 0&\text{for
}x>0.75;\end{cases}\end{aligned}\\\ Y_{l}(t)=0;\\\ h_{0}=2.\end{gathered}$

The exact entropy solution for this setup is given by (49), which gets
approximated well based on the profiles displayed in Fig. 12, although we
again notice the standing shock formations at $x=0.25$ and $x=0.75$ after a
certain time based on Fig. 12(b).

(a) Exact time plots for $Y(x,t)$

(b) Line plots for $Y(x,t)$ obtained using the dual scheme

Figure 12: The two plots have been drawn on same scale. Fig.(b) represents the
DtP mapping generated $Y$ profiles. The dots in various colors at $x=0.5$
correspond to the exact value of $Y$ derived from (49), corresponding to the
time indicated by the color of each respective dot, after the kink formation
takes place at $t=0.125$.

###  5.3 Viscous Burgers-HJ

Evident from the examples presented in the last section, the dual scheme does
not recover the entropy solutions in the case of the Fan 5.2.1, the Half
N-wave 5.2.4 and the N-wave 5.2.5. As mentioned previously, for any Riemann
problem of the form (26) with $U_{R}>U_{L}$, the equation has non-unique weak
solutions and the entropy condition, i.e. information along characteristics
travel towards the shock-curve in space-time and not away from it, acts to
select a solution, in this case a fan. This ambiguity disappears when
viscosity is introduced; indeed, entropy solutions for inviscid Burgers(-HJ)
are known to be limits of Burgers solutions as $\nu\to 0$. While our scheme
recovers entropy solutions when dealing with the Burgers equation examples in
Sec. 5.1, this does not happen with the computational scheme for Burgers-HJ,
and it appears to recover a greater variety of weak solutions (which, in a
sense, is what the scheme, in and of itself, is required to do without any
further conditions). The examples presented in the last section demonstrate
this point.

Considering the functional (19), it is evident that the derivative
$\partial_{x}\lambda$ is not being penalized. This allows for spatial
discontinuities to exist in $\lambda$ (that may arise from even inevitable
numerical round-off errors in $\lambda$), and correspondingly in $\hat{u}$ due
to the DtP mapping (16b). These discontinuous will be propagated in the
spatial domain due to the inherent characteristics of the inviscid Burgers
equation. To test this hypothesis, we work with a dual formulation of Burgers-
HJ (12) with $\nu\neq 0$.

The first order system associated with (12) can be written as

$\displaystyle\partial_{t}Y=-\frac{u^{2}}{2}+\nu\partial_{x}u\quad\text{in
}\Omega;$ (27a) $\displaystyle\partial_{x}Y=u\qquad\text{in }\Omega.$ (27b)

The formulation for this set of equations is analogous to the formulation
presented in Sec. 3. Consider the following pre-dual functional given by

$\hat{S}_{H}[\,Y,u,\lambda,\gamma\,]=\int_{\Omega}\,\mathcal{L}(Y,u,\mathcal{D},x,t)\,dt\,dx-\left.\int_{0}^{T}\,f_{i}(t)\,\lambda(x,t)\,dt\right|_{x=0}^{x=L}\\\
-\int_{0}^{L}\,Y_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,Y_{l}(t)\,\gamma(0,t)\,dt;$
(28)
$\mathcal{L}(Y,u,\mathcal{D},x,t)=-Y\,\partial_{t}\lambda\,+\frac{u\,^{2}\lambda}{2}-Y\,\partial_{x}\gamma+\nu\,u\,\partial_{x}\lambda-u\,\gamma+\frac{\beta_{Y}}{2}(Y-\bar{Y})^{2}+\frac{\beta_{u}}{2}(u-\bar{u})^{2},$

where $f_{i}$, $i=1,2$, are prescribed functions on the left and right
boundaries defining natural boundary conditions on those boundaries when
$\lambda$ is allowed to freely vary there. Employing the following conditions
yields the DtP map

$\displaystyle\frac{\partial\mathcal{L}}{\partial Y}=0:\quad
Y=Y^{(H)}(\mathcal{D},x,t)=\bar{Y}+\frac{\partial_{t}\lambda+\partial_{x}\gamma}{\beta_{Y}};$
(29a) $\displaystyle\frac{\partial\mathcal{L}}{\partial u}=0:\quad
u=u^{(H)}(\mathcal{D},x,t)=\bar{u}+\frac{\gamma-\lambda\bar{u}-\nu\,\partial_{x}\lambda}{\beta_{u}+\lambda},$
(29b)

where
$\mathcal{D}:=(\lambda,\gamma,\partial_{t}\lambda,\partial_{x}\lambda,\partial_{x}\gamma)$.
Substituting the DtP mapping back into (28) gives us the functional
$S_{H}[\lambda,\mu]$ which can be explicitly written as

$S_{H}[\lambda,\gamma]=\int_{\Omega}\,\left(-\frac{\mathbb{K}_{1}}{\beta_{Y}}\,(\partial_{t}\lambda\,+\,\partial_{x}\gamma)^{2}-\frac{\mathbb{K}_{2}}{\beta_{u}}\,(\nu\,\partial_{x}\lambda+\bar{u}\lambda-\gamma)^{2}\right)\,dt\,dx\\\
+\bigintsss_{\Omega}\,\bigg{(}\bar{Y}(\partial_{t}\lambda+\partial_{x}\gamma)+\bar{u}\left(\frac{\lambda\bar{u}}{2}-\gamma+\nu\,\partial_{x}\lambda\right)\bigg{)}\,dt\,dx\\\
-\int_{0}^{L}\,Y_{0}(x)\,\lambda(x,0)\,dx-\int_{0}^{T}\,Y_{l}(t)\,\gamma(0,t)\,dt,$

where

$\mathbb{K}_{1}=\frac{1}{2};\qquad\mathbb{K}_{2}=\frac{1}{2\left(1+\frac{\lambda}{\beta_{u}}\right)},$

and for $\beta_{u}\gg|\lambda|$, $\partial_{x}\lambda$ now gets penalized
(multiplied by $\nu$). Of course, the E-L equations of $S_{H}$ is (27) with
the standard substitution of $(Y,u)\to(\hat{Y},\hat{u})$.

Considering the ellipticity of E-L equations of $S_{H}$, i.e. (27) with the
above replacement, the matrix in (21) now has a modification:

$\mathbb{A}_{1212}=\frac{\nu^{2}}{\beta_{u}+\lambda},$

with the other terms unchanged. Correspondingly, we have

$c_{i}A_{ij}c_{j}=\frac{(c_{1}+c_{4})^{2}}{\beta_{Y}}+\frac{c_{2}^{2}\,\nu^{2}}{\beta_{u}+\lambda}=\frac{(c_{1}+c_{4})^{2}}{\beta_{Y}}+\frac{c_{2}^{2}\,\nu^{2}}{\beta_{u}\left(1+\frac{\lambda}{\beta_{u}}\right)}\qquad\forall
c\in\mathbb{R}^{4}.$

In case when $\lambda=0$, we always have $c_{i}A_{ij}c_{j}\geq 0$ for non-
trivial $c$. Additionally in a neighborhood of $\lambda=0$, where
$|\lambda|\ll\beta_{u}$,

$c_{i}A_{ij}c_{j}\geq 0\quad\forall c\in\mathbb{R}^{4}.$

This establishes that $\mathbb{A}_{ijkl}$ is positive semi-definite in this
neighborhood and the equation set (27) is locally degenerate elliptic. In
contrast to the inviscid Burgers-HJ, the only rank-one matrices along which
$C:\mathbb{A}C=0$ are of the form $(0,a_{2})\otimes(1,0)$ for any
$a_{2}\in\mathbb{R}$. Based on the discussion in Sec. 3, it follows that
ellipticity is lost only along the normal $(1,0)$, where the first argument
represents the component of the normal along the time direction. This is
consistent with solutions of Burgers equation $(\nu\neq 0)$ as spatial
discontinuities in $u$ are not allowed. The only discontinuity allowed is in
$\partial_{t}\gamma$ across space-like surfaces ($t=$ constant) but the DtP
mapping equations (29) shows that this does not affect $u^{(H)}$ and hence the
mapped primal solutions are continuous; a similar conclusion was also obtained
in our dual formulation for the heat equation (see-[1, Sec. 5]).

Similar to the inviscid case, we obtain the following side conditions:

$\hat{Y}(x,0)=Y_{0}(x);\qquad\hat{Y}(0,t)=Y_{l}(t)$ (30)

where $\hat{Y}=Y^{(H)}(\mathcal{D}(x,t),x,t)$ . Additionally, on each of the
boundaries, the following operation can be performed which lead to different
classes of boundary conditions and hence different problem setups:

  * •

$\lambda$ can be specified as a Dirichlet BC with a fixed arbitrary value. In
such a case, $\delta\lambda=0$ on the boundary under consideration and we
obtain (30) as the only natural b.cs.

  * •

If $\lambda$ is not specified and allowed to vary freely, then an extra
natural BC for $S_{H}$,

$\nu\,\hat{u}=f_{i}$

emerges, corresponding to the boundary under consideration.

Our focus is primarily on observing the effect of the $\nu$ term in the bulk
of the domain. Accordingly, we consider the second BC above, weakly enforcing
$u=0$ on both ends. The impact of this term on the boundaries becomes more
pronounced for higher values of $\nu$, and we want to operate with as small a
value of $\nu$ as possible. In the following, we show calculations for
$\nu=10^{-3}$ for the Fan, the Half N-wave and the N-wave problems.

The residual (22) and Jacobian (25) and their corresponding discrete versions
are now appropriately modified.

####  5.3.1 Fan

Considering 5.2.1 with an active $\nu$, the outcomes are presented in Fig. 13.
It is apparent from the figure that the results have exhibited enhancement and
closely resemble the exact outcomes of the inviscid case. Upon close
inspection, a slight deviation from the exact expression can be noticed close
to the right boundary. This discrepancy can be attributed to the influence of
an active $\nu$ and more specifically, to the natural imposition of $\nu u=0$
on that boundary.

(a) $Y(x,t)$

(b) Exact $Y(x,t)$ based on (39)

(c) Line plots for $Y(x,t)$

Figure 13: DtP mapping generated primal field $Y$ for the fan example with
$\nu=10^{-3}$. In Fig.(c), the dots in various colors at last spatial Gauss
point correspond to the exact value of $Y$ derived from (39), corresponding to
the time indicated by the color of each respective dot.

####  5.3.2 Half N-wave

Considering 5.2.4 again with an active $\nu$ term, the results as shown in
Fig. 14. The boundary effects in this problem are not very pronounced because
the Half N-wave itself has a value of $u=0$ on the right boundary for a
significant duration based on the exact expression. With a small viscosity
$\nu$, one would not anticipate the solution for this problem to be
significantly different from the case where $\nu=0$, which is also evident
from the results presented.

(a) $Y(x,t)$

(b) Exact $Y(x,t)$ based on (45)

(c) Line plots for $Y(x,t)$

Figure 14: DtP mapping generated primal field $Y$ for the Half N-wave with
$\nu=10^{-3}$. In Fig.(c), the coloured dots at $Y=0.25$ represent the
trajectory of the kink in the space-time domain based on the exact expressions
(43) for $\nu=0$.

####  5.3.3 N-wave

Considering the problem 5.2.5 again with an active $\nu$, the results are
presented in 15, from which it is evident that the two dips have been
resolved. Furthermore, the line plots for $Y$ derived from the $L^{2}$
projected data at different nodal times are displayed, revealing the existence
of a kink at the node located at $x=0.5$. A slight deviation in the speed at
which the kink travels down from the exact expression may again be attributed
to the presence of an active $\nu$.

(a)

(b)

(c)

(d)

Figure 15: DtP mapping generated primal field $Y$ for the N-wave with
$\nu=10^{-3}$. In Fig. 15(c) and Fig. 15(d), the colored dots at $X=0.5$
represent the trajectory of the kink in the space-time domain based on the
exact expressions (49) for $\nu=0$.

###  5.4 Inviscid Burgers-HJ with Viscous base states

In this section we will explore solutions to the inviscid Burgers-HJ equation
generated from the dual scheme with base states provided from exact solutions
to the (viscous) Burgers-HJ equation with small viscosity.

Given an initial condition $u_{0}(x)$ for the Burgers equation, for
$x\in\mathbb{R}$, the Hopf-Cole transformation can be used to derive the exact
expression for $u(x,t)$ [11] along with its antiderivative with respect to
$x$, $Y(x,t)$, such that for $t>0$ and $x\in\mathbb{R}$,

$\displaystyle Y(x,t)=-2\nu\log\left[\frac{1}{2\sqrt{\pi\nu
t}}\int_{-\infty}^{\infty}\,e^{-\frac{Y_{0}(y)}{2\nu}-\frac{(x-y)^{2}}{4\nu
t}}\,dy\right];$ (31a) $\displaystyle
u(x,t)=\frac{\int_{-\infty}^{\infty}\frac{x-y}{t}e^{-\frac{(x-y)^{2}}{4\nu
t}-\frac{Y_{0}(x)}{2\nu}}\,dy}{\int_{-\infty}^{\infty}e^{-\frac{(x-y)^{2}}{4\nu
t}-\frac{Y_{0}(x)}{2\nu}}dy},$ (31b)

where $Y_{0}$ can be generated via

$Y_{0}(x,t)=\int_{0}^{x}\,u_{0}(y,t)\,dy$

and we set $Y_{0}(0)=0$ without loss of generality. These expressions will be
denoted as the ‘viscous’ formulae.

For $\nu\rightarrow 0$, the formula (31) for $u(x,t)$ recovers the solution to
the inviscid Burgers equation. Motivated by this fact, we set a small value
$\nu=10^{-3}$ and generate the solutions to the Burgers equation using the
expressions (31). We consider the algorithm presented for the Burgers equation
(Sec. 5.1) and use the same settings as used for the Half N-wave (Sec. 5.1.4),
except that the base state for $u$ is now set based the expression obtained
using the formula (31b). The result for such a setup is shown in Fig. 16(a).
Similarly for the Burgers-HJ, we utilize the algorithm on the Half N-wave
(Sec. 5.2.4), with the base states for $Y$ and $u$ obtained from the formulae
31a and 31b, respectively. We use coarser meshes as compared to 5.2.4 and set
$N_{x}=50$, $N_{t}=400$ and $T^{(s)}=0.01$ (fine meshes work as well). The
result for this setup has been shown in Fig. 16(b).

Based on the results we observe the following:

  * •

In the case of the inviscid Burgers equation, the unique entropy solution is
again recovered, now with no smoothing operator for base state generation
(which is to be expected).

  * •

In the case of the dual inviscid Burgers-HJ, using the viscous base states is
one way of enforcing the entropy condition.

(a) $u(x,t)$

(b) $Y(x,t)$

Figure 16: (a) is obtained from the formulation based on the Burgers equation
Sec. 5.2. (b) is produced using the formulation based on Burgers-HJ Sec. 5.1.
Respective base states used for both the problems are based on the viscous
formulas (31).

## Acknowledgment

This work was supported by the grant NSF OIA-DMR #2021019 and was also
supported by a Simons Pivot Fellowship to AA. It was completed while he was on
sabbatical leave at a) the Max Planck Institute for Mathematics in the
Sciences in Leipzig, and b) the Hausdorff Institute for Mathematics at the
University of Bonn funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC-2047/1 –
390685813, as part of the Trimester Program on Mathematics for Complex
Materials. The support and hospitality of both institutions is acknowledged.

## Appendices

##  Appendix A Review of the general formalism

The following two subsections are excerpted from [8, 6] to make this paper
self-contained.

###  A.1 The essential idea: An optimization problem for an algebraic system
of equations

Consider a generally nonlinear system of algebraic equations in the variables
$x\in\mathbb{R}^{n}$ given by

$A_{\alpha}(x)=0,$ (32)

where $A:\mathbb{R}^{n}\to\mathbb{R}^{N}$ is a given function (a simple
example would be $A_{\alpha}(x)=\bar{A}_{\alpha i}\,x^{i}-b_{\alpha}$,
$\alpha=1\ldots N,i=1\ldots n$, where $\bar{A}$ is a constant matrix, not
necessarily symmetric (when $n=N$), and $b$ is a constant vector). We allow
for all possibilities $0<n\lesseqqgtr N>0$.

The goal is to construct an objective function whose critical points solve the
system (32) (when a solution exists) by defining an appropriate
$x^{*}\in\mathbb{R}^{n}$ satisfying $A_{\alpha}(x^{*})=0$.

For this, consider first the auxiliary function

$\widehat{S}_{H}(x,z)=z^{\alpha}A_{\alpha}(x)+H(x)$

(where $H$ belongs to a class of scalar-valued function to be defined shortly)
and define

$S_{H}(z)=z^{\alpha}A_{\alpha}(x_{H}(z))+H(x_{H}(z))$

with the requirement that the system of equations

$z^{\alpha}\frac{\partial A_{\alpha}}{\partial x^{i}}(x)+\frac{\partial
H}{\partial x^{i}}(x)=0$ (33)

be solvable for the function $x=x_{H}(z)$ through the choice of $H$, and any
function $H$ that facilitates such a solution qualifies for the proposed
scheme.

In other words, given a specific $H$, it should be possible to define a
function $x_{H}(z)$ that satisfies

$z^{\alpha}\partial_{x^{i}}A_{\alpha}(x_{H}(z))+\partial_{x^{i}}H(x_{H}(z))=0\quad\forall
z\in\mathbb{R}^{N}$

(the domain of the function $x_{H}$ may accommodate more intricacies, but for
now we stick to the simplest possibility). Note that (33) is a set of $n$
equations in $n$ unknowns regardless of $N$ ($z$ for this argument is a
parameter).

Assuming this is possible, we have

$\frac{\partial S_{H}}{\partial
z^{\beta}}(z)=A_{\beta}(x_{H}(z))+\left(z^{\alpha}\frac{\partial
A_{\alpha}}{\partial x^{i}}(x_{H}(z))+\frac{\partial H}{\partial
x^{i}}(x_{H}(z))\right)\frac{\partial x^{i}_{H}}{\partial
z^{\beta}}(z)=A_{\beta}(x_{H}(z)),$

using (33). Thus,

  * •

if $z_{0}$ is a critical point of the objective function $S_{H}$ satisfying
$\partial_{z^{\beta}}S_{H}(z_{0})=0$, then the system $A_{\alpha}(x)=0$ has a
solution defined by $x_{H}(z_{0})$;

  * •

if the system $A_{\alpha}(x)=0$ has a unique solution, say $y$, and if
$z^{H}_{0}$ is any critical point of $S_{H}$, then
$x_{H}\left(z^{H}_{0}\right)=y$, for all admissible $H$.

  * •

If $A_{\alpha}(x)=0$ has non-unique solutions, but
$\partial_{z^{\beta}}S(z)=0$ ($N$ equations in $N$ unknowns) has a unique
solution for a specific choice of the function $z\mapsto x_{H}(z)$ related to
a choice of $H$, then such a choice of $H$ may be considered a selection
criterion for imparting uniqueness to the problem $A_{\alpha}(x)=0$.

  * •

Finally, to see the difference of this approach with the Least-Squares (LS)
Method, we note that the optimality condition for the objective
$A_{\alpha}(x)A_{\alpha}(x)$ is
$A_{\alpha}(x)\partial_{x^{i}}A_{\alpha}(x)=0\centernot\implies
A_{\alpha}(x)=0$.

For a linear system $\bar{A}x=b$, the LS governing equations are given by

$\bar{A}^{T}\bar{A}z=\bar{A}^{T}b,$

with LS solution defined as $z$ even when the original problem $\bar{A}x=b$
does not have a solution (i.e., when $b$ is not in the column space of
$\bar{A}$). The LS problem always has a solution, of course. In contrast, in
the present duality-based approach with quadratic $H(x)=\frac{1}{2}x^{T}x$ the
governing equation is

$-\bar{A}\bar{A}^{T}z=b$

with solution to $\bar{A}x=b$ given by $x=-\bar{A}^{T}z$, and the problem has
a solution only when $\bar{A}x=b$ has a solution, since the column spaces of
the matrices $\bar{A}$ and $\bar{A}\bar{A}^{T}$ are identical.

###  A.2 The idea behind the general formalism

The proposed scheme for generating variational principles for nonlinear PDE
systems may be abstracted as follows: We first pose the given system of PDE as
a first-order system (introducing extra fields representing (higher-order)
space and time derivatives of the fields of the given system); as before let
us denote this collection of primal fields by $U$. ‘Multiplying’ the primal
equations by dual Lagrange multiplier fields, the collection denoted by $D$,
adding a function $H(U)$, solely in the variables $U$ (the purpose of which,
and associated requirements, will be clear shortly), and integrating by parts
over the space-time domain, we form a ‘mixed’ functional in the primal and
dual fields given by

$\widehat{S}_{H}[U,D]=\int_{[0,T]\times\Omega}\
\mathcal{L}_{H}(\mathcal{D},U)\,dtdx$

where $\mathcal{D}$ is a collection of local objects in $D$ and at most its
first order derivatives. We then require that the family of functions $H$ be
such that it allows the definition of a function $U_{H}(\mathcal{D})$ such
that

$\frac{\partial\mathcal{L}_{H}}{\partial U}(\mathcal{D},U_{H}(\mathcal{D}))=0$

so that the dual functional, defined solely on the space of the dual fields
$D$, given by

$S_{H}[D]=\int_{[0,T]\times\Omega}\mathcal{L}_{H}(\mathcal{D},U_{H}(\mathcal{D}))\,dtdx$

has the first variation

$\delta
S_{H}=\int_{[0,T]\times\Omega}\frac{\partial\mathcal{L}_{H}}{\partial\mathcal{D}}\delta\mathcal{D}\,dtdx.$

By the process of formation of the functional $\widehat{S}_{H}$, it can then
be seen that the (formal) E-L equations arising from $\delta S_{H}$ have to be
the original first-order primal system, with $U$ substituted by
$U_{H}(\mathcal{D})$, regardless of the $H$ employed.

Thus, the proposed scheme may be summarized as follows: we wish to pursue the
following (local-global) critical point problem

$\begin{smallmatrix}\mbox{extremize}\\\
D\end{smallmatrix}\int_{[0,T]\times\Omega}\begin{smallmatrix}\mbox{extremize}\\\
U\end{smallmatrix}\ \mathcal{L}_{H}(\mathcal{D}(t,x),U)\,dtdx,$

where the pointwise extremization of $\mathcal{L}_{H}$ over $U$, for fixed
$\mathcal{D}$, is made possible by the choice of $H$.

Furthermore, assume the Lagrangian $\mathcal{L}_{H}$ can be expressed in the
form

$\mathcal{L}_{H}(\mathcal{D},U):=-P(\mathcal{D})\cdot U+f(U,D)+H(U)$

for some function $P$ defined by the structure of the primal first-order
system ((linear terms in) first derivatives of $U$ after multiplication by the
dual fields and integration by parts always produce such terms), and for some
function $f$ which, when non-zero, does not contain any linear dependence in
$U$. Our scheme requires the existence of a function $U_{H}$ defined from
‘solving $\frac{\partial\mathcal{L}}{\partial U}(\mathcal{D},U)=0$ for $U$,’
i.e. $\exists\ U_{H}(P(\mathcal{D}),\mathcal{D})$ s.t. the equation

$-P(\mathcal{D})+\frac{\partial f}{\partial
U}(U_{H}(P(\mathcal{D}),\mathcal{D}),\mathcal{D})+\frac{\partial H}{\partial
U}\left(U_{H}(P(\mathcal{D}),\mathcal{D})\right)=0$

is satisfied. This requirement may be understood as follows: define

$f(U,\mathcal{D})+H(U)=:M(U,\mathcal{D})$

and assume that it is possible, through the choice of $H$, to make the
function $\frac{\partial M}{\partial U}(U,\mathcal{D})$ monotone in $U$ so
that a function $U_{H}(P,\mathcal{D})$ can be defined that satisfies

$\frac{\partial M}{\partial
U}(U_{H}(P,\mathcal{D}),\mathcal{D})=P,\quad\forall P.$

Then the Lagrangian is

$\mathcal{L}(\mathcal{D},U_{H}(P(\mathcal{D}),\mathcal{D}))=-P(\mathcal{D})\cdot
U_{H}(P(\mathcal{D}),\mathcal{D})+M(U_{H}(P(\mathcal{D}),\mathcal{D}),\mathcal{D})=:-M^{*}(P(\mathcal{D}),\mathcal{D})$

where $M^{*}(P,\mathcal{D})$ is the Legendre transform of the function $M$
w.r.t $U$, with $\mathcal{D}$ considered as a parameter.

Thus, our scheme may also be interpreted as designing a concrete realization
of abstract saddle point problems in optimization theory [12], where we
exploit the fact that, in the context of ‘solving’ PDE viewed as constraints
implemented by Lagrange multipliers to generate an unconstrained problem,
there is a good deal of freedom in choosing an objective function(al) to be
minimized. We exploit this freedom in choosing the function $H$ to develop
dual variational principles corresponding to general systems of PDE.

##  Appendix B Application of base states: a simple example

Figure 17: Schematic of Circle-line intersection.

Consider the following algebraic system of equations for
$(x,y)\in\mathbb{R}^{2},\alpha\in\mathbb{R}$:

$\displaystyle x^{2}+y^{2}$ $\displaystyle=1$ (34) $\displaystyle x$
$\displaystyle=\alpha.$

A schematic is shown in Fig. 17, and solutions are given by
$\\{(\alpha,\pm\sqrt{1-\alpha^{2}}):|\alpha|\leq 1\\}$. We use the logic of
Appendix A.1 to multiply each of the above equations by a dual multiplier and
add a quadratic auxiliary potential $H$:

$\widehat{S}(x,y,\lambda,\gamma)=\lambda(x^{2}+y^{2}-1)+\gamma(x-\alpha)+\frac{1}{2}(x-\bar{x})^{2}+\frac{1}{2}(y-\bar{y})^{2},$

where $\bar{x},\bar{y}\in\mathbb{R}$ are constants (in this algebraic
problem), and we will refer to the pair as a ‘base state.’ Next we need to
generate the analog of the mapping function $x_{H}$ of Appendix A.1, also
referred to as the DtP mapping in the text:

$\displaystyle\frac{\partial S}{\partial x}=0:\quad 2\lambda
x+\gamma+(x-\bar{x})=0$ $\displaystyle\Rightarrow\quad\mbox{for}\
\lambda\neq-\frac{1}{2},\quad
x_{H}(\lambda,\gamma)=\frac{\bar{x}-\gamma}{2\lambda+1};$ (35)
$\displaystyle\frac{\partial S}{\partial y}=0:\quad 2\lambda y+(y-\bar{y})=0$
$\displaystyle\Rightarrow\quad\mbox{for}\ \lambda\neq-\frac{1}{2},\quad
y_{H}(\lambda)=\frac{\bar{y}}{2\lambda+1}.$

$\lambda\neq-\frac{1}{2}$: Considering only the case $\lambda\neq-\frac{1}{2}$
for the moment, the dual objective function is now obtained by substituting
the DtP mapping into $\widehat{S}$:

$S(\lambda,\gamma)=\lambda\left(x_{H}^{2}(\lambda,\gamma)+y_{H}^{2}(\lambda,\gamma)-1\right)+\gamma\left(x_{H}(\lambda,\gamma)-\alpha\right)+\frac{1}{2}\left(x_{H}(\lambda,\gamma)-\bar{x}\right)^{2}+\frac{1}{2}\left(y_{H}(\lambda,\gamma)-\bar{y}\right)^{2}.$

The critical point equations for this objective, by design, are the equations
(34) with the substitution $(x\to x_{H},y\to y_{H})$:

$\displaystyle\frac{(\bar{x}-\gamma)^{2}}{(2\lambda+1)^{2}}+\frac{\bar{y}^{2}}{(2\lambda+1)^{2}}$
$\displaystyle=1$ $\displaystyle\frac{\bar{x}-\gamma}{2\lambda+1}$
$\displaystyle=\alpha.$

A necessary condition for solutions is

$(2\lambda+1)^{2}(1-\alpha^{2})-\bar{y}^{2}=0$

which implies that dual solutions exist only for $|\alpha|\leq 1$ and when
$|\alpha|=1$ only if $\bar{y}=0$.

Thus, for $|\alpha|<1$, $\bar{y}\neq 0$,

$\lambda=\frac{1}{2}\left(\pm\frac{|\bar{y}|}{\sqrt{1-\alpha^{2}}}-1\right);\qquad\gamma=\bar{x}-\alpha(2\lambda+1).$

are the extrema of $S$.

Since only $\lambda\neq-\frac{1}{2}$ is being considered, we do not consider
the case $|\alpha|<1$, $\bar{y}=0$ here.

For $|\alpha|=1$, $\bar{y}=0$, the pairs

$-\frac{1}{2}\neq\lambda\in\mathbb{R}\
\mbox{arbitrary};\qquad\gamma=\bar{x}-\alpha(2\lambda+1).$

are the extrema of $S$.

Putting these dual solutions back into the DtP mapping (35), we recover the
correct primal solutions as expected.

$\lambda=-\frac{1}{2}$: When $\lambda=-\frac{1}{2}$, extrema of $S$ exist only
for $\bar{y}=0$ and are given by $(\lambda=-\frac{1}{2},\gamma=\bar{x})$.
However, note that this class of dual solutions does not define solutions to
the primal problem in a non-vacuous manner (i.e. while all primal solutions
are admitted, it does not provide specific guidance for generating primal
solutions).

The following conclusion can be drawn from this simple, yet non-trivial,
example:

  * •

A ‘good’ choice of the auxiliary function can be crucial for the success of
the dual scheme in generating solutions to the primal problem. For example, if
‘no base states’ are invoked within this class of quadratic auxiliary
functions, i.e. $(\bar{x},\bar{y})=(0,0)$, while dual extrema exist, the
scheme essentially fails to define primal solutions, except for the case
$|\alpha|=1$. And when the latter is the case for the primal problem, then a
base state with non-zero $\bar{y}$ is not feasible for defining dual and
primal solutions to the problem.

##  Appendix C $L^{2}$ projection

Let $u$ represent the primal field under consideration, which may depend on
the derivatives of the dual fields (approximated using linear FE shape
functions). To achieve a $C^{0}$ continuous approximation of $u$ in the domain
$\Omega\subset\mathbb{R}\mbox{ or }\mathbb{R}^{2}$, we employ the following
method: Let $u_{h}$ represent the projection of $u$ onto a space $V_{h}$
formed by the linear span of globally continuous, piecewise smooth finite
element shape functions corresponding to a FE mesh for $\Omega$. We enforce
the following condition upon $u_{h}$:

$u_{h}=\operatorname*{arg\,min}_{v\in
V_{h}}\int_{\Omega}\frac{1}{2}|u-v|^{2}\,d\Omega.$

Let $N^{A}$ represent the basis functions associated with any node with an
index $A$. The discrete version of the optimality condition of the above
statement is

$\sum_{A=1}^{N}\sum_{B=1}^{N}\delta
u^{A}_{h}\biggl{(}\int_{\Omega}N^{A}\,N^{B}\,d\Omega\biggl{)}u_{h}^{B}=\sum_{A=1}^{N}\delta
u^{A}_{h}\int_{\Omega}N^{A}\,u\,d\Omega,$ (36)

where $u_{h}^{B}$ denotes the nodal value at node $B$ of the sought continuous
projection of $u$, and $\delta u_{h}:=\delta u^{A}_{h}N^{A}$ is a test
function. The primal data for the integral in the right hand side of (36) is
obtained from the Gauss points of 1-D/2-D element. In solving (36), we impose
any known data, e.g. initial and boundary conditions of the primal problem, as
known function values $u^{A}_{h}$, at corresponding nodes $A$ where the data
is known, with $\delta u^{A}_{h}=0$. Using the arbitrariness of the remaining
‘free’ $\delta u^{A}_{h}$, we obtain a system of linear equations to solve for
the unconstrained nodal values $u_{h}^{B},B\in\\{1,\ldots,N\\}$.

If the $u$ profile of the primal field is required only at a specific time,
then conducting a 1-D $L^{2}$ projection along that timeline is advisable
since it is computationally inexpensive. As a consequence of $C^{0}$ type FE
shape functions being employed, the time derivatives of the dual fields
experience discontinuities across nodal timelines. In such a case the value of
the time derivative on any nodal timeline can be approximated as an average of
the value obtained from the element above and below the spatial point of
interest (these values remain constant in the direction of time within any
element). With the inputs to the DtP mapping generated in this manner, the
value of primal fields at the Gauss points (of this nodal timeline) are
obtained.

##  Appendix D The base state smoothing operator

The smoothing operator $\mathcal{S}$ serves to mitigate the high-wave number
contributions in any function. Let $u:=\mathcal{S}[f]$ represent the smoothed
output produced for an input function $f$. Then we require $u$ to satisfy:

$u-\eta\,\partial_{xx}u=f,\qquad\eta>0.$ (37)

Here, $\eta$ represents a diffusion-like control parameter.

Let $t_{s}$ represent the time at which the smoothing operation occurs. When
boundary conditions (BC) are available at $x=0$ $(u_{l}(t_{s}))$ or (and)
$x=L$ $(u_{r}(t_{s}))$ from the primal problem, we employ them as the BCs for
$u$. However, when the boundary data is not available, we calculate the
average of the Gauss data values for the first element ($f_{1}^{g_{1}}$ and
$f_{1}^{g_{2}}$) and set it as the left BC, and the average of the Gauss data
values for the last element ($f_{l}^{g_{1}}$ and $f_{l}^{g_{2}}$) set it as
the right BC on $u$. Correspondingly, we generate a weak form for equation
(37). Let $\delta u$ be the space of continuous test functions on $(0,L)$ with
$\delta u(0)=\delta u(L)=0$. The we intend to find $u$ which satisfies, for
such test functions,

$\int_{0}^{L}dx\bigl{(}u\,\delta u+\eta\,\partial_{x}u\,\partial_{x}(\delta
u)-f\,\delta u\bigl{)}=0.$

We use a standard Galerkin FEM discretization to solve the above with an
$\eta=10^{-4}$.

##  Appendix E Exact entropy solutions for the provided examples

The exact entropy solutions presented in the following examples are based on
the concepts and examples outlined in [10]. Given an initial condition
$x\mapsto u_{0}(x)\in\mathbb{R}$, entropy solutions to the Burgers equation
($u$) are obtained using two important concepts, recalled below from [10]:

  1. 1.

Rankine-Hugoniot/jump condition: This is simply a consequence of the
underlying conservation statement for Burgers equation as well as looking for
weak solutions (in the sense of distributions) to it. If the solution to
$\partial_{t}u+\partial_{x}\left(\frac{u^{2}}{2}\right)=0$ has different
values across a shock curve $x=X(t)$, then the shock speed $s=\frac{dX}{dt}$
must satisfy

$s=\frac{1}{2}\left(u(X^{+}(t))+u(X^{-}(t))\right),$

where $X^{-}$ and $X^{+}$ denote the regions immediately to the left and right
of the shock, respectively.

  2. 2.

The Entropy condition: The characteristic lines must go into the shock as $t$
increases, so those on the left go faster (measured by $\frac{dX}{dt}$) than
the shock and those on the right go slower. With respect to the Burgers
equation this implies

$u(X^{-}(t))>s>u(X^{+}(t)).$

###  E.1 Expansion fan

Consider the initial condition

$u_{0}(x)=\begin{cases}0&\mbox{for }x<0.5\\\ 1&\mbox{for }x>0.5.\end{cases}$

The initial signals for $x>0.5$ travel along characteristics with slope $1$,
while the signals at $x<0.5$ travel along vertical characteristics (slope $0$,
i.e. do not move in space). This leads to an expansion fan filling the space-
time in between $x=0.5$ and $x=0.5+t$, which is a rarefaction wave traveling
along $x$. The corresponding expression for $u(x,t)$ can then be given as

$u(x,t)=\begin{cases}\quad\,0&\mbox{for }\quad\,x<0.5\\\\[2.0pt]
\dfrac{x-0.5}{t}&\mbox{for }\,0.5+t>x\geq 0.5\\\\[4.0pt] \quad\,1&\mbox{for
}\quad\,x\geq 0.5+t.\end{cases}$ (38)

The corresponding $Y$ profile for $Y_{0}(0)=0$ is given by

$Y(x,t)=\begin{cases}\qquad\quad\,0&\mbox{for }\qquad\,x<0.5\\\\[2.0pt]
\dfrac{x^{2}}{2t}-\dfrac{0.5\,x}{t}+\dfrac{1}{8t}&\mbox{for }\,0.5+t>x\geq
0.5\\\\[6.0pt] \quad\,x-0.5-\dfrac{t}{2}&\mbox{for }\quad\,x\geq
0.5+t.\end{cases}$ (39)

###  E.2 Shock

Consider the initial condition

$u_{0}(x)=\begin{cases}1&\mbox{for }x<0.5\\\ 0&\mbox{for }x>0.5.\end{cases}$

Given the presence of a jump in the initial condition from 1 to 0 at $x=0.5$,
the characteristics originating from the left will merge with the
characteristics from the right, resulting in the traveling shock wave with
velocity $s=\frac{1+0}{2}=\frac{1}{2}$. The exact entropy solution is given by

$u(x,t)=\begin{cases}1&\mbox{for }x<0.5+t/2\\\ 0&\mbox{for
}x>0.5+t/2.\end{cases}$ (40)

The corresponding $Y$ profile for $Y_{0}(0)=0$ is given by

$Y(x,t)=\begin{cases}x-\frac{t}{2}&\mbox{for }x<0.5+t/2\\\ 0.5&\mbox{for
}x>0.5+t/2.\end{cases}$ (41)

###  E.3 Double shock

Consider the initial condition

$u_{0}(x)=\begin{cases}\begin{aligned} 1&\quad\text{for }x<0.25\\\
0.5&\quad\text{for }0.25<x<0.5\\\ 0&\quad\text{for
}x>0.5.\end{aligned}\end{cases}$

Following the same ideas from example E.2, the shock starting at $x=0.25$ will
travel with a velocity of $s_{1}=\frac{1+0.5}{2}=0.75$, while the shock
starting at $x=0.5$ will travel at with a velocity
$s_{2}=\frac{0.5+0}{2}=0.25$. Since, $s_{1}>s_{2}$, the first shock will merge
with the second shock at a time $t_{m}$ and travel with a velocity
$s=\frac{1+0}{2}=0.5$. The time $t_{m}$ for the faster shock to catch up with
the slower shock is given by

$s_{1}\,t_{m}=s_{2}\,t_{m}+0.25\quad\Rightarrow\quad t_{m}=0.5.$

Consequently the shocks meet at $x=0.5+s_{2}t_{m}=0.625$. The solution
$u(x,t)$ is then be given by

$u(x,t)=\begin{cases}\begin{rcases}1&\,\text{for }\,\qquad x<0.25+0.75t\\\
0.5&\,\text{for }\,0.25+0.75t<x<0.5+0.25t\\\ 0&\,\text{for }\,\qquad
x>0.5+0.25t\end{rcases}t<t_{m}\\\ \\\ \begin{rcases}1&\quad\text{for }\,\qquad
x<0.625+0.5(t-t_{m})\\\ 0&\quad\text{for }\,\qquad
x>0.625+0.5(t-t_{m})\end{rcases}t>t_{m}.\end{cases}$ (42)

The corresponding $Y$ profile for $Y_{0}(0)=0$ is given by

$Y(x,t)=\begin{cases}\begin{rcases}x-t/2&\,\text{for }\,\qquad x<0.25+0.75t\\\
x/2+0.125-t/8&\,\text{for }\,0.25+0.75t<x<0.5+0.25t\\\ 0.375&\,\text{for
}\,\qquad x>0.5+0.25t\end{rcases}t<t_{m}\\\ \\\
\begin{rcases}x-t/2&\qquad\text{for }\,\qquad x<0.625+(t-t_{m})\\\
0.375&\qquad\text{for }\,\qquad
x>0.625+(t-t_{m})\end{rcases}t>t_{m}.\end{cases}$ (43)

Figure 18: The right-angled triangle with a height $h_{0}$ indicates the
initial profile for the Half N-wave. The right-angled triangle with the height
$h$ indicates the anticipated profile at some later time.

###  E.4 Half N-Wave

Consider the initial condition

$\displaystyle u_{0}(x)=\begin{cases}0&\text{for }x<x_{0}\\\\[2.0pt]
\dfrac{h_{0}}{l_{0}}(x-x_{0})&\text{for }x_{0}\leq x<x_{0}+l_{0}\\\\[4.0pt]
0&\text{for }x>x_{0}+l_{0}.\end{cases}$

This profile has been depicted in Fig. 18, and it is apparent that the point
at $x=x_{0}+l_{0}$ corresponds to a shock, as
$u_{0}((x_{0}+l_{0})^{-})>u_{0}((x_{0}+l_{0})^{+})$. Furthermore, the point at
$x_{0}$ serves as a source for an expansion fan due to the linear increase
observed from $x=x_{0}$ to $x=x_{0}+l_{0}$. For any time $t$, let the shock be
present at a point $X(t)=x_{0}+l(t)$, where $l(0)=l_{0}$ and the height of
shock be given by $h(t)$, where $h(0)=h_{0}$. To obtain the behavior of the
shock front, we utilize the following information:

  * •

Burgers equation on $\mathbb{R}$ is a conservation law for $u$ which is
conserved:

$\int_{x_{l}}^{x_{r}}\partial_{t}u\,dx=0\quad\Rightarrow\frac{d}{dt}\int_{x_{l}}^{x_{r}}u\,dx=0\quad\Rightarrow\frac{d}{dt}\int_{x_{l}}^{x_{r}}u\,dx=0$

for such time that the shock stays with $[x_{l},x_{r}]$. This implies

$\frac{1}{2}h_{0}\,l_{0}=\frac{1}{2}h(t)\,l(t).$

  * •

Using the Rankine-Hugoniot condition and the last obtained expression across
the shock, we establish that

$s(t)=\frac{dX(t)}{dt}=h(t)/2\quad\Rightarrow\frac{dl(t)}{dt}=\frac{1}{2}\frac{h_{0}\,l_{0}}{l(t)}$

The last expression leads to

$l(t)=\sqrt{h_{0}\,l_{0}\,t+l_{0}^{2}};\qquad X(t)=x_{0}+l(t);\quad
h(t)=\frac{h_{0}\,l_{0}}{l(t)}.$

The solution $u(x,t)$ is given by

$\displaystyle u(x,t)=\begin{cases}0&\text{for }x<x_{0}\\\\[2.0pt]
\dfrac{h(t)}{l(t)}(x-x_{0})&\text{for }x_{0}\leq x<x_{0}+l(t)\\\\[4.0pt]
0&\text{for }x>x_{0}+l(t).\end{cases}$ (44)

The corresponding $Y$ profile for $Y_{0}(0)=0$ is given by

$\displaystyle Y(x,t)=\begin{cases}0&\text{for }x<x_{0}\\\\[4.0pt]
\dfrac{h(t)}{l(t)}\left(\dfrac{x^{2}+x_{0}^{2}}{2}-x_{0}\,x\right)&\text{for
}x_{0}\leq x<x_{0}+l(t)\\\\[10.0pt] \dfrac{h_{0}\,l_{0}}{2}&\text{for
}x>x_{0}+l(t).\end{cases}$ (45)

For the example considered in this paper, we have used the values
$x_{0}=0.25$, $l_{0}=0.25$ and $h_{0}=2$.

###  E.5 N-wave

Consider the profile given by

$\displaystyle u_{0}(x)=\begin{cases}0&\text{for }x<0.25\\\
-4h_{0}(x-0.5)&\text{for }0.25<x<0.75\\\ 0&\text{for }x>0.75.\end{cases}$ (46)

For the example considered in this paper ($h_{0}=2$), the characteristic plane
has been shown in Fig. 19(a). Each characteristic ray from $x^{*}$ at time
$t=0$ is a straight line with slope $u(x^{*})$ (Strang, pg. 591 [10]).
Accordingly, for any $x_{1},x_{2}\in(0.25,0.75)$, the characteristics from
each of these points will be given by

$t=\frac{1}{u_{0}(x_{1})}(x-x_{1});\quad t=\frac{1}{u_{0}(x_{2})}(x-x_{2}).$
(47)

Since $u_{0}(x)$ is a linearly decreasing function in the considered domain,
if $(x_{m},t_{m})$ represents the point where two such characteristics meet,
substituting (46) into (47) and eliminating $t_{m}$ gives

$-\frac{x_{m}-x_{1}}{4h_{0}(x_{1}-0.5)}=-\frac{x_{m}-x_{2}}{4h_{0}(x_{2}-0.5)}\quad\Rightarrow
x_{m}=0.5.$

Using (47), $t_{m}=\frac{1}{4h_{0}}$. Hence, for linearly decreasing functions
$u_{0}(x)$, we find that $(x_{m},t_{m})$ is independent of $x_{1}$ and $x_{2}$
and all such characteristics converge at the single point in space-time. For
$h_{0}=2$, $(x_{m},t_{m})=(0.5,0.125)$.

Beyond $t=t_{m}$, shock formation takes place as a consequence of
characteristics traveling in the opposite directions. In the regions $x<0.25$
and $x>0.75$, the rays travel parallel to time-axis. Furthermore, the jumps at
$x=0.25$ and $x=0.75$ turn into expansion fans as shown in Fig. 19(a). The
solution for any point $\\{(x,t):t>t_{m},x\neq 0.5\\}$ can be found based on
the characteristics coming from the expansion fans.

(a) Characteristics

(b) $u(x,t)$

Figure 19: In Fig. 19(a), slope of any characteristic in the $x-t$ plane is
given by $1/u(x,t)$. Three distinct regions are delineated by the red-colored
lines. The red coloured, vertical x-patterned line beyond $t_{m}$ at $x=0.5$
represents the shock.

We define the following regions:

$\displaystyle\text{I}:\\{(x,t):-2t+x-0.25\geq
0\\}\quad\cap\quad\\{(x,t):2t+x-0.75\leq 0\\};$
$\displaystyle\text{II}:\\{(x,t):\quad 0.25\leq
x<0.5\\}\quad\cap\quad\\{(x,t):-2t+x-0.25\leq 0\\};$
$\displaystyle\text{III}:\\{(x,t):\quad 0.5<x\leq
0.75\\}\quad\cap\quad\\{(x,t)2t+x-0.75\geq 0\\}.$

The complete $u$ profile can thus be given by

$u(x,t)=\begin{cases}\dfrac{8(x-0.5)}{8\,t-1}&\,\text{for Region I}\\\\[8.0pt]
\dfrac{x-0.25}{t}&\,\text{for Region II}\\\\[8.0pt]
\dfrac{x-0.75}{t}&\,\text{for Region III}\\\ 0&\,\text{Otherwise.
}\par\end{cases}$ (48)

The continuous $Y$ profile corresponding to $Y(0,t)=0$ can be obtained as

$Y(x,t)=\begin{cases}\begin{rcases}\dfrac{2x^{2}-x}{4t}+\dfrac{1}{32t}&\,\text{for
}\,0.25\leq x<0.25+2t\\\\[8.0pt]
\dfrac{4(x^{2}-x)}{8t-1}+\dfrac{8t+3}{4(8t-1)}&\text{for }\,0.25+h_{0}\,t\leq
x<0.75-h_{0}\,t\\\\[8.0pt] \dfrac{2x^{2}-3x}{4t}+\dfrac{9}{32t}&\text{for
}\,0.75-h_{0}\,t\leq x<0.75\\\\[8.0pt] 0&\text{Otherwise }\\\
\end{rcases}0<t<t_{m}\\\\[10.0pt] \\\
\begin{rcases}\dfrac{2x^{2}-x}{4t}+\dfrac{1}{32t}&\,\,\text{for }\,0.25\leq
x<0.5\\\\[8.0pt] \dfrac{2x^{2}-3x}{4t}+\dfrac{9}{32t}&\,\,\text{for }\,0.5\leq
x<0.75\\\\[8.0pt] 0&\,\text{Otherwise }\end{rcases}t\geq t_{m}.\end{cases}$
(49)

## References

  * [1] U. Kouskiya and A. Acharya, “Hidden convexity in the heat, linear transport, and euler’s rigid body equations: A computational approach.” _To appear in Quarterly of Applied Mathematics_ , https://arxiv.org/abs/2304.09418, 2023.
  * [2] Y. Brenier, “The initial value problem for the Euler equations of incompressible fluids viewed as a concave maximization problem,” _Communications in Mathematical Physics_ , vol. 364, no. 2, pp. 579–605, 2018.
  * [3] ——, “Examples of hidden convexity in nonlinear PDEs,” https://hal.science/hal-02928398/document, 2020.
  * [4] J. M. Burgers, “A mathematical model illustrating the theory of turbulence,” _Advances in Applied Mechanics_ , vol. 1, pp. 171–199, 1948.
  * [5] H. Bateman, “Some recent researches on the motion of fluids,” _Monthly Weather Review_ , vol. 43, no. 4, p. 163–170, 1915.
  * [6] A. Acharya, “A dual variational principle for nonlinear dislocation dynamics,” _Journal of Elasticity_ , https://doi.org/10.1007/s10659-023-09998-5, 2023.
  * [7] ——, “A hidden convexity in continuum mechanics, with application to classical, continuous-time, rate-(in) dependent plasticity,” _arXiv preprint arXiv:2310.03201_ , 2023.
  * [8] ——, “Variational principles for nonlinear PDE systems via duality,” _Quarterly of Applied Mathematics_ , vol. LXXXI, pp. 127–140, 2023, Article electronically published on September 26, 2022.
  * [9] P. D. Lax, _Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves_. Society for Industrial and Applied Mathematics, 1973.
  * [10] G. Strang, _Introduction to Applied Mathematics_. Wellesley Cambridge Press, 1977.
  * [11] T. P. Liu, “Hopf-Cole Transformation,” _Academia Sinica (New Series), Bulletin of the Institute of Mathematics_ , vol. 12, no. 7, pp. 71–101, 2017.
  * [12] R. T. Rockafellar, _Conjugate duality and optimization_. SIAM, 1974.