\documentclass[reqno]{amsart}
\usepackage{graphicx}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
{\em Electronic Journal of Differential Equations},
Vol. 2007(2007), No. 114, pp. 1--22.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2007/22\hfil Nonlinear waves formation]
{Global solution to a Hopf equation and its application to
non-strictly hyperbolic systems of conservation laws}

\author[D. Mitrovic, J. Susic\hfil EJDE-2007/22\hfilneg]
{Darko Mitrovic, Jela Susic}  % in alphabetical order

\address{Darko Mitrovic \newline
Department of Mathematical Sciences, Norwegian Institute
of Science and Technology, Alfred Getz vei 1, NO-7491 Trondheim, Norway}
\email{mitrovic@math.ntnu.no}

\address{Jela Susic \newline
 Faculty of Mathematics and Natural Sciences, University of Montenegro, 81000
Podgorica, Montenegro}
\email{jela@rc.pmf.cg.ac.yu}

\thanks{Submitted December 5, 2006. Published August 22, 2007.}
\thanks{The first author is partially supported by the Research Council of Norway}
\subjclass[2000]{35L65}
\keywords{Weak asymptotic method; Hopf equation; shock wave formation; \hfill\break\indent
pressureless gas dynamics; system of conservation laws; multidimensional shocks}

\begin{abstract}
 From a Hopf equation we develop a recently introduced technique, the
 weak asymptotic method, for describing the shock wave formation and
 the interaction processes. Then, this technique is applied to a
 system of conservation laws arising from pressureless gas dynamics.
 As an example, we study the shock wave formation process in a
 two-dimensional scalar conservation laws arising in oil reservoir
 problems.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}

\section{Introduction}

The starting point of this paper is the Hopf equation
\begin{equation}
u_t+(u^2)_x=0, \label{avg715}
\end{equation}
with the initial condition
\begin{equation}
u_0(x)=\begin{cases}
U, & x\leq a_2\\
-Kx+b, & a_2< x < a_1\\
u_0^0, & a_1\leq x.
\end{cases}
\label{avg915}
\end{equation}
Here $U>u_0^0$ and $K$, $b$ are constants that satisfy
 $-Ka_1+b=u_0^0$ and $-Ka_2+b=U$.

Our aim is to find global approximating solutions for this problem;
more precisely,  to describe the shock wave formation process.
Although this problem sounds simple and well known,
this is a very interesting model for developing the
technique in this paper. Namely, if we understand the problem
of shock wave formation properly in this case, we can apply
the same procedure to the scalar conservation laws with arbitrary
nonlinearity \cite{DM5}, to  various systems of conservation laws,
and to multidimensional scalar conservation laws.
The latter two task will be presented here.

Global in time, $t\in \mathbb{R}^+$, approximating solutions to
\eqref{avg715}--\eqref{avg915} can be obtained by means of vanishing
viscosity regularization combined with the Florin-Hopf-Cole
transformation. Here, we shall use more general procedure - the
\emph{weak asymptotic method} (for more information about the method
see \cite{dan,DM,do,DSO,DSH}). Solutions obtained using this method
are called \emph{weak asymptotic solutions} and are defined as
follows.

