# A combination of Residual Distribution and the Active Flux formulations or a
new class of schemes that can combine several writings of the same hyperbolic
problem: application to the 1D Euler equations

R. Abgrall
Institute of Mathematics, University of Zürich
Winterhurerstrasse 190, CH 8057 Zürich
email<EMAIL_ADDRESS>
###### Abstract

We show how to combine in a natural way (i.e. without any test nor switch) the
conservative and non conservative formulations of an hyperbolic system that
has a conservative form. This is inspired from two different class of schemes:
the Residual Distribution one [1], and the Active Flux formulations [2, 3, 4,
5, 6]. The solution is globally continuous, and as in the active flux method,
described by a combination of point values and average values. Unlike the
”classical” active flux methods, the meaning of the point-wise and cell
average degrees of freedom is different, and hence follow different form of
PDEs: it is a conservative version of the cell average, and a possibly non
conservative one for the points. This new class of scheme is proved to satisfy
a Lax-Wendroff like theorem. We also develop a method to perform non linear
stability. We illustrate the behaviour on several benchmarks, some quite
challenging.

##  1 Introduction

The notion of conservation is essential in the numerical approximation of
hyperbolic systems of conservation: if it is violated, there is no chance, in
practice, to compute the right weak solution in the limit of mesh refinement.
This statement is known since the celebrated work of Lax and Wendroff [7], and
what happens when conservation is violated has been discussed by Hou and Le
Floch [8]. This conservation requirement imposes the use of the conservation
form of the system. However, in many practical situations, this is not really
the one one would like to deal with, since in addition to conservation
constraints, one also seeks for the preservation of additional features, like
contacts for fluid mechanics, or entropy decrease for shocks.

In this paper, we are interested in compressible fluid dynamics. Several
authors have already considered the problem of the correct discretisation of
the non conservative form of the system. In the purely Lagrangian framework,
when the system is described by the momentum equation and the Gibbs equality,
this has been done since decades: one can consider the seminal work of
Wilkins, to begin with, and the problem is still of interest: one can consider
[9, 10, 11] where high order is sought for. In the case of the Eulerian
formulation, there are less work. One can mention [12, 13, 14] where staggered
meshes are used, the thermodynamic variables are localised in the cells, while
the kinetic ones are localised at the grid points, or [15] where a non
conservative formulation with correction is used from scratch. The first two
references show how to construct at most second order scheme, while the last
one shows this for any order. All constructions are quite involved in term of
algebra, because one has to transfert information from the original grid and
the staggered one.

In this paper, we aim at showing how the notion of conservation introduced in
the residual distribution framework [16] is flexible enough to allow to deal
directly with the non conservative form of the system, while the correct
solutions are obtained in the limit of mesh refinement. More precisely, we
show how to to deal both with the conservative and non conservative form of
the PDE, without any switch, as it was the case in [17]. We illustrate our
strategy on several versions of the non conservative form, and provide first,
second order and third order accurate version of the scheme. More than a
particular example, we describe a general strategy which is quite simple. The
systems on which we will work are descriptions of the Euler equations for
fluid mechanics:

  * •

The conservation one:

$\dfrac{\partial}{\partial t}\begin{pmatrix}\rho\\\ \rho u\\\
E\end{pmatrix}+\dfrac{\partial}{\partial x}\begin{pmatrix}\rho u\\\ \rho
u^{2}+p\\\ u(E+p)\end{pmatrix}=0$ (1)

  * •

the primitive formulation:

$\dfrac{\partial}{\partial t}\begin{pmatrix}\rho\\\ u\\\
p\end{pmatrix}+\begin{pmatrix}\dfrac{\partial\rho u}{\partial x}\\\
u\dfrac{\partial u}{\partial x}+\frac{1}{\rho}\dfrac{\partial p}{\partial
x}\\\ u\dfrac{\partial p}{\partial x}+(e+p)\dfrac{\partial u}{\partial
x}\end{pmatrix}=0$ (2)

  * •

The ”entropy” formulation:

$\dfrac{\partial}{\partial t}\begin{pmatrix}p\\\ u\\\
s\end{pmatrix}+\begin{pmatrix}u\dfrac{\partial p}{\partial
x}+(e+p)\dfrac{\partial u}{\partial x}\\\ u\dfrac{\partial u}{\partial
x}+\frac{1}{\rho}\dfrac{\partial p}{\partial x}\\\ u\dfrac{\partial
s}{\partial x}\end{pmatrix}=0$ (3)

where as usual $\rho$ is the density, $u$ the velocity, $p$ the pressure,
$E=e+\tfrac{1}{2}\rho u^{2}$ is the total energy, $e=(\gamma-1)p$ and
$s=\log(p)-\gamma\log(\rho)$ is the entropy. The ratio of specific heats,
$\gamma$ is supposed to be constant here, mostly for simplicity.

This paper has several source of inspirations. The first one is the residual
distribution (RD) framework, and in particular [16]. The second one is the
family of active flux [2, 18, 3, 4, 5], where the solution is represented by a
cell average and point values. The conservation is recovered from how the
average is updated. Here the difference comes from the fact that in addition
several forms of the same system can be conserved, as (1), (2), (3) for the
point value update while a Lax Wendroff like result can still be shown. If the
same system were used, both for the cell average and the point values, this
would easily fit into the RD framework, using the structure of the polynomial
reconstruction. The difference with Active Flux is that we use only the
representation of the solution within one cell, and not a fancy flux
evaluation. Another difference is about the way the solution is evolved in
time: the AF method uses the method of characteristics to evolve the point
value, while here we rely on more standard Runge-Kutta methods.

The format of the paper is  as follows. In a first part, we explain the
general principles of our method, and justify why, under the assumptions made
on the numerical sequence for the Lax-Wendroff theorem (boundedness in
$L^{\infty}$ and strong convergence in a $L^{p}$, $p\geq 1$, of a subsequence
toward a $v\in L^{p}$, then this $v$ is a weak solution of the problem), we
can also show the convergence of a subsequence to a weak solution of the
problem, under the same assumptions. In the second part, we describe several
discretisations of the method, and in a third part we provide several
simulations to illustrate the method.

In this paper, the letter $C$ denotes a constant, and we uses the standard
”algebra”, for example $C\times C=C$, $C+C=C$, or $\alpha C=C$ for any
constant $\alpha\in\mathbb{R}$.

##  2 The method

###  2.1 Principle

We consider the problem

$\dfrac{\partial\mathbf{u}}{\partial
t}+\dfrac{\partial{\mathbf{f}}(\mathbf{u})}{\partial x}=0,\qquad
x\in\mathbb{R}$ (4a) with the initial condition
$\mathbf{u}(x,0)=\mathbf{u}_{0}(x),\qquad x\in\mathbb{R}$ (4b) Here
$\mathbf{u}\in\mathcal{D}_{\mathbf{u}}\subset\mathbb{R}^{p}$. For smooth
solutions, we also consider an equivalent formulation in the form
$\dfrac{\partial\mathbf{v}}{\partial t}+J\dfrac{\partial\mathbf{v}}{\partial
x}=0$ (4c)

where $\mathbf{v}=\Psi(\mathbf{u})\in\mathcal{D}_{\mathbf{v}}$ and
$\Psi:\mathcal{D}_{\mathbf{u}}\rightarrow\mathcal{D}_{\mathbf{v}}$ is assumed
to be one-to-one and $C^{1}$ (as well as the inverse function). For example,
if (4) corresponds to (1), then

$\mathcal{D}_{\mathbf{u}}=\\{(\rho,\rho u,E)\in\mathbb{R}^{3}\text{ such that
}\rho>0\text{ and }E-\frac{1}{2}\rho u^{2}>0\\}.$

If (4c) corresponds to (2), then

$\mathcal{D}_{\mathbf{v}}=\\{(\rho,u,p)\text{ such that }\rho>0\text{ and
}p>0\\}$

and (for a perfect gas) the mapping $\Psi$ corresponds to $(\rho,\rho
u,E)\mapsto\big{(}\rho,u,p=(\gamma-1)\big{(}E-\tfrac{1}{2}\rho
u^{2}\big{)}\big{)},$ while

$J=\begin{pmatrix}u&\rho&0\\\ 0&u&\frac{1}{\rho}\\\ 0&(e+p)&u\end{pmatrix}.$

For (3),

$\mathcal{D}_{\mathbf{u}}=\\{(p,u,s)\in\mathbb{R}^{3},p>0\\}.$

More generally, we have
$J=\big{[}\nabla_{\mathbf{u}}\big{(}\Psi^{-1}\big{)}\big{]}\nabla_{\mathbf{u}}{\mathbf{f}}$.

The idea is to discretise _simultaneously_ (4a) and (4c). Forgetting the
possible boundary conditions, $\mathbb{R}$ is divided into non overlapping
intervals $K_{j+1/2}=[x_{j},x_{j+1}]$ where $x_{j}<x_{j+1}$ for all
$j\in\mathbb{Z}$. We set $\Delta_{j+1/2}=x_{j+1}-x_{j}$ and
$\Delta=\max_{j}\Delta_{j+1/2}$. At the grid points, we will estimate
$\mathbf{v}_{j}$ in time, while in the cells we will estimate the average
value

$\bar{\mathbf{u}}_{j+1/2}=\frac{1}{\Delta_{j+1/2}}\int_{x_{j}}^{x_{j+1}}\mathbf{u}(x)\;dx$

When needed, we have $\mathbf{u}_{j}=\Psi^{-1}(\mathbf{v}_{j})$, however
$\bar{\mathbf{v}}_{j+1/2}=\Psi(\bar{\mathbf{u}}_{j+1/2})$ is meaningless since
the $\Psi$ does not commute with the average.

In $K_{j+1/2}$ any continuous function can be represented by
$\mathbf{u}_{j}=\mathbf{u}(x_{j})$, $\mathbf{u}_{j+1}=\mathbf{u}(x_{j+1})$ and
$\bar{\mathbf{u}}_{j+1/2}$: one can consider the polynomial $R_{\mathbf{u}}$
defined on $K_{j+1/2}$ by

$\big{(}R_{\mathbf{u}}\big{)}_{|K_{j+1/2}}(x)=\mathbf{u}_{j}L_{j+1/2}^{0}+\mathbf{u}_{j+1}L_{j+1/2}^{1}+\bar{\mathbf{u}}_{j+1/2}L_{j+1/2}^{1/2},$

with

$L_{j+1/2}^{\xi}(x)=\ell_{\xi}\big{(}\frac{x-x_{j}}{x_{j+1}-x_{j}}\big{)}$

and

$\ell_{0}(s)=(1-s)(1-3s),\quad\ell_{1}(s)=s(3s-2),\qquad\ell_{1/2}(x)=6s(1-s).$

We see that

$\begin{split}\ell_{0}(0)=1,&\quad\ell_{0}(1)=0,\quad\int_{0}^{1}\ell_{0}(s)ds=0\\\
\ell_{1}(1)=1,&\quad\ell_{1}(0)=0,\quad\int_{0}^{1}\ell_{1}(s)ds=0\\\
\ell_{1/2}(0)=0,&\quad\ell_{1/2}(1)=0,\quad\int_{0}^{1}\ell_{1/2}(s)\;ds=1.\end{split}$

How to evolve $\bar{\mathbf{u}}_{j+1/2}$ following (4a) and $v_{j}$ following
(4c) in time? The solution is simple for the average value: since

$\Delta_{j+1/2}\;\dfrac{d\bar{\mathbf{u}}_{j+1/2}}{dt}+{\mathbf{f}}(\mathbf{u}_{j+1}(t))-{\mathbf{f}}(\mathbf{u}_{j}(t))=0,$

we simply take

$\Delta_{j+1/2}\dfrac{d\bar{\mathbf{u}}_{j+1/2}}{dt}+\big{(}\hat{\mathbf{f}}_{j+1/2}-\hat{\mathbf{f}}_{j-1/2}\big{)}=0$
(5a) where $\hat{\mathbf{f}}_{j+1/2}$ is a consistent numerical flux that
depends continuously of its arguments. In practice, since the approximation is
continuous, we take
$\hat{\mathbf{f}}_{j+1/2}={\mathbf{f}}(\mathbf{u}_{j}){={\mathbf{f}}\big{(}\Psi^{-1}(\mathbf{v}_{j})\big{)}.}$
(5b) For $\mathbf{v}$, we assume a semi-discrete scheme of the following form:
$\dfrac{d\mathbf{v}_{j}}{dt}+\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}+\overrightarrow{\Phi}^{\mathbf{v}}_{j-1/2}=0$
(5c) such that
$\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}+\overrightarrow{\Phi}^{\mathbf{v}}_{j+1/2}$
is a consistent approximation of $J\dfrac{\partial\mathbf{v}}{\partial x}$ in
$K_{j+1/2}$. We will give examples later, for now we only describe the
principles. In general the residuals
$\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}$ and
$\overrightarrow{\Phi}^{\mathbf{v}}_{j-1/2}$ need to depend on some
$\mathbf{v}_{l}$ and $\mathbf{v}_{l+1/2}\approx\mathbf{v}(x_{l+1/2})$. We can
recover the missing informations at the half points in two steps:

  1. 1.

From $\mathbf{v}_{j}$, we can get $\mathbf{u}_{j}=\Psi(\mathbf{v}_{j})$,

  2. 2.

Then in $[x_{j},x_{j+1}]$ we approximate $\mathbf{u}$ by

$R_{\mathbf{u}}(x)=\mathbf{u}_{j}\ell_{0}\big{(}\frac{x-x_{j}}{\Delta_{j+1/2}}\big{)}+\bar{\mathbf{u}}_{j+1/2}\ell_{1/2}\big{(}\frac{x-x_{j}}{\Delta_{j+1/2}}\big{)}+\mathbf{u}_{j+1}\ell_{1}\big{(}\frac{x-x_{j}}{\Delta_{j+1/2}}\big{)},$

which enable to provide $\mathbf{u}_{j+1/2}:=R_{\mathbf{u}}(x_{j+1/2})$, i.e

$\mathbf{u}_{j+1/2}=\frac{3}{2}\bar{\mathbf{u}}_{j+1/2}-\frac{\mathbf{u}_{j}+\mathbf{u}_{j+1}}{4}.$
(5d)

Note that this relation is simply
$\bar{\mathbf{u}}_{j+1/2}=\tfrac{1}{6}\big{(}\mathbf{u}_{j}+\mathbf{u}_{j+1}+4\mathbf{u}_{j+1/2}\big{)},$
i.e. Simpson’s formula.

  3. 3.

Finally, we state

$\mathbf{v}_{j+1/2}=\Psi^{-1}(R_{\mathbf{u}}(x_{j+1/2}))$

In some situations, described later, we will also make the approximation:

$\mathbf{v}_{j+1/2}=\Psi^{-1}(\bar{\mathbf{u}}_{j+1/2})$

which is nevertheless consistent (but only first order accurate). As written
above, the fluctuations $\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}$ and
$\overrightarrow{\Phi}^{\mathbf{v}}_{j+1/2}$ are functionals of the form
$\Phi\big{(}\\{\mathbf{v}_{l},\mathbf{v}_{l+1/2}\\},j-p\leq l\leq j+p)$ for
some fixed value of $p$. We will make the following assumptions:

  1. 1.

Lipschiz continuity: There exists $C$ that depends only on $\mathbf{u}^{0}$
and $T$ such that for any $j\in\mathbb{Z}$

$\|\Phi\big{(}\\{\mathbf{v}_{l},\mathbf{v}_{l+1/2}\\},j-p\leq l\leq
j+p)\|\leq\frac{C}{\Delta_{j+1/2}}\bigg{(}\sum_{l=-p}^{p}\|\mathbf{v}_{l}-\mathbf{v}_{l+1/2}\|\bigg{)},$
(6a)

  2. 2.

Consistency. Setting $\mathbf{v}^{h}=R_{\mathbf{u}}$,

$\sum_{j\in\mathbb{Z}}\int_{K_{j+1/2}}\|\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}+\overrightarrow{\Phi}^{\mathbf{v}}_{j+1/2}-J\dfrac{\partial\mathbf{v}^{h}}{\partial
x}\|\;d\mathbf{x}\leq C\;\Delta.$ (6b)

  3. 3.

Regular mesh: the meshes are regular in the finite element sense.

The ODE systems (5) are integrated by a standard ODE solver. We will choose
the Euler forward method, and the second order and third order SSP Runge-Kutta
scheme.

