Delta shock waves with Dirac delta function in both components for systems of conservation laws

Available online at www.sciencedirect.com

ScienceDirect J. Differential Equations 257 (2014) 4369–4402 www.elsevier.com/locate/jde

Delta shock waves with Dirac delta function in both components for systems of conservation laws ✩ Hanchun Yang ∗ , Yanyan Zhang Department of Mathematics, Yunnan University, Kunming 650091, PR China Received 1 June 2014; revised 11 July 2014 Available online 10 September 2014

Abstract We study a class of non-strictly and weakly hyperbolic systems of conservation laws which contain the equations of geometrical optics as a prototype. The Riemann problems are constructively solved. The Riemann solutions include two kinds of interesting structures. One involves a cavitation where both state variables tend to zero forming a singularity, the other is a delta shock wave in which both state variables contain Dirac delta function simultaneously. The generalized Rankine–Hugoniot relation and entropy condition are proposed to solve the delta shock wave. Moreover, with the limiting viscosity approach, we show all of the existence, uniqueness and stability of solution involving the delta shock wave. The generalized Rankine–Hugoniot relation is also confirmed. Then our theory is successfully applied to two typical systems including the geometric optics equations. Finally, we present the numerical results coinciding with the theoretical analysis. © 2014 Elsevier Inc. All rights reserved. Keywords: Hyperbolic systems of conservation laws; Geometric optics; Cavitation; Delta shock wave; Generalized Rankine–Hugoniot relation; Limiting viscosity approach

1. Introduction This work is the twin of a recent Ref. [8]. In [8], we studied the Riemann problems of the hyperbolic systems of conservation laws ✩

Supported by the NSF of China 11361073.

* Corresponding author.

http://dx.doi.org/10.1016/j.jde.2014.08.009 0022-0396/© 2014 Elsevier Inc. All rights reserved.

4370

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402



  ut + φ(u, v)u x = 0,   vt + φ(u, v)v x = 0

(1.1)

with Riemann initial data (u, v)(0, x) = (u± , v± )

(±x > 0),

(1.2)

where φ(u, v) satisfies the assumption: (H1 ) φ(u, v) = φ(au + bv) is a given smooth function satisfying a 2 + b2 = 0, a and b are constants. A theory of delta shock wave with Dirac delta function developing in both state variables was established and applied to a large variety of systems. With regard to the delta shock wave and its related topics, we refer to [8] for a more detailed review. The theory of delta shock wave with Dirac delta function developing in both state variables is quite different from the aforementioned those, on which only one state variable contains the Dirac delta function. However, the theory does not work in some typical systems. For instance, the following equations of geometrical optics ⎧

⎪ u2 ⎪ ⎪ + = 0, u √ ⎨ t 2 + v2 u x (1.3)

⎪ uv ⎪ ⎪ vt + √ =0 ⎩ u2 + v 2 x were proposed by Engquist and Runborg [1] in 1996 when they considered multiple scale problems in numerical computations of propagating waves in vacuum medium, where (u, v) represents the direction and strength (particle density) of ray. The system (1.3) is only weakly hyperbolic. Generally, it means that (1.3) is not well-posed in the strongly hyperbolic sense. Besides, the system (1.3) possesses extreme degeneracy and strong nonlinearity in flux functions. These features may lead to a tough problem in the system (1.3). According to rough analysis and numerical evidence, they suggested that (1.3) can have delta function type solutions which may be used to describe the phenomenon where the phase should have split into two new phases. Since then, the existence of solutions to (1.3) has been an open question. The Riemann problem has not been solved so far. To our knowledge, some efforts have been done to explore the theoretical issues, for example [2,7] in which the authors tried to weaken the nonlinearity of flux functions by transformation of variables, but the results are far from the desired aim because the delta shock wave of (1.3) may contain Dirac delta function simultaneously in both state variables u and v. Fortunately, the problem has been fulfilled as below. Solving Riemann problem (1.3), (1.2) is interesting and exciting, but it is just one of the objectives of the present paper. To this end, we focus here on establishing a theory of the delta shock wave for a class of systems (1.1) with the assumption: (H2 ) φ(u, v) is a smooth function satisfying φ(u, v) = φ(αu, αv), α > 0 is constant. Obviously, in this situation, we can rewrite φ(u, v) = φ( uv ) = φ(r) (r = uv ) if v = 0. By taking φ(u, v) =  2u 2 , one can find that the system (1.3) is the very prototype of (1.1). u +v

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4371

For the systems (1.1) with the assumptions (H1 ) and (H2 ), respectively, they have substantial differences in spite of the homologous model. The system (1.1) with (H1 ) possesses two eigenvalues. One is linearly degenerate, the other is genuinely nonlinear except for certain degenerate points. So the classical elementary waves consist of rarefaction waves, contact discontinuities and shock waves. While the system (1.1) with (H2 ) has a double eigenvalue. It is linearly degenerate and weakly hyperbolic, thus the classical elementary waves only involve contact discontinuities. Besides, the latter has extreme degeneracy and strong nonlinearity in flux functions. It is exactly a stiff feature which causes an essential difficulty to be addressed. This is also just one of the main reasons why the theory of delta shock wave in [8] cannot be applied to the system (1.3). Here we use a new approach and skill to achieve our goal. In what follows, we outlook the context of each section in this paper. Here and after, we study the system (1.1) under the assumption (H2 ). In Section 2, with the phase plane and characteristic analysis, we construct the Riemann solutions of (1.1), (1.2) by two kinds. For the case r− < r+ , where r± = uv±± , the solution involves contact discontinuity and cavitation. The cavitation is a singular solution with a violated point (u, v) = (0, 0) which forms a singularity. In nonlinear geometrical optics, the cavitation may be used to indicate the fact that ray vanishes. It is a quite rare phenomenon different from the previous vacuum state in zero-pressure gas dynamics in which the density becomes zero and the velocity is a smooth function [5]. For the other case r− > r+ , we prove that u, v and ( uv )x blow up simultaneously in finite time even if the initial data are smooth enough. It means that no classical weak solutions to (1.1), (1.2) exist, and the delta shock wave with both components u and v simultaneously containing Dirac delta measures develops in solution. In Section 3, based on the idea and techniques in [8], by the definition of delta shock type solution to (1.1), for the delta shock wave, since vu has a bounded variation, we can arrange the function uv to take certain values on the discontinuity line so that we can overcome the difficulty of multiplying a Dirac delta function by a discontinuous function in both equations to be defined in the sense of distributions. Then we propose the generalized Rankine–Hugoniot relation to describe the relationship among the location and propagation speed of the discontinuity, rate of change of its weights (resp. strengths) and reassignment of u and v on its discontinuity. This relation applied especially to Riemann problem (1.1), (1.2) can be reduced to a function equation of one variable. Thus, the existence and uniqueness of solution involving the delta shock wave is equivalent to studying the existence and uniqueness of solution of the function equation, which is solved completely under the entropy condition. By the limiting viscosity approach, say [3,4], Sections 4–6 study the stability of the obtained delta shock wave solution. We consider the following viscous perturbations 

  ut + φ(u, v)u x = εtuxx ,   vt + φ(u, v)v x = εtvxx

(1.4)

with initial data (1.2). The process of proof with some special modifications and techniques is similar to that in [8]. Firstly, we show that the viscosity regularized problem (1.4), (1.2) has a smooth solution (uε (ξ ), v ε (ξ )) (ξ = x/t) on (−∞, +∞) for every ε > 0. Secondly, we prove that the solution of (1.4), (1.2) obtained in Section 4 generates the delta shock solution of (1.1), (1.2) which is weak star limit of (uε (ξ ), v ε (ξ )) (ξ = x/t) as ε → 0. Then from the results in Sections 4–6, one can find that r(t, x) = u(t, x)/v(t, x) is always a bounded monotone function of ξ = x/t, while u(t, x) and v(t, x) may be unbounded along the ray x = τ0 t simultaneously.

4372

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

In particular, when r− > r+ , both u(t, x) and v(t, x) are the sums of a step function and a Dirac δ-function matching respective strengths with the discontinuity line x = σ t as their support. Meanwhile, r(t, x) is a step function, u(t, x) and v(t, x) are required to take certain values satisfying φ(r(σ )) = σ on x = σ t, so that (1.1) holds in the sense of distributions. These facts show that the delta shock solution constructed is stable to the reasonable viscous perturbations (1.4), and confirm the generalized Rankine–Hugoniot relation (3.5), (3.6) proposed in Section 3. In Section 7, we give two typical examples including the equations of geometrical optics (1.3) to present an application for results of the paper. Thus the Riemann problem (1.3), (1.2) is solved completely. It is also noticed, in [1], that the model of geometrical optics provides a physical interpretation for this type of delta shock wave which shows the research of significance in actual applications. They will be further explored in the future. In addition, the work together with that in [8] composes a relatively complete theory on the delta shock waves possessing the generality and practicability. Finally, by employing the weighted essentially non-oscillatory scheme (WENO) [6], we simulate the cavitation and delta shock wave for the equations of geometrical optics (1.3), as presented in Section 8. The numerical results confirm the theoretical analysis. 2. Classical elementary waves and Riemann solutions Consider the Riemann problem (1.1), (1.2). The eigenvalue and corresponding right eigenvector are λ = φ(u, v),