\begin{definition} \label{def1} \rm
By $O_{\mathcal{D}'}(\varepsilon^\alpha)$, $\varepsilon\in (0,1)$,  we denote a family of
distributions depending $t\in \mathbb{R}^+$ such that for any test
function $\eta(x)\in C^1_0(\mathbb{R})$, the estimate
$$
\langle O_{\mathcal{D}'}(\varepsilon^\alpha),\eta(x)\rangle=O(\varepsilon^\alpha)
$$
holds, where the estimate on the right-hand side is understood in
the usual sense and is locally uniform in $t$; i.e.,
$|O(\varepsilon^\alpha)|\leq C_T\varepsilon^\alpha$ for $t\in[0,T]$.
\end{definition}

\begin{definition} \label{def2} \rm
A family of functions $u_\varepsilon=u_\varepsilon(x,t)$,
$\varepsilon>0$, is called a weak asymptotic solution of problem
\eqref{avg715}, \eqref{avg915} if
\begin{equation*}%(6)
\frac{\partial u_\varepsilon}{\partial t}+\frac{\partial u^2_\varepsilon}{\partial
x}=O_{\mathcal{D}'}(\varepsilon^{\alpha}), \quad
u_\varepsilon\Big|_{t=0}-u\Big|_{t=0}=O_{\mathcal{D}'}(\varepsilon^{\alpha}), \quad
\alpha>0.
\end{equation*}
\end{definition}

As we can see from these definitions, in the framework of the weak
asymptotic method, the discrepancy is assumed to be small in the
sense of space of functionals $\mathcal{D}'(\mathbb{R})$ over test
functions depending only on the ``space'' variable $x$. Such
approach allows us to reduce the problem of nonlinear wave
interaction (in this case interaction of weak discontinuities) to
the problem of solving a system of ordinary differential equations
(see \eqref{avg735}, \eqref{avg745}).

A more general situation than the one considered here and in
 \cite{dan} was analyzed in \cite{DM}.
There the passage from continuous to discontinuous state of
the solution (including the uniform,
 $t\in \mathbb{R}^+$, description of interaction of weak
discontinuities) was described for scalar conservation laws with arbitrary
convex nonlinearity.

In this paper we propose another procedure for describing the
shock wave formation process and show how to apply the obtained
results to a non-strictly hyperbolic system of conservation laws
as well as to a multidimensional scalar conservation law.

A step forward from this paper is  the form of the ansatz of the
solution. In \cite{dan} and \cite{DM} very special form of
ansatz is used, which can be an obstacle for describing the
passage from continuous to discontinuous state of the solution in
general situations such as systems of conservation laws.
Furthermore, in \cite{DM} rather sophisticated
(complicated) mathematical tools are used (such as complex germ
lemma, asymptotic linear independence and some nontrivial
estimates).  In our approach, the ansatz has a rather general form
and the method used  can be generalized for scalar
conservation laws with arbitrary nonlinearity, for systems
of conservation laws and for arbitrary multidimensional scalar
conservation laws, almost without any changes.
The approaches in \cite{dan} and \cite{DM} are rather different,
although the problem \cite{dan} is special case of the problem
considered in \cite{DM}.


Concerning the shock wave formation process, as in \cite{DM},
standard characteristics are replaced by \emph{new characteristics}
(of course, new characteristics here and in \cite{DM} have different
form) which, unlike standard characteristics, never intersect
(otherwise they bear the same information). However, along the new
characteristics, the solution to the problem remains constant, and,
as $\varepsilon\to 0$ the new characteristics coincide with the
standard characteristics
 (with a discontinuity line in the appropriate domain).
Accordingly, the solution of our problem is
found along the characteristics and it is defined as long as new
characteristics do not intersect; and that is along the entire time
axis.

The usage of new characteristics is rather important  since they may
represent the first effective attempt to generalize the notion of
standard characteristics on the realm of weak solutions (see
\cite{bre, daf}). After we  develop the technique for the Hopf
equation, we use the  approximating solution to problem
\eqref{avg715}-\eqref{avg915} to solve the Riemann problem
\eqref{jul165}-\eqref{jul166}. A good introduction to the study of
the Riemann problem \eqref{jul165}-\eqref{jul166} can be found in
\cite{DSH}. For a complete study of this problem and its variants,
see for example \cite{chen,DSH2,DSH,DSH3,She,PSH,ind,mn1,mn} and
references therein.

At the end of this paper, as an example, we apply our technique to
the problem of shock wave formation  in two dimensional
scalar conservation laws (more precisely, for an oil
reservoir problem). As far as we know, the problem of shock wave
formation in the multidimensional case was mainly treated with the
techniques based on  geometry \cite{kos, nak}. Here, from a
simple example, we propose principles for an analytical approach.

\section{Main result}

First we introduce an auxiliary statement proved in \cite{dan,
DSO} which is called nonlinear superposition law in the
quadratic case.


\begin{theorem}\label{th**}
Let $\omega_i\in C^\infty_0(\mathbb{R})$, $i=1,2$, where $\lim_{z\to
+\infty}\omega_i=1$, $\lim_{z\to -\infty}\omega_i=0$ and
$\frac{d\omega_i(z)}{dz}\in \mathcal{S}(\mathbb{R})$, $i=1,2$, where
$\mathcal{S}(\mathbb{R})$ is the space of rapidly decreasing
functions. For $\varphi_i\in \mathbb{R}$, $i=1,2$, we have
\begin{equation}
\begin{aligned}
&\omega(\frac{\varphi_1-x}{\varepsilon})\omega(\frac{\varphi_2-x}{\varepsilon}) \\
&=B_1\Big(\frac{\varphi_2-\varphi_1}{\varepsilon}\Big)H(\varphi_1-x)
+B_2\Big(\frac{\varphi_2-\varphi_1}{\varepsilon}\Big)H(\varphi_2-x)+
\mathcal{O}_{\mathcal{D}'}(\varepsilon), \quad x\in \mathbb{R},
\end{aligned}\label{an0}
\end{equation}
where $H$ is the Heaviside function and for $\rho\in \mathbb{R}$,
\begin{equation}
B_1(\rho)=
\int\dot{\omega}_1(z)\omega_2(z+\rho)dz, \quad
B_2(\rho)=
\int\dot{\omega}_2(z)\omega_1(z-\rho)dz,
\label{111}
\end{equation} and
$B_1(\rho)+B_2(\rho)=1$.
\end{theorem}

In the sequel we use the following notation (as usual $x\in
\mathbb{R}$, $t\in \mathbb{R}^+$):
\begin{gather*}
u_1=u_1(x,t,\varepsilon), \quad B_i=B_i(\rho), \quad
 \varphi_i=\varphi_i(t,\varepsilon), \\
\theta_{i}=\theta(\varphi_i-x), \quad
 \delta_{i}=\delta(\varphi_i-x), \qquad i=1,2,
 \\
\rho=\frac{\varphi_2(t,\varepsilon)-\varphi_1(t,\varepsilon)}{\varepsilon},
\end{gather*}where $H$ is the Heaviside function and $\delta$ is the Dirac
distribution.

The following theorem is analogue to a special case of the main
result in \cite{DM}. We use a much simpler approach and propose two
possible solutions. The first solution is given in Theorem
\ref{t**}, and the second in Theorem \ref{tjul861}. The approach
used in Theorem \ref{t**} can be used for arbitrary continuous
initial data \cite{DM4}. Also, Theorem \ref{t**} represents
motivation for Theorem \ref{tjul861}. The difference between Theorem
\ref{t**} and Theorem \ref{tjul861} is explained at the end of the
section.

\begin{theorem} \label{t**}
The weak asymptotic solution of problem \eqref{avg715}, \eqref{avg915}
has the form
\begin{equation} \label{avg725}
\begin{aligned}
u_{\varepsilon}(x,t)&=u^0_0+\left(u_1(x,t,\varepsilon)-u^0_0\right)\omega_{1}
(\frac{\varphi_1(t,\varepsilon)-x}{\varepsilon})\\
&\quad + \left(U-u_1(x,t,\varepsilon)\right)
\omega_{2}(\frac{\varphi_2(t,\varepsilon)-x}{\varepsilon}),
\end{aligned}
\end{equation}
where $\omega_i \in C^\infty_0(\mathbb{R})$, $i=1,2$, satisfy the
conditions from Theorem \ref{th**}


The functions $\varphi_{i}(t,\varepsilon)$, $t\in \mathbb{R}^+$, $i=1,2$,
are given by
\begin{gather}
\varphi_{1}(t,\varepsilon)=\int_0^t(2u_0^0B_2(\rho)+2UB_1(\rho))dt'+a_1+\varepsilon
A \frac{a_1+a_2}{2},
\label{avg1835}\\
\varphi_2(t,\varepsilon)=\int_0^t
(2u_0^0B_1(\rho)+2UB_2(\rho))dt'+a_2-\varepsilon A \frac{a_1+a_2}{2},
\label{avg1845}
\end{gather}
for  constant $A$ which is large enough.

The function $\rho=\rho(\tau)=\rho(\tau(t,\varepsilon))$ appearing in the
previous formulas is the (global) solution of Cauchy problem
\begin{equation}
\rho_{\tau}=1-2B_1(\rho), \quad
 \lim_{\tau\to -\infty}\frac{\rho}{\tau}=1, \label{no3oct}
\end{equation}
and
$$
\tau=\frac{2Ut+a_2-2u_0^0t-a_1}{\varepsilon}.
$$
For each $\varepsilon>0$, the function $u_1(x,t,\varepsilon)$ is
defined as
$$
u_1(x,t,\varepsilon)=u_0(x_0(x,t,\varepsilon))
$$
where $x_0$ is the inverse
function to the function $x=x(x_0,t,\varepsilon)$, $t>0$, $\varepsilon>0$, of
the ``new characteristics'' defined trough the Cauchy problem
\begin{equation}
\begin{gathered}
\dot x=2u_1(x,t,\varepsilon)(B_2(\rho)-B_1(\rho))+\left(2U+2u_0^0\right)B_1(\rho),\\
\dot u_1=0, \\
u_1(0)=u_0(x_0),  \quad x(0)=x_0+\varepsilon A(x_0-\frac{a_1+a_2}{2}), \quad
x_0 \in [a_2,a_1].
\end{gathered} \label{avg45}
\end{equation}
\end{theorem}

\begin{proof}On the beginning, note that the distributional limit of
$\omega_{i}(\frac{\varphi_i-x}{\varepsilon})$ is the Heaviside
function $H_i=H_i(\varphi_i-x)$, $i=1,2$. Having this in mind, we
have after substituting (\ref{avg725}) into \eqref{avg715} and using
formula (\ref{an0}):
\begin{align*}
&\Big[u^0_0+\left(u_1-u^0_0\right)H_1+
\left(U-u_1\right)H_2\Big]_t\\
&+\Big\{(u_0^0)^2+\left[
(u_1-u_0^0)^2+2u_0^0(u_1-u_0^0)+2(u_1-u_0^0)(U-u_1)B_1(\rho)\right]H_1
 \\
&\qquad\qquad+\left[(U-u_1)^2+2u_0^0(U-u_1)+2(u_1-u_0^0)(U-u_1)B_2(\rho)
\right]H_2\Big\}_x\\
&=\mathcal{O}_{\mathcal{D}'}(\varepsilon), \quad \text{as } \varepsilon\to 0,
\end{align*}
where (see Theorem \ref{th**})
\begin{equation}
\rho=\frac{\varphi_2-\varphi_1}{\varepsilon}. \label{no2oct}
\end{equation}
After finding derivative in the previous expression (recall that
$\delta_i=-\frac{d}{dx}H_i$, $i=1,2$) and collecting terms
multiplying $H_i$ and $\delta_i$, we have (below we also use
$B_2+B_1=1$):
\begin{equation}
\begin{aligned}
&\Big[ \frac{\partial u_1}{\partial t}+\left(2u_1(B_2-B_1)
 + \left(2U+2u_0^0\right)B_1\right) \frac{\partial u_1}{\partial x}\Big]H_1 \\
&+\Big[-\frac{\partial u_1}{\partial t}-\left(2u_1(B_2-B_1)+\left(2u^0_0+2U_0 \right)
 B_1 \right)\frac{\partial u_1}{\partial x}]H_2 \\
&+\left[(u_1-u^0_0)\varphi_{1t}-2u^0_0(u_1-u^0_0)-(u_1-u^0_0)^2-2(u_1-u^0_0)
 (U-u_1)B_1 \right]\delta_1 \\
&+\left[(U-u_1)\varphi_{2t}-2u_0^0(U-u_1)-(U-u_1)^2-2(u_1-u^0_0)(U-u_1)B_2
\right]\delta_2\\
&=\mathcal{O}_{\mathcal{D}'}(\varepsilon).
\end{aligned}\label{avg765}
\end{equation}
 We rewrite the previous expression in the following manner
(we use
$M\theta_{1\varepsilon}+N\theta_{2\varepsilon}=(M+N)\theta_{1\varepsilon}
+N(\theta_{2\varepsilon}-\theta_{1\varepsilon})$):
\begin{equation}\label{sep2615}
\begin{aligned}
&\Big[ \frac{\partial u_1}{\partial t}+
\left(2u_1(B_2-B_1)+\left(2U+2u_0^0\right)B_1\right)\frac{\partial
u_1}{\partial x} \Big]\left(H_1-H_2\right)\\
&+\left[(u_1-u^0_0)\varphi_{1t}-2u^0_0(u_1-u^0_0)-(u_1-u^0_0)^2-2
 (u_1-u^0_0)(U-u_1)B_1 \right]\delta_1\\
&+\left[(U-u_1)\varphi_{2t}-2u_0^0(U-u_1)-(U-u_1)^2-2(u_1-u^0_0)(U-u_1)B_2
\right]\delta_2\\
&=\mathcal{O}_{\mathcal{D}'}(\varepsilon).
\end{aligned}
\end{equation}
We equate with zero coefficient multiplying $H_1-H_2$. We put
\begin{equation}
\frac{\partial u_1}{\partial t}+
\left(2u_1(B_2-B_1)+\left(2U+2u_0^0\right)B_1\right)\frac{\partial
u_1}{\partial x}=0. \label{1}
\end{equation}
We will prove that the last equation has continuous piecewise smooth
solution on $\mathbb{R}^+ \times \mathbb{R}$ for the initial condition:
\begin{equation}
u_{1}(x,0,\varepsilon)=-Kx+b, \ \ x\in [a_2,a_1]. \label{2}
\end{equation}
System of characteristics to equation (\ref{1}) has the form (those
are ``almost'' equations of ``new characteristics'' from (\ref{avg45});
see (\ref{no10oct*})):
\begin{equation} \label{no10oct}
\begin{gathered}
\frac{dx}{dt}=2u_1(B_2-B_1)+(2U+2u_0^0)B_1, \quad x(0)=x_0\in [a_2,a_1],
\\
\dot{u}_1=0, \quad u_1(0)=u_0(x_0)=-Kx_0+b. \quad
\text{(since $x_0\in[a_2,a_1]$)}
\end{gathered}
\end{equation}
To prove global solvability of (\ref{1}), (\ref{2}) it is
sufficient to prove the global existence of the inverse function
$x_0=x_0(x,t,\varepsilon)$ to the function $x$ defined by previous
equations of characteristics (\ref{no10oct}). It appears that it is
much easier to accomplish this if we perturb initial data for $x$ in
the previous system for a parameter of order $\varepsilon$. More precisely,
instead of (\ref{no10oct}) we shall consider the following system
(the same is done in \cite{DM}):
\begin{equation}
\begin{gathered}
\frac{dx}{dt}=2u_1(B_2-B_1)+(2U+2u_0^0)B_1, \quad
x(0)=x_0+\varepsilon A \big(x_0-\frac{a_1+a_2}{2}\big),\\
\dot{u}_1=0, \quad u_1(0)=u_0(x_0)=-Kx_0+b, \quad
 x_0\in [a_2,a_1].
\end{gathered} \label{no10oct*}
\end{equation}
Since our initial data are continuous, such perturbation will
change exact solution of (\ref{1}), (\ref{2}) for
$\mathcal{O}_{\mathcal{D}'}(\varepsilon)$.

However, before we are able to prove existence of the inverse function
$x_0$ to the function $x$ defined by (\ref{no10oct*}), we need to
define equations for the functions $\varphi_i$, $i=1,2$, and prove
that $\rho$ given by (\ref{no2oct}) satisfies (\ref{no3oct}).


As the characteristics emanating from $a_2$ and $a_1$ we have for
$\varphi_2$ and $\varphi_1$:
\begin{gather}
\varphi_{1t}=2u_0^0B_2+2UB_1, \quad \varphi_1(0,\varepsilon)=a_1+\varepsilon A
\frac{a_1-a_2}{2},
\label{avg735}
\\
\varphi_{2t}=2u_0^0B_1+2UB_2, \quad \varphi_2(0,\varepsilon)=a_2-\varepsilon A
\frac{a_1-a_2}{2}.
\label{avg745}
\end{gather}
Now, we are interested in the behavior of
$\varphi_2-\varphi_1$. As usual \cite{dan,DM,DSH}, we introduce the fast
variable
$$
\tau=\frac{\varphi_{20}(t)-\varphi_{10}(t)}{\varepsilon},
$$
where $\varphi_{10}$ and $\varphi_{20}$ are standard characteristics
emanating from $a_1$ and $a_2$ respectively:
\begin{gather*}
\varphi_{10}(t)=2u_0^0 t+a_1,\\
\varphi_{20}(t)=2U t+a_2.
\end{gather*}
Note that $\tau$ can be considered independent on $t$ thanks to
small parameter $\varepsilon$. Also, from the equation
$\varphi_{10}(t)=\varphi_{20}(t)$ we can compute the moment of
blowing up of the classical solution (since the choice of initial
data provides admissible weak solution of problem \eqref{avg715},
\eqref{avg915} to lie in algebra $\mathcal{L}\{ 1, \theta(x-ct)\}$
where $c=U+u_0^0$ (Rankine-Hugoniot condition); see also Figure
\ref{fig*trnd1}). We denote the moment of blowing up of the
classical solution by $t^*$ and appropriate space point by $x^*$.
We easily infer that
\begin{equation}
t^*=\frac{a_1-a_2}{U-u_0^0}, \quad x^*=2Ut^*+a_2=2u_0^0t^*+a_1.
 \label{trnd1}
\end{equation}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig1}
\end{center}
\caption{Standard characteristics for \eqref{avg715}, \eqref{avg915}. Dotted point in
$(t,x)$ plane is $(t^*,x^*)$.}\label{fig*trnd1}
\end{figure}


Note that before $t^*$ we have $\tau\to -\infty$ as $\varepsilon\to 0$
and for $t>t^*$ we have $\tau\to \infty$ as $\varepsilon\to 0$. So,
variable $\tau$ can be understood as indicator of state of the
solution. When it is large toward $-\infty$ we have classical
solution of the problem (since $t<t^*$), and when $\tau$ is large
toward $+\infty$ the classical solution blew up and we have only
weak solution to the problem (since $t>t^*$).

Subtracting \eqref{avg735} from \eqref{avg745}. We
have
\begin{equation}
(\varphi_2-\varphi_1)_t=2(U-u_0^0)(1-2B_1(\rho)). \label{no1oct}
\end{equation}
Since $\tau=\frac{2Ut+a_2-2u_0^0t-a_1}{\varepsilon}$ we have
\begin{gather*}
(\varphi_2-\varphi_1)_t=(\varepsilon\rho)_t=\varepsilon\rho_{\tau}\tau_t=2(U-u_0^0)\rho_{\tau},
\end{gather*} combining this with (\ref{no1oct}) we have:
\begin{equation*}
\rho_{\tau}=1-2B_1(\rho), \quad \lim_{\tau\to -\infty}\frac{\rho}{\tau}=1.
\end{equation*}
We explain the condition $\lim_{\tau\to-\infty}\frac{\rho}{\tau}=1$.
We have from \eqref{avg735} and \eqref{avg745}
$$
\frac{\rho}{\tau}=\frac{\int_0^t2(U-u_0^0)(B_2-B_1)dt'
+a_2-a_1}{2(U-u_0^0)t+a_2-a_1}.
$$
Putting $t=0$ in the previous relation we see that
\begin{equation}
\frac{\rho}{\tau}\Big|_{t=0}=1.
\label{no12oct}
\end{equation}
When we let $\varepsilon\to 0$ when $t=0$ we have $\tau\to -\infty$.
Therefore, from (\ref{no12oct}),
$$
\frac{\rho}{\tau}\Big|_{\tau\to-\infty}=1.
$$
This relation practically means that new characteristics
emanating from $a_i$, $i=1,2$, coincides at least in the initial
moment with standard characteristics up to some small parameter
$\varepsilon$. Still, since $\tau\to -\infty$, i.e. $B_1\to 0$, for every
$t<t^*$ we see from \eqref{avg735} and \eqref{avg745} that new
characteristics coincides with standard ones for every $t<t^*$ up to
some small parameter $\varepsilon$.

Thus, we have proved that $\rho$ given in (\ref{no2oct}) indeed
satisfies (\ref{no3oct}). From the classical ODE theory one sees
that problem (\ref{no3oct}) has global solution $\rho$ such that
$\rho\to \rho_0$ as $\tau \to +\infty$ where $\rho_0$ is constant
such that $B_1(\rho_0)=B_2(\rho_0)=1/2$ (more precisely, $\rho_0$ is
stationary solution to equation from (\ref{no3oct})).

Replacing $\rho=\rho_0$ in the expressions for $\varphi_{it}$
(expressions \eqref{avg735}, \eqref{avg745}) and using
$B_i(\rho_0)=1/2$, $i=1,2$, we obtain that in the limit:
\begin{equation}
\varphi_{1t}=\varphi_{2t}=U+u_0^0, \label{oslo1}
\end{equation}
which means that after the interaction the points $a_1$ and
$a_2$ continue to move with the same velocity (which, as expected,
coincides with the velocity given by Rankine-Hugoniot condition).

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig2}
\end{center}
\caption{System of characteristics for $u_{\varepsilon}$. The points
$a_2-\varepsilon A\frac{a_1+a_2}{2}$ and $a_1+\varepsilon
A\frac{a_1+a_2}{2}$ are dotted. ``New characteristics'' emanate from
the interval $[a_2-\varepsilon A\frac{a_1+a_2}{2},a_1+\varepsilon
A\frac{a_1+a_2}{2}]$.}\label{fig1}
\end{figure}


As we have mentioned earlier, problem (\ref{1}), (\ref{2}) is
globally solvable if the (new) characteristics defined trough
(\ref{no10oct*}) do not intersect. To prove that we will use the
inverse function theorem. We will prove that for every $t$ we have
$\frac{\partial x}{\partial x_0}>0$ which means that for every $x=x(x_0,t)$,
$x_0\in [a_2,a_1]$, we have unique $x_0=x_0(x,t,\varepsilon)$ and we can
write $u_1(x,t)=u_0(x_0(x,t,\varepsilon))$.


Since $u_1(x_0,0,\varepsilon)=-Kx_0+b$ (see (\ref{2})), from
(\ref{no10oct*}) we conclude:
\begin{equation}\label{123}
x=\int_0^t \left(-2Kx_0+b)(B_2-B_1)+(2U+2u_0^0)B_1
\right)dt'+x_0+\varepsilon A \big(x_0-\frac{a_1+a_2}{2}\big).
\end{equation}
Finding derivative of (\ref{123}) in $x_0$ we obtain (we use
$B_2+B_1=1$):
\begin{equation}
\frac{\partial x}{\partial x_0}=1+\varepsilon A-2K\int_0^t(1-2B_1)dt'. \label{avg835}
\end{equation}
For $t\in [0,t^*]$ we have (below we use $1-2Kt^*=0$)
\begin{align*}
\frac{\partial x}{\partial x_0}
&=1+\varepsilon A-2K\int_0^t(1-2B_1)dt'\\
&\geq 1+\varepsilon A-2K\int_0^{t^*}(1-2B_1)dt'\\
&=\varepsilon A+4\int_0^{t^*}B_1dt'>0
\end{align*}
since $B_1>0$.
So, everything is correct for $t\in [0,t^*]$. To see what is
happening for $t>t^*$, initially we estimate $\rho_{\tau}$ when
$\tau\to \infty$.


 From (\ref{no3oct}) we have (we use Taylor expansion):
\begin{equation}
\rho_{\tau}=1-2B_1(\rho)=-2(\rho-\rho_0)B_1'(\tilde{\rho}),
\label{trnd2}
\end{equation}
for some $\tilde{\rho}$ belonging to the interval with
ends in $\rho$ and $\rho_0$. From here we see:
$$
\rho-\rho_0=C{\rm exp}(\int_{\tau_0}^{\tau}-
2B_1'(\tilde{\rho})d\tau')=C{\rm
exp}((-\tau+\tau_0)2B_1'(\tilde{\rho}_1))
$$
for some fixed
$\rho_0\in \mathbb{R}$ and
$\tilde{\rho}_1\in(\rho(\tau_0),\rho(\tau))\subset [\rho(\tau_0),\rho_0]$.
We remind that  $B_1'(\tilde{\rho}_1))\geq c >0$, for some constant $c$,
since $B_1$ is increasing function and $\tilde{\rho}_1$ belongs to
the compact interval $[\rho(\tau_0),\rho_0] $, letting
$\tau\to \infty$ we conclude that for any $N\in {\bf N}$
$$
\rho-\rho_0=\mathcal{O}(1/\tau^N), \quad  \tau\to \infty.
$$
This in turn means that for $t>t^*$ we have
\begin{equation} \rho-\rho_0=\mathcal{O}(\varepsilon^N), \quad \varepsilon\to 0.
\label{trnd3}
\end{equation}
Now, we can prove resoluteness of problem (\ref{1}), (\ref{2}) for
$t>t^*$. We have
\begin{equation} \label{no60oct}
\begin{aligned}
\frac{\partial x}{\partial x_0}
&= 1+\varepsilon A-2K\int_0^t(1-2B_1)dt'\\
&= 1+\varepsilon A-2K\int_0^{t^*}(1-2B_1)dt'- 2K\int_{t^*}^{t}(1-2B_1)dt'\\
&= \varepsilon A+4\int_0^{t^*}B_1dt'-2K\int_{t^*}^{t}(1-2B_1)dt'>\varepsilon
A-2K\int_{t^*}^{t}(1-2B_1)dt'.
\end{aligned}
\end{equation}
Recall that
$$
B_1=B_1(\rho(\tau))=B_1\Big(\rho\big(\frac{\psi_0(t)}{\varepsilon}\big)\Big),
$$
where
$\psi_0(t)=2(U-u_0^0)t+(a_2-a_1)$.
Consider the last term in expression (\ref{no60oct}):
\begin{align*}
2K\int_{t^*}^{t}(1-2B_1)dt'
&=2K\int_{t^*}^{t}\Big(1-2B_1\Big(\rho\big(\frac{\psi_0(t')}{\varepsilon}
\big)\Big)\Big)dt'\\
&=2K\varepsilon\int_0^{\frac{\psi_0(t)}{\varepsilon}}(1-2B_1(\rho(z)))dz<\varepsilon 2KC,\\
&\left(
\begin{array}{cc} \frac{\psi_0(t')}{\varepsilon}= z \implies&
(u-u_0^0)dt'=\varepsilon dz\\
t^*<t'<t\implies& 0<z<\frac{\psi_0(t)}{\varepsilon}
\end{array}
\right)
\end{align*}
where
$$
C=\int_0^{\infty}(1-2B_1(\rho(z)))dz<\infty,
$$
since $1-2B_1(\rho(z))=\mathcal{O}(z^{-N})$, $z\to \infty$ and
$N\in {\bf N}$ arbitrary (see (\ref{trnd2}) and (\ref{trnd3}) for
this).
Therefore, for $A$ large enough (more precisely for
$A>2K \int_0^{\infty}(1-2B_1(\rho(z)))dz$) we have
$\frac{\partial x}{\partial x_0}>0$ what we wanted to prove.

Now, we return to (\ref{sep2615}). Taking into account (\ref{1}),
from (\ref{sep2615}) we have
\begin{equation}\label{*}
\begin{aligned}
&\left[(u_1-u^0_0)\varphi_{1t}-u^0_0-u_1-2(U-u_1)B_1 \right]\delta_1 \\
&+\left[(U-u_1)\varphi_{2t}-2u_0^0-U+u_1-2(u_1-u^0_0)B_2
\right]\delta_2\\
&=\mathcal{O}_{\mathcal{D}'}(\varepsilon).
\end{aligned}
\end{equation}
After substituting values for $\varphi_{it}$, $i=1,2$, into the
last expression we have
\begin{equation}\label{sep2715}
(B_2-B_1)(u_1-u^0_0)(u_1+u^0_0)\delta_1+
(B_2-B_1)(U-u_1)(U+u_1)\delta_2
=\mathcal{O}_{\mathcal{D}'}(\varepsilon).
\end{equation} We have from the definition of the Dirac distribution, after multiplying (\ref{sep2715}) by $\eta\in
C^1_0(\mathbb{R})$ and integrating over $\mathbb{R}$,
\begin{align*}
&(B_2-B_1)(U-u_1(\varphi_2,t))(U+u_1(\varphi_2,t))
\eta(\varphi_2)\\
&+(B_2-B_1)(u_1(\varphi_1,t)-u^0_0)(u_1(\varphi_1,t)+u^0_0)
\eta(\varphi_1)dx=\mathcal{O}(\varepsilon),
\end{align*}which is correct
 since $u_1\equiv U$ for $x\in (-\infty, \varphi_2]$ and $u_1\equiv u_0^0$
for $x\in [\varphi_1,\infty)$. This proves (\ref{sep2715}) and
finishes the proof of the theorem.
\end{proof}

The following corollary is obvious. It claims that the weak
asymptotic solution defined in arbitrary of the previous theorems
tends to the shock wave with the states $U$ on the left and $u^0_0$
on the right (see (\ref{oslo1}) to remove ambiguities).

\begin{corollary} \label{coro5}
With the notation from the previous theorems,
for $t>t^*$ the weak asymptotic solution $u_{\varepsilon}$ to problem
\eqref{avg715}, \eqref{avg915} we have for every fixed $t>0$:
\begin{equation}
u_{\varepsilon}(x,t)\rightharpoonup \begin{cases}
U,   & x<(U+u_0^0)(t-t^*)+x^*,\\
U_0, & x>(U+u_0^0)(t-t^*)+x^*,
\end{cases}
\label{ostr71}
\end{equation}
where $\rightharpoonup$ means convergence in the weak sense with
respect to the real variable.
\end{corollary}

The following theorem is motivated by the previous one and based on
the following observation. Once the shock wave is formed, it
continuous to move according to Rankine-Hugoniot conditions and it
does not change its shape along entire time axis. Therefore, the
linear equation
\begin{equation}
\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0, \quad
 c=U+u_0^0, \label{jul161}
\end{equation}
and equation \eqref{avg715} with the same initial condition
\begin{equation*}
u|_{t=0}=\begin{cases}
U,  & x<0,\\
u_0^0, & x\geq 0,
\end{cases}
\end{equation*}
will have the same solutions.
The question is: If we do not have Riemann initial conditions, how
to replace \eqref{avg715} by (\ref{jul161}) in domains where we
can do it (i.e. after shock wave formation) without loosing
properties of the solution of original problem. As we will see,
Theorem \ref{tjul861} will prove that one of the possible answer
is to describe passage from \eqref{avg715} to (\ref{jul161})
smoothly in $t\in \mathbb{R}^+$.
In the following theorem the notions and notation are the same
as in the previous theorem.

\begin{theorem} \label{tjul861}
The weak asymptotic solution $u_{\varepsilon}$, $\varepsilon>0$,
to Cauchy problem \eqref{avg715}, \eqref{avg915} is given by
\begin{equation}
u_{\varepsilon}(x,t)=\hat{u}(x_0(x,t,\varepsilon)), \label{avg2956}
\end{equation}
where $x_0$ is inverse function to the function $x=x(x_0,t,\varepsilon)$, $t>0$,
$\varepsilon>0$, of 'new characteristics' defined trough the Cauchy problem
\begin{equation}
\begin{gathered}
\dot{x}=f'(u_{\varepsilon})(B_2(\rho)-B_1(\rho))+cB_1(\rho), \\
 x(0)=x_0+\varepsilon A \Big(x_0-\frac{a_1+a_2}{\varepsilon} \Big),\\
\dot{u_{\varepsilon}}=0, \quad u_{\varepsilon}(0)=\hat{u}(x_0), \quad x_0\in \mathbb{R},
\end{gathered} \label{jul661}
\end{equation}
where $A$ is large enough, the functions $B_1$ and $B_2$ are defined
in Theorem \ref{th**}, and constant $c$ such that
$$
\frac{c}{2}=U+u_0^0,
$$
and
$\rho=\rho(\psi_0(t)/\varepsilon)$ is the solution of the Cauchy
problem
\begin{gather*}
\rho_{\tau}=1-2B_1(\rho), \quad \lim_{\tau\to-\infty}\frac{\rho}{\tau}=1.
\end{gather*}
\end{theorem}

\begin{proof}
Consider the family of Cauchy problems (recall that $B_i=B_i(\rho)$, $i=1,2$):
\begin{equation}
\frac{\partial u_{\varepsilon}}{\partial t}+\left(2u_{\varepsilon}(B_2-B_1)+2B_1
\left(U+u_0^0\right)\right)\frac{\partial u_{\varepsilon}}{\partial x}=0, \quad
 x\in \mathbb{R}, \quad t\in \mathbb{R}^+,
\label{jul162}
\end{equation}
Note that the ``new characteristics'' given by (\ref{jul661})
correspond to Cauchy problem (\ref{jul162}), \eqref{avg915} up to
$\mathcal{O}(\varepsilon)$ (since we have perturbed initial data for the
characteristic $x$ in (\ref{jul661})). Since initial conditions to
equations \eqref{avg715} and (\ref{jul162}) are the same, it is
enough to prove that the solution to Cauchy problem
(\ref{jul162}), \eqref{avg915} (possibly perturbed by term of
order $\varepsilon$), represents the weak asymptotic solution to
\eqref{avg715}, \eqref{avg915}.


First, we have to solve Cauchy problem (\ref{jul162}),
\eqref{avg915}. We use standard method of characteristics. The
characteristics of given Cauchy problem are
\begin{gather*}
\dot{x}=2u_{\varepsilon}(x,t)(B_2-B_1)+ 2B_1(U+u_0^0),\quad  x(0)=x_0,\\
\dot{u}_{\varepsilon}=0, \quad  u_{\varepsilon}(0)=u_0(x_0).
\end{gather*}
We perturb initial data for $x$ i.e. we put
\begin{equation}
x(0)=x_0+\varepsilon A\big(x_0 -\frac{a_1+a_2}{2}\big)
\label{trnd4}
\end{equation}
and use $\dot{u}_{\varepsilon}=0\Rightarrow u_{\varepsilon}(x,t)=u_0(x_0)$:
\begin{equation}
\dot{x}=\begin{cases}
B_2 U+B_1u_0^0, & x_0<a_2,\\
(-2Kx_0+2b)(B_2-B_1)+ B_1\left(2U+2u_0^0 \right), & x_0\in [a_1,a_1]\\
B_1u_0^0+B_2U, & x_0>a_1.
\end{cases}
\label{jul862}
\end{equation}
After integrating from $0$ to $t$ and finding derivative in $x_0$
we have (using (\ref{trnd4})),
\begin{equation}
\frac{\partial x}{\partial x_0}=\begin{cases}
1, & x_0<a_2,\\
1+\varepsilon A-2K\int_0^t(B_2-B_1)dt'=\frac{\varphi_1-\varphi_2}{a_1-a_2},
&\ x_0\in [a_2,a_1],\\
1, & x_0>a_1.
\end{cases}
\label{avg2626}
\end{equation}
According to the part of the previous theorem between formulas
(\ref{avg835}) and (\ref{*}),
we see that $\frac{\partial x}{\partial x_0}>0$ for every $\varepsilon>0$.
According to inverse function theorem, this means that characteristics
(\ref{jul862}) never intersects, i.e.
we can define solution of (\ref{jul162}), \eqref{avg915} along
characteristics for every $\varepsilon>0$:
\begin{equation}
u_{\varepsilon}(x,t)=u_0(x_0(x,t,\varepsilon)),
\label{jul961}
\end{equation}
where $x_0$ is inverse function to the function $x$ defined
trough (\ref{jul862}).

Now, we have to prove that family $u_{\varepsilon}$, $\varepsilon>0$, of solutions
to (\ref{jul162}), \eqref{avg915} defines weak asymptotic solution
to \eqref{avg715}, \eqref{avg915}. More precisely, we have to prove
that that for the solution $u_{\varepsilon}$ of problem (\ref{jul162}),
\eqref{avg915} it holds
\begin{equation}
\frac{\partial u_{\varepsilon}}{\partial t}+\frac{\partial u^2_{\varepsilon}}{\partial x}
=\mathcal{O}_{\mathcal{D}'}(\varepsilon).
\label{no40oct}
\end{equation}
We have
\begin{equation}
\begin{aligned}
\frac{\partial u_{\varepsilon}}{\partial t}+\frac{\partial u^2_{\varepsilon}}{\partial x}
&=\frac{\partial u_{\varepsilon}}{\partial t}+2u_{\varepsilon}\frac{\partial u_{\varepsilon}}{\partial x}\\
&=\frac{\partial u_{\varepsilon}}{\partial t}+\left(2u_{\varepsilon}(B_2-B_1)+ B_1
\left(2U+2u_0^0\right) \right)\frac{\partial u_{\varepsilon}}{\partial x}\\
&\quad - B_1 \left(2U+2u_0^0-4u_{\varepsilon}\right)\frac{\partial u_{\varepsilon}}{\partial x}\\
&= \mathcal{ O}_{\mathcal{D}'}(\varepsilon).
\end{aligned}\label{jul164}
\end{equation}
Since we assumed that $u_{\varepsilon}$ is the solution
(\ref{jul162}), \eqref{avg915},   from (\ref{jul164}) we have
\begin{align}
 B_1\left(2U+2u_0^0-4u_{\varepsilon}\right)\frac{\partial u_{\varepsilon}}{\partial x}
 =\mathcal{O}_{\mathcal{D}'}(\varepsilon).
 \label{oslo2}
\end{align}
Note that we have $|\rho B_1|<\infty$ for every
$\tau\in \mathbb{R}$. Namely,
\begin{equation}\label{oslo3}
\begin{gathered}
|\rho B_1(\rho)|\to 0 \text{ as $\tau\to -\infty $  since in this case  }
B_1(\rho(\tau))\sim B_1(\tau)\sim\frac{1}{\tau^N}\sim\frac{1}{\rho^N},
\\
|\rho B_1(\rho)|\to \rho_0B_1(\rho_0) \text{ as $\tau\to \infty$  since
in this case }
 \rho\to \rho_0.
\end{gathered}
\end{equation}
Knowing this, we multiply (\ref{oslo2}) with $\eta\in C^1_0(\mathbb{R})$,
integrate and apply partial integration (we bear in mind that
$u_{\varepsilon}\equiv U$ for $x\leq \varphi_2$ and $u_{\varepsilon}\equiv u_0^0$
for $x\geq \varphi_1$):
\begin{align*}
&\int B_1 \left(2U+2u_0^0-2u_{\varepsilon} \right)u_{\varepsilon}\eta'(x)dx\\
&=B_1\Big(\int_{-\infty}^{\varphi_2} \left(2U+2u_0^0-2u_{\varepsilon} \right)
 u_{\varepsilon}\eta'(x)dx+\int_{\varphi_2}^{\varphi_1}
 \left(2U+2u_0^0-2u_{\varepsilon} \right)u_{\varepsilon}\eta'(x)dx\\
&\quad + \int_{\varphi_1}^{\infty} \left(2U+2u_0^0-2u_{\varepsilon}
\right)u_{\varepsilon}\eta'(x)dx\Big)\\
&= 2u_0^0UB_1\eta(\varphi_2)
 + \varepsilon \rho B_1(\rho)\frac{1}{\varphi_2-\varphi_1}\int_{\varphi_2}^{\varphi_1}
  \left(2U+2u_0^0-2u_{\varepsilon} \right)u_{\varepsilon} dx
 -2u_0^0 U B_1 \eta(\varphi_1)\\
&=2U u_0^0 \varepsilon \rho B_1\frac{\eta(\varphi_2)
 -\eta(\varphi_1)}{\varphi_2-\varphi_1} +\mathcal{O}(\varepsilon)\\
&=\mathcal{O}(\varepsilon).
\end{align*}
which proves (\ref{jul164}) and concludes the proof of the theorem.
\end{proof}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig3}
\end{center}
\caption{System of characteristics for $u_{\varepsilon}$ defined in Theorem
\ref{tjul861}. The points $a_1+\varepsilon A\frac{a_1+a_2}{2}$ and
$a_2-\varepsilon A\frac{a_1+a_2}{2}$ are dotted on the $x$
axis.}\label{fig3}
\end{figure}

Before we consider the system of equations we will explain
difference between weak asymptotic solution of problem
\eqref{avg715}, \eqref{avg915} we have constructed in Theorem
\ref{t**} and exact solution of (\ref{jul162}), \eqref{avg915},
perturbed possibly for term of order $\varepsilon$, which is, at the same
time, weak asymptotic solution to \eqref{avg715}, \eqref{avg915}.
The solution of (\ref{jul162}), \eqref{avg915} is constructed by
standard method of characteristics. The characteristics are well
defined, i.e. they do not mutually intersect. In other words,
solution of (\ref{jul162}), \eqref{avg915} forms continuous group of
transformations (see Figure 2) while solution constructed in Theorem
\ref{t**} forms only continuous semigroup of transformations since
in that case characteristics intersects along lines $\varphi_i$,
$i=1,2$ (see Figure 1).

\section{Application of the method to a system of PDEs}

In this section we consider the system
\begin{equation}
\begin{gathered}
u_t+(\frac{1}{2}u^2)_x=0,\\
v_t+(uv)_x=0,
\end{gathered}
\label{jul165}
\end{equation}
with Riemann initial data
\begin{equation}
\begin{gathered}
u|_{t=0}=u_0(x)=\begin{cases}
U, & x<0\\
u_0^0, & x\geq 0
\end{cases}\\
v|_{t=0}=v_0(x)=\begin{cases}
v_0, & x<0\\
v_1, & x\geq 0.
\end{cases}
\end{gathered}
\label{jul166}
\end{equation}
This non-strictly hyperbolic system of conservation laws arises
from pressureless gas dynamics and it is intensively investigated
in many papers (see the Introduction).
Here, we will demonstrate how delta shock wave
naturally arises if we ``smooth'' a little bit our Riemann initial
data.

In the sequel, all the notions and notation are the same as in
the previous section.
To solve problem \eqref{jul165}, \eqref{jul166} we use the
following procedure.

On the first step we replace initial data \eqref{jul166} by
perturbed continuous initial data:
\begin{equation}
\begin{gathered}
u_{\varepsilon}|_{t=0}=u_{0\varepsilon}(x)=\begin{cases}
U, & x\leq a_2=-\varepsilon^{1/2}\\
-\frac{U-u_0^0}{2\varepsilon^{1/2}}x+\frac{U+u_0^0}{2},
 & -\varepsilon^{1/2}=a_2 < x< a_1=\varepsilon^{1/2},\\
u_0^0, & x\geq \varepsilon^{1/2}
\end{cases} \\
v_{\varepsilon}|_{t=0}=v_{0\varepsilon}(x,t)\begin{cases}
v_0, & x\leq a_2= -\varepsilon^{1/2}\\
-\frac{v_1-v_0}{2\varepsilon^{1/2}}x+\frac{v_1+v_0}{2},
  & -\varepsilon^{1/2}=a_2< x< a_1=\varepsilon^{1/2},\\
v_1, & x\geq a_1=\varepsilon^{1/2}.
\end{cases}
\end{gathered}
\label{jul167}
\end{equation}
Note that in this case gradient catastrophe (blow up of classical solution)
will happen in the moment
$$
t^*=\frac{2\varepsilon^{1/2}}{U-u_0^0}.
$$
Next, as in the previous section we put
\begin{gather*}
\varphi_1(t,\varepsilon)=\int_0^t(u_0^0B_2+UB_1)dt'+\varepsilon^{1/2}+\varepsilon^{3/2} A ,\\
\varphi_2(t,\varepsilon)=\int_0^t(U B_2+u_0^0 B_1)dt'-\varepsilon^{1/2}-\varepsilon^{3/2}A,
\end{gather*}
while $B_i=B_i(\rho)$, $i=1,2$, are defined by Theorem \ref{th**},
$\rho=\rho(\tau)$ is defined by Cauchy problem
(\ref{no3oct}), and
$$
\tau=\frac{2(U-u_0^0)t-2\varepsilon^{1/2}}{\varepsilon}.
$$

On the second step we replace system \eqref{jul165} by the system
\begin{equation}
\begin{gathered}
u_{\varepsilon t}+\Big(\frac{1}{2}u_{\varepsilon}^2(B_2-B_1)+B_1
\left( u_0^0+ U \right)u_{\varepsilon}\Big)_x=0,\\
v_{\varepsilon t}+\left(u_{\varepsilon} v_{\varepsilon}(B_2-B_1) +B_1
v_{\varepsilon}(u_0^0+U)\right)_x=F(x,t,\varepsilon),
\end{gathered}
\label{jul168}
\end{equation}
where $u_{\varepsilon}=u_{\varepsilon}(x,t)$ and $v_{\varepsilon}=v_{\varepsilon}(x,t)$, and $F$
is function to be determined
from the condition of equivalence in the weak asymptotic sense of
systems \eqref{jul165} and (\ref{jul168}).

Since we have proved in Theorem \ref{tjul861} that the first
equation of system (\ref{jul168}) is equivalent in the weak
asymptotic sense to the first equation of \eqref{jul165} we
investigate relation between the second equation from (\ref{jul168})
and the second equation from \eqref{jul165}.

We have to determine $F$ so that for arbitrary weak asymptotic
solution $(u_\varepsilon,v_\varepsilon)$ of (\ref{jul168}) we have:
$$
v_{\varepsilon t}+(u_{\varepsilon}v_{\varepsilon})_x=\mathcal{O}_{\mathcal{
D}'}(\varepsilon^{1/2}).
$$
 From here we  have after adding and subtracting appropriate terms
and using $B_2+B_1=1$:
\begin{align*}
&v_{\varepsilon t}+\left(u_{\varepsilon} v_{\varepsilon}(B_2-B_1)+B_1 v_{\varepsilon} (u_0^0+U)
 \right)_x-F(x,t,\varepsilon)-\\
&B_1(v_{\varepsilon} u_0^0+U
v_{\varepsilon}-2u_{\varepsilon}v_{\varepsilon})_x+F(x,t,\varepsilon)=\mathcal{O}_{\mathcal{
D}'}(\varepsilon^{1/2}).
\end{align*}
 We use (\ref{jul168}) to deduce
$$
B_1(v_{\varepsilon} u_0^0+ v_{\varepsilon} U-2u_{\varepsilon}v_{\varepsilon})_x
=F(x,t,\varepsilon)+\mathcal{O}_{\mathcal{D}'}(\varepsilon^{1/2}).
$$
Now, we multiply the last expression by $\eta\in
C^1_0(\mathbb{R})$, integrate over $\mathbb{R}$ and use partial
integration
\begin{equation}
-B_1\int(u_0^0 v_{\varepsilon}+U
v_{\varepsilon}-2u_{\varepsilon}v_{\varepsilon})\eta'(x)dx=\int F\eta
dx+\mathcal{O}(\varepsilon^{1/2}). \label{jul261}
\end{equation}
Here, as usual, $F=F(x,t,\varepsilon)$.
Clearly, for $x<\varphi_2$ we have $v_{\varepsilon}\equiv v_1$ and for
$x>\varphi_1$ we have $v_{\varepsilon}\equiv v_0$.
Therefore, from (\ref{jul261}) we have
\begin{equation} \label{jul264}
\begin{aligned}
&-B_1\Big(\int_{-\infty}^{\varphi_2}(v_1 u_0^0-v_1 U)\eta'(x)dx
 +\int_{\varphi_2}^{\varphi_1}(v_{\varepsilon} u_0^0+U v_{\varepsilon}
 -2u_{\varepsilon}v_{\varepsilon})\eta'(x)dx\\
&+ \int_{\varphi_1}^{\infty} (U v_0-u_0^0 v_0) \eta'(x)dx\Big) \\
&=B_1\left((v_0 U -v_0 u_0^0)\eta(\varphi_1)+(v_1 U-v_1
u_0^0)\eta(\varphi_2)\right) \\
&\quad -B_1\int_{\varphi_2}^{\varphi_1}(u_0^0 v_{\varepsilon}+U
v_{\varepsilon}-2u_{\varepsilon}v_{\varepsilon})\eta'(x)dx\\
&=\int F\eta dx+\mathcal{O}(\varepsilon^{1/2}).
\end{aligned}
\end{equation}
We will see later (\ref{**}) that
\begin{equation}
B_1\int_{\varphi_2}^{\varphi_1}(u_0^0 v_{\varepsilon}+U
v_{\varepsilon}-2u_{\varepsilon}v_{\varepsilon})\eta'(x)dx=\mathcal{O}(\varepsilon^{1/2}).
\label{jul2610}
\end{equation}
Therefore, it follows from (\ref{jul264}) that
\begin{equation*}
B_1\left((v_0 U -v_0 u_0^0)\eta(\varphi_1)+(v_1 U-v_1
u_0^0)\eta(\varphi_2)\right)=\int F\eta
dx+\mathcal{O}(\varepsilon^{1/2}).
\end{equation*}
We rewrite the above expression as
$$
B_1(v_0+v_1)(U-u_0^0)\eta(\varphi_1)+\varepsilon\rho B_1 (v_1 U-v_1 u_0^0)
 \frac{\eta(\varphi_2)-\eta(\varphi_1)}{\varphi_2-\varphi_1}
 =\int F\eta dx+\mathcal{O}(\varepsilon^{1/2}).
$$
Recalling (\ref{oslo3}) we know
$\varepsilon\rho B_1 (v_1 U-v_1 u_0^0)\frac{\eta(\varphi_2)
-\eta(\varphi_1)}{\varphi_2-\varphi_1}=\mathcal{O}(\varepsilon)$
 which implies that unknown function $F$ should satisfy
\begin{equation}
B_1(v_0+v_1)(U-u_0^0)\eta(\varphi_1)=\int F\eta
dx+\mathcal{O}(\varepsilon^{1/2}). \label{jul265}
\end{equation}
It is obvious from here that the function $F$ should represent
regularization of the Dirac $\delta$ distribution supported in $x=\varphi_1$.
 We will choose a regularization  which will make further
computations easier. Accordingly, let
$$
F(x,t,\varepsilon)=B_1(v_0+v_1)(U-u_0^0)\frac{\kappa((\varphi_2,\varphi_1))}
{\varphi_1-\varphi_2},
$$
where $\kappa((a,b))=\kappa((a,b))(x)$ is
characteristic function of the interval $(a,b)$. We prove that
(\ref{jul265}) is satisfied for such choice of $F$:
\begin{align*}
&B_1(v_0+v_1)(U-u_0^0)\eta(\varphi_1) \\
&=B_1(v_0+v_1)(U-u_0^0)\int\frac{\kappa((\varphi_2,\varphi_1))}
 {\varphi_1-\varphi_2}\eta(x) dx+\mathcal{O}(\varepsilon^{1/2})\\
&=B_1(v_0+v_1)(U-u_0^0)\frac{1}{\varphi_1-\varphi_2}
  \int_{\varphi_2}^{\varphi_1}\eta(x)dx+\mathcal{O}(\varepsilon^{1/2})\\
&=B_1(v_0+v_1)(U-u_0^0)\frac{1}{\varphi_1-\varphi_2}
  \int_{\varphi_2}^{\varphi_1}\left(\eta(\varphi_1)
  +(x-\varphi_1)\eta'(\hat{x}))\right)dx+\mathcal{O}(\varepsilon^{1/2})\\
&=B_1(v_0+v_1)(U-u_0^0)\eta(\varphi_1)\\
&\quad +B_1(v_0+v_1)(U-u_0^0)
  \frac{1}{\varphi_1-\varphi_2}\int_{\varphi_2}^{\varphi_1}(x-\varphi_1)
  \eta'(\hat{x})dx + \mathcal{O}(\varepsilon^{1/2}).
\end{align*}
 From here it follows that
$$
B_1(v_0+v_1)(U-u_0^0)\frac{1}{\varphi_1-\varphi_2}
 \int_{\varphi_2}^{\varphi_1}(x-\varphi_1)\eta'(\hat{x})dx
=\mathcal{O}(\varepsilon^{1/2}),
$$
which is true due to (\ref{oslo3}) and since
\begin{align*}
&|B_1(v_0+v_1)(U-u_0^0)\frac{1}{\varphi_1-\varphi_2}
 \int_{\varphi_2}^{\varphi_1}(x-\varphi_1)\eta'(\hat{x})dx|\\
&\leq \varepsilon \rho B_1 (v_0+v_1)(U-u_0^0){\rm sup}_{x\in
(\varphi_2,\varphi_1)}|\eta'(x)|=\mathcal{O}(\varepsilon).
\end{align*}
This implies that we have chosen the function $F$ correctly and
we have to solve the system
\begin{equation}
\begin{gathered}
u_{\varepsilon t}+\Big(\frac{1}{2}u_{\varepsilon}^2(B_2-B_1)+B_1\left( u_0^0+ U
\right)u_{\varepsilon}\Big)_x=0\\
 v_{\varepsilon t}+\left(u_{\varepsilon} v_{\varepsilon}
+B_1(v_{\varepsilon} u_0^0+U v_{\varepsilon}-2u_{\varepsilon}v_{\varepsilon})
\right)_x=B_1(v_0+v_1)(U-u_0^0)
\frac{\kappa((\varphi_2,\varphi_1))}{\varphi_1-\varphi_2},
\end{gathered}
\label{jul266}
\end{equation}
with initial conditions (\ref{jul167}).
We remind that family $u_{\varepsilon}$, $\varepsilon>0$, of exact solutions of
problem (\ref{jul266}), (\ref{jul167}), perturbed possibly for term
of order $\varepsilon$, represents the weak asymptotic solution to problem
\eqref{jul166}, (\ref{jul167}).

So, we can pass on the third step of our procedure. At this instance
we want to solve Cauchy problem (\ref{jul266}), (\ref{jul167}). In
the previous section we found smooth solution to the first equation
from (\ref{jul266}) with initial data \eqref{avg915} and we pass to
the second one. After applying Leibnitz rule to the second equation
it becomes
\begin{align*}
&v_{\varepsilon t}+\left(u_{\varepsilon}(B_2-B_1)  +B_1( u_0^0+U) \right)v_{\varepsilon x}\\
&= -v_{\varepsilon}u_{\varepsilon x}(1-2B_1)+
B_1(v_0+v_1)(U-u_0^0)\frac{\kappa((\varphi_2(t),\varphi_1(t)))}{\varphi_1(t)
-\varphi_2(t)}.
\end{align*}
To simplify notation, in the sequel we will not use
perturbations as in e.g. (\ref{trnd4}). Clearly, this does not
affect on the generality of our considerations.

System of characteristics for this equation is
\begin{equation}
\begin{gathered}
\dot{x}=u_{\varepsilon}(B_2-B_1)+B_1(u_0^0+U),\\
\dot{v}_{\varepsilon}= -v_{\varepsilon}u_{\varepsilon x}(1-2B_1)+
B_1(v_0+v_1)(U-u_0^0)\frac{\kappa((\varphi_2,\varphi_1))}{\varphi_1-\varphi_2},
\\
v_{\varepsilon}(0)=v_{\varepsilon}|_{t=0}(x_0), \quad  x(0)=x_0
\end{gathered}
\label{no50oct}
\end{equation}
The first equation of the system is the same as the
first equation from (\ref{no10oct}). Therefore,
$\varphi_i(t,\varepsilon)=x(a_i,t,\varepsilon)$ where $x$ represents the solution
to the first equation in (\ref{no50oct}). Using the fact that the
characteristics are non-intersecting we know that for $x_0<a_2$ we
have $x<\varphi_2$ and for $x_0>a_1$ we have $x>\varphi_1$.
Accordingly, we can rewrite (\ref{no50oct}) as
\begin{gather*}
\dot{x}=u_{\varepsilon}(B_2-B_1)+B_1(u_0^0+U),\\
\dot{v_{\varepsilon}}=
\begin{cases}
-v_{\varepsilon}u_{\varepsilon x}(1-2B_1), & x_0<a_2,\\
-v_{\varepsilon}u_{\varepsilon x}(1-2B_1)+ B_1(v_0+v_1)(U-u_0^0)
 \frac{1}{\varphi_1-\varphi_2}, & x_0\in[a_2,a_1],\\
-v_{\varepsilon}u_{\varepsilon x}(1-2B_1), & x_0>a_1,
\end{cases},\\
v_{\varepsilon}(0)=v_{\varepsilon}|_{t=0}(x_0), \quad x(0)=x_0, \quad
x_0\in[a_2,a_1],
\end{gather*}
This is linear system of ODEs and it is not difficult to integrate it.
For the function $v_{\varepsilon}$, we have
\begin{equation}
v_{\varepsilon}=\begin{cases}
\frac{v_{0\varepsilon}(x_0)}{\frac{\partial x}{\partial x_0}}, & x_0<a_2,\\
\frac{v_{0\varepsilon}(x_0)}{\frac{\partial x}{\partial x_0}}+
\frac{(v_0+v_1)(U-u_0^0)}{\frac{\partial x}{\partial x_0}}
\int_0^t B_1 \frac{\frac{\partial x}{\partial x_0}}{(\varphi_1(t',\varepsilon)
-\varphi_2(t',\varepsilon))}dt', & x_0\in[a_2,a_1],\\
\frac{v_{0\varepsilon}(x_0)}{\frac{\partial x}{\partial x_0}}, & x_0>a_1.
\end{cases}
\label{avg2616}
\end{equation}
Recalling (\ref{avg2626}), from (\ref{avg2616}) it follows that
the solution of (\ref{jul266}), (\ref{jul167}) has the form
\begin{equation}
v_{\varepsilon}(x,t)=\begin{cases}
v_{0\varepsilon}(x_0(x,t,\varepsilon)), & x<\varphi_2,\\
\frac{v_{0\varepsilon}(x_0(x,t,\varepsilon))}{\frac{\partial x}{\partial x_0}}+
\frac{(v_0+v_1)(U-u_0^0)}{\frac{\partial x}{\partial x_0}}\int_0^t B_1
 \frac{\frac{\partial x}{\partial x_0}}{(\varphi_1(t',\varepsilon)-\varphi_2(t',\varepsilon))}dt',
 & x\in [\varphi_2,\varphi_1],\\
v_{0\varepsilon}(x_0(x,t,\varepsilon)), & x>\varphi_1,\\
\end{cases}
\label{jul267}
\end{equation}
where $x_0(x,t,\varepsilon)$ is inverse function of the function $x$ defined
as the solution to (\ref{no50oct}).
The existence of the function $x_0$ is proved in Theorem \ref{t**}.

Now we  return to (\ref{jul2610}). It remains to check if:
$$
B_1\int_{\varphi_2}^{\varphi_1}(u_0^0 v_{\varepsilon}
+U v_{\varepsilon}-2u_{\varepsilon}v_{\varepsilon})\eta'(x)dx=\mathcal{O}(\varepsilon).
$$
We substitute here expressions for $v_{\varepsilon}$ and $u_{\varepsilon}$. After
recalling (\ref{avg2626}) we have
\begin{equation} \label{**}
\begin{aligned}
&B_1\int_{\varphi_2}^{\varphi_1}
\Big( (u_0^0 +U )\frac{v_{0}(x_0(x,t,\varepsilon))}{\frac{\partial x}{\partial x_0}}
-2u_{0}(x_0(x,t,\varepsilon))
 \frac{v_{0}(x_0(x,t,\varepsilon))}{\frac{\partial x}{\partial x_0}}\Big)\eta'(x)dx\\
&+B_1\int_0^t B_1(\rho(\tau(t'))) dt'\cdot
\int_{\varphi_2}^{\varphi_1}\Big((u_0^0+U)\frac{(v_0+v_1)
(U-u_0^0)}{\varphi_1-\varphi_2}\\
&-2u_{0}(x_0(x,t,\varepsilon))\frac{(v_0+v_1)(U-u_0^0)}{\varphi_1-\varphi_2}\Big)
 \eta'(x)dx=\mathcal{O}(\varepsilon).
\end{aligned}
\end{equation}
 We change variable here passing from $x$ to $x_0$, i.e. we put
$x=x(x_0,t,\varepsilon)$ which implies
$dx=\frac{\partial x}{\partial x_0}dx_0$. Recall that we also have
$\varphi_i=x(a_i,t,\varepsilon)$, $i=1,2$, $a_1=\varepsilon^{1/2}$,
$a_2=-\varepsilon^{1/2}$. So, the above expression becomes
\begin{equation} \label{oslo4}
\begin{aligned}
&B_1\int_{-\varepsilon^{1/2}}^{\varepsilon^{1/2}}\left(
(u_0^0+U)v_0(x_0)-2u_0(x_0)v_0(x_0)
\right)\eta'(x(x_0,t,\varepsilon)dx_0\\
&+B_1\int_0^t B_1(\rho(\tau(t')))dt'
\int_{-\varepsilon^{1/2}}^{\varepsilon^{1/2}}(v_0+v_1)(U-u_0^0)
(u_0^0 +U \\
&\quad -2u_0(x_0))\eta'(x(x_0,t,\varepsilon))dx_0\\
&=\mathcal{O}(\varepsilon^{1/2}),
\end{aligned}
\end{equation}
and this is obviously true since $u_0$ and $v_0$ are bounded functions.
In that way, we have proved that the functions $u_{\varepsilon}$ and
$v_{\varepsilon}$ given by (\ref{jul961}), (\ref{jul267}), respectively,
represent weak asymptotic solution of problem \eqref{jul165},
\eqref{jul166}.

Finally, we let $\varepsilon \to 0$ to see what we obtain as a weak limit
of the weak asymptotic solution of problem \eqref{jul165},
\eqref{jul166}. For $u_{\varepsilon}$ we know that (Corollary \ref{coro5}):
\begin{equation*}
w-\lim_{\varepsilon\to 0}u_{\varepsilon}=
\begin{cases}
U, & x<(U+u_0^0)t/2,\\
u_0^0, & x\geq (U+u_0^0)t/2.
\end{cases}
\end{equation*}
We inspect weak limit of $v_{\varepsilon}$. We multiply $v_{\varepsilon}$ by
arbitrary $\eta\in C^1_0(\mathbb{R})$ and integrate over
$\mathbb{R}$,
\begin{align*}
\int v_{\varepsilon}(x,t) \eta(x) dx
&=\int_{-\infty}^{\varphi_2}v_{0\varepsilon}(x_0(x,t,\varepsilon))\eta(x)dx
  +\int_{\varphi_1}^{\infty} v_{0\varepsilon}(x_0(x,t,\varepsilon))\eta(x)dx\\
&\quad+ \int_{\varphi_2}^{\varphi_1}\Big(\frac{v_{0\varepsilon}
 (x_0(x,t,\varepsilon))}{\frac{\partial x}{\partial x_0}}\eta(x)dx \\
&\quad + \frac{(v_0+v_1)(U-u_0^0)}{\frac{\partial x}{\partial
x_0}}\int_0^t B_1 \frac{\frac{\partial x}{\partial
x_0}}{\varphi_1(t',\varepsilon)-\varphi_2(t',\varepsilon)}dt'\Big)\eta(x)dx.
\end{align*}
Passing from variable $x$ to variable $x_0$ in the second line of
the previous expression (as in (\ref{**}-\ref{oslo4})) and using
(\ref{avg2626}) we have
\begin{equation} \label{jul1061}
\begin{aligned}
&\int v_{\varepsilon}(x,t) \eta(x) dx\\
&=\int_{-\infty}^{\varphi_2} v_{0\varepsilon}(x_0(x,t,\varepsilon))\eta(x)dx
+\int_{\varphi_1}^{\infty} v_{0\varepsilon}(x_0(x,t,\varepsilon))\eta(x)dx
+\int_{-\varepsilon^{1/2}}^{\varepsilon^{1/2}}v_0(x_0)dx_0\\
&\quad +\frac{1}{\varepsilon^{1/2}-(-\varepsilon^{1/2})}\int_{-\varepsilon^{1/2}}^{\varepsilon^{1/2}}
(v_0+v_1)(U-u_0^0)\eta(x(x_0,t,\varepsilon))dx_0 \int_0^t B_1dt'.
\end{aligned}
\end{equation}
Letting $\varepsilon \to 0$ here we conclude that (see explanation below)
\begin{equation}
w-\lim_{\varepsilon\to 0}v_{\varepsilon}(x,t)\to
\frac{1}{2}t(v_0+v_1)(U-u_0^0)\delta(x-(U+u_0^0)t'/2)
+\begin{cases}
v_1, & x<(U+u_0^0)t/2\\
v_0, & x\geq (U+u_0^0)t/2.
\end{cases}
\label{no51oct}
\end{equation}
Now, we explain this passage in  detail. Recall that for every fixed $t>0$
we have $\tau=\frac{(U-u_0^0)t-2\varepsilon^{1/2}}{\varepsilon}\to \infty$ as $\varepsilon\to 0$.
Therefore, $B_1\to 1/2$ for every fixed $t$ and (\ref{no51oct})
follows. Similarly, $x(x_0,t,\varepsilon))\to \frac{(U+u_0^0)t}{2}$ as
$\varepsilon\to 0$ according to (\ref{ostr71}) and the fact that $t^*\to 0$
and $x^*\to 0$ as $\varepsilon\to 0$.

We collect the previous considerations in the next theorem. The result
coincides with the results in \cite{DSH,ind,mn}.

\begin{theorem} \label{thm7}
The Riemann problem \eqref{jul165}, \eqref{jul166} has
weak asymptotic solution $(u_{\varepsilon}, v_{\varepsilon})$ given by
\begin{align*}
u_{\varepsilon}(x,t)&=u_{0\varepsilon}(x_0(x,t,\varepsilon),\\
v_{\varepsilon}(x,t)&=\begin{cases}
v_{0\varepsilon}(x_0(x,t,\varepsilon)), &  x<\varphi_2,\\[5pt]
\frac{v_{0\varepsilon}(x_0(x,t,\varepsilon))}{\frac{\partial x}{\partial x_0}}+
\frac{(v_0+v_1)(U-u_0^0)}{\frac{\partial x}{\partial x_0}}\int_0^t B_1
\frac{\frac{\partial x}{\partial x_0}}{(\varphi_1(t',\varepsilon)-\varphi_2(t',\varepsilon))}dt',
 & x\in [\varphi_2,\varphi_1],\\[5pt]
v_{0\varepsilon}(x_0(x,t,\varepsilon)), & x>\varphi_1.
\end{cases}
\end{align*}
where $\varphi_i=\varphi_i(t,\varepsilon)$, $i=1,2$. Weak limit of the
weak asymptotic solution to \eqref{jul165}, \eqref{jul166} is
\begin{gather*}
w-\lim_{\varepsilon\to 0}u_{\varepsilon}=
\begin{cases}
U, & x<(U+u_0^0)t/2,\\
u_0^0, & x\geq (U+u_0^0)t/2.
\end{cases}
\\
\begin{aligned}
&w-\lim_{\varepsilon\to 0}v_{\varepsilon}(x,t)\\
&\to \frac{1}{2}t(v_0+v_1)(U-u_0^0)\delta(x-(U+u_0^0)t/2)
+\begin{cases}
v_1, & x<(U+u_0^0)t/2\\
v_0, & x\geq (U+u_0^0)t/2.
\end{cases}
\end{aligned}
\end{gather*}
\end{theorem}

Finally, as an example we show how the method can be applied
to the shock wave formation process in the case of multidimensional
scalar conservation law.

\section{Example}
We consider Cauchy problem (\ref{fri2}), (\ref{nonov1}) which is
special case of one appearing in the oil reservoir problems. As we
will see, geometrically, this problem is very simple, but if we
perturb geometry of our problem only a little bit, geometrical
approach becomes very complicated (see \cite{kos,nak}). On the other
hand, our approach is almost the same for large class of different
geometries. On this simple example we demonstrate basic principles
of the method in more then one space dimension. Complete treatment
will be done elsewhere.
\begin{gather}
L(u)=\partial_t u+\partial_{x_1}u^2+\partial_{x_2}u^2=0,
 \label{fri2}\\
u|_{t=0}=\hat{u}_0(x_1,x_2)
=\begin{cases} 1, & x_1<-2x_2-1\\
u_0(x_1,x_2), & -2x_2-1<x_1<-2x_2+1\\
-1, & x_1>-2x_2+1
\end{cases}
\label{nonov1}
\end{gather}
where the function $u_0$ we determine from the
continuity condition i.e.,
it has to be
\begin{gather*}
u_0\equiv 1 \quad  \text{on the line } x_1=-2x_2-1, \\
u_0\equiv -1 \quad \text{on the line } x_1=-2x_2+1,
\end{gather*}
and from the condition
\begin{equation}
\begin{gathered}
2\frac{\partial u_0}{\partial x_1}+2\frac{\partial u_0}{\partial x_2}+K=0, \\
u_0|_{x_1=-2x_2-1}=1, \quad u_0|_{x_1=-2x_2+1}=-1
\end{gathered} \label{nonov3}
\end{equation}
for some $K=K(s)$ where $s$ is a parameter of parametrization of the line
$x_1=-2x_2-1$. We take so, since the characteristics of problem
(\ref{nonov3}) start from
$$
\Gamma_1=\{(x_1,x_2):x_1=-2x_2-1\}
$$
 and end on
$$
\Gamma_2=\{(x_1,x_2):x_1=-2x_2+1\}.
$$
We explain this condition more closely. We begin with the remark
that it is analogical to the one dimensional condition which is
satisfied by initial data \eqref{avg915}.
Namely, system of characteristics for problem (\ref{fri2}),
(\ref{nonov1}) has the form
\begin{equation}
\begin{gathered}
\dot{x}_1=2u, \quad x_1(0)=x_{10}\\
\dot{x}_2=2u, \quad x_2(0)=x_{20}\\
\dot{u}=0, \quad u(0)=\hat{u}_0(x_{10},x_{20})
\end{gathered}
\label{nonov2}
\end{equation}
As is well known, our problem has classical solution as long as there
exists inverse function $(x_{10},x_{20})$ of the function
$(x_1,x_2)$ defined by (\ref{nonov2}) for $(x_{10},x_{20})\in
\{(x_1,x_2): -2x_2-1<x_1<-2x_2+1\}$ (since characteristics emanating
out of that interval are parallel), i.e., according to the inverse
function theorem, as long as (see e.g. \cite{prpDan} or \cite{nak}):
\begin{equation}
J=\det\Big|\frac{\partial x}{\partial x_0}\Big|
=t\big(2\frac{u_0}{\partial x_1}+2\frac{u_0}{\partial x_2} \big)+1\neq 0.
\label{ostr73}
\end{equation}
The point $(t^*,x^*_1,x^*_2)$ such that $J=0$, where $t^*$ is minimal such
that $J=0$, is usually called the point of the gradient catastrophe.
It appears that it is much easier to describe the shock wave
formation when we have 'the line of the gradient catastrophe' (see
\cite{DM} and compare with matching method \cite{Il}), i.e. the
curve $(x_1(\tau),x_2(\tau))$ such that $J=0$ for fixed minimal
$t^*$ (that is, for every $t<t^*$ we have $J\neq 0$) and every
$(x_1,x_2)\in (x_1(\tau),x_2(\tau))$. Of course, $\tau$ appearing
here is such that the point $(x_1(\tau),x_2(\tau))$ always lies
between $\Gamma_1$ and $\Gamma_2$.

Therefore, we look for the initial condition which will generate
curves of gradient catastrophe. According to all said above (compare
(\ref{nonov3}) and (\ref{ostr73})), such initial condition is given
exactly by boundary problem (\ref{nonov3}).

It is not difficult to integrate (\ref{nonov3}) and to determine
$K$. We have $K\equiv 6$ (in this case it does not depend on $s$)
and
\begin{gather*}
u_0(\tau)=1-\frac{1}{3}\tau,\\
x_1(\tau)=x_{10}+2\tau,\\
x_2(\tau)=x_{20}+2\tau.
\end{gather*}
 From here it is easy to find the function $u_0$. We
have
\begin{align*}
u_0(x_1,x_2)=1-\frac{x_1+2x_2+1}{18},
\end{align*}

After determining the function $u_0$ we continue as follows.
Since we have two dimensional problem we have to modify a little bit
the method we have used in one dimensional case. Here, it is not
convenient to write $x=\varphi(t)$ since $x$ has two dimensions and
we do not have appropriate relation of order in this case (that
means that it is very difficult to describe mutual position of the
point; compare to $\varphi_{i0}$, $i=1,2$ from Section 1).
Therefore, we write $t=\psi(x)$ (in \cite{nak} it was used
$x_2=\psi(t,x_1)$) which, roughly speaking, renders our problem on
one dimension.

In the sequel, by $\psi_{i0}(x)$, $i=1,2$, $x\in \mathbb{R}^2$, we
denote time necessary a point $x_{i0}\in \Gamma_i$ to reach to the
point $x$.
One can verify  that (see \cite[page 6]{prpDan})
$$
\Psi_0(x)=(\psi_{20}-\psi_{10})(x)=-2(t-1/6),
$$
for $x\in \psi_{10}(\Gamma_1)\cap \psi_{20}(\Gamma_2)$.

Now, as in the one dimensional case we replace  (\ref{fri2}) by
its weak asymptotic analogue
\begin{align}
L_{\varepsilon}(u)=u_{\varepsilon t}+{\rm div}\left(u_\varepsilon
^2(B_2(\rho)-B_1(\rho))+c\cdot u_\varepsilon B_1(\rho)\right)=0,
\label{fri1}
\end{align}where $\rho=\rho(\tau)$ is solution of the Cauchy problem:
\[
\dot\rho=1-2B_1(\rho), \quad  \lim_{\tau\to -\infty}\frac{\rho}{\tau}=1,
\]
and
$$
\tau=\frac{\Psi_0(x)}{\varepsilon}.
$$
Furthermore, $c=(c_1,c_2)=(0,0)$ since
for such $c$ we have
$$
L(u_{\varepsilon})=\mathcal{O}_{\mathcal{D}'}(\varepsilon)
$$
where $u_{\varepsilon}$ is global solution to problem (\ref{fri1}), (\ref{nonov1}).
Number 2 appears here since $B_1\to 1/2$ as $\varepsilon\to 0$ and $t>t^*$.
More precisely,
$$
\rho=\frac{\varphi_2(t,s,\varepsilon)-\varphi_1(t,s,\varepsilon)}{\varepsilon}, \quad
t\in \mathbb{R}^+, \; s\in \mathbb{R},
$$
and
$\varphi_2(0,s,\varepsilon)\in \Gamma_2$ and $\varphi_1(0,s,\varepsilon)\in
\Gamma_1$ are connected by the characteristics of equation
(\ref{nonov3}). More precisely,
$$
\varphi_2(0,s,\varepsilon)=\varphi_1(0,s,\varepsilon)+(2\tau,2\tau)\in \Gamma_2,
$$
for some $\tau>0$ and, as in the one dimensional case, for every fixed $s$,
\begin{gather*}
\frac{d}{dt}\varphi_1(t,s,\varepsilon)=-2(B_2(\rho)-B_1(\rho)),\\
\frac{d}{dt}\varphi_2(t,s,\varepsilon)= 2(B_2(\rho)-B_1(\rho)).
\end{gather*}
Note that if the initial function $u_0$ is not constant on the
lines $\Gamma_{i}$, $i=1,2$, the righthand side of the latter
equations (defining $\varphi_i$, $i=1,2$), will depend explicitly on $s$.
Also, note that we can look for the (asymptotic) solution along
lines $(x_1(\tau), x_2(\tau))$ since solution is globally smooth
everywhere.

As $\varepsilon\to 0$ we see that for $t<1/6$ we have classical solution to
the problem and for $t>1/6$ the solution is stationary shock shock
concentrated on the straight line
$$
x_1=-2x_2.
$$
Details of the construction will be done elsewhere for the general
case of multidimensional scalar conservation law and more general
situation of initial data.


\begin{thebibliography}{99}

\bibitem{bre}Y. Brenier,
\emph{Averaged Multivalued Solutions for Scalar Conservation Laws},
SIAM Journal of Numerical Analysis, Vol. 21, No. 6. (Dec., 1984),
pp. 1013-1037.

\bibitem{chen} G-Q. Chen, H. Liu, Formation of $\delta$-shocks and
vacuum states in the vanishing pressure limit of solutions to the
Euler equations for isentropic fluids,  SIAM J. Math. Anal.
{\textbf{34}} (2003), no. 4, 925--938.

\bibitem{daf} C. M. Dafermos
\emph{Hyperbolic Conservation Laws in Continuum Physics},
Berlin; Heidelberg; New York; Barcelona; Hong Kong; London; Milan;
Paris; Singapore; Tokyo: Springer, 2000.

\bibitem{dan} V. G. Danilov,
\emph{Generalized Solution Describing Singularity Interaction},
International Journal of Mathematics and Mathematical Sciences,
Volume 29, No. 22. February 2002, pp. 481-494.

\bibitem{prpDan} V. G. Danilov
\emph{Remarks on the formation and decay of multidimensional shock
waves}, preprint available on http://www.math.ntnu.no/conservation,
2004-033

\bibitem{DM} V. G. Danilov, D. Mitrovic,
\emph{Weak asymptotic of shock wave formation process}, Nonlinear
Analysis, 61(2005) 613-635.

\bibitem{DM5} V. G. Danilov, D. Mitrovic,
\emph{Delta shock wave formation in the case of triangular hyperbolic
system of conservation laws}, preprint available at
http://www.math.ntnu.no /conservation/2006/057.html

\bibitem{DM4} V. G. Danilov, D. Mitrovic,
\emph{Smooth approximations of global in time solutions to scalar
conservation laws}, preprint.

\bibitem{do} V. G.  Danilov, G. A. Omelianov \emph{Weak asymptotic
method for the study of propagation and interaction of infinitely
narrow $\delta$-solitons}, Electron. J. Differential Equations 2003,
No. 90, 27 pp. (electronic).

\bibitem{DSO} V. G. Danilov, G. A. Omelianov, V. M. Shelkovich,
\emph{Weak Asymptotic Method and Interaction of Nonlinear Waves} in:
M.Karasev (Ed.), Asymptotic Methods for Wave and Quantum Problems,
American Mathematical Society Translation Series, vol. 208, 2003,
pp. 33-165.

\bibitem{DSH1} V. G. Danilov,  V. M. Shelkovich,
\emph{Propagation and interaction of nonlinear waves}, in: \emph{
Proceedings of Eight International Conference on Hyperbolic
Problems. Theory-Numerics-Applications}, Univ. Magdeburg,
Magdeburg, 2000, pp. 326--328.

\bibitem{DSH2} V. G. Danilov,  V. M. Shelkovich,
\emph{Propagation and interaction of $\delta$-shock waves of
hyperbolic systems of conservation laws}, Dokl. Akad. Nauk 394
(2004), no. 1, 10-14.

\bibitem{DSH} V. G. Danilov, V. M. Shelkovich,
\emph{Dynamics of propagation and interaction of $\delta$-shock waves
in conservation law system}, Journal of Differential Equations,
211(2005) 333-381.

\bibitem{DSH3} V. G. Danilov, V. M. Shelkovich,
\emph{Delta-shock wave type solution of hyperbolic systems of
conservation laws}, Quart. Appl. Math. 63 (2005), no. 3, 401-427.

\bibitem{PSH} E. Yu. Panov, V. M. Shelkovich, \emph{$\delta'$-shock
waves as a new type of solutions to systems of conservation laws},
J.Differential Equations 228 (2006), no. 1, 49-86.

\bibitem{She} V. M. Shelkovich,
\emph{The Riemann problem admitting $\delta$-, $\delta'$-shocks, and
vacuum states (the vanishing viscosity approach)}, J. Differential
Equations 231 (2006), no. 2, 459-500.

\bibitem{Il}A. M. Il'in,
\emph{Matching of Asymptotic Expansions of Solutions of Boundary
Value Problems}, Nauka, Moscow, 1989; English transl., AMS,
Providence, RI, 1992.

\bibitem{kos}S. Izumiya, G. Kossioris,
\emph{Geometric Singularities for Solutions of Single Conservation
Laws}, Arch. Rational Mech. Anal. 139 (1997) 255-290.

\bibitem{ind} K. T. Joseph,
\emph{A Rieman problem whose viscosity solution contain $\delta$
measures}, Asymptotic Analysis 7 (1993), 105-120.

\bibitem{mn1}D. Mitrovic, M. Nedeljkov,
\emph{Delta shock waves as a limit of shock waves}, J. of Hyperbolic
Differential Equations, to appear.

\bibitem{nak} S. Nakane,
\emph{Formation of shocks for a single conservation law}, SIAM J.
Math. Anal., Vol. 19, No. 6, November 1988

\bibitem{mn}M. Nedeljkov,
\emph{Delta and singular delta locus for one-dimensional systems of
conservation laws}, Math. Meth. Appl. Sci. {\textbf{27}} (2004),
931--955.

\bibitem{chi} H. Yang,
\emph{Riemann problems for class of coupled hyperbolic system of
conservation laws}, Journal of Differential Equations, 159(1999)
447-484.

\end{thebibliography}

\end{document}