###  2.2 Analysis of the method

In order to explain why the method can work, we will choose the simplest ODE
integrator, namely the Euler forward method. The general case can be done in
the same way, with more technical details. So we integrate (5) by:

$\bar{\mathbf{u}}_{j+1/2}^{n+1}=\bar{\mathbf{u}}_{j+1/2}^{n}-\frac{\Delta
t_{n}}{\Delta_{j+1/2}}\big{(}\underbrace{{\mathbf{f}}(\mathbf{u}_{j+1}^{n})-{\mathbf{f}}(\mathbf{u}_{j}^{n})}_{:=\delta_{j+1/2}{\mathbf{f}}}\big{)},$
(7)

and

$\mathbf{v}_{j}^{n+1}=\mathbf{v}_{j}^{n}-{\Delta
t_{n}}\big{(}\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}+\overrightarrow{\Phi}^{\mathbf{v}}_{j-1/2}\big{)}$
(8)

Setting $\Delta_{j}$ as the average of $\Delta_{j+1/2}$ and $\Delta_{j-1/2}$,
we rewrite (8) as

$\mathbf{v}_{j}^{n+1}=\mathbf{v}_{j}^{n}-\frac{\Delta
t_{n}}{\Delta_{j}}\delta_{x}\mathbf{v}_{j}$ (9)

and we note that, using the assumption (6a) as well as the fact that the mesh
is shape regular, that there exists $C>0$ depending only on $\mathbf{u}^{0}$
and $T$ such that

${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\|\delta_{x}\mathbf{v}_{j}}\|\leq
C\sum_{j=p-1}^{p+1}\|\mathbf{v}_{j}-\mathbf{v}_{j+1/2}\|.$

Using the transformation (5d), from (8), we can evaluate
$\mathbf{u}_{j}^{n+1}=\Psi(\mathbf{v}_{j}^{n+1})$, and then write the update
of $\mathbf{u}$ as

${\Delta_{j}\big{(}\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n}\big{)}+\Delta
t_{n}\delta_{x}\mathbf{u}_{j}=0}$ (10)

where

$\delta_{x}\mathbf{u}_{j}=\dfrac{\Delta_{j}}{\Delta
t_{n}}\bigg{(}\Psi(\mathbf{v}_{j}^{n}-\frac{\Delta
t_{n}}{\Delta_{j}}\delta\mathbf{v}_{j})-\Psi(\mathbf{v}_{j}^{n})\bigg{)},$

which, thanks to the assumptions we have made on $\Psi$ satisfies

$\|\delta_{x}\mathbf{u}_{j+1/2}\|\leq C\|\delta_{x}\mathbf{v}_{j}\|\leq
C\sum_{j=-p}^{p}\|\mathbf{v}_{j+l}-\mathbf{v}_{j+l+1}\|\leq
C\sum_{l=-p}^{p}\|\mathbf{u}_{j+l}-\mathbf{u}_{j+1+l}\|$

for some constants that depends on the gradient of $\Psi$ and the maximum of
the $\mathbf{v}_{i}^{n}$ for $i\in\mathbb{Z}$.

To explain the validity of the approximation, we start by the Simpson formula,
which is exact for quadratic polynomials:

$\int_{x_{j}}^{x_{j+1}}f(x)\;dx\approx\dfrac{\Delta_{j+1/2}}{6}\big{(}f(x_{j})+4f(x_{j+1/2})+f(x_{j+1})\big{)}.$

From the point values $\mathbf{u}_{j}$, $\mathbf{u}_{j+1}$ and
$\mathbf{u}_{j+1/2}$ at times $t_{n}$ and $t_{n+1}$, we define the quadratic
Lagrange interpolant $R_{\mathbf{u}^{n}}$ and $R_{\mathbf{u}^{n+1}}$ and then
write

$\int_{x_{j}}^{x_{j+1}}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\varphi(x,t)}\big{(}R_{\mathbf{u}^{n+1}}-R_{\mathbf{u}^{n}}\big{)}\;dx\approx\dfrac{\Delta_{j+1/2}}{6}\big{(}\varphi_{j+1}(\mathbf{u}_{j+1}^{n+1}-\mathbf{u}_{j+1}^{n})+4\varphi_{j+1/2}(\mathbf{u}_{j+1/2}^{n+1}-\mathbf{u}_{j+1/2}^{n})+\varphi_{j}(\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n})\bigg{)}.$

Accuracy is not an issue here. Using (8) and (10), setting
$\delta_{j}^{n+1/2}\mathbf{u}=\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n}$, we get

$\begin{split}\Sigma:=\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\dfrac{\Delta_{j+1/2}}{6}\bigg{(}&\varphi_{j+1}^{n}\delta_{j}^{n+1/2}\mathbf{u}+4\varphi_{j+1/2}^{n}\delta_{j+1/2}^{n+1/2}\mathbf{u}+\varphi_{j}^{n}\delta_{j}^{n+1/2}\mathbf{u}\bigg{)}\\\
&=\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\dfrac{\Delta_{j+1/2}}{6}\bigg{(}\varphi_{j+1}^{n}\delta_{j+1}^{n+1/2}\mathbf{u}\\\
&\qquad\qquad\qquad\qquad+4\varphi_{j+1/2}^{n}\big{(}\frac{3}{2}\delta_{j+1/2}^{n+1/2}\bar{\mathbf{u}}-\frac{\delta_{j+1}^{n+1/2}\mathbf{u}+\delta_{j}^{n+1/2}\mathbf{u}}{4}\big{)}+\varphi_{j}^{n}\delta_{j}^{n+1/2}\mathbf{u}\bigg{)}\\\
&=\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\Delta_{j+1/2}\varphi_{j+1/2}^{n}\delta_{j+1/2}^{n+1/2}\bar{\mathbf{u}}\\\
&\qquad\qquad\qquad\qquad+\underbrace{\sum\limits_{j\in\mathbb{Z}}\frac{\delta_{j}^{n+1/2}\mathbf{u}}{6}\bigg{\\{}\Delta_{j+1/2}\big{(}\varphi_{j}-\varphi_{j+1/2}\big{)}+\Delta_{j-1/2}\big{(}\varphi_{j}-\varphi_{j-1/2}\big{)}\bigg{\\}}}_{S_{n}}.\end{split}$

so that we get, using (10)

$\begin{split}\sum\limits_{n\in\mathbb{N}}\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\dfrac{\Delta_{j+1/2}}{6}\bigg{(}&\varphi_{j+1}^{n}\delta_{j}^{n+1/2}\mathbf{u}+4\varphi_{j+1/2}^{n}\delta_{j+1/2}^{n+1/2}\mathbf{u}+\varphi_{j}^{n}\delta_{j}^{n+1/2}\mathbf{u}\bigg{)}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}-}\sum\limits_{n\in\mathbb{N}}\Delta
t_{n}\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\varphi_{j+1/2}^{n}\delta_{j+1/2}{\mathbf{f}}\\\
&-\sum\limits_{n\in\mathbb{N}}S_{n}=0\end{split}$ (11)

Then we use again (10), use the fact that the mesh is regular, and observe
that

$S_{n}=\Delta
t_{n}\sum_{j}\Delta_{j}\;\delta_{x}\mathbf{u}_{j}+O(\Delta^{3}).$

In appendix A, we will show that in the limit, the contribution of the $S_{n}$
term will converges towards $0$, while the first term of (11) will converge to

$\int_{0}^{+\infty}\int_{\mathbb{R}}\dfrac{\partial\varphi}{\partial
t}\mathbf{u}\;dxdt-\int_{\mathbb{R}}\mathbf{u}_{0}\;dx$

while the second term will converge towards

$\int_{0}^{+\infty}\int_{\mathbb{R}}\dfrac{\partial\varphi}{\partial
x}{\mathbf{f}}(\mathbf{u})\;d\mathbf{x}.$

This will be shown, using classical arguments, in the appendix A, so that we
have

######  Proposition 2.1.

We assume that the mesh is regular: there exists $\alpha$ and $\beta$ such
that $\alpha\leq\Delta_{j+1/2}/\Delta_{j-1/2}\leq\beta$. If
$\max\limits_{j\in\mathbb{Z}}\|\mathbf{u}_{j}^{n}\|_{\infty}$ and
$\max\limits_{j\in\mathbb{Z}}\|v_{j+1/2}^{n}\|_{\infty}$ are bounded, and if a
subsequence of $\mathbf{u}_{\Delta}$ converges in $L^{1}$ towards
$\mathbf{u}$, then $\mathbf{u}$ is a weak solution of the problem.

######  Remark 2.2.

Indeed, the definition of a precise $\Delta_{j}$ is not really needed, and we
come back to this in the next section. What is needed is a spatial scale that
relates the updates in $\mathbf{v}$ and $\mathbf{u}$ in an incremental form of
the finite difference type. This is why the asssumption of mesh regularity is
fundemental.

##  3 Some examples of discretisation

We list possible choices: for $\dfrac{\partial\mathbf{v}}{\partial
t}+J\dfrac{\partial\mathbf{v}}{\partial x}=0,$ where $J$ is the Jacobian of
${\mathbf{f}}$ with respect to $\mathbf{u}$; they have been used in the
numerical tests. The question here is to define $\overleftarrow{\Phi}_{j+1/2}$
and $\overrightarrow{\Phi}_{j+1/2}$ that are the contributions of $K_{j\pm
1/2}$ to $J\dfrac{\partial\mathbf{v}}{\partial x}$ so that

$J\dfrac{\partial\mathbf{v}}{\partial
x}(x_{i})\approx\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}+{\overrightarrow{\Phi}^{\mathbf{v}}_{j-1/2}}.$

We follow the work of Iserle [19] who gives all the possible schemes that
guaranty a stable (in $L^{2}$) semi-discretisation of the convection equation,
for a regular grid which we assume. The only difference in his notations and
ours is that the grid on which are defined the approximation of the derivative
is made of the mesh points $x_{j}$ and the half points $x_{j+1/2}$.

The first list of examples have an upwind flavour:

$\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}=\big{(}{J(\mathbf{v}_{j})}\big{)}^{-}\frac{\delta^{-}_{j}\mathbf{v}}{\Delta_{j+1/2}/2}\text{
and
}\overrightarrow{\Phi}^{\mathbf{v}}_{j+1/2}=\big{(}J(\mathbf{v}_{j+1})\big{)}^{+}\frac{\delta_{j+1}^{+}\mathbf{v}}{\Delta_{j+1/2}/2}$
(12)

where $\delta^{\pm}_{j}$ is an approximation of $\Delta_{j+1/2}\dfrac{\partial
v}{\partial x}$ obtained from [19]111The author works on $\dfrac{\partial
u}{\partial t}=-\dfrac{\partial u}{\partial x}$ which is a bit confusing
w.r.t. to ”modern” habits. It is true that British drive left.:

  * •

First order approximation: we take

$\delta_{j}^{+}\mathbf{v}=\mathbf{v}_{j}-\mathbf{v}_{j-1/2},\qquad\delta_{j}^{-}\mathbf{v}=\mathbf{v}_{j+1/2}-\mathbf{v}_{j}.$
(13)

  * •

Second order: we take

$\begin{split}\delta_{j}^{-}\mathbf{v}&=-\frac{3}{2}\mathbf{v}_{j}+2\mathbf{v}_{j+1/2}-\frac{\mathbf{v}_{j+1}}{2}\\\
\delta_{j}^{+}\mathbf{v}&=\frac{\mathbf{v}_{j-1}}{2}-2\mathbf{v}_{j-1/2}+\frac{3}{2}\mathbf{v}_{j}\end{split}$
(14)

  * •

Third order: We take

$\begin{split}{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\delta_{j}^{-}=-\frac{v_{i+1}}{6}+v_{i+1/2}-\frac{v_{i}}{2}-\frac{v_{i-1/2}}{3}},\\\
{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\delta_{j}^{+}=\frac{v_{i-1}}{6}-v_{i-1/2}+\frac{v_{i}}{2}+\frac{v_{i+1/2}}{3}}\end{split}$
(15)

  * •

Fourth order: The fully centered scheme would be

$\delta_{j}^{\pm}\mathbf{v}={\dfrac{\mathbf{v}_{j+1}-\mathbf{v}_{j-1}}{12}+2\dfrac{\mathbf{v}_{j+1/2}-\mathbf{v}_{j-1/2}}{3}}$

but we prefer

$\begin{split}\delta_{j}^{-}\mathbf{v}&=\frac{\mathbf{v}_{j-1/2}}{4}+\frac{5}{6}\mathbf{v}_{j}-\frac{3}{2}\mathbf{v}_{j+1/2}+\frac{\mathbf{v}_{j+1}}{2}-\frac{\mathbf{v}_{j+3/2}}{12}\\\
\delta_{j}^{+}\mathbf{v}&=\frac{\mathbf{v}_{j+1/2}}{4}+\frac{5}{6}\mathbf{v}_{j}-\frac{3}{2}\mathbf{v}_{j-1/2}+\frac{1}{2}\mathbf{v}_{j-1}-\frac{\mathbf{v}_{j-3/2}}{12}\end{split}$
(16)

  * •

Etc…

It can be useful to have more dissipative versions of a first order scheme. We
take:

$\bigg{(}J\dfrac{\partial v}{\partial
x}\bigg{)}_{j}=\overleftarrow{\Phi}_{j+1/2}+\overrightarrow{\Phi}_{j-1/2}$

with

$\begin{split}\frac{\Delta_{j+1/2}}{2}\overleftarrow{\Phi}_{j+1/2}&=\frac{1}{2}\widehat{J\dfrac{\partial\mathbf{v}}{\partial
x}}_{j}+\alpha\big{(}\mathbf{v}_{j}-\frac{\mathbf{v}_{j}+\mathbf{v}_{j+1/2}}{2}\big{)}\\\
\frac{\Delta_{j+1/2}}{2}\overrightarrow{\Phi}_{j+1/2}&=\frac{1}{2}\widehat{J\dfrac{\partial\mathbf{v}}{\partial
x}}_{j+1}+\alpha\big{(}\mathbf{v}_{j+1}-\frac{\mathbf{v}_{j+1}+\mathbf{v}_{j+1/2}}{2}\big{)}\end{split}$

where $\widehat{J\dfrac{\partial v}{\partial x}}_{l}$ is a consistent
approximation of $J\dfrac{\partial u}{\partial x}$ at $x_{l}$ and $\alpha$ is
an upper-bound of the spectral radius of $J(\mathbf{v}_{j})$,
$J(\mathbf{v}_{j+1/2})$ and $J(\mathbf{v}_{j+1})$. We take, for simplicity,
$\mathbf{v}_{j+1/2}=\Psi^{-1}(\bar{\mathbf{u}}_{j+1/2})$. For the model (2),
we take

$\frac{\Delta_{j+1/2}}{2}\widehat{J\dfrac{\partial v}{\partial
x}}_{j}=\begin{pmatrix}(\rho u)_{j+1/2}-(\rho u)_{j}\\\
\frac{1}{2}\big{(}u_{j+1/2}^{2}-u_{j}^{2})+\frac{1}{\tilde{\rho}_{j+1/2}}\big{(}p_{j+1/2}-p_{j}\big{)}\\\
\tilde{u}_{j+1/2}\big{(}p_{j+1/2}-p_{j})+\tilde{\rho
c^{2}}\big{(}u_{j+1/2}-u_{j}\big{)}\end{pmatrix}$

where $\tilde{\rho}_{j+1/2}$ is the geometric average of $\rho_{j}$ and
$\rho_{j+1/2}$, $\tilde{u}_{j+1/2}$ is the arithmetic average of $u_{j}$ and
$u_{j+1/2}$, while $\widetilde{\rho
c^{2}}_{j+1/2}=\gamma\frac{p_{j}+p_{j+1/2}}{2}$. For the model (3), we take:

$\frac{\Delta_{j+1/2}}{2}\widehat{J\dfrac{\partial v}{\partial
x}}_{j}=\begin{pmatrix}\tilde{u}_{j+1/2}\big{(}s_{j+1/2}-s_{j}\big{)}\\\
\frac{1}{2}\big{(}u_{j+1/2}^{2}-u_{j}^{2})+\frac{1}{\tilde{\rho}_{j+1/2}}\big{(}p_{j+1/2}-p_{j}\big{)}\\\
\tilde{u}_{j+1/2}\big{(}p_{j+1/2}-p_{j})+\widetilde{\rho
c^{2}}\big{(}u_{j+1/2}-u_{j}\big{)}\end{pmatrix}.$