r = (φv , −φu )T ,

and a simple calculation shows that ∇λ · r ≡ 0. Thus (1.1) is non-strictly hyperbolic and λ is linearly degenerate. For convenience, we assume φ(0) = 0,

φr > 0,

v > 0,

(2.1)

where the subscript r denotes the derivative with respect to r, then φ(r) is strictly monotone and increasing with respect to r. The remaining cases can be treated in a same way. As usual, we seek the self-similar solution (u, v)(t, x) = (u, v)(ξ ), then (1.1), (1.2) can be reduced to 

ξ=

  −ξ uξ + φ(u, v)u ξ = 0,   −ξ vξ + φ(u, v)v ξ = 0

x , t

(2.2)

and the boundary values at infinity, i.e., (u, v)(±∞) = (u± , v± ).

(2.3)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4373

For smooth solution, (2.2) is equivalent to

−ξ + uφu + φ(u, v) uφv vφu −ξ + vφv + φ(u, v)



u = 0, v ξ

which provides either the constant state (u, v)(ξ ) = constant,

(2.4)

⎧ ⎨ ξ = φ(u, v), u(φu du + φv dv) = 0, ⎩ v(φu du + φv dv) = 0,

(2.5)

or the singular solution

which, by integrating, yields

ξ = φ(u, v), u = kv,

(2.6)

where k is an arbitrary constant. The singular solution (2.5) contains a violated point (u, v) = (0, 0) which does not satisfy the system (1.1). In terms of the vacuum state in fluid dynamics, it is called a cavitation or generalized vacuum state, denoted by Cav, which may be used to interpret the phenomenon where ray vanishes. If ξ = σ is a bounded discontinuity, the Rankine–Hugoniot condition 

 −σ [u] + φ(u, v)u = 0,

 −σ [v] + φ(u, v)v = 0

(2.7)

holds, where [p] = p− − p+ denotes the jump of p across the discontinuity. By solving (2.7), we obtain the contact discontinuity J : ξ = σ = φ(u− , v− ) = φ(u+ , v+ ),

(2.8)

which is just the characteristic line for both sides in (x, t)-plane. Thanks to the monotonicity of φ(r), it follows that two states (u− , v− ) and (u+ , v+ ) can be connected by a contact discontinuity J if and only if φ(u− , v− ) = φ(u+ , v+ ), namely, uv−− = uv++ . That is, the two states (u− , v− ) and (u+ , v+ ) are on the line uv = uv−− = uv++ emitting from the origin in the (u, v)-phase plane. Now with constants, cavitation and contact discontinuity, one can construct solutions of the Riemann problem (1.1), (1.2) by two cases. For the case r− < r+ , no characteristic lines from the x-axis pass through the domain {(x, t)|φ(r− )t  x  φ(r+ )t, t > 0} in the (x, t)-plane, therefore the cavitation should appear in this domain. In the (u, v)-phase plane, we also can draw lines uv = uv−− and uv = uv++ from (u− , v− ) and (u+ , v+ ), respectively. The two lines intersect only at (0, 0), an isolated point. Thus we obtain the solution which consists of two contact discontinuities and a cavitation besides two constant states. As shown in Fig. 1, the solution can be expressed as

4374

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Fig. 1. Riemann solution of case r− < r+ .

Fig. 2. Characteristic analysis of delta shock wave.

⎧ ⎨ (u− , v− ), (u, v)(ξ ) = Cav, ⎩ (u+ , v+ ),

−∞ < ξ < φ(r− ), φ(r− )  ξ  φ(r+ ), φ(r+ ) < ξ < +∞.

(2.9)

In the case r− > r+ , the characteristic lines from the x-axis will overlap in the domain Ω = {(x, t)|φ(r+ )t  x  φ(r− )t, t > 0} in the (x, t)-plane, thus the singularity of solutions must happen in the domain Ω shown in Fig. 2. We assert that there exist no solutions in bounded variation space, that is, the solution cannot be constructed by classical waves. Indeed, from (1.1), by the characteristic equations dx = φ(r), dt

du = −uφ(r)x , dt

dv = −vφ(r)x , dt

(2.10)

one can find that rt + φ(r)rx = 0.

(2.11)

Then we consider Eq. (2.11) with sufficiently smooth initial data r(0, x) = r0 (x) satisfying r0 (x) < 0. From (2.10), we get that the characteristic line passing through any given point (0, α) on the x-axis is

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4375

  x = α + φ r0 (α) t,

(2.12)

on which r takes the constant value r = r0 (α). Differentiating Eq. (2.11) with respect to x, we have rtx + φ (r)rx2 + φ(r)rxx = 0,

(2.13)

φ (r)rtx + φ (r)φ(r)rxx = −φ 2 (r)rx2 .

(2.14)

namely,

From (2.14), we can deduce that 2  d(φ (r)rx ) = − φ (r)rx , dt

(2.15)

which is a type of the Riccati equation. Hence, along the characteristic (2.12), by (2.10) and (2.15), we can obtain that rx =

r0 (α) . 1 + φ (r0 (α))r0 (α)t

(2.16)

Combining with (2.10), we have u(t, x) =

u0 (α) , 1 + φ (r0 (α))r0 (α)t

v(t, x) =

v0 (α) . 1 + φ (r0 (α))r0 (α)t

(2.17)

Noting that φ (r0 (α)) > 0, r0 (α) < 0, it follows that lim

t→−(φ (r0 (α))r0 (α))−1

(u, v, rx ) = (∞, ∞, ∞)

(2.18)

along the characteristic (2.12). It shows that the solution u, v and rx must blow up simultaneously at a finite time. Define        −1  . E = (x, t)x = α + φ r0 (α) t, t = − φ r0 (α) r0 (α)

(2.19)

Then the set E is the envelope of all characteristics, on which all of u, v and rx become singular. Therefore the smooth solutions can be defined by (u, v)(t, x) =

u(α) v(α) ,



1 + φ (r0 (α))r0 (α)t 1 + φ (r0 (α))r0 (α)t

(2.20)

only for t < −(φ (r0 (α))r0 (α))−1 . The above discussion shows the mechanism of occurrence of delta shock wave with both variables u and v simultaneously containing Dirac delta functions. Thus, motivated by [8], for the case r− > r+ , we can construct the solution using a delta shock wave illustrated in Fig. 2.

4376

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

3. Generalized Rankine–Hugoniot relation of delta shock waves In view of discussion in the above section, for the case r− > r+ , we will seek the delta-shock solution of (1.1), (1.2). We will also clarify the generalized Rankine–Hugoniot relation and entropy condition for the delta shock wave, and then apply them to solving the Riemann problem (1.1), (1.2) again. In order to accomplish this, we introduce the definitions of two-dimensional weighted delta function and delta shock wave type solution as follows. Definition 3.1. A two-dimensional weighted delta function w(s)δS supported on a smooth curve S = {(x(s), t (s)) : a  s  b} can be defined by 

 w(s)δS , ϕ(t, x) =

b

  w(s)ϕ t (s), x(s) t (s)2 + x (s)2 ds

(3.1)

a

for all the test functions ϕ ∈ C0∞ ([0, +∞) × (−∞, +∞)). Definition 3.2. A pair (u, v) is called a delta shock wave type solution of (1.1) in the sense of distributions if there exist a smooth curve S and a weight w such that u and v are represented in the following form u = U (t, x) + rδ wδS ,

v = V (t, x) + wδS ,

(3.2)

U, V ∈ L∞ ((0, +∞) × R; R), w ∈ C 1 (S), φ(r)|S = σ = φ(rδ ), σ is the tangential derivative of curve S, and it satisfies   u, ϕt + φ(u, v)u, ϕx = 0,   v, ϕt + φ(u, v)v, ϕx = 0

(3.3)

for all test functions ϕ ∈ C0∞ ([0, +∞) × (−∞, +∞)), where +∞ +∞ u, ϕ = U ϕdxdt + rδ wδS , ϕ , 0 −∞

+∞ +∞ φ(u, v)u, ϕ = φ(U, V )U ϕdxdt + σ rδ wδS , ϕ ,





