\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 $tt^*$). 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 $t0$ 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^*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_0a_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_0a_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_0a_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_0a_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_0a_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-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-1t^*$. 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}