All this has a Local Lax-Friedrichs’ flavour, and seems to be positivity
preserving for the velocity and the pressure.

Using this, the method is :

$\begin{split}\dfrac{d\mathbf{v}_{j}}{dt}+\overleftarrow{\Phi}_{j+1/2}^{\mathbf{v}}+\overrightarrow{\Phi}_{j-1/2}^{\mathbf{v}}=0\end{split}$
(17a) combined with $\Delta
x\dfrac{d\bar{\mathbf{u}}_{j+1/2}}{dt}+{{\mathbf{f}}(\mathbf{u}_{j+1})-{\mathbf{f}}(\mathbf{u}_{j})}=0.$
(17b)

We see in (17a) that the time derivative of $\mathbf{v}$ is obtained by adding
two fluctuations, one computed for the interval $K_{j+1/2}=[x_{j},x_{j+1}]$
and one for the interval $K_{j-1/2}=[x_{j-1},x_{j}]$. These fluctuations are
obtained from (12) with the increments in $\mathbf{v}$ defined by (13), (15),
(16), etc. In the sequel, we denote the scheme applied on the interval
$K_{j+1/2}$ by $S_{j+1/2}(k)$ where the average are integrated by (17b) and
$\mathbf{v}$ by (17a) with the fluctuations (13) for $k=1$, (14) for $k=2$ and
(15) for $k=3$, etc. To make sure that the first order scheme is positivity
preserving (at least experimentaly), we may also consider the case denoted by
$k=0$ where $S_{j+1/2}(0)$ is the local Lax Friedrichs scheme defined above.
Both fluctuation (13) and the local Lax Friedrichs scheme are first order
accurate, but the second one is quite dissipative but positivity preserving
while the scheme (13) is not (experimentaly) positivity preserving. The system
(17) is integrated in time by a Runge-Kutta solver: RK1, RK SSP2 and RK SSP3.

###  3.1 Error analysis in the scalar case

Here, the mesh is uniform, so that $\Delta_{j+1/2}=\Delta$ for any
$j\in\mathbb{Z}$. It is easy to check the consistency, and on figure 1 we show
the $L^{1}$ error on $u$ and $\bar{u}$ for (17) with SSPKR2 and SSPRK3
(CFL=$0.4$) for a convection problem

$\dfrac{\partial u}{\partial t}+\dfrac{\partial u}{\partial x}=0$

with periodic boundary conditions and the initial condition $u_{0}=\cos(2\pi
x)$.

Figure 1:  Error plot for $u$ and $\bar{u}$ for (17) with SSPKR2 and SSPRK3
(cfl=$0.4$). Here $f(u)=u$. The second order results are obtained with SSPRK2
with (17a)-(17b), the third order results is obtained by (17a)-(17b).

######  Remark 3.1 (Linear stability).

In the appendix B, we perform the $L^{2}$ linear stability and we get, with
$\lambda=\tfrac{\Delta t}{\Delta}$,

  * •

First order scheme, $|\lambda|\leq 0.92$,

  * •

Second order scheme, $|\lambda|\leq 0.6$,

  * •

Third order scheme, $|\lambda|\leq 0.5$.

We also have run this scheme for the Burgers equation, and compared it with a
standard finite volume (with local Lax-Friedrichs). The conservative form of
the PDE is used for the average, and the non conservative one for the point
values: $J=u$ and $\psi(u)=u$. This is an experimental check of conservation.
The initial condition is

$\mathbf{u}_{0}(x)=\sin(2\pi x)+\frac{1}{2}$

on $[0,1]$, so that there is a moving shock.

Figure 2:  Solution of Burgers with $10\,000$ points, $t_{fin}=0.4$, $CFL=0.4$
for the second order (a and b) ((17a)-(17b) with SSPRK2) and third order (c
and d) ((17a)-(17b) with SSPRK3). The global solution is represented in (a)
and (c), and a zoom around the discontinuity is shown in (b) and (c).

We can see that the agreement is excellent and that the numerical solution
behaves as expected.

###  3.2 Non linear stability

As such, the scheme is at most linearly stable, with a CFL condition based on
the fine grid. However, in case of discontinuities or the occurence of
gradients that are not resolved by the grid, we have to face oscillations, as
usual.

In order to get high order oscillation free results, a natural option would be
to extend the MUSCL approach to the present context. However, it is not very
clear how to proceed, so we have relied on the MOOD paradigm [20, 21]. The
idea is to work with several schemes ranging from order $p$ to $1$, with the
lowest order one able to provide results with positive density and pressure.
These schemes are the $S_{j+1/2}(k),k=1,\ldots,3$ scheme defined above. They
are assumed to work for a given CFL range, and the algorithm is as follows:
For each Runge-Kutta sub-step, starting from
$U^{n}=\\{\bar{\mathbf{u}}_{j+1/2}^{n},\bar{\mathbf{v}}_{j}^{n}\\}_{j\in\mathbb{Z}}$,
we compute

$\begin{split}\tilde{\bar{\mathbf{u}}}_{j+1/2}^{n+1}&=\bar{\mathbf{u}}_{j+1/2}^{n}-\lambda_{n}\big{(}{\mathbf{f}}(\mathbf{u}_{j+1}^{n})-{\mathbf{f}}(\mathbf{u}_{j}^{n})\big{)},\qquad\lambda_{n}=\frac{\Delta
t_{n}}{\Delta_{j+1/2}}\\\
\tilde{\mathbf{v}}_{j}^{n+1}&=\tilde{\mathbf{v}}_{j}^{n}-2\Delta
t_{n}\overleftarrow{\Phi}^{\mathbf{v}}_{j+1/2}\\\
\tilde{\mathbf{v}}_{j+1}^{n+1}&=\tilde{\mathbf{v}}_{j+1}^{n}-2\Delta
t_{n}\overrightarrow{\Phi}^{\mathbf{v}}_{j+1/2}\end{split}$ (18)

by the scheme $S_{j+1/2}(p)$. Then we test the validity of these results in
the interval $[x_{j},x_{j+1}]$ for the density (and possibly the pressure).
This is described a little bit later. The variable $\mathbf{v}$ is updated as
in (18), because at $t_{n+1}$, the true update of $\mathbf{v}_{j}$ is the half
sum of $\tilde{\mathbf{v}}_{j}^{n+1}$ and $\tilde{\mathbf{v}}_{j+1}^{n+1}$.

If the test is positive, then we keep the scheme $S_{j+1/2}(p)$ in that
interval, else we start again with $S_{j+1/2}(p-1)$, and repeat the procedure
unless all the intervals $K_{j+1/2}$ have successfully passed the test. This
is described in Algorithm 1 where $\mathcal{S}_{j+1/2}$ is the stencil used in
$K_{j+1/2}$.

Algorithm 1 Description of the MOOD loop. The algorithm stops because
$\SS_{j+1/2}=0$ corresponds to the local Lax Friedrichs scheme for which the
test is always true.

0: $U^{n}=\\{\bar{u}_{j+1/2}^{n},\bar{v}_{j}^{n}\\}_{j\in\mathbb{Z}}$

0: Allocate $\\{\SS_{j+1/2}\\}_{j\in\mathbb{Z}}$ an array of integers. It is
initialized with $\SS_{j+1/2}=S_{j+1/2}(p)$, the maximum order.

for $k=p,\ldots,2$ do

for all For all $K_{j+1/2}$ do

Define $\tilde{\bar{\mathbf{u}}}_{j+1/2}^{n+1}$,
$\tilde{\mathbf{v}}_{j}^{n+1}$ and $\tilde{\mathbf{v}}_{j+1}^{n+1}$ as in (18)

Apply the test on $\tilde{\bar{\mathbf{u}}}_{j+1/2}^{n+1}$,
$\tilde{\mathbf{v}}_{j}^{n+1}$ and $\tilde{\mathbf{v}}_{j+1}^{n+1}$ :

if test=.true. then

$\SS_{j+1/2}=S_{j+1/2}(k-1)$

end if

end for

end for

Now, we describe the tests. We do, in the following order, for each element
$K_{j+1/2}$, at the iteration $k>0$ of the loop of 1: the tests are performed
on variables evaluated from $\mathbf{u}$ and $\mathbf{v}$. For the scalar
case, they are simply the point values at $x_{j},x_{j+1/2}$ and $x_{j+1}$. For
the Euler equations they are the density, and possibly the pressure

  1. 1.

We check if all the variables are numbers (i.e. not NaN). If not, we state
$\SS_{j+1/2}=S_{j+1/2}(k-1)$,

  2. 2.

(Only for the Euler equations) We check if the density is positive. We can
also request to check if the pressure is also positive. If the variable is
negative, the we state that $\SS_{j+1/2}=S_{j+1/2}(k-1)$.

  3. 3.

Then we check if at $t_{n}$, the solution was not constant in the numerical
stencils of the degrees of freedom in $K_{j+1}$, this in order to avoid to
detect a fake maximum principle. We follow the procedure of [21]. if we
observe that the solution was locally constant, the $\SS_{j+1/2}$ is not
modified.

  4. 4.

Then we apply a discrete maximum principle, even for systems though it is not
very rigorous. For the variable $\xi$ (in practice the density, and we may
request to do the same on the pressure), we compute $\min_{j+1/2}\xi$ (resp.
$\max_{j+1/2}\xi$) the minimum (resp. maximum) of the values of $\xi$ on
$K_{j+1/2}$, $K_{j-1/2}$ and $K_{j+3/2}$. We say we have a potential maximum
if
$\tilde{\xi}^{n+1}\not\in[\min_{j+1/2}\xi^{n}+\varepsilon_{j+1/2},\max_{j+1/2}\xi^{n}-\varepsilon_{j+1/2}]$,
with $\epsilon_{j+1/2}$ estimated as in [20]. Then:

     * •

If
$\tilde{\xi}^{n+1}\in[\min_{j+1/2}\xi^{n}+\varepsilon_{j+1/2},\max_{j+1/2}\xi^{n}-\varepsilon_{j+1/2}]$,
$\SS_{j+1/2}$ is not modified

     * •

Else we use the following procedure introduced in [21]. In each $K_{l+1/2}$,
we can evaluate a quadratic polynomial $p_{l+1/2}$ that interpolates $\xi$.
Note that its derivative is linear in $\xi$. We compute

$p^{\prime}_{j-1/2}(x_{j}),p^{\prime}_{j+3/2}(x_{j+1}),p^{\prime}_{j+1/2}(x_{j})\text{
and }p^{\prime}_{j+1/2}(x_{j+1}).$

       * –

If