0 −∞

and v has the similar integral identities as above. With the definition, we propose to find a solution with the discontinuity x = x(t) for system (1.1) of the form ⎧ ⎨ (u− , v− )(t, x), (u, v)(t, x) = (rδ w(t), w(t))δ(x − x(t)), ⎩ (u+ , v+ )(t, x),

x < x(t), x = x(t), x > x(t),

(3.4)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4377

where (u− , v− )(t, x) and (u+ , v+ )(t, x) are piecewise smooth solutions of system (1.4), δ(·) is the standard Dirac measure, x(t) ∈ C 1 , and rδ is to be determined later. If (3.4) satisfies the relations ⎧ dx ⎪ ⎪ ⎪ dt = σ, ⎪ ⎪ ⎪ ⎨ √

 d 1 + σ 2 rδ w(t) = −σ [u] + φ(u, v)u , ⎪ dt ⎪ ⎪ √ ⎪ ⎪ 2 ⎪ ⎩ d 1 + σ w(t) = −σ [v] + φ(u, v)v  dt

(3.5)

φ(r)|x=x(t) = φ(rδ ) = σ,

(3.6)

and

then the solution (u, v)(t, x) defined in (3.4) satisfies (1.1) in the sense of distributions. In particular, here it is noticed, from (2.11) and (2.18), that r = uv has a bounded variation. Thus we arrange φ(r) on the discontinuous line x = σ t as (3.6), so that we can multiply a delta function by a discontinuity function in a similar way to that in [5,8], etc. Eqs. (3.5)–(3.6) are called the generalized Rankine–Hugoniot relation. It reflects the exact relationship among the limit states on both sides of the discontinuity, location, propagation speed, weights and reassignment of uv on the discontinuity. In fact, if (3.5) and (3.6) hold, for any test function ϕ ∈ C0∞ ([0, +∞) × (−∞, +∞)), by Green’s formulation and integrating by parts, it follows that   I = u, ϕt + φ(u, v)u, ϕx +∞ σ t +∞ +∞     ϕt + φ(u− , v− )ϕx u− dxdt + ϕt + φ(u+ , v+ )ϕx u+ dxdt = 0 −∞ +∞

 ϕt (t, σ t) + σ ϕx (t, σ t) 1 + σ 2 rδ w(t)dt

0



ϕφ(u− , v− )u− dt − ϕu− dx −

= +

σt



+ 

0

ϕφ(u+ , v+ )u+ dt − ϕu+ dx

+∞ 

 ϕt (t, σ t) + σ ϕx (t, σ t) 1 + σ 2 rδ w(t)dt

0

=

+∞ 



φ(u− , v− ) − σ u− ϕdt −

0

+∞ 

+∞ √ d 1 + σ 2 rδ w(t) ϕdt dt

0

0

 φ(u+ , v+ ) − σ u+ ϕdt −

√

+∞

 d 1 + σ 2 rδ w(t) = ϕdt −σ [u] + φ(u, v)u − dt 0

= 0,

4378

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

which implies the first equation of (3.3). A completely similar argument leads to the second equation of (3.3). This claim is true. To guarantee uniqueness, we should supplement the entropy condition (3.7)

φ(r+ ) < φ(rδ ) < φ(r− ),

which means that all characteristics on both sides of the discontinuity are in-coming. Under condition φ (r) > 0, it is equivalent to (3.8)

r+ < rδ < r− .

A discontinuity satisfying (3.5)–(3.8) is called a delta shock wave of the system (1.1) with (H2 ), symbolized by δ. Next, we apply the generalized Rankine–Hugoniot relation to solve Riemann problem (1.1), (1.2) for the case r− > r+ . At this moment, the Riemann problem can be reduced to solve the ordinary differential equations (3.5)–(3.6) with initial data t = 0: x(0) = 0,

w(0) = 0.

(3.9)

In view of knowledge concerning delta shock waves, say [5,8], we find that σ is constant and w(t) is a linear function of t . So a delta shock solution of (1.1), (1.2) can be assumed to take the following form δ:

x(t) = σ t,

w(t) = w0 t,

σ = φ(rδ )

(3.10)

which satisfies (3.9), where σ , w0 and rδ are to be determined constants. Substituting (3.10) into (3.5)–(3.6), one can rewrite the generalized Rankine–Hugoniot relation as ⎧

 2 ⎪ ⎨ 1 + σ rδ w0 = −σ [u] + φ(r)u , 

 2 ⎪ 1 + σ w0 = −σ [v] + φ(r)v , ⎩ σ = φ(rδ ).

(3.11)

From (3.11), we obtain the function equation of one variable

  rδ φ(rδ )[v] − φ(rδ )[u] − rδ φ(r)v + φ(r)u = 0.

(3.12)

Under the entropy condition (3.8), we now consider a solution of (3.12). Setting

  f (rδ ) = rδ φ(rδ )[v] − φ(rδ )[u] − rδ φ(r)v + φ(r)u ,

(3.13)

we can calculate that   f (r+ ) = v− φ(r− ) − φ(r+ ) (r− − r+ ) > 0,   f (r− ) = −v+ φ(r− ) − φ(r+ ) (r− − r+ ) < 0, because of r− > r+ and φ (r) > 0. Again, differentiating f (rδ ) with respect to rδ gives

(3.14)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

    f (rδ ) = v− φ(rδ ) − φ(r− ) + v+ φ(r+ ) − φ(rδ )   + v− (rδ − r− ) + v+ (r+ − rδ ) φ (rδ ) < 0.

4379

(3.15)

According to the zero point theorem in mathematical analysis, the function f (rδ ) has one and only one zero point in (r+ , r− ). That is, Eq. (3.12) has a unique solution, denoted by rδ , which satisfies r+ < rδ < r− . Returning to Eqs. (3.11), we solve the σ and w0 uniquely. Thus we conclude with the following result. Theorem 3.1. Under assumption (2.1), let φ(r− ) > φ(r+ ). Then the Riemann problem (1.1), (1.2) admits one and only one delta shock wave solution of the form ⎧ ⎨ (u− , v− ), (u, v)(t, x) = (rδ w(t), w(t))δ(x − x(t)), ⎩ (u+ , v+ ),

x < σ t, x = σ t, x > σ t,

(3.16)

where w(t) = w0 t and all of the three constants σ , w0 , and rδ are determined uniquely by (3.11) and (3.8). Remark. Theorem 3.1 is true for the case φ (r) > 0, v < 0. Furthermore, we find that the result is similar to Theorem 3.1 for the case φ (r) < 0, v > 0 and φ (r) < 0, v < 0, which can be obtained by using the same way as the above proof with only few modifications. These facts can also be observed from the below proofs in Sections 4–6. In addition, we also discuss signs of w0 and rδ in order to present numerical results in Section 8 below. From (3.11), under entropy condition (3.8), one can get that w0 = √ =√

1

−φ(rδ )(v− − v+ ) + φ(r− )v− − φ(r+ )v+

1 + σ2 1 1 + σ2



    v+ φ(rδ ) − φ(r+ ) + v− φ(r− ) − φ(rδ ) > 0.

(3.17)

This fact means that the weight of particles in the sticking process is getting heavier as time increases. For the sign of rδ , we consider the following three cases. Case 1.

u− v−

>

u+ v+

> 0.

By (3.8), we get Case 2. 0 >

u− v−

>

u− v−

> rδ >

u+ v+

> 0. At this moment, weights of u and v are positive.

u+ v+ .

It is easy to see 0 > uv−− > rδ > weight of u is negative. Case 3.

u− v−

>0>

u+ v+ .

u+ v+ ,

so we can find that the weight of v is positive, while the

4380

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

In this case, we have u− > 0 > u+ . The sign of rδ is changing. From (3.12), when φ(r− )u− = φ(r+ )u+ , we have   φ(rδ )(u+ − u− ) − rδ φ(r− )v− − φ(r+ )v+ + rδ φ(rδ )(v− − v+ ) = 0, which solves rδ = 0 because f (rδ ) has only one zero point in ( uv++ , uv−− ). When φ(r− )u− > φ(r+ )u+ , we have

u− f (0) > 0, f < 0, f (rδ ) < 0, v−

(3.18)

and we can deduce that rδ > 0. Similarly, when φ(r− )u− < φ(r+ )u+ , we get rδ < 0. In summary, for the third case, the sign of rδ depends on the jump [φ(r)u]. 4. Existence of solutions to viscous problem (4.1), (4.2) In this section, we show that the viscosity regularized problem (1.4), (1.2) has a smooth solution. Equivalently, we consider the boundary value problem    −ξ uξ + φ(u, v)u ξ = εuξ ξ ,   (4.1) −ξ vξ + φ(u, v)v ξ = εvξ ξ with boundary conditions (u, v)(±∞) = (u± , v± ). If (u(ξ ), v(ξ )) is a solution of (4.1), (4.2), then v(ξ ) and r(ξ ) = 

(4.2) u(ξ ) v(ξ )

satisfy

  −ξ(rv)ξ + φ(r)rv ξ = ε(rv)ξ ξ ,   −ξ vξ + φ(r)v ξ = εvξ ξ ,

(4.3)

and (r, v)(±∞) = (r± , v± ).

(4.4)

Because our focus is the delta shock wave, we mainly discuss the situation r− > r+ . In order to show the existence of a smooth solution to (4.1), (4.2) on (−∞, +∞) for every fixed ε > 0, we consider the two-parameter boundary value problem 

  −ξ uξ + μ φ(u, v)u ξ = εuξ ξ ,   −ξ vξ + μ φ(u, v)v ξ = εvξ ξ ,

(4.5)

and (u, v)(±L) = (μu± , μv± ),

(4.6)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4381

where parameters μ ∈ [0, 1], L  1. Corresponding to (4.3), (4.4), we have    −ξ(rv)ξ + μ φ(r)rv ξ = ε(rv)ξ ξ ,   −ξ vξ + μ φ(r)v ξ = εvξ ξ ,

(4.7)

and (r, v)(±L) = (μr± , μv± ).

(4.8)

When μ = 1 and L → +∞, (4.5)–(4.8) are reduced to (4.1)–(4.4). As an application of Theorem 2.1 in [4], we have Theorem 4.1. Assume that there are positive constants M and δ depending on u+ , u− , v+ , v− and ε (and thus independent of μ and L) such that any solution (u(ξ ), v(ξ )) of (4.5), (4.6), corresponding to any 0  μ  1, L  1, satisfies sup

ξ ∈[−L,+L]

   u(ξ ) + v(ξ )  M, inf

ξ ∈[−L,+L]

v(ξ )  δ.

(4.9)

Then there exists a solution of (4.1), (4.2). The following lemma is crucial in establishing the estimates (4.9). Lemma 4.2. Let (r(ξ ), v(ξ )) be a nonconstant solution of (4.7), (4.8) in [−L, L] for some 0 < μ  1 and L  1. Assume r− > r+ . Then r(ξ ) is strictly decreasing in [−L, L] while v(ξ ) satisfies one of the following: (1) v(ξ ) is monotone in [−L, L], (2) v(ξ ) has only maximum points at some τ in [−L, L]. Proof. Let (r(ξ ), v(ξ )) be a nonconstant solution of (4.7), (4.8). From (4.7), we have

  vξ rξ = εrξ ξ . −ξ + μφ r(ξ ) − 2ε v

(4.10)

Using (4.8), from (4.10), we obtain that  ξ  

1 r(ξ ) = λ exp A(s)ds , ε



(4.11)

−L

v

where A(ξ ) = −ξ + μφ(r(ξ )) − 2ε v(ξξ ) and λ is given by L r(L) − r(−L) = λ −L

 t  1 exp A(s)ds dt. ε

(4.12)

−L

From (4.11), (4.12) and r− > r+ , we deduce that r(ξ ) is monotone decreasing in [−L, L].

4382

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Next we suppose v(ξ ) has a critical point τ in [−L, L]. Then v (τ ) = 0. By the second equation of (4.7), we have εv

(τ ) = μv(τ )φ (r(τ ))r (τ ) < 0 since r (τ ) < 0, v(τ ) > 0 and φ (r(τ )) > 0. Thus v(ξ ) must take a maximum value at τ . The proof is complete. 2 From Lemma 4.2, one can get that r(ξ ) is uniformly bounded in [−L, L] with respect to 0 < μ  1, L  1 and ε > 0. Apparently, v(ξ ) is uniformly bounded in [−L, L] with respect to μ, L for case (1) in Lemma 4.2. It remains to estimate v(ξ ) for case (2) in Lemma 4.2. Lemma 4.3. For case (2) in Lemma 4.2, there exist positive constants M0 and δ, independent of 0 < μ  1 and L  1, such that δ  v(ξ )  M0 . Proof. Motivated by [3], we first prove β v(ξ )dξ  N + v(β ¯ − α)

(4.13)

α

for every interval (α, β) ⊂ (−L, L), where v¯ = max{v+ , v− } and N = v(φ(r ¯ − ) − φ(r+ )). In fact, we set    θ1 = sup ξ ∈ [−L, α)v(ξ )  v¯ if v(α) > v, ¯ while    θ1 = inf ξ ∈ (α, β)v(ξ )  v¯ if v(α)  v¯ (if this set is empty, (4.13) is automatically satisfied). Similarly, we set    θ2 = inf ξ ∈ (β, L]v(ξ )  v¯ if v(β) > v, ¯ while    θ2 = sup ξ ∈ (α, β)v(ξ )  v¯ if v(β)  v¯ (if this set is empty, (4.13) is automatically satisfied). Since v(θ1 ) = v(θ2 ) = v, ¯ we have β α



 v(ξ ) − v¯ dξ 

θ2 



θ2

v(ξ ) − v¯ dξ = −

θ1

ξ vξ dξ.

(4.14)

θ1

Noticing that vξ (θ1 )  0 and vξ (θ2 )  0, integrating the second equation in (4.7) over (θ1 , θ2 ), we obtain that

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

β

4383



      v(ξ ) − v¯ dξ  εvξ (θ2 ) − εvξ (θ1 ) + μv¯ φ r(θ1 ) − φ r(θ2 )

α

   v¯ φ(r− ) − φ(r+ ) = N,

(4.15)

which implies (4.13). We now apply (4.13) to estimate v(ξ ) from above for case (2) in Lemma 4.2. By (4.13), it follows that v  v(ξ )  v¯ +

N , |τ − ξ |

ξ ∈ [−L, L] \ {τ },

(4.16)

where v = min{v− , v+ }. Without loss of generality we assume that v(τ ) > v. ¯ We fix β ∈ [−L, τ ) such that v(β) = v. ¯ For any ξ1 ∈ [β, τ ) we let ξ2 be a point in (τ, L] with the property v(ξ2 ) = ¯ Integrating the second equation in (4.7) over (ξ1 , ξ2 ), v(ξ1 ) (such a point exists since v(L)  v). we obtain ξ2 εvξ (ξ2 ) − εvξ (ξ1 ) = −

     ξ vξ dξ + μv(ξ1 ) φ r(ξ2 ) − φ r(ξ1 ) .

(4.17)

ξ1

We note that vξ (ξ2 )  0 and ξ2 ξ vξ dξ = ξ1

ξ2 

 v(ξ1 ) − v(ξ ) dξ  0.

ξ1

Therefore, (4.17) gives that        εvξ (ξ1 )  μv(ξ1 ) φ r(ξ1 ) − φ r(ξ2 )  v(ξ1 ) φ(r− ) − φ(r+ ) .

(4.18)

Integrating (4.18) with respect to ξ1 over (ξ0 , τ ), we deduce that   v(τ )  v(ξ0 ) exp A(τ − ξ0 ) ,

(4.19)