$p^{\prime}_{j+1/2}(x_{j})\in[\min(p^{\prime}_{j-1/2}(x_{j}),p^{\prime}_{j+3/2}(x_{j+1})]\text{
and
}p^{\prime}_{j+1/2}(x_{j+1})\in[\min(p^{\prime}_{j-1/2}(x_{j}),p^{\prime}_{j+3/2}(x_{j+1})]$

we say it is a true regular extrema and $\SS_{j+1/2}$ will not be modified,

       * –

Else the extrema is declared not to be regular, and
$\SS_{j+1/2}=S_{j+1/2}(k-1)$

As a first application, to show that the oscillations are well controlled
without sacrificing the accuracy, we consider the advection problem (with
constant speed unity) on $[0,1]$, periodic boundary conditions with initial
condition:

$u_{0}(x)=\left\\{\begin{array}[]{ll}0&\text{ if }y\in[-1,-0.8[\\\
\frac{1}{6}\big{(}G(y,\beta,z-\delta)+G(y,\beta,z+\delta)+4G(y,\beta,z)&\text{
if }y\in[-0.8,-0.6]\\\ 1&\text{ if }y\in[-0.4,-0.2]\\\ 1-|10y-1|&\text{ if
}y\in[0,0.2]\\\
\frac{1}{6}\big{(}F(y,\beta,z-\delta)+G(y,\beta,z+\delta)+4F(y,\beta,z)&\text{
else,}\end{array}\right.\text{ with }y=2x-1$

Here $a=0.5$, $z=-0.7$, $\delta=0.005$, $\alpha=10$,

$\beta=\dfrac{\log 2}{36\delta^{2}}$

and

$G(t,\beta,z)=\exp\big{(}-\beta(t-z)^{2}\big{)},\qquad
F(t,a,\alpha)=\sqrt{\max\big{(}0,1-\alpha(t-a)^{2}\big{)}}.$

Using the MOOD procedure with the third order scheme, the results obtained for
$300$ points for $T=10$ are displayed in figure 3. They look very reasonable.

Figure 3:  Shu Jiang problem, CFL=0.4, third order scheme with MOOD, 300
points, periodic conditions, 10 periods. The point values and cell average are
almost undistinguishable.

##  4 Numerical results for the Euler equations

In this section, we show the flexibility of the approach, where conservation
is recovered only by the equation (17a), and so lots of flexibility is
possible with the relations on the $\mathbf{u}_{i}$. To illustrate this, we
consider the Euler equations. We will consider the conservative formulation
(1) for the average value, so $\mathbf{u}=(\rho,\rho u,E)^{T}$ and either the
form (2), i.e. $\mathbf{v}=(\rho,u,p)$ or the form (3) with
$\mathbf{v}=(p,u,s)^{T}$.

###  4.1 Sod test case

The Sod case is defined for $[0,1]$, the initial condition is

$(\rho,u,p)^{T}=\left\\{\begin{array}[]{ll}(1,0,1)^{T}&\text{ for }x<0.5\\\
(0.125,0,0.1)^{T}\text{ else.}\end{array}\right.$

The final time is $T=0.16$. The problem is solved with (1)-(2) and displayed
in figures 4, 5, 6 and 7, while the solution obtained with the combination
(1)-(3) is shown on figure 8 and (9). When the MOOD procedure is on, it is
applied with $\rho$ and $p$ and all the test are performed.

Figure 4:  100 grid points, and the second order SSPRK2 scheme with CFL=0.1.
(a): density, (b): velocity, (c): pressure.

Figure 5:  100 grid points, and the second order SSPRK2 scheme with CFL=0.1.
(a): density, (b): velocity, (c): pressure. MOOD test made on $\rho$ and $p$

The exact solution is also shown every time. Different order in time/space are
tested. The results are good, eventhough the MOOD procedure is not perfect.

Figure 6:  100 grid points, and the third order SSPRK3 scheme with CFL=0.1.
(a): density, (b): velocity, (c): pressure.

Figure 7:  100 grid points, and the third order SSPRK3 scheme with CFL=0.1.
(a): density, (b): velocity, (c): pressure. Mood Test made on $\rho$ and $p$

This the use of the combination (1)-(3) seems more challenging, we have
performed a convergence study (with 10000 points). This is shown on figure
(9), and a zoom around the contact discontinuity is also shown.

(a) $\rho$

(b) $p$

(c) $u$

(d) $s$

Figure 8:  Solution with the variables (s,u,p) for 100 points, comparison with
the exact solution, third order in time/space with Mood and non Mood. Mood is
done on $\rho$ and $p$. Cfl=0.2

(a) $\rho$ (b) $p$

(c) $u$ (d) $p$ zoom

(e) $u$ zoom

Figure 9:  Solution with the variables (s,u,p) for 10000 points, comparison
with the exact solution. Cfl=0.1, no mood The zoomed figures are for
$x\in[0.6,0.7]$ and the ticks are for $10^{-7}$. We plot $u$ and $p$ across
the contact

We can observe a numerical convergence to the exact one in all cases. In the
appendix C we show some results on irregular meshes, with the same
conclusions.

###  4.2 A smooth case

We consider a fluid with $\gamma=3$: the characteristics are straight lines.
The initial condition is inspired from Toro: in $[-1,1]$,

$\begin{split}\rho_{0}(x)&=1+\alpha\sin(2\pi x)\\\ u_{0}(x)&=0\\\
p_{0}(x0&=\rho_{0}(x)^{\gamma}\end{split}$ (19)

The classical case is is for $\alpha=0.999995$ where vaccum is almost reached.
Here, since we do not want to test the robustness of the method, we take
$\alpha=\frac{3}{4}$. The final time is set to $T=0.1$.

The exact density and velocity in this case can be obtained by the method of
characteristics and is explicitly given by

$\rho(x,t)=\dfrac{1}{2}\big{(}\rho_{0}(x_{1})+\rho_{0}(x_{2})\big{)},\quad
u(x,t)=\sqrt{3}\big{(}\rho(x,t)-\rho_{0}(x_{1})\big{)},$

where for each coordinate $x$ and time $t$ the values $x_{1}$ and $x_{2}$ are
solutions of the nonlinear equations

$\displaystyle x+\sqrt{3}\rho_{0}(x_{1})t-x_{1}=0,$ $\displaystyle
x-\sqrt{3}\rho_{0}(x_{2})t-x_{2}=0.$

A example of numerical solution, superimposed with the exact one, is shown on
figure 10. It is obtained with the third order (time and space) scheme, and
here we have used the model $(\rho,u,p)$. The CFL number is set to $0.2$.

(a) $\rho$

(b) $u$

(c) $p$

Figure 10:  Solution (numerical and exact) for the conditions (19). The number
of grid points is set to $80$, with periodic boundary conditions.

The errors are shown in table 1.

$h=1/N$ | $L^{1}$ |  | $L^{2}$ |  | $L^{\infty}$ |
---|---|---|---|---|---|---
20 | $2.136\;10^{-4}$ | $-$ | $2.968\;10^{-4}$ | $-$ | $6.596\;10^{-4}$ | $-$
40 | $1.912\;10^{-5}$ | $-3.48$ | $2.702\;10^{-5}$ | $-3.45$ | $5.750\;10^{-5}$ | $-3.52$
80 | $1.398\;10^{-6}$ | $-3.77$ | $2.138\;10^{-6}$ | $-3.65$ | $4.673\;10^{-6}$ | $-3.62$
160 | $1.934\;10^{-7}$ | $-2.85$ | $2.595\;10^{-7}$ | $-3.04$ | $5.753\;10^{-7}$ | $-3.02$
320 | $3.641\;10^{-8}$ | $-2.40$ | $5.523\;10^{-8}$ | $-2.23$ | $1.276\;10^{-7}$ | $-2.17$
Table 1:  $L^{1}$, $L^{2}$ and $L^{\infty}$ error for the initial conditions
(19) with the third order scheme.

The errors, computed in $[-1,1]$ are in reasonable agreement with the $-3$
expected slopes. We also have done the same test with the non linear
stabilisation procedure described in section 3.2. _Exactly_ the same errors
are obtained: the order reduction test are never activated.

###  4.3 Shu-Osher case

The initial condition are:

$(\rho,u,p)=\left\\{\begin{array}[]{ll}(3.857143,2.629369,10.3333333)&\text{if
}x<-4\\\ (1+0.2\sin(5x),0,1)&\text{else}\end{array}\right.$

on the domain $[-5,5]$ until $T=1.8$. We have used the combination (1)-(2),
since the other one seems less robust. The density is compared to a reference
solution (obtained with a standard finite volume scheme with $20\;000$ points,
and the solution obtained with the third order scheme with $CFL=0.3$ and
$200$, $400$, $800$ and $1600$ points. The mood procedure uses the first order
upwind scheme if a PAD, a NaN or DMP is detected, the other cases use the 3rd
order scheme.The solutions are displayed in 11. With little resolution, the
results are very close of the reference one.

Figure 11:  (a): Solution of the Shu Osher problem , (b): zoom of the solution
around the shock.

For figure 11, the second order scheme is used as a rescue scheme.

###  4.4 Le Blanc case

The initial conditions are

$(\rho,u,e)=\left\\{\begin{array}[]{ll}(1,0,0.1)&\text{ if }x\in[-3,3]\\\
(0.001,0.10^{-7})&\text{ if }x\in[3,6]\end{array}\right.$

where $e=(\gamma-1)p$ and $\gamma=\tfrac{5}{3}$. The final time is $t=6$. This
is a very strong shock tube. The combination (1)-(2). It is not possible to
run higher that first order without the MOOD procedure. We show the second and
third order results are shown on figure 12, and zooms around the shocks and
the fan are showed in 13.

(a) $\rho$ (b) $\rho$

(c) $p$ (d) $p$

(e) $u$ (f) $u$

Figure 12:  Le Blanc test case, $CFL=0.1$, from 400 to 800 points. Left
column: MOOD test on $\rho$ and $p$, second order, right column: MOOD test on
$\rho$ and $p$, third order

Figure 13:  Le Blanc test case, zooms, comparison on the pressure between
second order and third order with 400 points

At time $t=6$ the shock wave should be at $x=8$: in addition to the extreme
conditions, it is generally difficult to get a a correct position of the shock
wave; this is why a convergence study is shown in figure 14. It is performed
with 400, 800, 10000 grid points, and the third order SSPRK3 scheme with
CFL=0.1. It is compared to the exact solution, and the results are good, see
for example [22] for a comparison with other methods, or [23] for a comparison
with Lagrangian methods.

Figure 14:  Convergence study on the density for the LeBlanc test case.

##  5 Conclusion

This study is preliminary and should be seen as a proof of concept. We show
how to combine, without any test, several formulations of the same problem,
one conservative and the other ones in non conservative form, in order to
compute the solution of hyperbolic systems. The emphasis is mostly put on the
Euler equations.

We explain why the formulation leads to a method that satisfies a Lax-Wendroff
like theorem. We also propose a way to provide non linearity stability, this
method works well but is not yet completely satisfactory.

Besides the theoretical results, we also show numerically that we get the
convergence to the correct weak solution. This is done on standard benchmark
problems, some being very challenging.

We intend to extend the method to several space dimensions, and improve the
limiting strategy. Different systems, such as the shallow water system, will
also be considered.

## Acknowledgements.

This work was done while the author was partially funded by SNF project
200020$\\_$175784\. The support of Inria via the International Chair of the
author at Inria Bordeaux-Sud Ouest is also acknowledged. Discussions with Dr.
Wasilij Barsukow are acknowledged, as well as the encouragements of Anne
Burbeau (CEA DAM, France). Last, I would like to thank, warmly, the two
anonymous referees: their critical comments have led to big improvements.

## References

  * [1] R. Abgrall.  The notion of conservation for residual distribution schemes (or fluctuation splitting schemes), with some applications.  Commun. Appl. Math. Comput., 2(3):341–368, 2020.
  * [2] T.A. Eyman and P.L. Roe.  Active flux.  49th AIAA Aerospace Science Meeting, 2011.
  * [3] T.A. Eyman.  Active flux.  PhD thesis, University of Michigan, 2013.
  * [4] C. Helzel, D. Kerkmann, and L. Scandurra.  A new ADER method inspired by the active flux method.  Journal of Scientific Computing, 80(3):35–61, 2019.
  * [5] W. Barsukow.  The active flux scheme for nonlinear problems.  J. Sci. Comput., 86(1):Paper No. 3, 34, 2021.
  * [6] P.L. Roe.  Is discontinuous reconstruction really a good idea?  Journal of Scientific Computing, 73:1094–1114, 2017.
  * [7] P. Lax and B. Wendroff.  Systems of conservation laws.  Comm. Pure Appl. Math., 13:381–394, 1960.
  * [8] T. Y. Hou and P. G. Le Floch.  Why nonconservative schemes converge to wrong solutions : error analysis.  Math. Comp., 62(206):497–530, 1994.
  * [9] V. A. Dobrev, T. Kolev, and R. N. Rieben.  High-order curvilinear finite element methods for Lagrangian hydrodynamics.  SIAM J. Sci. Comput., 34(5):B606–B641, 2012.
  * [10] R. Abgrall and S. Tokareva.  Staggered grid residual distribution scheme for Lagrangian hydrodynamics.  SIAM J. Sci. Comput., 39(5):A2317–A2344, 2017.
  * [11] R. Abgrall, K. Lipnikov, N. Morgan, and S. Tokareva.  Multidimensional staggered grid residual distribution scheme for Lagrangian hydrodynamics.  SIAM J. Sci. Comput., 42(1):A343–A370, 2020.
  * [12] R. Herbin, J.-C. Latché, and T.T. Nguyen.  Consistent segregated staggered schemes with explicit steps for the isentropic and full Euler equations.  ESAIM Math. Model. Numer. Anal., 52(3):893–944, 2018.
  * [13] G. Dakin, B. Després, and S. Jaouen.  High-order staggered schemes for compressible hydrodynamics. Weak consistency and numerical validation.  J. Comput. Phys., 376:339–364, 2019.
  * [14] R. Abgrall and K. Ivanova.  High order schemes for compressible flow problems with staggered grids.  in preparation, 2021.
  * [15] R. Abgrall, P. Bacigaluppi, and S. Tokareva.  High-order residual distribution scheme for the time-dependent Euler equations of fluid dynamics.  Comput. Math. Appl., 78(2):274–297, 2019.
  * [16] R. Abgrall.  Some remarks about conservation for residual distribution schemes.  Comput. Methods Appl. Math., 18(3):327–351, 2018.
  * [17] S. Karni.  Multicomponent flow calculations by a consistent primitive algorithm.  J. Comput. Phys., 112(1):31–43, 1994.
  * [18] T.A. Eyman and P.L. Roe.  Active flux for systems.  20 th AIAA Computationa Fluid Dynamics Conference, 2011.
  * [19] A. Iserles.  Order stars and saturation theorem for first-order hyperbolics.  IMA J. Numer. Anal., 2:49–61, 1982.
  * [20] S. Clain, S. Diot, and R. Loubère.  A high-order finite volume method for systems of conservation laws—Multi-dimensional Optimal Order Detection (MOOD).  J. Comput. Phys., 230(10):4028–4050, 2011.
  * [21] F. Vilar.  A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction.  J. Comput. Phys., 387:245–279, 2019.
  * [22] R. Ramani, J. Reisner, and S. Shkoller.  A space-time smooth artificial viscosity method with wavelet noise indicator and shock collision scheme, part 2: The 2-d case.  Journal of Computational Physics, 387:45–80, Jun 2019.
  * [23] R. Loubère.  Validation test case suite for compressible hydrodynamics computation.  http://loubere.free.fr/images/test_suite.PDF, 2005.
  * [24] E. Godlewski and P.-A. Raviart.  Hyperbolic systems of conservation laws, volume 3/4 of Mathématiques & Applications (Paris) [Mathematics and Applications].  Ellipses, Paris, 1991.

##  Appendix A Proof of proposition 2.1.

We show the proposition 2.1 in the scalar case, the system case is identical.

We start by some notations: $\mathbb{R}$ is subdivided into intervals
$K_{j+1/2}=[x_{j},x_{j+1}]$ with $x_{j}<x_{j+1}$, and $h$ will be the maximum
of the length of the $K_{j+1/2}$. On each interval, from the point values
$u_{i}$ and $u_{i+1}$, as well as the average $\bar{u}_{j+1/2}$ we can
construct a quadratic polynomial . From this and as above, we can construct a
globally continuous piecewise quadratic function, that for simplicity of
notations we will denote by $R_{u_{\Delta}}$.

Let $T>0$ and a time discretisation $0<t_{1}<\ldots<t_{n}<\ldots<t_{N}\leq T$
of $[0,T]$. We define $\Delta t_{n}=t_{n+1}-t_{n}$ and $\Delta
t=\max\limits_{n}\Delta t_{n}$. We are given the sequences
$\\{u_{j}^{p}\\}_{j\in\mathbb{Z}}^{p=0\ldots N}$ and
$\\{{\bar{u}}_{j+1/2}^{n}\\}_{j\in\mathbb{Z}}^{p=0,\ldots,N}$. We can define a
function $u_{\Delta}$ by:

$\text{ if }(x,t)\in[x_{j},x_{j+1}]\times[t_{n},t_{n+1}[,\text{ then
}u_{\Delta}(x,t)=R_{u^{n}_{\Delta}}(x).$

The set of these functions is denoted by $X_{\Delta}$ and is equipped with the
$L^{\infty}$ and $L^{2}$ norms.

We have the following lemma

######  Lemma A.1.

Let $T>0$ , $\\{t_{n}\\}_{n=0,\ldots,N}$ an increasing subdivision of $[0,T]$,
$[a,b]$ a compact of $\mathbb{R}$. Let $(u_{\Delta})_{h}$ a sequence of
functions of $X_{\Delta}$ defined on $\mathbb{R}\times\mathbb{R}^{+}$. We
assume that there exists $C\in\mathbb{R}$ independent of $\Delta$ and $\Delta
t$, and $\mathbf{u}\in L^{2}_{loc}([a,b]\times[0,T])$ such that

$\sup\limits_{\Delta}\sup\limits_{x,t}|u_{\Delta}(x,t)|\leq C\quad\text{ and
}\lim\limits_{\Delta,\Delta t\rightarrow
0}|u_{\Delta}-u|_{L^{2}([a,b]\times[0,T])}=0.$

Then

$\lim\limits_{\Delta,\Delta t\rightarrow 0}\sum_{n=0}^{N}\Delta
t_{n}\bigg{[}\sum_{j\in\mathbb{Z}}\Delta_{j+1/2}\bigg{(}|u_{j}^{n}-\bar{u}_{j+1/2}^{n}|+|u_{j+1}^{n}-\bar{u}_{j+1/2}^{n}|+|u_{j}^{n}-u_{j+1}^{n}|\bigg{)}\bigg{]}=0.$
(20)

###### Proof.

First, because the vector space of polynomials of degree 3 on
$[x_{j},x_{j+1}]$ is finite dimensional and with a dimension independent of
$j$, there exists $C_{1}$ and $C_{2}$ such that

$\begin{split}C_{1}\Delta_{j+1/2}\bigg{(}|u_{j}^{n}-\bar{u}_{j+1/2}^{n}|+|u_{j+1}^{n}-\bar{u}_{j+1/2}^{n}|\bigg{)}\leq\int_{x_{j}}^{x_{j+1}}&|u_{\Delta}(x,t_{n})-{\bar{u}}_{j+1/2}^{n}|\;dx\\\
&\leq
C_{2}\Delta_{j+1/2}\bigg{(}|u_{j}^{n}-\bar{u}_{j+1/2}^{n}|+|u_{j+1}^{n}-\bar{u}_{j+1/2}^{n}|\bigg{)}\end{split}$

so that

$\sum_{n=0}^{K}\Delta
t_{n}\sum\limits_{j,K_{j+1/2}\subset[a,b]}\Delta_{j+1/2}\bigg{(}|u_{j}^{n}-\bar{u}_{j+1/2}^{n}|+|u_{j+1}^{n}-\bar{u}_{j+1/2}^{n}|\bigg{)}\leq
C_{1}^{-1}\int_{0}^{T}\int_{a}^{b}|u_{\Delta}-\bar{u}_{\Delta}|\;dx,$

where for simplicity we have introduced $\bar{u}_{\Delta}$ the function
defined by:

$\text{ if
}(x,t)\in[x_{j},x_{j+1}[\times[t_{n},t_{n+1}[,\bar{u}_{\Delta}(x,t)={\bar{u}}_{j+1/2}^{n}.$

Then we rely on classical arguments of functional analysis: since
$(u_{\Delta})$ is bounded, and since $L^{\infty}([a,b]\times[0,T])\subset
L^{1}([a,b]\times[0,T])$, there exists $u^{\prime}\in
L^{\infty}([a,b]\times[0,T]$ such that $u_{\Delta}\rightarrow u^{\prime}$ in
the weak-$\star$ topology. Similarly, there exists $\bar{u}\in
L^{\infty}([a,b]\times[0,T])$ such that ${\bar{u}}_{\Delta}\rightarrow\bar{u}$
for the weak-$\star$ topology.

Since $u_{\Delta}\rightarrow u$ in $L^{2}_{loc}$, we have $u^{\prime}=u$
because $[a,b]\times[0,T]$ is bounded and $C_{0}^{\infty}([a,b]\times[0,T])$
is dense in $L^{1}([a,b]\times[0,T])$. Let us show that $\bar{u}=u$. let
$\varphi\in C_{0}^{\infty}(\mathbb{R}\times\mathbb{R}^{+})$. We have, setting

$\bar{\varphi}_{j+1/2}^{n}=\dfrac{1}{\Delta_{j+1/2}\Delta
t_{n}}\int_{t_{n}}^{t_{n+1}}\int_{x_{j}}^{x_{j+1}}\varphi(x,t)\;dxdt,$
$\begin{split}\int_{0}^{T}\int_{a}^{b}\big{(}\bar{u}_{\Delta}-u_{\Delta}\big{)}\varphi\;dxdt&=\sum_{n}\sum\limits_{a\leq
x_{j}<x_{j+1}\leq
b}\int_{t_{n}}^{t_{n+1}}\int_{x_{j}}^{x_{j+1}}\big{(}\bar{u}_{\Delta}-u_{\Delta}\big{)}\varphi\;dxdt\\\
&=\sum_{n}\sum\limits_{a\leq x_{j}<x_{j+1}\leq
b}\Bigg{(}\int_{t_{n}}^{t_{n+1}}\int_{x_{j}}^{x_{j+1}}\big{(}\bar{u}_{\Delta}-u_{\Delta}\big{)}\varphi\;dxdt-\int_{t_{n}}^{t_{n+1}}\int_{x_{j}}^{x_{j+1}}\big{(}\bar{u}_{\Delta}-u_{\Delta}\big{)}\bar{\varphi}_{j+1/2}^{n}\;dxdt\Bigg{)}\\\
&=\sum_{n}\sum\limits_{a\leq x_{j}<x_{j+1}\leq
b}\int_{t_{n}}^{t_{n+1}}\int_{x_{j}}^{x_{j+1}}\big{(}\bar{u}_{\Delta}-u_{\Delta}\big{)}\big{(}\varphi-\bar{\varphi}_{j+1/2}^{n}\big{)}\;dxdt\end{split}$

using the fact that for any $[x_{j},x_{j+1}]\times[t_{n},t_{n+1}]$, we have
$\int_{t_{n}}^{t_{n+1}}\int_{x_{j}}^{x_{j+1}}\big{(}\bar{u}_{\Delta}-u_{\Delta}\big{)}\;dx\;dt=0$.

Since $\varphi\in C_{0}^{\infty}(\mathbb{R}\times\mathbb{R}^{+})$, there
exists $C$ that depends only on
$\|\dfrac{d\varphi}{dx}\|_{L^{\infty}(\mathbb{R}\times\mathbb{R}+)}$ such that

$\bigg{|}\int_{x_{j}}^{x_{j+1}}\big{(}\bar{u}_{\Delta}-u_{\Delta}\big{)}\big{(}\varphi-\bar{\varphi}\big{)}\;dxdt\bigg{|}\leq\Delta
t\Delta_{j+1/2}\;\Delta\max\limits_{j\in\mathbb{Z},n\leq
N}\big{(}|u_{j}^{n}|,|{\bar{u}}_{j+1/2}^{n}|\big{)}\leq\Delta
t\Delta^{2}\max\limits_{j\in\mathbb{Z},n\leq
N}\big{(}|u_{j}^{n}|,|{\bar{u}}_{j+1/2}^{n}|\big{)}$

and then,

$\bigg{|}\int_{0}^{T}\int_{a}^{b}\big{(}{\bar{u}}_{\Delta}-u_{\Delta}\big{)}\varphi\;dxdt\bigg{|}\leq
C\;\Delta$

and passing to the limit, $\bar{u}=u^{\prime}$. Since a subsequence of
$u_{\Delta}$ converges to $u$ in $L^{2}$, we have $\bar{u}=u^{\prime}=u$.

The same method shows that $(u_{\Delta}^{2})$ and $({\bar{u}}_{\Delta}^{2})$
have the same weak-$\star$ limit. Let us show it is $u^{2}$. Since
$C_{0}^{\infty}([a,b]\times[0,T]$ is dense in $L^{1}([a,b]\times[0,T]$), and
since $u_{\Delta}^{2}$ is bounded independently of $\Delta$ and $\Delta t$, we
can choose functions $\phi$ in $C_{0}^{\infty}([a,b]\times[0,T])$. This test
function is bounded in $[a,b]\times[0,T]$ and then, we have, at least for a
subsequence,

$\int_{a}^{b}\int_{0}^{T}|u-u_{\Delta}|^{2}\;dxdt\rightarrow 0,$

and then

$\int_{a}^{b}\int_{0}^{T}u_{\Delta}^{2}\;dxdt-2\int_{a}^{b}\int_{0}^{T}u_{\Delta}\,u\;dxdt+\int_{a}^{b}\int_{0}^{T}u^{2}\;dxdt\rightarrow
0.$

By the Cauchy-Schwarz inequality, $u\phi\in L^{1}([a,b]\times[0,T])$: the
second term tends towards

$\int_{a}^{b}\int_{0}^{T}u^{2}\;dxdt,$

so that

$\int_{a}^{b}\int_{0}^{T}u_{\Delta}^{2}\;dxdt-\int_{a}^{b}\int_{0}^{T}u^{2}\;dxdt\rightarrow
0.$

and $u_{\Delta}^{2}\rightarrow u^{2}$ in $L^{\infty}$ weak-$\star$.

Last, again by the same argument for $\phi=1$, since
$u_{\Delta}^{2}\rightarrow u^{2}$ in $L^{\infty}$ weak-$\star$, we get

$\int_{a}^{b}\int_{0}^{T}|\bar{u}_{\Delta}-u|^{2}\;dxdt\rightarrow 0,$

and finally

$\int_{a}^{b}\int_{0}^{T}|\bar{u}_{\Delta}-u_{\Delta}|^{2}\;dxdt\rightarrow
0.$

Since $[a,b]\times[0,T]$ is bounded, $L^{1}([a,b]\times[0,T])\subset
L^{2}([a,b]\times[0,T])$, we obtain

$\lim\limits_{\Delta,\Delta t\rightarrow 0}\sum_{n=0}^{N}\Delta
t_{n}\bigg{[}\sum_{j\in\mathbb{Z}}\Delta_{j+1/2}\bigg{(}|u_{j}^{n}-\bar{u}_{j+1/2}^{n}|+|u_{j+1}^{n}-\bar{u}_{j+1/2}^{n}|\bigg{)}\bigg{]}=0$

From this we get (20) because

$|u_{j}^{n}-u_{j+1}^{n}|\leq|u_{j}^{n}-\bar{u}_{j+1/2}^{n}|+|u_{j+1}^{n}-\bar{u}_{j+1/2}^{n}|.$

∎

Then we can proof proposition 2.1. We proceed the proof in several lemma.

######  Lemma A.2.

Under the conditions of proposition 2.1, for any $\varphi\in
C_{0}^{\infty}(\mathbb{R}\times\mathbb{R}^{+})$ we have

$\begin{split}\lim\limits_{\Delta t\rightarrow 0,\Delta\rightarrow
0}\sum\limits_{n=0}^{\infty}\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}&\dfrac{\Delta_{j+1/2}}{6}\bigg{(}\varphi_{j+1}^{n}(\mathbf{u}_{j+1}^{n+1}-\mathbf{u}_{j+1}^{n})+4\varphi_{j+1/2}^{n}(\mathbf{u}_{j+1/2}^{n+1}-\mathbf{u}_{j+1/2}^{n})+\varphi_{j}^{n}(\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n})\bigg{)}\\\
&=-\int_{\mathbb{R}\times\mathbb{R}^{+}}\dfrac{\partial\varphi}{\partial
t}u\;dxdt+\int_{\mathbb{R}}\varphi(x,0)u_{0}\;dxdt.\end{split}$

###### Proof.

This is a simple adaptation of the classical proof, see for example [24]. We
have, using that

$\delta u_{j+1/2}=\frac{3}{2}\overline{\delta u}_{j+1/2}-\frac{\delta
u_{j}+\delta u_{j+1}}{4}$

and the compactness of the support of $\varphi$,

$\begin{split}\sum\limits_{n=0}^{\infty}\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}&\dfrac{\Delta_{j+1/2}}{6}\bigg{(}\varphi_{j+1}^{n}(\mathbf{u}_{j+1}^{n+1}-\mathbf{u}_{j+1}^{n})+4\varphi_{j+1/2}^{n}(\mathbf{u}_{j+1/2}^{n+1}-\mathbf{u}_{j+1/2}^{n})+\varphi_{j}^{n}(\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n})\bigg{)}\\\
&=\underbrace{\sum\limits_{n=0}^{\infty}\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\Delta_{j+1/2}\varphi_{j+1/2}\overline{\delta
u}_{j+1/2}}_{(I)}\\\
&\qquad+\underbrace{\sum\limits_{n=0}^{\infty}\sum\limits_{j\in\mathbb{Z}}\bigg{(}\frac{\Delta_{j+1/2}}{6}\big{(}\varphi_{j}-\varphi_{j+1/2}\big{)}+\frac{\Delta_{j-1/2}}{6}\big{(}\varphi_{j}-\varphi_{j-1/2}\big{)}\bigg{)}\delta
u_{j}}_{(II)}\end{split}$

where
$\overline{\delta\mathbf{u}}_{j+1/2}=\bar{\mathbf{u}}_{j+1/2}^{n+1}-\bar{\mathbf{u}}_{j+1/2}^{n}$
and $\delta\mathbf{u}_{j}=\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n}$. The first
part, $(I)$, is the classical term, and under the condition of the lemma,
converges to

$-\int_{\mathbb{R}\times\mathbb{R}^{+}}\dfrac{\partial\varphi}{\partial
t}u\;dxdt+\int_{\mathbb{R}}\varphi(x,0)u_{0}\;dxdt.$

Since $\varphi\in C_{0}^{\infty}(\mathbb{R}\times\mathbb{R}^{+})$, there
exists $C$ that depends only on the $L^{\infty}$ norm of the first derivative
of $\varphi$ such that the term $(II)$ can be bounded by

$\begin{split}\bigg{|}\sum\limits_{n=0}^{N}\Delta
t_{n}\sum\limits_{j\in\mathbb{Z}}\bigg{(}\frac{\Delta_{j+1/2}}{6}\big{(}\varphi_{j}-\varphi_{j+1/2}\big{)}&+\frac{\Delta_{j-1/2}}{6}\big{(}\varphi_{j}-\varphi_{j-1/2}\big{)}\bigg{)}\delta
u_{j}\bigg{|}\leq C\sum\limits_{n=0}^{N}\Delta
t\;\Delta\sum\limits_{j\in\mathbb{Z}}\Delta_{j+1/2}|\delta u_{j}|\\\
&\quad\leq
CT(b-a)\Delta\max\limits_{j,p\in\mathbb{N}}|\mathbf{u}_{j}^{p}|.\end{split}$

This tends to zero because
$\max\limits_{j\in\mathbb{Z},p\in\mathbb{N}}|\mathbf{u}_{j}^{p}|$ is finite. ∎

######  Lemma A.3.

Under the assumptions of proposition 2.1,

$\lim\limits_{\Delta t,\Delta\rightarrow
0}\sum\limits_{n=0}^{\infty}\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\Delta
t_{n}\varphi_{j+1/2}\delta_{j+1/2}{\mathbf{f}}=-\int_{\mathbb{R}\times\mathbb{R}^{+}}\dfrac{\partial\varphi}{\partial
x}f(u)\;dxdt.$

###### Proof.

This is again a simple adaptation of the classical proof since
$\delta_{j+1/2}f=f(u_{j+1})-f(u_{j})$. We have

$\sum\limits_{[x_{j},x_{j+1}],j\in\mathbb{Z}}\varphi_{j+1/2}\delta_{j+1/2}{\mathbf{f}}=-\sum\limits_{j\in\mathbb{Z}}{\mathbf{f}}(\mathbf{u}_{j})\big{(}\varphi_{j+1/2}-\varphi_{j-1/2}\big{)},$

Then using the boundedness of the solution and Lebesgue dominated convergence
theorem, we get the result. ∎

Then we have

######  Lemma A.4.

Under the conditions of proposition 2.1, we have

$\lim\limits_{\Delta t,\Delta\rightarrow
0}\bigg{(}\sum\limits_{n\in\mathbb{N}}\sum\limits_{j\in\mathbb{Z}}\delta_{j}^{n+1/2}\mathbf{u}\bigg{\\{}\Delta_{j+1/2}\big{(}\varphi_{j}-\varphi_{j+1/2}\big{)}+\Delta_{j-1/2}\big{(}\varphi_{j}-\varphi_{j-1/2}\big{)}\bigg{\\}}\Bigg{)}=0$

###### Proof.

Since $\varphi\in C_{0}^{\infty}(\mathbb{R}\times\mathbb{R}^{+})$, there
exists $C$ that depends only on the first derivative of $\varphi$ such that

$|\Delta_{j+1/2}\big{(}\varphi_{j}-\varphi_{j+1/2}\big{)}+\Delta_{j-1/2}\big{(}\varphi_{j}-\varphi_{j-1/2}\big{)}|\leq
C\Delta\;(\Delta_{j+1/2}+\Delta_{j-1/2}).$

Then using (10), we see that

$\delta_{j}^{n+1/2}:=\mathbf{u}_{j}^{n+1}-\mathbf{u}_{j}^{n}=-\dfrac{\Delta
t_{n}}{\Delta_{j}}\delta_{x}\mathbf{u}_{j},$

so that

$\begin{split}\|\sum\limits_{j\in\mathbb{Z}}\delta_{j}^{n+1/2}\mathbf{u}\bigg{\\{}\Delta_{j+1/2}\big{(}\varphi_{j}-\varphi_{j+1/2}\big{)}&+\Delta_{j-1/2}\big{(}\varphi_{j}-\varphi_{j-1/2}\big{)}\bigg{\\}}\Bigg{)}\|\leq
C\Delta
t_{n}\sum\limits_{j\in\mathbb{Z}}\|\delta_{x}^{n+1/2}\mathbf{u}\|\Delta_{j}\\\
&=C\Delta\Delta
t_{n}\sum\limits_{j\in\mathbb{Z}}\|\delta_{x}\mathbf{u}_{j}\|=C\Delta\sum\limits_{j\in
Z}\sum\limits_{l=-p}^{l=p}|\mathbf{u}_{j+l}-\bar{\mathbf{u}}_{j+l+1/2}|\end{split}$

using the Lipschitz continuity of the fluctuations and the regularity of the
transformation $\mathbf{v}\mapsto\mathbf{u}$ together with the boundedness of
the solution. Then, from lemma A.1, we see that the results holds true. ∎

Then ends the proof of proposition 2.1.

##  Appendix B Linear stability analysis

The scheme writes, setting $\lambda=a\tfrac{\Delta t}{\Delta x}$ and assuming
$a>0$,

$\begin{split}u_{j}^{n+1}&=u^{n}_{j}-2\lambda\delta_{j}u^{n}\\\
\bar{u}_{j+1/2}^{n+1}&=\bar{u}_{j+1/2}^{n}-\lambda\big{(}u_{j+1}^{n}-u_{j}^{n}\big{)}\end{split}$

on with periodicity $1$. We set $\Delta x=\tfrac{1}{N}$ It is more convenient
to work with point values only, and we will use the form:

$\begin{split}u_{j}^{n+1}&=u^{n}_{j}-2\lambda\delta_{j}u^{n}\\\
u_{j+1/2}^{n+1}&=u_{j+1/2}^{n}-\lambda\bigg{(}\frac{3}{2}\big{(}u_{j+1}^{n}-u_{j}^{n}\big{)}-\frac{1}{4}\big{(}\delta_{j}u^{n}+\delta_{j+1}u^{n}\big{)}\bigg{)}.\end{split}$

We perform a linear stability analysis by Fourier analysis. What is not
completely standard is that the grid points do not play the same role. For
ease of notations, we double the indices, this avoids to have half integer in
the Fourier analysis. In other points, the quantities $u_{j}$ associated to
the grid points $x_{j}$ are denoted by $u_{2j}$: this will be the even terms.
Those associated to the intervals $[x_{j},x_{j+1}]$, i.e. $\bar{u}_{j+1/2}$
and $u_{j+1/2}$ will be denoted as $\bar{u}_{2j+1}$ and $u_{2j+1}$, they are
the odd terms, so that we use

$\begin{split}u_{2j}^{n+1}&=u^{n}_{j}-2\lambda\delta_{2j}u^{n}\\\
u_{2j+1}^{n+1}&=u_{2j+1}^{n}-\lambda\bigg{(}\frac{3}{2}\big{(}u_{2j+2}^{n}-u_{2j}^{n}\big{)}-\frac{1}{4}\big{(}\delta_{2j}u^{n}+\delta_{2j+2}u^{n}\big{)}\bigg{)}.\end{split}$
(21)

We have the Parseval equality,

$\frac{1}{2N}\sum\limits_{0}^{2N-1}u_{j}^{2}=\sum_{k=0}^{2N-1}|\hat{u}(k)|^{2},$

with

$\hat{u}(k)=\frac{1}{2N}\sum\limits_{j=0}^{2N-1}u_{j}e^{2i\pi\frac{kj}{2N}}=\frac{1}{2}\big{(}\hat{u}_{o}(k)+\hat{u}_{e}(k)\big{)}$

with, setting $\omega=e^{\frac{i\pi}{N}}$

$\hat{u}_{o}(k)=\frac{1}{N}\sum\limits_{j=0}^{N-1}u_{2j+1}\omega^{(2j+1)k},\quad\hat{u}_{e}(k)=\frac{1}{N}\sum\limits_{j=0}^{N-1}u_{2j}\omega^{(2j)k}$

The usual shift operator $[S(u)]_{j}=u_{j+1}$ gives:

$\hat{S(u)}_{o}=\omega^{-k}\hat{u}_{e},\quad\hat{S(u)}_{e}=\omega^{-k}\hat{u}_{o}.$

Using this, we see that the Euler forward method (21) gives

$\begin{pmatrix}\hat{u}_{o}^{n+1}\\\
\hat{u}_{e}^{n+1}\end{pmatrix}=\bigg{(}\text{Id}-\lambda
H\bigg{)}\begin{pmatrix}\hat{u}_{o}^{n}\\\ \hat{u}_{e}^{n}\end{pmatrix}$ (22)

with

$H_{1}(k)=\begin{pmatrix}\frac{1+\omega^{2k}}{4}&\frac{5\omega^{-k}+7\omega^{k}}{4}\\\
0&2\big{(}1-\omega^{k}\big{)}\end{pmatrix}$

for the first order in space scheme,

$H_{2}(k)=\begin{pmatrix}\frac{1+\omega^{2k}}{2}&\frac{9}{8}\omega^{-k}-2\omega^{-k}-\frac{\omega^{3k}}{8}\\\
0&2\big{(}1-\omega^{k}\big{)}\end{pmatrix}$

for the second order scheme and

$H_{3}(k)=\begin{pmatrix}\frac{\omega^{-2k}}{12}+\frac{1}{24}+\frac{\omega^{2k}}{12}-\frac{\omega^{4k}}{24}&\frac{31}{24}\omega^{-k}-\frac{3}{2}\omega^{k}+\frac{\omega^{2k}}{12}+\frac{5}{24}\omega^{3k}\\\
0&2\big{(}1-\omega^{k}\big{)}\end{pmatrix}$

for the third order in time space.

Combined with the RK time stepping we end up with an update of the form

$\begin{pmatrix}\hat{u}_{o}^{n+1}\\\
\hat{u}_{e}^{n+1}\end{pmatrix}=G_{k}\begin{pmatrix}\hat{u}_{o}^{n}\\\
\hat{u}_{e}^{n}\end{pmatrix}$

and writing

$G_{k}=\begin{pmatrix}\alpha_{k}&\beta_{k}\\\
\gamma_{k}&\delta_{k}\end{pmatrix}$

we end up with

$\begin{split}\hat{u}_{o}^{n+1}&=\alpha_{k}\hat{u}_{o}^{n}+\beta_{k}\hat{u}_{e}^{n}\\\
\hat{u}_{e}^{n+1}&=\gamma_{k}\hat{u}_{o}^{n}+\delta_{k}\hat{u}_{e}^{n}\\\
\end{split}$

so that

$\hat{u}^{n+1}=\frac{\alpha_{k}+\gamma_{k}}{2}\hat{u}_{o}^{n}+\frac{\beta_{k}+\delta_{k}}{2}\hat{u}_{e}^{n},$

from which we get

$|\hat{u}^{n+1}|^{2}=\frac{1}{4}\overline{\begin{pmatrix}\hat{u}_{o}^{n}&\hat{u}_{e}^{n}\end{pmatrix}}M_{k}\begin{pmatrix}\hat{u}_{o}^{n}\\\
\hat{u}_{e}^{n}\end{pmatrix}$

with

$M_{k}=\frac{1}{4}\begin{pmatrix}|\alpha_{k}+\gamma_{k}|^{2}&\big{(}\alpha_{k}+\gamma_{k}\big{)}\overline{\big{(}\beta_{k}+\delta_{k}\big{)}}\\\
\overline{\big{(}\alpha_{k}+\gamma_{k}\big{)}}{\big{(}\beta_{k}+\delta_{k}\big{)}}&|\beta_{k}+\delta_{k}|^{2}\end{pmatrix}$

We have stability if the spectral radius of these matrices is always $\leq 1$,
and we immediately see that

$\rho(M_{k})=\frac{1}{4}\bigg{(}|\alpha_{k}+\gamma_{k}|^{2}+|\beta_{k}+\delta_{k}|^{2}\bigg{)}.$

After calculations, we see that the stability limits are:

  * •

First order scheme, $|\lambda|\leq 0.92$,

  * •

Second order scheme, $|\lambda|\leq 0.6$,

  * •

Third order scheme, $|\lambda|\leq 0.5$.

##  Appendix C Some numerical results on irregular meshes

In order to support the theoretical analysis of the method, we have applied it
on irregular meshes. The goal is to show that even here, one gets convergence
of the solution to a weak solution that appears to be the right one. Since we
use the same schemes, there is no hope to get anything but first order
accuracy. Accuracy on irregular meshes will be the topic of future work. The
mesh is defined by: for $0\leq i\leq N$,

$y_{0}=0,\quad y_{i+1}=y_{i}+\Delta y,\quad\Delta
y=\dfrac{1+\epsilon_{i}/2}{N}$

and

$\epsilon_{0}=-1,\epsilon_{i+1}=-\epsilon_{i}.$

Then we define the actual mesh by

$x_{i}=\dfrac{y_{i}}{y_{N}}.$

On the Sod problem, with $N=10000$, we get the results of the figure 15.

(a) $\rho$ (b) $p$ (c) $u$ (d) $s$

Figure 15:  Plot of the density, pressure, velocity and the pressure. This is
obtained with the ”third order” with MOOD test on the density and the
pressure.

On figure 16, we show a zoom of the density around the shock wave. The
discretisation points as well as the numerical and exact solution are shown.

Figure 16:  Zoom of the density around the shock wave.


# An elementary proof of Acerbi Fusco minimizer existence theorem

Tomáš Roskovec  and  Filip Soudský  T. G. Roskovec: Faculty of Economics,
University of South Bohemia, Studentská 13, České Budějovice, Czech Republic;
<EMAIL_ADDRESS>F. Soudský: Faculty of Science, Humanities and Education,
Technical University of Liberec, Studentská 1402/2, Liberec, Czech Republic;
<EMAIL_ADDRESS>
###### Abstract.

The weak lower semicontinuity of the functional

$F(u)=\int_{\Omega}f(x,u,\nabla u)\operatorname{d}x$

is a classical topic that was studied thoroughly. It was shown that if the
function $f$ is continuous and convex in the last variable, the functional is
sequentially weakly lower semicontinuous on $W^{1,p}(\Omega)$. However, the
known proofs use advanced instruments of real and functional analysis. Our aim
here is to present proof that can be easily understood by students familiar
only with the elementary measure theory.

##  1\. Introduction and the main result

Many optimization problems in physics may be formulated in terms of
minimization of the functional in the form of

(M) $F(v)=\int_{\Omega}f(x,v,\nabla v)\operatorname{d}x,$

over some set $\mathcal{W}$ of admissible functions. This problem is a very
classical one and its roots can be traced to the 17th century. The systematic
study was launched later by Lagrange, du Boys-Reymond, and others (see for
instance [21, 7]). They formulated the Euler-Lagrange equations as necessary
conditions to hold for the minimizer. However, the existence of the minimizer
has to be proved first (otherwise even solving Euler-Lagrange equations does
not guarantee to find the solution). As many examples show, the existence of
the minimizer cannot be expected in general. Therefore we have to add some
additional assumptions on the function $f$ in (M) to assure the existence of a
minimizer. The existence problem, however, is typically approached with
advanced methods including the knowledge of certain parts of functional
analysis, measure theory, and function spaces. We present an alternative proof
of the existence of a minimizer under the convexity condition of $f$ with a
given boundary data. Unlike the classical ones, our proof is using only
elementary methods of measure theory and real analysis and its direction is
easy to follow. Our calculation is done in detail and we intentionally avoid
shortening it with advanced tools.

Now we specify the problem. Let us consider a domain
$\Omega\subset\operatorname{\mathbb{R}}^{n}$ with a Lipschitz boundary and let
$f\in\mathcal{C}(\overline{\Omega}\times\operatorname{\mathbb{R}}\times\operatorname{\mathbb{R}}^{n})$.
Note that the condition of continuity of $f$ can be relaxed and one may
consider the Caratheodory function instead (see [6]), however, this is done
easily using the Scorza-Dragoni theorem 4.1, thus we restrict ourselves to
continuous functions. Let $u_{0}\in
W^{1,p}(\Omega,\operatorname{\mathbb{R}}^{d})$, we denote the set of all
functions in given Sobolev space with the same boundary data as $u_{0}$ by

(1) $\mathcal{W}:=u_{0}+W^{1,p}_{0}(\Omega,\operatorname{\mathbb{R}}^{d}).$

Let $p\in(1,\infty)$, by p-growth condition of $f$ we mean, if there exists
$1\leq q<p$ such that

(p-G) $\displaystyle(\exists c_{0}>0)(\exists
c_{1}\in\operatorname{\mathbb{R}},c_{2}\in
L^{1}(\Omega))(\forall(x,u,\xi)\in\overline{\Omega}\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}^{z})$
$\displaystyle f(x,u,\xi)\geq c_{0}|\xi|^{p}+c_{1}|u|^{q}+c_{2}(x).$

The proof can be slightly simplified by considering $c_{1}$ and $c_{2}$
positive constants, as some authors do. Also, $d$ is the dimension of the
image of functions in (1), which can be $d=1$ but also higher for vector-
valued functions, we leave $d$ general in statements, but we focus on the case
$d=1$. Note that the case of $p=1$ follows immediately from Lemma 3.1.
However, since the Sobolev space $W^{1,1}$ is not reflexive, even the
sequential weak lower semi-continuity does not guarantee the existence of the
minimizer. The case $p=\infty$ is also covered in some books (for instance see
[6]), however, this would require a different technique. We consider the
convexity property of $f$ in the last variable

(conv) $\xi\mapsto f(x,u,\xi)\quad\textup{is convex for
all}\quad(x,u)\in\Omega\times\operatorname{\mathbb{R}}^{d}.$

In the following text we consider only functional $F$ and set $\mathcal{W}$
satisfying

(fin) $(\exists u\in\mathcal{W})\quad F(u)<\infty\quad\wedge\quad(\forall
u\in\mathcal{W})\quad F(u)>-\infty.$

Our problem is to prove that the minimizer of (M) on $\mathcal{W}$ exists
under conditions (conv), (fin), and (p-G). Note that the existence follows
immediately from Theorems 1.1 and 1.2. Our contribution is the new student-
friendly proof of Theorem 1.1. The reader is encouraged to compare our proof
with the classical proofs in [1, 6, 24] or [14], in these proofs useful tools
are both approximations by the functions with controlled growth and regularity
or replacing and extending the functions on the irregular parts of the domain,
also the maximal operator is used. The shorter proof was done later by
Kałamajska [17] with help of the theory of Young measures. For the sake of
completeness and convenience of the reader, we also include proof of Theorem
1.2 although this one is similar to the known ones.

######  Theorem 1.1 (Weak sequential lower semicontinuity).

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a domain with a Lipchitz
boundary and let $u_{0}\in W^{1,p}(\Omega,\operatorname{\mathbb{R}}^{d})$. Let
$f\in\mathcal{C}(\overline{\Omega}\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}^{nd})$
be a function such that (conv) holds. Then the functional $F$ given by (M) is
sequentially weakly lower semicontinuous.

######  Theorem 1.2 (Coercivity).

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a domain with a Lipchitz
boundary and let $u_{0}\in W^{1,p}(\Omega,\operatorname{\mathbb{R}}^{d})$. Let
$f\in\mathcal{C}(\overline{\Omega}\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}^{nd})$
be a function satisfying (p-G). Let $u_{k}\in\mathcal{W}$ given by (1) be a
minimizing sequence of functional $F$ given by (M). Moreover, let condition
(fin) hold. Then $u_{k}$ is bounded in
$W^{1,p}(\Omega,\operatorname{\mathbb{R}}^{d})$.

######  Corollary 1.3 (Existence of minimizer).

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a domain with a Lipchitz
boundary. Let
$f\in\mathcal{C}(\overline{\Omega}\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}^{nd})$
be a function satisfying (conv) and (p-G). Let $F$ be the functional given by
(M) satisfying (fin). Then $F$ has a minimizer in $\mathcal{W}$ defined by
(1).

The paper is structured into Section 2, where we prepare some observations and
lemmata and prove Theorem 1.2, Section 3, where we prove the main result, and
Appendix 4, where we recall some well-known theorems for convenience of less
experienced readers.

###  1.1. Development and connection to literature

As for the development in the literature, the result we reprove is of Acerbi
and Fusco in 1984 [1]. The methods they used include maximal operators, an
extension of functions, and a lower semicontinuous envelope, the paper also
contains the proof even for quasi-convex vector valued $f$. The result
succeeded the previous results by Tonelli [32], where the Theorem 1.1 was
proved with much stricter conditions. Later the conditions were relaxed in the
regularity of $f$, for development of this see [30, 10] or [25]. For vector-
valued $u$, the assumption of convexity is no longer necessary and can be
relaxed to various weaker properties. Originally, the result was given by
Morrey [28] and generalized by Meyers [26] under stronger regularity requiring
only quasi-convexity. More results were published by Ball, Liu, Dacorogna,
Marcellini, and others with finer properties under more relaxed assumptions
such as polyconvexity or rank-one convexity. As the definitions of particular
generalizations of convexity differ, we are led to study differences and
inclusions between these classes, the topic was addressed by Šverák [34] and
Alibert and Dacorogna [2] or recent results by Sil [31] and Grabovsky [15].
The ongoing trend is to ensure an improvement of convexity type by studying
additional properties [16]. Note that even though (p-G) conditions can not be
left out, they can be relaxed into finer scales than powers [33]. An important
tool is controlling $f$ by the convex envelope, see [3] for studying its
properties and dependence on growth conditions. Let us emphasise that the
variation problems without convexity are also approached, see [8, 9].

Let us briefly overview some milestones and interesting papers in the field
surrendering or following the Acerbi-Fusco theorem [1]. In 1982, an important
result about the regularity of the minimizer was done by Giaquinta and Giusti
[13], for a wide survey on this topic we recommend the paper by Mingione [27].
Approximation methods by Marcellini [24] improve the Acerbi-Fusco theorem in
1985. Significant relaxation result was presented in 1997 by Fonseca and Malý
[11]. A version of convexity based on curl called A-convexity initially
considered by Dacorogna [5] is shown to be the optimal variant of definition
for the lower semicontinuity property by Fonseca and Müller in 1999 [12]. The
Young measure theory is used by Kałamajska in 1997 to shorten the proof of the
Acerbi-Fusco result. Also by the theory of Young measure, Kristensen in 1999
provided new approximation results in [18] and in 2015 he redefine growth
conditions for gradient Young functions and characterise the lower
semicontinuity in this setting [19]. The case of the lower semicontinuity
among $W^{1,1}$ and BV functions was studied in this setting by Kristensen and
Rindler in 2010 [20]. Recent result by Prinari covers the lower semicontinuity
and approximation properties for $L^{\infty}$ functionals [29]. The splendid
result by Bourdin, Francfort, and Gilles put the variational approach to
address Griffith fracture models [4] in 2007. The proper list of references
would be overwhelming thus we had to omit a lot.

##  2\. Some auxiliary results and observations

In this and the following text we use measure theory, but even though the
formulation of lemmata and observations are more general, we apply them only
for the Lebesgue measure in spaces $\mathbb{R}$ and $\mathbb{R}^{n}$.

######  Observation 2.1.

Let $(\Omega,\mu)$ be a finite measure space and let $f\in L^{0}(\Omega)$ and
$g\in L^{1}(\Omega)$ be functions for which $g(x)\leq f(x)$ holds for a.e.
$x\in\Omega$. Let $\Omega_{k}\subset\Omega$ be sets such that
$|\Omega\setminus\Omega_{k}|\to 0$. Then

$\lim_{k}\left(\int_{\Omega_{k}}f\operatorname{d}\mu\right)=\int_{\Omega}f\operatorname{d}\mu.$

Note, that partition of a set $M$ means the pair-wise disjoint family $M_{m}$
such that

$\bigcup_{m}M_{m}=M.$

Given a partition $\mathcal{P}$ we denote the norm of the partition by

$\nu(\mathcal{P}):=\max_{j\in m}\left\\{\textup{diam}(P_{j})\right\\}.$

######  Lemma 2.2.

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a bounded domain, $u\in
L^{1}(\Omega)$ and let $\mathcal{P}=(P_{j})_{j=1}^{m}$ be a partition of
$\Omega$. Define

$u_{\mathcal{P}}:=\sum_{j=1}^{m}\chi_{P_{j}}\fint_{P_{j}}u.$

Then

$\lim_{\nu(\mathcal{P})\to
0}|\\{|u_{\mathcal{P}}-u|>\varepsilon\\}|=0\quad(\forall\varepsilon>0).$

###### Proof.

Let us first show that the statement is valid for continuous functions. Pick
$\varepsilon,\gamma>0$. Given a positive $\eta$ let us denote

$\Omega_{\eta}:=\\{x\in\Omega:\operatorname*{dist}(x,\Omega^{c})>\eta\\}.$

Now we choose $\eta$ so small that

$|\Omega\setminus\Omega_{2\eta}|<\gamma.$

Since $u$ is uniformly continuous on $\Omega_{\eta}$, there exists $\delta>0$
such that

$|x-y|<\delta\Rightarrow|u(x)-u(y)|<\varepsilon\quad\forall
x,y\in\Omega_{\eta}.$

Now choose a partition $\mathcal{P}$ with

(2) $\nu({\mathcal{P}})<\min\\{\eta,\delta\\}.$

For all $x\in\Omega_{2\eta}$ we have that

$|u(x)-u_{\mathcal{P}}(x)|<\varepsilon.$

Hence

$|\\{|u-u_{\mathcal{P}}|>\varepsilon\\}|<\gamma$

as soon as (2) holds. The convergence in measure for the continuous function
is proved.

Given a $\alpha,\varepsilon>0$ and $u\in L^{1}(\Omega)$ we first choose $v$
continuous such that

$\|u-v\|_{L^{1}}\leq\varepsilon\alpha,$

then by the triangle inequality

$\begin{aligned}
\left|\\{|u-u_{\mathcal{P}}|>3\varepsilon\\}\right|&\leq|\\{|u-v|>\varepsilon\\}|+|\\{|v-v_{\mathcal{P}}|>\varepsilon\\}|+|\\{|u_{\mathcal{P}}-v_{\mathcal{P}}|>\varepsilon\\}|\end{aligned}\\\
=:I+II+III.$

Note that $u\mapsto u_{\mathcal{P}}$ is non-expansive mapping in $L^{1}$, so
we estimate

$\|u_{\mathcal{P}}-v_{\mathcal{P}}\|_{L^{1}}\leq\varepsilon\alpha.$

Hence, by Chebyschev inequality 4.8 we have $I+III<2\alpha$. Since $v$ is
continuous it is enough to choose the norm of the partition small to obtain
$II<\alpha$. ∎

###### Proof of Theorem 1.2.

Let $\tilde{u}_{k}\in\mathcal{W}$ be a minimizing sequence of $F$. Denote

(3) $u_{k}:=\tilde{u}_{k}-u_{0}.$

First, we prove that the sequence is bounded in $W^{1,p}(\Omega)$. Let us
suppose the contrary that $u_{k}$ is unbounded. Without loss of generality
(otherwise pass to a proper subsequence) suppose

$0<\|u_{k}\|_{W^{1,p}}\to\infty,$

the Friedrich inequality thm. 4.2 implies

$\|\nabla u_{k}\|_{p}\to\infty.$

This implication does not have to be true in general, but thanks to
$u_{k}\equiv 0$ on the boundary it has to hold. Now, by (p-G) we estimate

$\displaystyle F(\tilde{u}_{k})$ $\displaystyle\geq
c_{0}\|\nabla\tilde{u}_{k}\|_{p}^{p}-|c_{1}|\|\tilde{u}_{k}\|_{q}^{q}+\|c_{2}\|_{1}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:oddel}}}{{\geq}}c_{0}(\|\nabla
u_{k}\|_{p}^{p}-\|\nabla
u_{0}\|_{p}^{p})-|c_{1}|(\|u_{k}\|_{q}^{q}+\|u_{0}\|_{q}^{q})+\|c_{2}\|_{1}$
$\displaystyle\stackrel{{\scriptstyle\text{Thm.
\ref{thmFr}}}}{{\geq}}c_{0}\|\nabla u_{k}\|_{p}^{p}-C(\Omega)c_{1}\|\nabla
u_{k}\|_{q}^{q}+C(u_{0},c_{1},c_{2},\Omega)$
$\displaystyle\stackrel{{\scriptstyle\text{H\"{o}lder
in.}}}{{\geq}}c_{0}\|\nabla u_{k}\|_{p}^{p}-C(c_{1},\Omega)\|\nabla
u_{k}\|_{p}^{q}+C(u_{0},c_{1},c_{2},\Omega)$ $\displaystyle=\|\nabla
u_{k}\|_{p}^{p}\left(c_{0}-C(c_{1},\Omega)\|\nabla u_{k}\|^{q-p}-\|\nabla
u_{k}\|_{p}^{-p}C(u_{0},c_{1},c_{2},\Omega)\right).$

As $p>q$, the second and the third term in the bracket tend to zero and
$c_{0}>0$, so the remaining expression diverges to infinity as $k$ tends to
infinity. But this is a contradiction to the fact that $u_{k}$ is a minimizing
sequence.

∎

##  3\. Proof of Theorem 1.1 and Corollary 1.3

Let $u\in L^{0}(\Omega,\operatorname{\mathbb{R}}^{d})$ and $v\in
L^{p}(\Omega,\operatorname{\mathbb{R}}^{z})$. Given such functions let us
consider

(4) $J_{\Omega}(u,v):=\int_{\Omega}f(x,u(x),v(x))\operatorname{d}x.$

Let us denote the norm topology in $L^{p}$ by $\mathcal{L}^{p}$ and the
topology of weak convergence by $\mathcal{L}^{p}_{w}$. We also denote the
convergence in the Lebesgue measure by
$\stackrel{{\scriptstyle\operatorname{d}x}}{{\to}}$. Let us define the
topology

$\tau=\mathcal{L}^{0}\times\mathcal{L}^{p}_{w}$

for $(u_{k},v_{k})\in L^{0}(\Omega,\operatorname{\mathbb{R}}^{d})\times
L^{p}(\Omega,\operatorname{\mathbb{R}}^{z})$ by

$(u_{k},v_{k})\stackrel{{\scriptstyle\tau}}{{\to}}(u,v)\quad\textup{if
}\left(u_{k}\stackrel{{\scriptstyle\operatorname{d}x}}{{\to}}u\wedge
v_{k}\stackrel{{\scriptstyle L^{p}}}{{\rightharpoonup}}v\right).$

######  Lemma 3.1.

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a bounded domain,
$p\in[1,\infty)$,
$f\in\mathcal{C}(\overline{\Omega}\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}^{z})$
be function satisfying (conv) and let $b\in L^{1}(\Omega)$ be a function such
that

(5) $f(x,u,v)\geq
b(x)\quad\forall(x,u,v)\in\Omega\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}^{z}.$

Then the functional $J_{\Omega}$ defined in (4) is sequentially lower
semicontinuous with respect to topology $\tau$.

######  Lemma 3.2.

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a bounded domain,
$p\in(1,\infty)$ and let
$f\in\mathcal{C}(\overline{\Omega}\times\operatorname{\mathbb{R}}^{d}\times\operatorname{\mathbb{R}}^{z})$
be a function satisfying (p-G) and (conv). Then $J_{\Omega}$ defined in (4) is
sequentially lower semicontinuous with respect to topology
$\mathcal{L}^{q}\times\mathcal{L}^{p}_{w}$.

###### Proof of Lemma 3.1.

Choose $\varepsilon>0$. Let $(u_{k},v_{k})$ be a sequence convergent to
$(u,v)$ in the topology $\tau$. Since $u_{k}$ converges to $u$ in measure, by
the Riesz thm. 4.4 one may assume that $u_{k}(x)\to u(x)$ for a.e.
$x\in\Omega$. Hence by the Yegorov thm. 4.6 and the Luzin thm. 4.7 for every
$l$ natural there exists a compact set $\Omega_{l}\subset\Omega$ such that

$\left(|\Omega\setminus\Omega_{l}|\leq\frac{1}{l}\right)\wedge\left(u_{k}\rightrightarrows
u\quad\textup{on
}\Omega_{l}\right)\wedge\left(u\in\mathcal{C}(\Omega_{l})\right).$

Define

$S:=\sup_{k}\sup_{x\in\Omega_{l}}|u_{k}(x)|.$

We may assume that $S<\infty$. Note that the function $(x,u)\mapsto
f(x,u,\xi)$ is uniformly continuous on
$\Omega_{l}\times[-S,S]\times\overline{B(0,2j)}$ for arbitrary $j$ natural.
Now, let us denote by $\eta_{j}$ the strictly positive number such that for
all $x,\overline{x}\in\Omega_{l},u,\overline{u}\in[-S,S]$

(6) $|(\overline{x},\overline{u})-(x,u)|<6\eta_{j}$

implies

$|f(\overline{x},\overline{u},\xi)-f(x,u,\xi)|<\varepsilon\quad(\forall|\xi|\leq
2j).$

Without loss of generality one may suppose that $\eta_{j}$ is decreasing and

(7) $|u(x)-u_{k}(x)|\leq\eta_{k}\quad(\forall x\in\overline{\Omega_{l}}).$

Moreover, since the function $u$ is uniformly continuous, by the triangle
inequality and double usage of (7), for each $j$ there exists positive
$\delta_{j}<\eta_{j}$ ($\delta_{j}$ decreasing) such that for all
$x,y\in\overline{\Omega}$, one has

$|x-y|<\delta_{j}\quad\textup{implies}\quad|(x,u_{k}(x))-(y,u_{k}(y))|<3\eta_{j}\quad(\forall
k\geq j).$

Use the previous estimates to observe

(8)
$|x-y|<\delta_{j}\quad\textup{implies}\quad|(x,u_{k}(x))-(y,u(y))|<6\eta_{j}\quad(\forall
k\geq j).$

For $j,k$ naturals denote

(9) $E_{k,j}:=\\{|v_{k}|\leq j\\}\cap\Omega_{l},\quad G_{j}:=\\{|v|\leq
j\\}\cap\Omega_{l}.$

Note that, using the Chebychev inequality thm. 4.8 and boundedness of $v_{k}$
in $L^{p}$ we obtain that there exists $C>0$ independent of $j,k$ such that

(10) $|\Omega_{l}\setminus G_{j}|\leq\frac{C}{j^{p}}\quad\text{ and
}\quad|\Omega_{l}\setminus E_{k,j}|\leq\frac{C}{j^{p}}.$

By uniform continuity of $(x,u,\xi)\mapsto f(x,u,\xi)$ on
$\Omega_{l}\times[-S,S]\times\overline{B(0,2j)}$, for every natural $j$ there
exists $\gamma_{j}$ such that

(11)
$|\xi-\overline{\xi}|<\gamma_{j}\Rightarrow|f(x,u,\xi)-f(x,u,\overline{\xi})|<\varepsilon\quad(\forall(x,u)\in\overline{\Omega_{l}}\times[-S,S]\quad\textup{and}\quad\forall\xi,\overline{\xi}\in\overline{B(0,2j)}$

By Lemma 2.2 we obtain that for each $j$ natural there exists a finite
partition of $G_{j}$, consider such a partition
$\mathcal{P}=(K^{j}_{m})_{m=1}^{M_{j}}$ satisfying

  1. (i)

$\nu(\mathcal{P})<\delta_{j}$.

  2. (ii)

$\sum_{m}\left|\\{x\in
K^{j}_{m}:|v-\fint_{K^{j}_{m}}v|>\gamma_{j}\\}\right|<1/j.$

By the weak convergence of $v_{k}$ to $v$ in $L^{p}$ one has

$\lim_{k}\fint_{K^{j}_{m}}(v-v_{k})=0.$

Therefore, one may assume without the loss of generality

(12) $\left|\fint_{K^{j}_{m}}(v-v_{k})\right|\leq\gamma_{k}\quad(\forall k\geq
j)$

(if this doesn’t hold one can pass to the appropriate sub-sequence). Set

(13) $Q^{j}_{m}:=\left\\{x\in
K^{j}_{m}:|v(x)-\fint_{K^{j}_{m}}v|\leq\gamma_{j}\right\\}$

and denote

(14) $Q_{j}:=\bigcup_{m}Q^{j}_{m}.$

Note that by (9) and the property (ii) of the parts $K^{j}_{m}$ one has
$|G_{j}\setminus Q_{j}|\leq\frac{1}{j}.$ Pick $a_{j}>j$ such that

(15)
$\left|\fint_{K^{j}_{m}}v_{j}\chi_{E^{c}_{j,a_{j}}}\operatorname{d}x\right|<\gamma_{j}.$

Moreover, let $Q_{m,t}^{j}$ be a partition of $Q^{j}_{m}$ such that

$\textup{diam}(Q_{m,t}^{j})<\delta_{a_{j}},$

and let us choose arbitrary $x^{j}_{m,t}\in Q_{m,t}^{j}$, creating
$2\delta_{a_{j}}$ net. For the following estimate denote

(16) $\displaystyle\xi_{j,m}:=$
$\displaystyle\fint_{K^{j}_{m}}v\operatorname{d}x,$
$\displaystyle\xi^{\spadesuit}_{j,m}:=$
$\displaystyle\fint_{K^{j}_{m}}v_{j}\operatorname{d}x,$
$\displaystyle\xi^{\heartsuit}_{j,m}:=$
$\displaystyle\fint_{K^{j}_{m}}v_{j}\chi_{E_{j,a_{j}}}\operatorname{d}x,$
$\displaystyle\tilde{v}_{j}:=$ $\displaystyle v_{j}\chi_{E_{j,a_{j}}},$
$\displaystyle T:=$
$\displaystyle\max_{(x,u)\in\overline{\Omega}_{l}\times[-S,S]}|f(x,u,0)|.$

Note that (15) estimate the difference of functions above as

(17)
$\left|\xi^{\spadesuit}_{j,m}-\xi^{\heartsuit}_{j,m}\right|=\left|\fint_{K^{j}_{m}}v_{j}\chi_{E^{c}_{j,a_{j}}}\operatorname{d}x\right|<\gamma_{j}.$

Let us estimate

$\displaystyle\int_{Q_{j}}f(x,u,v)\operatorname{d}x$
$\displaystyle\stackrel{{\scriptstyle\eqref{sit2},\eqref{DUL}}}{{\leq}}$
$\displaystyle\sum_{m}\int_{Q^{j}_{m}}f(x^{j}_{m},u_{j}(x^{j}_{m}),v)\operatorname{d}x+|\Omega_{l}|\varepsilon$
$\displaystyle\stackrel{{\scriptstyle\eqref{BytostQ},\eqref{PANTAU}}}{{\leq}}$
$\displaystyle\sum_{m}\int_{Q^{j}_{m}}f(x^{j}_{m},u_{j}(x^{j}_{m}),\xi_{j,m})\operatorname{d}x+2|\Omega_{l}|\varepsilon$
$\displaystyle\stackrel{{\scriptstyle\eqref{konvprum},\eqref{PANTAU}}}{{\leq}}$
$\displaystyle\sum_{m}\int_{Q^{j}_{m}}f(x^{j}_{m},u_{j}(x^{j}_{m}),\xi^{\spadesuit}_{j,m})\operatorname{d}x+3|\Omega_{l}|\varepsilon$
$\displaystyle\stackrel{{\scriptstyle\eqref{sit2},\eqref{DUL}}}{{\leq}}$
$\displaystyle\sum_{m}\sum_{t}\int_{Q_{m,t}^{j}}f(x^{j}_{m,t},u_{j}(x^{j}_{m,t}),\xi^{\spadesuit}_{j,m})\operatorname{d}x+4|\Omega_{l}|\varepsilon$
$\displaystyle\stackrel{{\scriptstyle\eqref{ajajaj},\eqref{PANTAU}}}{{\leq}}$
$\displaystyle\sum_{m}\sum_{t}\int_{Q_{m,t}^{j}}f(x^{j}_{m,t},u_{j}(x^{j}_{m,t}),\xi^{\heartsuit}_{j,m})\operatorname{d}x+5|\Omega_{l}|\varepsilon$
$\displaystyle\stackrel{{\scriptstyle\text{Thm }\ref{thmJe}}}{{\leq}}$
$\displaystyle\sum_{m}\sum_{t}\int_{Q_{m,t}^{j}}f(x^{j}_{m,t},u_{j}(x^{j}_{m,t}),\tilde{v}_{j})\operatorname{d}x+5|\Omega_{l}|\varepsilon.$

Note that the usage of the Jensen inequality in the last estimate is the
crucial step, where convexity (conv) comes to play. Note that in this step, we
may question if more general versions such as quasi-convexity can be used.
Following the fact that $x^{j}_{m,t}$ is a 2$\delta_{a_{j}}$ net and the
inequalities (6),(8) and (7), we estimate

$\displaystyle\leq$
$\displaystyle\sum_{m}\sum_{t}\int_{Q_{m,t}^{j}}f(x,u_{j},\tilde{v}_{j})\operatorname{d}x+6|\Omega_{l}|\varepsilon$
$\displaystyle\stackrel{{\scriptstyle\eqref{BytostQ naporcovana},\text{part.
of }Q_{m}^{j}}}{{=}}$
$\displaystyle\int_{Q_{j}}f(x,u_{j},\tilde{v}_{j})\operatorname{d}x+6|\Omega_{l}|\varepsilon.$

To estimate we split th domain into $Q_{j}\cap E_{j,a_{j}}$ where
$\tilde{v}_{j}=v_{j}$ and $Q_{j}\cap E_{j,a_{j}}^{c}$, where we estimate
$|f|\leq T$. In the last estimate, we add to the right-hand side integral of
non-negative $f-b$ over the set $\Omega_{l}\setminus(E_{j,a_{j}}\cap Q_{j})$,
so we estimate

$\displaystyle\stackrel{{\scriptstyle\eqref{znackaT}}}{{\leq}}$
$\displaystyle\int_{E_{j,a_{j}}\cap
Q_{j}}f(x,u_{j},v_{j})\operatorname{d}x+T|E_{j,a_{j}}^{c}|+6|\Omega_{l}|\varepsilon$
$\displaystyle\stackrel{{\scriptstyle\eqref{mir}}}{{\leq}}$
$\displaystyle\int_{\Omega_{l}}f(x,u_{j},v_{j})\operatorname{d}x+\frac{TC}{a_{j}^{p}}+6|\Omega_{l}|\varepsilon-\int_{\Omega_{l}\setminus(E_{j,a_{j}}\cap
Q_{j})}b(x)\operatorname{d}x.$

Passing to the limit with $\varepsilon\to 0_{+}$ yields

$J_{Q_{j}}(u,v)\leq
J_{\Omega_{l}}(u_{j},v_{j})+\frac{TC}{a_{j}^{p}}-\int_{\Omega_{l}\setminus(E_{j,a_{j}}\cap
Q_{j})}b(x)\operatorname{d}x.$

Using Observation 2.1 and (5), passing to limit for $j\to\infty$ we obtain

$J_{\Omega}(u,v)\leq\limsup_{j}J_{\Omega_{l}}(u_{j},v_{j}),$

which is equivalent to sequential lower semicontinuity of $J_{\Omega}$ on the
topological space $(L^{0}(\Omega)\times L^{p}(\Omega),\tau)$. ∎

###### Proof of the Lemma 3.2.

Given a function satisfying (p-G), define the auxiliary function

$\tilde{f}(x,u,\xi):=f(x,u,\xi)-c_{1}|u|^{q}$

and functional

$\tilde{J}(u,v):=\int_{\Omega}\tilde{f}(x,u,v)\operatorname{d}x.$

By the Lemma 3.1 the functional $\tilde{J}$ is sequentially lower
semicontinuous with respect to topology $\tau$. Hence it is sequentially lower
semicontinuous with respect to $\mathcal{L}^{q}\times\mathcal{L}^{p}_{w}$ and
thus

$J(u,v)=\tilde{J}(u,v)+c_{1}\int_{\Omega}|u|^{q}\operatorname{d}x\leq\liminf_{k}\tilde{J}(u_{k},v_{k})+c_{1}\lim_{k}\int_{\Omega}|u_{k}|^{q}\operatorname{d}x=\liminf_{k}J(u_{k},v_{k}).$

∎

###### Proof of Theorem 1.1.

Let $u_{k}\rightharpoonup u$ in $W^{1,p}$. Using the Rellich-Kondrachev
theorem 4.5 without loss of generality, we may suppose that $u_{k}\to u$
strongly in $L^{p}$ while $\nabla u_{k}\stackrel{{\scriptstyle
L^{p}}}{{\rightharpoonup}}\nabla u$. Hence by the preceding Lemma 3.2, we
obtain that

$F(u)=J(u,\nabla u)\leq\liminf_{k}J(u_{k},\nabla u_{k})=\liminf_{k}F(u_{k}).$

∎

###### Proof of Corollary 1.3.

Let
$f\in\mathcal{C}(\Omega\times\operatorname{\mathbb{R}}\times\operatorname{\mathbb{R}}^{n})$
be a function satisfying (p-G) and (conv). Let $u_{k}$ be a minimizing
sequence. By Theorem 1.2 we have that $u_{k}$ is bounded in $W^{1,p}(\Omega)$.
Alaoglu theorem implies that the closed ball in $W^{1,p}(\Omega)$ is compact,
therefore we may assume that $u_{k}$ is weakly convergent to $u\in
W^{1,p}(\Omega)$. Let us emphasise that by the Rellich-Kondrachev theorem 4.5
we also may assume that $u_{k}\to u$ in $L^{q}(\Omega)$. We get

$(u_{k},\nabla u_{k})\to(u,\nabla u)\quad\textup{in the
topology}\quad\mathcal{L}^{q}\times\mathcal{L}^{p}_{w}.$

Using the Theorem 1.1 we obtain that

$F(u)\leq\liminf_{k\to\infty}F(u_{k})=\inf\\{F(u)|u\in\mathcal{W}\\}.$

Hence $u$ is the minimizer. ∎

##  4\. Appendix

The following results are standard knowledge, the proofs can be found in [23],
[6], [22] etc.

######  Theorem 4.1 (Scorza-Dragoni).

Let $\varepsilon>0$ and
$f:\Omega\times\operatorname{\mathbb{R}}\times\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}$
be measurable in the first variable and continuous in the second and third
one. Then there exists compact set $K\subset\Omega$ such that
$|\Omega\setminus K|<\varepsilon$ and
$f|_{K\times\operatorname{\mathbb{R}}\times\operatorname{\mathbb{R}}^{n}}$ is
continuous.

######  Theorem 4.2 (Friedrichs inequality).

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a domain with Lipchitz
boundary, $p\in\left<1,\infty\right>$. There exists constant $C_{\Omega}$ such
that

$\|u\|_{p}\leq C_{\Omega}\|\nabla u\|_{p}$

holds for all $u\in W^{1,p}_{0}(\Omega)$.

######  Theorem 4.3 (Jensen inequality).

Let $N\subset\operatorname{\mathbb{R}}^{n}$ be a convex set. Let $\mu$ be a
measure with a finite total variation on $N$. Let
$\varphi:N\to\operatorname{\mathbb{R}}$ be a convex function. Let $E\subset N$
be a measurable set. Then the following inequality holds

$\varphi\left(\fint_{E}x\operatorname{d}\mu(x)\right)\leq\fint_{E}\varphi(x)\operatorname{d}\mu(x).$

######  Theorem 4.4 (Riesz).

Let $(\Omega,\mu)$ be a space with measure and let $u_{k}$ be a sequence of
measurable functions such that

$u_{k}\stackrel{{\scriptstyle\mu}}{{\to}}u\quad u_{k}\textup{ converge in
measure to }u.$

then there exists an increasing sequence of indices $l_{k}$ such that

$u_{l_{k}}(x)\to u(x)\quad\textup{for }\mu-\textup{ a.e. }x\in\Omega.$

######  Theorem 4.5 (Rellich-Kondrachev).

Let $p\in(1,\infty)$ and let $1\leq q<p^{*}$ where

$p^{*}:=\begin{cases}\frac{np}{n-p}\quad&\textup{for }p<n,\\\
\infty&\textup{otherwise.}\end{cases}$

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ be a bounded domain with
Lipchitz boundary. And let $u_{k}\in W^{1,p}(\Omega)$ be a bounded sequence in
$W^{1,p}(\Omega)$, then there exists an increasing sequence of indices $l_{k}$
and $u\in L^{q}(\Omega)$ such that

$\|u-u_{l_{k}}\|_{L^{q}(\Omega)}\to 0.$

######  Theorem 4.6 (Yegorov).

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ measurable, let $u_{k}$ be a
sequence of measurable functions such that

$u_{k}(x)\to u(x)\quad\textup{for }\mu-\textup{a.e. }x\in\Omega.$

Then there exist the sequence of sets $N_{j}$ and the strictly increasing
sequence of indices $k_{j}$ such that

  1. (Y1)

$N_{j}\subset N_{j+1},$

  2. (Y2)

$\mu(\Omega\setminus N_{j})\leq\frac{1}{j},$

  3. (Y3)

$|u-u_{k_{j}}|<\frac{1}{j}\quad\textup{on }N_{j}.$

######  Theorem 4.7 (Luzin).

Let $\mu$ be a complete Radon measure on locally compact space and measurable
function $u$ finite $\mu$-almost everywhere. Then for open $K$ and
$\varepsilon>0$ there exist open $G\subset K$ such that $u$ is continuous on
$K\setminus G$ and $\mu(G)<\varepsilon.$

######  Theorem 4.8 (Chebysev inequality).

Let $\Omega\subset\operatorname{\mathbb{R}}^{n}$ measurable, $u\in
L^{0}(\Omega)$, $p\in\left<1,\infty\right)$ and let $t\in(0,\infty)$ then one
has

$t^{p}\left|\left\\{|u|>t\right\\}\right|\leq\int_{\Omega}|u(x)|^{p}\operatorname{d}x.$

## Acknowledgments

Both named authors were supported by grant number GJ20-19018Y of the Grant
Agency of the Czech Republic.

## References

  * [1] Emilio Acerbi and Nicola Fusco.  Semicontinuity problems in the calculus of variations.  Archive for Rational Mechanics and Analysis, 86(2):125–145, 1984\.
  * [2] Jean-Jacques Alibert and Bernard Dacorogna.  An example of a quasiconvex function that is not polyconvex in two dimensions.  Arch. Rational Mech. Anal., 117(2):155–166, 1992.
  * [3] John M. Ball, Bernd Kirchheim, and Jan Kristensen.  Regularity of quasiconvex envelopes.  Calc. Var. Partial Differential Equations, 11(4):333–359, 2000\.
  * [4] Blaise Bourdin, Gilles A. Francfort, and Jean-Jacques Marigo.  The variational approach to fracture.  J. Elasticity, 91(1-3):5–148, 2008.
  * [5] Bernard Dacorogna.  Weak continuity and weak lower semicontinuity of nonlinear functionals, volume 922 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin-New York, 1982.
  * [6] Bernard Dacorogna.  Direct methods in the calculus of variations, volume 78.  Springer Science & Business Media, 2007.
  * [7] Paul du Bois-Reymond.  Erläuterungen zu den anfangsgründen der variationsrechnung.  Mathematische Annalen, 15(2):283–314, 1879.
  * [8] Ivar Ekeland.  On the variational principle.  J. Math. Anal. Appl., 47:324–353, 1974.
  * [9] Ivar Ekeland.  Nonconvex minimization problems.  Bull. Amer. Math. Soc. (N.S.), 1(3):443–474, 1979.
  * [10] Ivar Ekeland and Roger Temam.  Analyse convexe et problèmes variationnels.  Collection Études Mathématiques. Dunod, Paris; Gauthier-Villars, Paris-Brussels-Montreal, Que., 1974.
  * [11] Irene Fonseca and Jan Malý.  Relaxation of multiple integrals below the growth exponent.  Ann. Inst. H. Poincaré C Anal. Non Linéaire, 14(3):309–338, 1997.
  * [12] Irene Fonseca and Stefan Müller.  A-quasiconvexity, lower semicontinuity, and Young measures.  SIAM J. Math. Anal., 30(6):1355–1390, 1999.
  * [13] Mariano Giaquinta and Enrico Giusti.  On the regularity of the minima of variational integrals.  Acta Math., 148:31–46, 1982.
  * [14] Enrico Giusti.  Direct methods in the calculus of variations.  World Scientific, 2003.
  * [15] Yury Grabovsky.  From microstructure-independent formulas for composite materials to rank-one convex, non-quasiconvex functions.  Arch. Ration. Mech. Anal., 227(2):607–636, 2018.
  * [16] André Guerra and Jan Kristensen.  Automatic quasiconvexity of homogeneous isotropic rank-one convex integrands.  Arch. Ration. Mech. Anal., 245(1):479–500, 2022.
  * [17] Agnieszka Kałamajska.  On lower semicontinuity of multiple integrals.  Colloq. Math., 74(1):71–78, 1997.
  * [18] Jan Kristensen.  Lower semicontinuity in spaces of weakly differentiable functions.  Math. Ann., 313(4):653–710, 1999.
  * [19] Jan Kristensen.  A necessary and sufficient condition for lower semicontinuity.  Nonlinear Anal., 120:43–56, 2015.
  * [20] Jan Kristensen and Filip Rindler.  Characterization of generalized gradient Young measures generated by sequences in $W^{1,1}$ and BV.  Arch. Ration. Mech. Anal., 197(2):539–598, 2010.
  * [21] Joseph Louis Lagrange.  Mécanique analytique, volume 1.  Mallet-Bachelier, 1853.
  * [22] Giovanni Leoni.  A first course in Sobolev spaces, volume 181 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, second edition, 2017.
  * [23] Jaroslav Lukeš and Jan Malý.  Measure and integral.  Matfyzpress, Prague, second edition, 2005.
  * [24] Paolo Marcellini.  Approximation of quasiconvex functions, and lower semicontinuity of multiple integrals.  Manuscripta Math., 51(1-3):1–28, 1985.
  * [25] Paolo Marcellini and Carlo Sbordone.  Semicontinuity problems in the calculus of variations.  Nonlinear Anal., 4(2):241–257, 1980.
  * [26] Norman G. Meyers.  Quasi-convexity and lower semi-continuity of multiple variational integrals of any order.  Trans. Amer. Math. Soc., 119:125–149, 1965.
  * [27] Giuseppe Mingione.  Regularity of minima: an invitation to the dark side of the calculus of variations.  Appl. Math., 51(4):355–426, 2006.
  * [28] Charles B. Morrey, Jr.  Quasi-convexity and the lower semicontinuity of multiple integrals.  Pacific J. Math., 2:25–53, 1952.
  * [29] Francesca Prinari.  On the lower semicontinuity and approximation of $L^{\infty}$-functionals.  NoDEA Nonlinear Differential Equations Appl., 22(6):1591–1605, 2015\.
  * [30] James Serrin.  On the definition and properties of certain variational integrals.  Trans. Amer. Math. Soc., 101:139–167, 1961.
  * [31] Swarnendu Sil.  Calculus of variations: a differential form approach.  Adv. Calc. Var., 12(1):57–84, 2019.
  * [32] Leonida Tonelli.  La semicontinuitĂ nel calcolo delle variazioni.  Rendiconti del Circolo Matematico di Palermo (1884-1940), 44(1):167–249, 1920.
  * [33] Anna Verde and Gabriella Zecca.  Lower semicontinuity of certain quasiconvex functionals in Orlicz-Sobolev spaces.  Nonlinear Anal., 71(10):4515–4524, 2009.
  * [34] Vladimír Šverák.  Rank-one convexity does not imply quasiconvexity.  Proc. Roy. Soc. Edinburgh Sect. A, 120(1-2):185–189, 1992.