+) where A = φ(r− )−φ(r , β  ξ0  τ . ε If τ − β < 1, we choose ξ0 = β in (4.19) to get

v(τ )  ve ¯ A.

(4.20)

If τ − β  1, by choosing ξ0 = τ − 1, from (4.16) and (4.19), we obtain v(τ )  v¯ + The lemma is proved. 2

  N exp A(τ − ξ0 ) = eA (v¯ + N ). τ − ξ0

(4.21)

4384

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

From the above conclusions, we have that u(ξ ) = r(ξ )v(ξ ) is uniformly bounded on [−L, L] with respect to μ ∈ (0, 1] and L  1. Then the estimates (4.9) hold and the existence of a smooth solution (uε (ξ ), v ε (ξ )) of (4.1), (4.2) on (−∞, +∞) for every fixed ε > 0 can be established. 5. Existence of solutions to (1.1), (1.2) when uε (ξ ) and v ε (ξ ) are uniformly bounded In Section 4, we have obtained the existence of a smooth solution (uε (ξ ), v ε (ξ )) of viscous problem (4.1), (4.2) on (−∞, +∞) for every fixed ε > 0. Moreover, the same results in Lemma 4.2 are also true for the solution (r ε (ξ ), v ε (ξ )) of (4.3), (4.4). If {(uε (ξ ), v ε (ξ ))|0 < ε < 1} is uniformly bounded in ε, by means of the following result induced directly from Theorem 3.2 in [3], one can get the existence of a solution to Riemann problem (1.1), (1.2). Theorem 5.1. For every fixed ε > 0, let (uε (ξ ), v ε (ξ )) be a solution of (4.1), (4.2). Assume that {(uε (ξ ), v ε (ξ ))|0 < ε < 1} is of uniformly bounded variation. Then {(uε (ξ ), v ε (ξ ))|0 < ε < 1} possesses a subsequence which converges a.e. on (−∞, +∞) to function (u(ξ ), v(ξ )) of bounded variation. And (u(ξ ), v(ξ )) provides a weak solution to (1.1), (1.2). When r− > r+ , from Lemma 4.2, we know that r ε (ξ ) is always decreasing, v ε (ξ ) is a monotone function on [−L, L] for case (1). Furthermore, uε (ξ ) has a bounded variation since uε (ξ ) = r ε (ξ )v ε (ξ ). Thus, from Theorem 5.1, we obtain that Riemann problem (1.1), (1.2) has a weak solution in this case. It now remains to discuss case (2) in Lemma 4.2. At this time, v ε (ξ ) may tend to infinity as ε → 0. Let τ ε be a maximum point of v ε (ξ ) on (−∞, +∞) and τ ε → τ0 , |τ0 |  ∞ as ε → 0 (pass to a further subsequence if necessary). Theorem 5.2. If |τ0 | = ∞, then {(uε (ξ ), v ε (ξ ))|0 < ε < 1} is uniformly bounded. Proof. Assume firstly that τ0 = +∞. By (4.16), we obtain that v  v ε (ξ )  v¯ + N

(5.1)

for any ξ ∈ (−∞, a) and ε small, where a is any fixed real number. We take a = 2 and choose ξ0 ∈ [1, 2] such that  

0  v ε (ξ0 ) = v ε (2) − v ε (1)  v¯ + N − v.

(5.2)

We can take ε to be so small that (ξ0 − φ(r ε (ξ0 )))v ε (ξ0 ) + ε(v ε (ξ0 ))

 1, τ ε − φ(r ε (τ ε )) (τ ε − ξ0 )v¯ + N  2v¯ + 1, τ ε − φ(r ε (τ ε )) because φ(r+ )  φ(r ε (τ ε ))  φ(r− ) and τ ε → +∞ as ε → 0.

(5.3)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4385

Integrating the second equation in (4.3) over (ξ0 , τ ε ), one can get that   (ξ0 − φ(r ε (ξ0 )))v ε (ξ0 ) + ε(v ε (ξ0 )) + v τε = τ ε − φ(r ε (τ ε )) ε

 τε ξ0

v ε (ξ )dξ

(ξ0 − φ(r ε (ξ0 )))v ε (ξ0 ) + ε(v ε (ξ0 )) + (τ ε − ξ0 )v¯ + N τ ε − φ(r ε (τ ε ))  2 + 2v¯ 

(5.4)

by virtue of (5.2), (5.3) and (4.13). Thus, the theorem is proved for τ0 = +∞. The case τ0 = −∞ can be treated similarly. Moreover, noticing uε (ξ ) = r ε (ξ )v ε (ξ ), we immediately deduce that {uε (ξ )|0 < ε < 1} is uniformly bounded in ε. 2 Now we turn to consider the case |τ0 | < ∞. The singularity point of (4.1) is given by the unique solution of ξ = φ(r ε (ξ )) and denoted by ξαε . We set ξα = lim ξαε = limε→0 φ(r ε (ξαε )). ε→0

Then φ(r+ )  ξα  φ(r− ) since φ(r+ )  φ(r ε (ξ ))  φ(r− ). The point ξα plays an important role in the following argument. Theorem 5.3. ξα is defined as above. Then {(uε (ξ ), v ε (ξ ))|0 < ε < 1} is uniformly bounded in ε if τ0 = ξα . Proof. Suppose τ0 > ξα . By integrating the second equation in (4.3) over (τ ε , τ0 + 1), it follows that         ε v ε (τ0 + 1) = τ ε − φ r ε τ ε v ε τ ε 





τ0 +1

− τ0 + 1 − φ r (τ0 + 1) v (τ0 + 1) + ε

ε

v ε (ξ )dξ.

(5.5)

τε

Notice that (v ε (τ0 + 1))  0 and τ ε − φ(r ε (τ ε )) = τ ε − ξαε + (φ(r ε (ξαε )) − φ(r ε (τ ε ))) = (1 − φ (r ε (θ ε ))(r ε (θ ε )) )(τ ε − ξαε )  12 (τ0 − ξα ) > 0 for ε small since φ (r ε (θ ε ))(r ε (θ ε ))  0 and τ ε − ξαε → τ0 − ξα , where θ ε is between τ ε and ξαε . Therefore, (5.5) gives   (τ0 + 1 − φ(r ε (τ0 + 1)))v ε (τ0 + 1) − v τε  τ ε − φ(r ε (τ ε ))  2  τ0 + 1 − φ(r+ ) (v¯ + 2N )  τ 0 − ξα ε

 τ0 +1 τε

v ε (ξ )dξ

by virtue of (4.16). Thus v ε (ξ ) is uniformly bounded in ε. The case τ0 < ξα can be addressed similarly. Moreover, we can obtain that {uε (ξ )|0 < ε < 1} is uniformly bounded in ε if τ0 = ξα by uε (ξ ) = r ε (ξ )v ε (ξ ). The proof is finished. 2 Similarly to (4.16), we also have v  v ε (ξ )  v¯ +

|τ ε

N , − ξ|

  ξ ∈ (−∞, +∞) \ τ ε .

(5.6)

4386

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Given η > 0, from (5.6) we know v ε (ξ )  v¯ + 2N η , ξ ∈ (−∞, ξα − η) ∪ (ξα + η, +∞) for ε small ε if τ → τ0 = ξα as ε → 0. Combining this and Theorems 5.2–5.3, we easily get Theorem 5.4. When r− > r+ , {(uε (ξ ), v ε (ξ ))|0 < ε < 1} is uniformly bounded and of uniformly bounded variation over the interval (−∞, ξα − η) ∪ (ξα + η, +∞) for any fixed η > 0. Theorem 5.5. ξα is defined as above. Then for any η > 0, lim rξε (ξ ) = 0,

for |ξ − ξα |  η,

lim r ε (ξ ) = r+ ,

for ξ  ξα + η,

lim r ε (ξ ) = r− ,

for ξ  ξα − η

ε→0

ε→0

ε→0

(5.7)

uniformly in the above intervals. Proof. By Theorem 5.4, it is easily seen that there exists a positive M1 (independent of ε) such that



η η 0 < v  v ε (ξ )  M1 , ξ ∈ −∞, ξα − ∪ ξα + , +∞ (5.8) 2 2 for ε small. We set ξ0 to be a point in [ξα + η2 , ξα +

3η 4 ]

such that





 4 ε 3η η ε r (ξ0 ) = r ξα + − r ξα + η 4 2



ε

which yields that  

4 − (r− − r+ )  r ε (ξ0 )  0. η

(5.9)

From (4.3), we have   rξ vξ −ξ rξ + rξ φ r(ξ ) = εrξ ξ + 2ε , v then  ξ    ε  (v ε (ξ0 ))2 −s + φ(r ε (s)) r (ξ ) = r (ξ0 ) exp ds , ε (v ε (ξ ))2



ε

ξ  ξ0 .

(5.10)

ξ0

Observe that           −s + φ r ε (s) = φ r ε (s) − φ r ε ξαε + ξαε − s      

  η = φ r ε θ ε r ε θ ε − 1 s − ξαε  − , 4

(5.11)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4387

since s − ξαε → s − ξα  ξ0 − ξα  12 η as ε → 0, and φ (r ε (θ ε ))(r ε (θ ε ))  0, here θ ε is between s and ξαε . Combining (5.10) with (5.8), (5.9) and (5.11), one can easily get that

2  ε   4  r (ξ )   (r− − r+ ) M1 exp − η (ξ − ξ0 ) for ξ  ξ0 , η v 4ε

(5.12)

which implies that rξε (ξ ) → 0 uniformly on ξ  ξα + η. Letting ξ  ξα + η and noting

r+ − r (ξ ) = ε

+∞ 



r ε (s) ds,

ξ

we immediately have lim r ε (ξ ) = r+

ε→0

(5.13)

uniformly by virtue of (5.12). The results for ξ  ξα − η can be obtained in a completely similar way. 2 6. Existence of solutions to (1.1), (1.2) when uε (ξ ) and v ε (ξ ) may tend to infinity as τ0 = ξα Now we discuss the case when v ε (ξ ) may tend to infinity as τ0 = ξα . We prove that the family {(uε (ξ ), v ε (ξ ))|0 < ε < 1} of (4.1), (4.2) can also generate solution of Riemann problem (1.1), (1.2). When r− > r+ , the Riemann solution of (1.1), (1.2) contains a delta shock wave. Theorem 6.1. Let (uε (ξ ), v ε (ξ )) be a solution of (4.1), (4.2) and let v ε (ξ ) have its maximum at τ ε , τ ε → τ0 = ξα as ε → 0. Then the sequence {(uε (ξ ), v ε (ξ ))|0 < ε < 1} possesses a subsequence that converges a.e. to a pair of functions (u(ξ ), v(ξ )) on (−∞, ∞). Proof. By Theorem 5.4, v ε (ξ ) is uniformly bounded in ε on the finite domain I2 = [−2, τ0 − 1 1 2 ] ∪ [τ0 + 2 , 2] (without loss of generality, we assume that |τ0 | < 1). By Helly’s theorem, there exists a convergent subsequence of {v ε (ξ )} (still denoted by original one). Similarly, we can extract a convergent subsequence of {v ε (ξ )} on I3 = [−3, τ0 − 13 ] ∪ [τ0 + 13 , 3]. Continue the procedure on each In = [−n, τ0 − n1 ] ∪ [τ0 + n1 , n], n = 4, 5, · · · . Finally, extract the diagonal element at each enumerated sequence. The sequence of diagonal elements is convergent at each ξ = τ0 to a function v(ξ ) defined on (−∞, τ0 ) ∪ (τ0 , +∞). From Theorem 5.4, we observe that uε (ξ ) has the same property as v ε (ξ ). In a similar way as the above, we can also obtain a sequence which, at each ξ = τ0 , converges to a function u(ξ ) defined on (−∞, τ0 ) ∪ (τ0 , +∞). This theorem is verified. 2 Theorem 6.2. Assume that τ0 = ξα and r− > r+ . Then for any δ > 0, v(ξ ) = Hv (ξ − σ ) + w0 δ(ξ − σ ),

(6.1)

4388

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

where Hv (x) =

v− , x < 0, v+ , x > 0,

δ(x) is the Dirac function supported at x = 0, w0 = −σ [v] + [φ(r)v] is the strength of δ(x), and σ = ξα = limε→0 ξαε = limε→0 φ(r ε (ξαε )). Proof. We set ξ1 , ξ2 with the property −∞ < ξ1 < ξ2 < σ − δ. For each ψ ∈ C01 (ξ1 , ξ2 ), from the second equation in (4.3), it follows that ξ2



ξ2

v ψ dξ = ε

ε ξ1

    v ε ψ + ξ ψ − φ r ε ψ dξ.

(6.2)

ξ1

By (4.16), we have v  v ε (ξ )  v¯ +

N , δ

for ξ ∈ (ξ1 , ξ2 ), ε small.

(6.3)

Letting ε → 0 in (6.2), using (4.13), (6.3), (5.7) and the Lebesgue dominated convergence theorem, we get that ξ2

  v ψ + ξ ψ − φ(r− )ψ dξ = 0,

for each ψ(ξ ) ∈ C01 (ξ1 , ξ2 ),

(6.4)

ξ1

which leads to ξ2

v(ξ )ϕ (ξ )dξ = 0,

for each ϕ(ξ ) ∈ C01 (ξ1 , ξ2 )

(6.5)

ξ1

by setting ϕ(ξ ) = ψ(ξ )(ξ − φ(r− )) since ξ − φ(r− ) = 0 for ξ ∈ (ξ1 , ξ2 ). Therefore v(ξ ) is a constant on (−∞, σ − δ). Again v(−∞) = v− , thus v(ξ ) = v− . Analogously, for ξ1 , ξ2 satisfying σ + δ < ξ1 < ξ2 < +∞, we have also the similar formulae (6.3) and (6.4), in which we only need to replace r− by r+ . Furthermore, v(ξ ) = v+ for ξ > σ +δ. Combining this with Theorem 5.5, we find that (r(ξ ), v(ξ )) share the same discontinuity point ξ = ξα on (−∞, +∞). However, the Rankine–Hugoniot condition no longer holds for ξ = ξα since v− v+ (r− − r+ )(φ(r− ) − φ(r+ )) = 0. Thus we have to study in more detail the limiting behavior of v ε (ξ ) in the neighborhood of ξ = ξα = σ as ε → 0 in what follows. For the purpose, we take ψ(ξ ) ∈ C0∞ (ξ1 , ξ2 ) (ξ1 < σ < ξ2 ) with the property ψ(ξ ) = ψ(σ ) for ξ in a small neighborhood of σ . From (6.2), using (4.13), we have that ξ2 lim

ε→0 ξ1

    v ε ψ + ξ ψ − φ r ε ψ dξ = 0.

(6.6)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4389

For any α1 and α2 near σ with α1 < σ < α2 , from (6.6) and Theorem 5.5, we obtain that ξ2 lim

ε→0

   v ε ξ − φ r ε ψ dξ

ξ1

α1 = lim

  v ξ − φ r ε ψ dξ + lim ε

ε→0



ε→0

α1 v−

   v ε ξ − φ r ε ψ dξ

α2

ξ1

=

ξ2



 ξ − φ(r− ) ψ dξ +

ξ2

  v+ ξ − φ(r+ ) ψ dξ

α2

ξ1

 = v+ φ(r+ ) − v− φ(r− ) − (v+ α2 − v− α1 ) ψ(σ ) − 

α1

ξ2 v− ψ(ξ )dξ −

v+ ψ(ξ )dξ. α2

ξ1

Letting α1 → σ − and α2 → σ + , we conclude that ξ2

  

 v ξ − φ r ε ψ dξ = σ [v] − φ(r)v ψ(σ ) − ε

lim

ε→0



ξ1

ξ2 Hv (ξ − σ )ψ(ξ )dξ.

(6.7)

ξ1

Thus, it follows from (6.6) that

lim

ξ2 

 

 v ε − Hv (ξ − σ ) ψ(ξ )dξ = −σ [v] + φ(r)v ψ(σ )

ε→0

(6.8)

ξ1

holds for every test function ψ ∈ C0∞ (ξ1 , ξ2 ). Therefore we can define v(ξ ) = Hv (ξ − σ ) + w0 δ(ξ − σ ) on (−∞, +∞). This completes the proof. 2 Theorem 6.3. Assume that τ0 = ξα and r− > r+ . Then r = Hr (ξ − σ ) and φ(r(σ )) = σ is required so that (1.1) holds in the sense of distributions, where r , x < 0, Hr (x) = − r+ , x > 0. Proof. From Theorem 5.5, we can easily conclude that r = Hr (ξ − σ ). Again, by Theorems 6.1–6.2, we have shown the existence of a solution (r(ξ ), v(ξ )) of the boundary value problem 

and (4.4).

  −ξ(rv)ξ + φ(r)rv ξ = 0,   −ξ vξ + φ(r)v ξ = 0

4390

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Now we proceed to determine the value of φ(r(ξ )) at ξ = σ . Letting ξ1 < σ < ξ2 , ψ ∈ C0∞ (ξ1 , ξ2 ), from the second equation in the above, it follows that ξ2

  v ψ + ξ ψ − φ(r)ψ dξ = 0.

(6.9)

 

  Hv (ξ − σ ) + −σ [v] + φ(r)v δ(ξ − σ ) ψ + ξ ψ − φ(r)ψ dξ = 0.

(6.10)

ξ1

That is, ξ2  ξ1

The integral in (6.10) equals to σ

  Hv (ξ − σ ) ψ + ξ ψ − φ(r)ψ dξ +

ξ2

  Hv (ξ − σ ) ψ + ξ ψ − φ(r)ψ dξ

σ

ξ1

+

ξ2 

  −σ [v] + φ(r)v δ(ξ − σ ) ψdξ

ξ1

+

ξ2 

   −σ [v] + φ(r)v δ(ξ − σ ) ξ ψ − φ(r)ψ dξ

ξ1

    σ ξ 

 = v− ψ ξ − φ(r− ) ξ + v+ ψ ξ − φ(r+ ) σ2 + −σ [v] + φ(r)v ψ(σ ) 1 

   + −σ [v] + φ(r)v σ − φ r(σ ) ψ (σ ) 

   = −σ [v] + φ(r)v σ − φ r(σ ) ψ (σ ), therefore 

   −σ [v] + φ(r)v σ − φ r(σ ) ψ (σ ) = 0.

(6.11)

  rδ = lim r ε ξαε = r(σ ),

(6.12)

Denoting

ε→0

we can easily get that r+ < rδ < r− . Noticing    σ = ξα = lim ξαε = lim φ r ε ξαε = φ(rδ ), ε→0

ε→0

(6.13)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4391

we have (6.14)

φ(r+ ) < σ < φ(r− ).

From (6.11) and (6.14), we have φ(r(σ )) = σ since −σ [v] + [φ(r)v] = 0 and ψ is arbitrary. The proof is finished. 2 Theorem 6.4. Assume that τ0 = ξα . Then for any δ > 0, u(ξ ) = H (ξ − σ ) + rδ w0 δ(ξ − σ )

(6.15)

for r− > r+ , where H (x) =

r− v− , r+ v+ ,

x < 0, x > 0.

(6.16)

Proof. By Theorems 6.2–6.3 and uε (ξ ) = r ε (ξ )v ε (ξ ), we get that u(ξ ) = r− v−

for ξ < σ − δ;

u(ξ ) = r+ v+

for ξ > σ + δ.

(6.17)

Now we take ψ(ξ ) ∈ C0∞ (ξ1 , ξ2 ), ξ1 < σ < ξ2 , and ψ(ξ ) = ψ(σ ) for ξ in a small neighborhood Ω˜ of point σ . From the first equation of (4.3), we have ξ2



ξ2

r v ψ dξ = ε ε

ε ξ1

    r ε v ε ψ + ξ ψ − φ r ε ψ dξ.

(6.18)

ξ1

In a completely analogous way as the proof of Theorem 6.2, from (6.18) we can obtain that ξ2 

 

 r ε v ε − H (ξ − σ ) ψ(ξ )dξ = −σ [rv] + rvφ(r) ψ(σ )

lim

ε→0

(6.19)

ξ1

for all ψ ∈ C0∞ [ξ1 , ξ2 ]. Thus r ε v ε converges weakly star to the function r(ξ )v(ξ ) given by (6.16) plus a Dirac delta function with strength −σ [rv] + [rvφ(r)]. If we take the test function as ψ/(˜r ε + ν) in (6.18), where r˜ ε is a modified function satisfying ε ˜ and let ν → 0, then we can find the other formula r (σ ) in Ω˜ and r ε outside Ω,

lim

ξ2 

 

 v ε − Hv (ξ − σ ) ψ(ξ )dξ · rδ = −σ [rv] + rvφ(r) ψ(σ ),

ε→0

(6.20)

ξ1

where ψ ∈ C0∞ (ξ1 , ξ2 ). Comparing (6.20) with (6.8), we have

 w0 rδ = −σ [rv] + rvφ(r) .

(6.21)

4392

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Noting u(ξ ) = r(ξ )v(ξ ), thus we can define u(ξ ) = H (ξ − σ ) + rδ w0 δ(ξ − σ ) on (−∞, +∞). The proof of Theorem 6.4 is completed. 2 Then, from (6.13), (6.8) and (6.21), it follows that ⎧ ⎨ σ = φ(rδ ),

 w0 = −σ [v] + vφ(r) ,

 ⎩ w0 rδ = −σ [u] + uφ(r) .

(6.22)

From all the above results, we arrive at the following theorem. Theorem 6.5. When r− > r+ , let (uε , v ε )(t, x) be the solution of (1.4), (1.2). Then the limit functions u(t, x) and v(t, x) of uε (t, x) and v ε (t, x) exist in the sense of distributions, and (u, v)(t, x) solves (1.1), (1.2). The solution (u, v)(t, x) can be explicitly expressed as ⎧ ⎨ (u− , v− ), (u, v)(t, x) = (rδ w(t), w(t))δ(x − x(t)), ⎩ (u+ , v+ ),

x < σ t, x = σ t, x > σ t,

(6.23)

where w(t) = w0 t , σ , w0 and rδ are determined uniquely by (6.22), (6.14). This theorem shows the stability of the delta shock wave for (1.1) under viscous perturbations. In addition, from Theorems 6.2–6.4, we see that the limit functions u(ξ ) and v(ξ ) of uε (ξ ) and v ε (ξ ) are the sums of a step function and a Dirac δ-function with strengths rδ w0 and w0 , respectively; the limit function r(ξ ) of r ε (ξ ) is a step function and φ(r(σ )) = σ is reduced. All these coincide with the delta shock solution obtained in Section 3, and confirm the generalized Rankine–Hugoniot relation (3.5), (3.6). Meanwhile, just as remarked in [8], from Definitions 3.1–3.2 and (3.11), as well as (6.22), one can observe that the factor  1 2 arising 1+σ

in the line integrals is a part of the weights, but it has no contribution to the strengths of the delta shock wave. 7. Two typical examples In order to illustrate the application of our results and proofs in the above sections, we present here two typical examples. In these two models, we are more interested in the delta shock waves. Example 1. Engquist and Runborg’s model ⎧

⎪ u2 ⎪ ⎪ + = 0, u √ ⎨ t 2 + v2 u x

⎪ uv ⎪ ⎪ + = 0, v √ ⎩ t u2 + v 2 x which is deduced by taking φ(u, v) =

 u u2 +v 2

(7.1)

in system (1.1). In their work of [1,2], the Riemann

problem of (7.1) has not been solved since 1996.

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

In this situation, λ =

 u u2 +v 2

=



u v

1+( uv )2

When the Riemann initial data satisfy cavitation one, which can be expressed as

= u− v−

⎧ (u− , v− ), ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ (u, v)(ξ ) = Cav, ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (u+ , v+ ),

 r 1+r 2 < uv++ ,

4393 3

 φ(r), φ (r) = (1 + r 2 )− 2 > 0. as indicated in Section 2, the solution is a

−∞ < ξ <  u− 2 u2− +v−  u+ 2 u2+ +v+

 u− , 2 u2− +v−  u+

ξ 

2 u2+ +v+

,

(7.2)

< ξ < +∞.

When uv++ < uv−− , we can get from (3.5)–(3.6) that the generalized Rankine–Hugoniot relation of the delta shock wave for the system (7.1) reads ⎧ dx ⎪ ⎪ = σ, ⎪ ⎪ dt ⎪ ⎪ √   ⎪ ⎨ d 1 + σ 2 rδ w(t) u2 = −σ [u] + √ , dt ⎪ u2 + v 2 ⎪ ⎪ √   ⎪ ⎪ ⎪ d 1 + σ 2 w(t) uv ⎪ ⎩ = −σ [v] + √ , dt u2 + v 2

(7.3)

and √

r

   2

rδ = = σ, 1 + r x=x(t) 1 + rδ2

(7.4)

which is equivalent to ⎧   u2 ⎪ 2 r w = −σ [u] + √ ⎪ 1 + σ , ⎪ δ 0 ⎪ ⎪ u2 + v 2 ⎪ ⎪   ⎪ ⎨ uv 1 + σ 2 w0 = −σ [v] + √ , ⎪ u2 + v 2 ⎪ ⎪ ⎪ rδ ⎪ ⎪ , ⎪σ =  ⎪ ⎩ 1 + rδ2

(7.5)

where w(t) = w0 t , σ , w0 and rδ are constants to be determined later. From (7.5), we have     r2 rδ uv u2  δ [v] −  [u] − rδ √ + √ = 0, u2 + v 2 u2 + v 2 1 + rδ2 1 + rδ2

(7.6)

which is a fourth-order equation of one variable. Under the entropy condition (3.8), by setting

4394

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

f (rδ ) = 

rδ2 1 + rδ2

[v] − 

    uv u2 [u] − rδ √ + √ , u2 + v 2 u2 + v 2 1 + r2 rδ

(7.7)

δ

we obtain that f

u+ v+

u− f v−



= v−



u− u+ − v− v+

= −v+



u− u+ − v− v+





u− 2 u2− + v−



−

u− 2 u2− + v−



u+

> 0,

2 u2+ + v+

−

u+

< 0.

(7.8)

2 u2+ + v+

Differentiating f (rδ ) with respect to rδ yields



u− rδ rδ u+ f (rδ ) = v−  − − + v+  2 2 u2− + v− u2+ + v+ 1 + rδ2 1 + rδ2



u− u+ 1  − rδ < 0. + v− rδ − + v+ v− v+ ( 1 + r 2 )3

(7.9)

δ

With the help of zero point theorem, there exists one and only one zero point of f (rδ ) in ( uv++ , uv−− ). Thus Eq. (7.6) possesses a unique solution rδ ∈ ( uv++ , uv−− ). Returning to (7.5) we solve the σ and w0 uniquely. Therefore, the unique delta shock solution of (7.1), (1.2) can be expressed as (3.16), where w(t) = w0 t and constants σ , w0 , and rδ are determined uniquely by (7.5) and (3.8). Example 2. Consider the system ⎧ ⎨



u2 = 0, ut + v x ⎩ vt + ux = 0,

(7.10)

which is just a simple version of system (1.1) as φ(u, v) = uv . At this moment, λ = uv = r  φ(r), φ (r) = 1 > 0. Similarly, when uv++ < uv−− , from (3.5)–(3.6), one can directly write the generalized Rankine– Hugoniot relation of the delta shock wave as follows ⎧ dx ⎪ ⎪ = σ, ⎪ ⎪ dt ⎪ ⎪ √ ⎪  2 ⎨ d 1 + σ 2 r w(t) u δ = −σ [u] + , ⎪ dt v ⎪ ⎪ √ ⎪ ⎪ 2 ⎪ ⎪ ⎩ d 1 + σ w(t) = −σ [v] + [u], dt and

(7.11)

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4395

r|x=x(t) = rδ = σ,

(7.12)

⎧  2  u ⎪ ⎪ 2 r w = −σ [u] + ⎪ 1 + σ , δ 0 ⎨ v  ⎪ 1 + σ 2 w0 = −σ [v] + [u], ⎪ ⎪ ⎩ σ = rδ ,

(7.13)

which can be reduced to

where w(t) = w0 t , σ , w0 and rδ are constants to be determined later. From (7.13), it follows that  rδ2 [v] − 2rδ [u] +

 u2 = 0. v

(7.14)

When [v] = 0, (7.14) is a quadratic equation of one variable. Setting  f (rδ ) = rδ2 [v] − 2rδ [u] +

 u2 , v

(7.15)

an easy calculation gives f (r+ ) = v− (r− − r+ )2 > 0,

f (r− ) = −v+ (r− − r+ )2 < 0,

f (rδ ) = 2v− (rδ − r− ) + 2v+ (r+ − rδ ) < 0. So Eq. (7.14) has a unique solution rδ ∈ ( uv++ , uv−− ). Of course, one can also solve (7.14) to obtain σ = rδ =

[u] − ( uv−− −

u+ √ v+ ) v − v +

[v]

.

(7.16)

Then

1 u− u+ √ w(t) = w0 t = √ − v− v+ t. v+ 1 + σ 2 v−

(7.17)

When [v] = 0, one can easily get from (7.14) that σ = rδ =

u− + u+ , 2v−

1 w(t) = √ (u− − u+ )t. 1 + σ2

(7.18)

Under condition uv−− > uv++ , obviously, rδ satisfies the entropy condition (3.8). Therefore, the unique delta shock solution of (7.10), (1.2) is given in (3.16), where σ = rδ and w(t) = w0 t are shown as (7.16) and (7.17), or (7.18).

4396

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Fig. 3. Numerical results of u (left) and v (right) for Data 1 at t = 11.

8. Numerical simulations of Riemann solutions for (7.1), (1.2) In this section, to complement the theoretical analysis established in this paper, we calculate the Riemann solutions of (7.1), (1.2) by employing the weighted essentially non-oscillatory (WENO) scheme [6] with the mesh 100. Here, many more numerical tests have been performed to make sure that what we present are not numerical artifacts. In what follows, we exhibit the numerical results by four cases. Case 1. r− < r+ . We take three groups of initial data as follows Data 1:

u− = 0.20,

v− = 0.55,

u+ = 0.50,

v+ = 0.45,

Data 2:

u− = −0.45,

v− = 0.50,

u+ = −0.10,

Data 3:

u− = −0.10,

v− = 0.55,

u+ = 0.20,

v+ = 1.65, v+ = 0.50.

The numerical results are presented in Figs. 3–7. Figs. 3–7 clearly show the structure of a two-contact-discontinuity solution including a cavitation state. Along with the time increasing, both state variables u and v in the intermediate state tend to zero, and a cavitation is forming inside, while uv is a linear function. Case 2. r− > r+ > 0. We choose the following initial data Data 4: u− = 0.30,

v− = 0.15,

u+ = 0.20,

v+ = 0.40.

Case 3. 0 > r− > r+ . The initial data are Data 5:

u− = −0.10,

v− = 0.60,

The numerical results are presented by Figs. 8–10.

u+ = −0.30,

v+ = 0.50.

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Fig. 4. Numerical results of u (left) and v (right) for Data 2 at t = 20.

Fig. 5. Numerical results of uv for Data 1 (left) and Data 2 (right).

Fig. 6. Numerical results of u (left) and v (right) for Data 3 at t = 10.

4397

4398

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Fig. 7. Numerical result of uv for Data 3 at t = 10.

Fig. 8. Numerical results of u (left) and v (right) for Data 4 at t = 3.

Fig. 9. Numerical results of u (left) and v (right) for Data 5 at t = 2.5.

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4399

Fig. 10. Numerical results of uv for Data 4 (left) and Data 5 (right).

Fig. 11. Numerical results of u (left) and v (right) for Data 6 at t = 3. u2 ], u2 +v 2

Case 4. r− > 0 > r+ . According to the sign of [  be

we select three groups of initial data to

Data 6: u− = 0.10,

v− = 0.20,

u+ = −0.60,

v+ = 0.50,

Data 7: u− = 0.50,

v− = 0.40,

u+ = −0.10,

v+ = 0.20,

Data 8: u− = 0.10,

v− = 0.25,

u+ = −0.10,

v+ = 0.25.

For Data 6, we can easily get that know rδ < 0. Similarly,



u2− 2 2 u− +v−

u2  − 2 u2− +v−

−



−

u2+ 2 2 u+ +v+



u2+ 2 u2+ +v+

< 0. From the arguments in Section 3, we

> 0, rδ > 0 for Data 7, and

rδ = 0 for Data 8. The numerical results are presented by Figs. 11–15.



u2− 2 2 u− +v−

−



u2+ 2 2 u+ +v+

= 0,

4400

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

Fig. 12. Numerical results of u (left) and v (right) for Data 7 at t = 3.

Fig. 13. Numerical results of uv for Data 6 (left) and Data 7 (right).

Fig. 14. Numerical results of u (left) and v (right) for Data 8 at t = 5.

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

4401

Fig. 15. Numerical result of uv for Data 8 at t = 5.

For the case r− > r+ , one can clearly observe the structure of a delta shock solution. As the time increases, both state variables u and v develop extreme concentrations which indicate that a delta shock wave forms. Moreover, uv is a step function. From Section 3, we have w(t) > 0 for t > 0, thus the narrow region of v is always up sharp-angled. The narrow region of u is up or down sharp-angled dependent on the signs of rδ . In this situation, we have rδ > 0 when 2 r− > r+ > 0 and r− > 0 > r+ with [  u2 2 ] > 0. Fig. 8 and Fig. 12 confirm this fact. While u +v

rδ < 0 when 0 > r− > r+ and r− > 0 > r+ with [ 

u2 ] < 0. u2 +v 2

Fig. 9 and Fig. 11 correspond to

this result. When rδ = 0, at this moment, u may degenerate into a classical shock which is shown in Fig. 14. To sum up, all of the above numerical results clearly reveal the two kinds of geometric structures of analytical solutions obtained in Section 2. Acknowledgments Hanchun Yang would like to thank Professor Tong Zhang, his advisor, for constant encouragement and assistance. Hanchun Yang also thanks Professor Zhouping Xin for his kind assistance in the Chinese University of Hong Kong. The authors really appreciate the referees for their valuable comments and suggestions to revise this paper. References [1] B. Engquist, O. Runborg, Multiphase computations in geometrical optics, J. Comput. Appl. Math. 74 (1996) 175–192. [2] B. Engquist, Multiscale and multiphase methods for wave propagation, Doctoral dissertation, Royal Institute of Technology, Stockholm, Sweden, 1998. [3] C.M. Dafermos, Solutions of the Riemann problem for a class of hyperbolic systems of conservation laws by the viscosity method, Arch. Ration. Mech. Anal. 52 (1973) 1–9. [4] M. Slemrod, A. Tzavaras, A limiting viscosity approach for the Riemann problem in isentropic gas dynamics, Indiana Univ. Math. J. 38 (1989) 1047–1074. [5] W. Sheng, T. Zhang, The Riemann problem for transportation equations in gas dynamics, Mem. Amer. Math. Soc. 137 (564) (1999).

4402

H. Yang, Y. Zhang / J. Differential Equations 257 (2014) 4369–4402

[6] C. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, ICASE Report, 1998, pp. 7–65. [7] H. Yang, H. Cheng, Riemann problem for a geometrical optics system, Acta Math. Sin., Engl. Series (2014), http://dx.doi.org/10.1007/s10114-014-3561-1, to appear. [8] H. Yang, Y. Zhang, New developments of delta shock waves and its applications in systems of conservation laws, J. Differential Equations 252 (2012) 5951–5993.