- Email: [email protected]
Vanishing viscosity approximation to hyperbolic conservation laws
J. Differential Equations 244 (2008) 1692–1711 www.elsevier.com/locate/jde
Vanishing viscosity approximation to hyperbolic conservation laws ✩ Wen Shen ∗ , Zhengfu Xu Department of Mathematics, Penn State University, University Park, PA 16802, USA Received 7 December 2006 Available online 8 February 2008
Abstract We study high order convergence of vanishing viscosity approximation to scalar hyperbolic conservation laws in one space dimension. We prove that, under suitable assumptions, in the region where the solution is smooth, the viscous solution admits an expansion in powers of the viscosity parameter ε. This allows an extrapolation procedure that yields high order approximation to the non-viscous limit as ε → 0. Furthermore, an integral across a shock also admits a power expansion of ε, which allows us to construct high order approximation to the location of the shock. Numerical experiments are presented to justify our theoretical findings. Published by Elsevier Inc. Keywords: Hyperbolic conservation; Vanishing viscosity; High order convergence; Numerical simulation
1. Introduction Consider a system of equations of conservation laws with viscosity in one space dimension ut + f (u)x = εuxx ,
u(0, x) = u(x). ¯
(1.1)
We want to study qualitatively how the solution of this system approximates the one without ✩
This work is supported in part by NSF grant No. 0619587.
* Corresponding author.
E-mail address: [email protected] (W. Shen). 0022-0396/$ – see front matter Published by Elsevier Inc. doi:10.1016/j.jde.2008.01.005
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1693
viscosity ut + f (u)x = 0,
u(0, x) = u(x), ¯
(1.2)
as the viscosity parameter ε tends to zero. It is well known that, even for a smooth initial data, the solution of (1.2) can develop shocks in finite time. This accounts for most of the difficulties in the theoretical and numerical studies of the solutions of (1.2). For a short list of some classical literatures on hyperbolic conservation laws and the recent progresses, see e.g. [2,9,10,18,22,24]. For general solutions with small total variation, various approximating algorithms are known to converge to the unique entropy weak solution. Namely: the Glimm scheme [12], wave-front tracking [2,7,15], and vanishing viscosity approximations [1]. In all cases, the rate of convergence is rather slow. For a fixed time τ > 0, as the step sizes x, t of the grid approach zero at the same rate, the analysis in [3] proved that uGlimm (τ ) − uexact (τ )L1 = 0. √ x→0 x |ln x| lim
A similar bound was established in [4] for viscous approximations. Namely, let A(u) = Df (u) be the Jacobian matrix of the flux function f , and uε be the solution to the parabolic equation uεt + A uε uεx = εuεxx
(1.3)
with same initial data (1.2). Then, for a fixed time τ > 0 one has ε u (τ ) − u(τ )
L1
= O(1) ·
√ ε|ln ε|.
(1.4)
Higher order numerical approximations to (1.2) has a large literature. Many methods have been proposed, see [5,6,19,20,23,26] and reference therein. Most of high order schemes achieve high order convergence rates for smooth solutions. However, in the presence of shocks, the convergence is often degraded to fractional order. According to Glimm’s theorem, solutions to hyperbolic conservation laws are functions with bounded variation. Their set of discontinuities may well be everywhere dense. In this generality, it is unlikely that a high order approximation methods will ever be found. However, in most of the relevant applications, the solutions are much better behaved, with a finite number of shock curves and a finite set of points where new shocks form or two shocks interact. Throughout the paper, we assume the solution is piecewise smooth. In order to achieve high convergence rates also in the presence of shocks, one possibility is to use high order front tracing algorithms, developed by Glimm and his co-authors [13,14]. In a neighborhood of each shock line t → x(t), the system of quasi-linear equations ut + A(u)ux = 0
(1.5)
must be supplemented by the Rankine–Hugoniot conditions. For a shock connecting the left and right states u− , u+ , respectively, these take the form f (u+ ) − f (u− ) = x(u ˙ + − u− ).
(1.6)
1694
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Assuming that the solution is smooth to the right and to the left of the shock, one can use a Taylor expansion of (1.6) and derive a high order numerical algorithm. For a solution containing various shocks, it is essential that all shock curves are traced with high order. In the end, one has to solve a family of N ordinary differential equations, one for each shock curve, coupled with the quasi-linear hyperbolic system (1.5) in the various domains bounded by the shock curves, with appropriate boundary conditions derived from (1.6). Despite its complexities, the feasibility of this approach has been demonstrated in [13]. In the present paper we pursue a different approach. At the present stage, our study is only a theoretical one, and we consider only the scalar equations in one space dimension. The idea is to construct viscous approximate solutions uε of (1.1), with different viscosity coefficients ε1 > ε2 > · · · > 0. A high order approximation to the non-viscous solution u of (1.2) is then sought by means of Richardson’s extrapolation technique. More precisely, fix a time τ > 0. Assuming that the limit solution u(τ, ·) is piecewise smooth, with shocks located at x1 < · · · < xN . On regions where u is smooth, under suitable assumptions, we prove that the viscous solution uε has an expansion of the form uε (x) = u(x) + v1 (x)ε + v2 (x)ε 2 + · · · .
(1.7)
Knowing the values of uε (x) for different values of ε yields a high order approximation to u(x). Furthermore, the location of the j th shocks xj can be recovered by computing a scalar functional of the form bj Ij (u) =
ej · u(τ, x) dx,
(1.8)
aj
where [aj , bj ] is an interval containing one shock xj in its interior. We prove that Ij also satisfies an power expansion of the form I (ε) = I0 + εI1 + ε 2 I2 + · · · + ε k Ik + O ε k+1 .
(1.9)
Combining this with an expansion for the values of u(τ, ·) to the right and to the left of the shock, we can construct a high order approximation of the location of the shock xj . In the literature, similar high order convergence results were obtained in [11], for vanishing viscosity approximations to Hamilton–Jacobi equations, restricted to open sets where the limit solution is smooth. A high order numerical method based on [11] is presented in [25]. Our paper provides a counterpart for those two papers. In Section 2 we present our main results, and give the proof. Numerical experiments supporting our theoretical findings are included in Section 3, which shows that high order accuracy is obtained. In Appendix A we give the proof to a technical lemma, which provides an essential ingredient in the proof of our main theorem. 2. The main result We consider here a scalar conservation law ut + f (u)x = 0,
x ∈ R, t ∈ [0, T ],
(2.1)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1695
and a bounded initial data u(0, x) = u(x). ¯
(2.2)
We assume that the flux function f is smooth and strictly convex, so that fuu (u) κ > 0 for all u ∈ R.
(2.3)
In addition, we consider the viscous Cauchy problem with the same initial data uεt + f uε x = εuεxx ,
uε (0, x) = u(x). ¯
(2.4)
It is well known [17] that under these hypotheses the Cauchy problem (2.1)–(2.2) admits a unique entropy solution u = u(t, x). Moreover, the assumption (2.3) implies the Oleinik inequality [21] u(t, y) − u(t, x)
y −x κτ
∀x < y, t > 0.
(2.5)
In the following, the right and left limits of u at a point (t, x) are written . u+ (t, x) = lim u(t, y), y→x+
. u− (t, x) = lim u(t, y). y→x−
On regions where the limit solution is smooth, the main result on high order convergence of viscous approximations is the following. Theorem 1. Assume that the solution u of (2.1)–(2.2) is uniformly Lipschitz continuous and infinitely many times continuously differentiable restricted to a domain of the form Ω = (t, x); t ∈ [0, T ], a(t) < x < b(t) , where a(·), b(·) are Lipschitz continuous curves whose derivatives satisfy a(t) ˙ f u+ t, a(t) ,
˙ f u− t, b(t) , b(t)
(2.6)
for almost every t ∈ [0, T ]. For any δ > 0, consider the subdomain Ωδ = (t, x); t ∈ [0, T − δ], a(t) + δ x b(t) − δ .
(2.7)
Then, for any k 1, the family of viscous solutions uε admits an expansion at ε = 0: uε (t, x) = u(t, x) + εu1 (t, x) + · · · + ε k uk (t, x) + Rk,ε (t, x),
(2.8)
where u1 , . . . , uk are smooth functions and the remainder satisfies lim ε −k · Rk,ε (t, x) = 0
ε→0+
uniformly as (t, x) ∈ Ωδ .
(2.9)
1696
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Fig. 1. Illustration of α(t) and β(t).
Remark 1. A typical situation where (2.6) holds is when Ω is the region bounded by a maximal backward characteristic a(·) and a minimal backward characteristic b(·), as defined in [8]. Proof of Theorem 1. The result will be proved in several steps. 1. Let δ > 0 be given. We begin by constructing smooth curves t → α(t) and t → β(t) such that α(t), α(t) + η ⊂ a(t), a(t) + δ ,
β(t) − η , β(t) ⊂ b(t) − δ, b(t) ,
(2.10)
and
f u(t, x) α(t) ˙ − 2η, x ∈ α(t), α(t) + η , ˙ + 2η, x ∈ β(t) − η , β(t) , f u(t, x) β(t)
(2.11)
for some constants η, η > 0. To construct the curve α(·) we proceed as follows. Consider the backward Cauchy problem x(t) ˙ = min a(t), ˙ f u t, x(t) − 4η ,
x(T ) = a(T ) + η.
Since the right-hand side of the ODE is Lipschitz continuous, for every small η > 0 we obtain a unique solution x η (·). Clearly, x η (t) > a(t). Moreover, as η → 0, we have x η (t) → a(t) uniformly for all t ∈ [0, T ]. We can thus choose η > 0 small enough so that x η (t) ∈ a(t), a(t) + δ/2 ,
t ∈ [0, T ].
By a mollification, we now replace the Lipschitz function x η (·) with a smooth function t → α(t) > a(t). Using (2.6) and the continuity of u, we can achieve the bound α(t) ˙ f u t, α(t) − 3η,
t ∈ [0, T ].
Still by continuity, we can choose η > 0 such that the first inequality in (2.11) holds. The construction of the curve β(·) is similar. See Fig. 1 for an illustration.
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1697
2. Next, recalling that a(t) < α(t) < β(t) < b(t) for all t ∈ [0, T ] and that u is smooth on the domain Ω, we construct a smooth function u˜ : [0, T − δ/2] × R → R such that u(t, ˜ x) = u(t, x),
x ∈ α(t), β(t) ,
(2.12)
and u(t, ˜ x) = 0,
|x| M,
(2.13)
for all t ∈ [0, T − δ/2] and some constant M large enough. Defining the smooth function . ˜ x) x , φ(t, x) = u˜ t (t, x) + f u(t, it is trivial to check that u˜ provides a solution to the balance law ˜ x = φ. u˜ t + f (u)
(2.14)
Notice that, by construction, we have φ(t, x) = 0 whenever x ∈ [a(t) + η, b(t) − η], and also when |x| M. 3. We now consider the viscous approximations u˜ ε , defined as the solutions to the viscous Cauchy problems with source term u˜ ε + f u˜ ε x = φ + ε u˜ xx ,
u˜ ε (0, x) = u(0, ˜ x).
(2.15)
Since the limit solution u˜ is globally smooth, it is well known that the functions u˜ ε admit an asymptotic expansion in terms of powers of ε: k,ε . u˜ ε = u0 + εu1 + ε 2 u2 + · · · + ε k uk + R
(2.16)
We recall here the basic construction. The functions uj can be inductively determined by requiring that u0 + εu1 + · · · + ε j uj t + f u0 + εu1 + · · · + ε j uj x = φ + ε u0 + εu1 + · · · + ε j uj xx + O ε j +1
(2.17)
with initial data ˜ x), u0 (0, x) = u(0,
uj (0, x) = 0,
j 1.
In particular, we have . ˜ u0 = u, while u1 is found by solving the linear equation ˜ 1 )x + f (u) ˜ u˜ x u1 = u˜ xx . (u1 )t + f (u)(u
(2.18)
1698
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
In general, uj satisfies a linear equation of the form ˜ j )x = A(t, x)uj + B(t, x), (uj )t + f (u)(u
(2.19)
where the coefficients A, B are smooth and depend on the previous functions u0 , . . . , uj −1 . The linear equation (2.19) can of course be solved by the method of characteristics. Recalling (2.11), we conclude that the restriction of all functions uj to the region . Ω α,β = (t, x); t ∈ [0, T − δ/2], x ∈ α(t), β(t)
(2.20)
is independent of the choice of the smooth function u. ˜ k,ε in (2.16), we introduce the functions To estimate the remainder term R . u˜ ε,k = u0 + εu1 + ε 2 u2 + · · · + ε k uk . k,ε satisfies a nonlinear equation of the form, with w as the For each k 1, the remainder R unknown, (u˜ k,ε + w)t + f (u˜ k,ε + w)(u˜ k,ε + w)x = φ + ε(u˜ k,ε + w)xx
(2.21)
with zero initial data. Recalling (2.17), one can write (2.21) in the equivalent form wt + f (u˜ k,ε )wx + Φk,ε (w)w · (u˜ k,ε + w)x = εwxx + ψk,ε ,
(2.22)
where Φk,ε is the smooth function defined by . Φk,ε (w) =
1
f θ u˜ k,ε + (1 − θ )w dθ.
0
Moreover, by the smoothness of the functions u0 , . . . , uk , for every ν 0 we have lim sup ε −k−1 ψk,ε C ν < ∞.
(2.23)
ε→0+
Therefore, for every k, ν 1 there exists a constant Ck,ν such that k,ε C ν Ck,ν ε k+1 R
(2.24)
for every ε > 0 sufficiently small. 4. In the remainder of the proof we compare the solution uε of (2.4) with the solution u˜ ε . Uniformly on the domain Ω α,β defined in (2.20), for any k 1 we claim that |u˜ ε − uε | = 0. ε→0 εk lim
(2.25)
Together, (2.24) and (2.25) will establish the theorem. The estimate (2.25) will be established by a homotopy method. For θ ∈ [0, 1], define w ε,θ as the solution to the following Cauchy problem
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1699
wt + f (w)x = θ φ + εwxx ,
(2.26)
w(0, x) = θ u(0, ˜ x) + (1 − θ )u(0, x).
(2.27)
Observe that w ε,0 = uε while w ε,1 = u˜ ε . To prove the estimate (2.25), it suffices to prove an estimate on the partial derivative ∂ ε,θ w (t, x) Ck ε k+1 ∂θ . uniformly on Ω α,β . Toward a proof of (2.28), we observe that the function z = a solution to the linear, non-homogeneous Cauchy problem
∂ ε,θ ∂θ w
(2.28) provides
zt + f w ε,θ z x = φ + εzxx ,
(2.29)
z(0, x) = u(0, ˜ x) − u(0, x).
(2.30)
Before working out details, we explain here the heart of the matter. According to (2.29)–(2.30), the initial data vanishes in a neighborhood of the interval [α(0), β(0)]. Moreover, the source term φ vanishes on a whole neighborhood of Ω α,β . If ε = 0, Eq. (2.29) reduces to a first order linear equation, and the method of characteristics immediately yields z = 0 on a neighborhood of the domain Ω α,β . On the other hand, for ε > 0, the small diffusion, which possibly takes mass inside the domain, is overwhelmed by convection, which transports mass with a speed that points strictly outward across the boundaries of Ω α,β . The total mass which is found inside the domain is thus o(ε k ), for every k 1. It is convenient here to rescale time and space coordinates. In terms of the new coordinates (t , x ) = (t/ε, x/ε), the linear Cauchy problem (2.19)–(2.30) takes the form zt + f w ε,θ z x = εφ + zx x
(2.31)
z(0, x ) = u(0, ˜ εx ) − u(0, εx ).
(2.32)
with
Notice that in this case the solution is uniformly Lipschitz continuous on the domain . Ωα,β = (t, x); t ∈ [0, T /ε], x ∈ α(t)/ε, β(t)/ε .
(2.33)
Moreover, for all ε > 0 sufficiently small, the drift λ = f (w ε,θ ) ≈ f (u) points strictly outward, on the regions (t , x ); t ∈ [0, T /ε], α(εt ) εx α(εt ) + η , (t , x ); t ∈ [0, T /ε], β(εt ) − η εx β(εt ) . The key fact is that, as ε → 0, the width of these regions grows as η /ε. By linearity, the solution of the Cauchy problem (2.31)–(2.32) is the sum of two solutions: one with initial data and source term supported in the region where x < α(t), and one with initial data and source term supported
1700
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
in the region where x > β(t). We can thus apply Lemma 1 (see Appendix A) to each of these solutions, obtaining the bound
z(t , x ) 2Ce−0 η /ε , x ∈ α(εt ) + η , β(εt ) − η . (2.34) ε ε Returning to the original variables, for x ∈ [α(t) + η , β(t) − η ] using (2.34) we conclude ε u (t, x) − u˜ ε (t, x)
1 ∂ ε,θ w (t, x) dθ 2Ce−0 η /ε . ∂θ 0
Since 0 , η > 0, this implies (2.25), completing the proof of Theorem 1.
2
Next corollary gives an estimate to the remainder term in the expansion (2.8). Corollary 1. In the same setting as Theorem 1, in the expansion (2.8) the remainder term satisfies the estimate lim ε −k · Rk,ε C ν (Ωδ ) = 0
ε→0+
(2.35)
for every k, ν 1. Proof. In view of (2.14), it suffices to show that the difference u˜ ε − uε converges to zero in C ν (Ωδ ) faster than any power of ε. In the rescaled variables (t , x ) = (t/ε, x/ε), the functions uε , u˜ ε satisfy the uniformly parabolic equations uεt + f uε x = uεx x , u˜ εt + f u˜ ε x = εθ φ + u˜ εx x . Moreover, their initial data are smooth restricted to the intervals [a(0)/ε, b(0)/ε]. By standard α,β parabolic estimates, all their derivatives remain uniformly bounded on the rescaled domains Ωε at (2.33), namely m n ε m n ε ∂ ∂ u (t , x ) Cm,n , ∂ ∂ u (t , x ) Cm,n , (t , x ) ∈ Ω α,β . ε t x t x In the original variables t, x, for every ν 0 this yields the estimates ε ε u ν+1 α,β Cν ε −ν−1 , u˜ ν+1 α,β Cν ε −ν−1 C (Ω ) C (Ω )
(2.36)
for a suitable constant Cν . The bounds (2.36) are far from optimal, but suffice for our purposes. Indeed, let integers ν, 1 be given. Using (2.9) with k = ν + + 1 one obtains ε u − u˜ ε 0 α,β = Ck ε ν++2 . (2.37) C (Ω ) By interpolation, (2.36) and (2.37) yield ε u − u˜ ε ν α,β = O ε . C (Ω ) Since ν, are arbitrary, the corollary is proved.
2
(2.38)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1701
Next corollary shows that an integral across a shock also admits a similar expansion. Corollary 2. Fix a time τ > 0 and assume that u is smooth on a neighborhood of two backward ¯ Then the integral characteristics passing through the points (τ, a) ¯ and (τ, b). . I (ε) =
b¯ uε (τ, x) dx a¯
admits an expansion of the form I (ε) = I0 + εI1 + ε 2 I2 + · · · + ε k Ik + O ε k+1 .
(2.39)
Proof. Call t → a(t), t → b(t) the two characteristic curves. Since Eq. (1.3) is in conservation form, we can write b(0) τ ε I (ε) = f u t, a(t) − f uε t, b(t) dt u(0, x) dx + 0
a(0)
τ +
ε ux t, b(t) − ux t, a(t) dt.
(2.40)
0
Even if the limit solution u contains shocks inside the region bounded by the characteristics a(·) and b(·), applying Corollary 1 with ν = 1 we obtain that the last integrands in (2.40) admit an asymptotic expansion in powers of ε. This concludes the proof. 2 The previous results are obtained assuming that all viscous solutions have the same initial data (2.2). One can show that the same conclusions hold provided that the initial data uε (0, ·) converge to u¯ in a suitable way. Corollary 3. Assume that, as ε → 0+, the initial data uε (0, ·) remain uniformly bounded in L∞ and converge to u¯ in L1loc . Moreover, assume that an expansion of the form ¯ + εu1 (0, x) + · · · + ε k uk (0, x) + O ε k+1 uε (0, x) = u(x)
(2.41)
holds as ε → 0+, in neighborhood of the point y. Then the conclusion of Theorem 1 remains valid. If expansions of the form (2.7) hold in a neighborhood of the two points a = a − τf (u(τ, a)) and b = b − τf (u(τ, b)), then the conclusion of Corollary 1 remains valid. The proof is the same as in Theorem 1, except that the initial data for z now is not ≡ 0 but satisfies an estimate z(0, x) C0 ε k . The conclusion does not change.
1702
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
3. Numerical experiments In this section we present some numerical experiments to verify our theoretical findings in Section 2. Assume some quantity u(ε) admits a power expansion in ε such that u(ε) = u¯ + εu1 + ε 2 u2 + ε 3 u3 + · · · . Knowing the values of u for different values of ε allows us to obtain high order approximation to the zero limit u¯ by Richardson extrapolation technique. Say if we computed u1 = u(ε) and u2 = u(ε/2), then 2u2 − u1 is a second order approximation to u, ¯ 2u2 − u1 = u¯ + ε 2 u˜ 2 + · · · . By choosing a sequence of ε, we can generate a triangular table. If the rates of convergence are as expected, it is an indication that the expansion is valid. In our numerical experiments we consider the viscous Burgers’ equation ut +
u2 2
= εuxx
(3.1)
x
with periodic initial condition u(x, 0) = 1 − sin(2πx). For any positive ε, solutions of (3.1) are smooth for all time. Theoretically viscous shocks will form at t = 1. Several experiments are carried out. 3.1. Extrapolation in smooth region; Theorem 1 In this computation we want to verify Theorem 1, i.e., an expansion as in (2.8) is valid in the smooth region that satisfies the assumptions in Theorem 1. Our numerical approximation is computed on the interval x ∈ [0, 1], and we impose corresponding periodic boundary conditions at x = 0 and x = 1. We consider three cases: at t = 0.5 (before shock formation), at t = 1 (right at shock formation) and at t = π2 (after shock formation). To compute the solutions to the viscous Burgers’ equation, we use high order numerical approximation: we use a fifth order finite difference WENO reconstruction in space direction and third order TVD Runge–Kutta in time discretization. We refer to [16] for the detail of high order finite difference WENO schemes. Grid size is chosen very small to give an accurate approximation. We compute solutions to viscous Burgers’ equation with a sequence of ε, chosen in a way such that εi+1 = 12 εi . Provided that the expansion (2.8) is valid, one can design a Richardson extrapolation technique to obtain higher order approximation to the zero limit. Let u¯ ε be the (numerical) solutions to the viscous Burgers’ equation, and u the (exact) solution of the invicid Burgers’ equation. At a point x = xs where the invicid solution is smooth, the errors at different levels of extrapolations are computed as
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1703
Table 1 Errors of extrapolation in smooth region e2
r
e3
r
t = 1/2, before viscous shock is formed 1/20 3.9922e–2 1/40 2.0330e–2 0.9735 1/80 1.0255e–2 0.9873 1/160 5.1498e–3 0.9937 1/320 2.5804e–3 0.9972 1/640 1.2915e–3 0.9985
7.3812e–4 1.8040e–4 4.4387e–5 1.0994e–5 2.7349e–6
2.0395 2.0000 2.0324 2.0065
5.4983e–6 9.5282e–7 1.3672e–7 1.8176e–8
2.5300 2.7999 2.9109
t = 1, the time when viscous shock is formed 1/20 2.6790e–2 1/40 1.3030e–2 1.0398 1/80 6.4292e–3 1.0191 1/160 3.1937e–3 1.0093 1/320 1.5917e–3 1.0046 1/640 7.9459e–4 1.0023
7.2933e–4 1.7187e–4 4.1733e–5 1.0283e–5 2.5525e–6
2.0852 2.0421 2.0208 2.0103
1.3944e–5 1.6477e–6 1.9936e–7 2.4479e–8
3.0811 3.0470 3.0257
t = π/2, after viscous shock is formed 1/20 9.5957e–3 1/40 4.6714e–3 1.0385 1/80 2.3060e–3 1.0184 1/160 1.1458e–3 1.0090 1/320 5.7115e–4 1.0044 1/640 2.8513e–4 1.0022
2.5282e–4 5.9290e–5 1.4380e–5 3.5425e–6 8.7922e–7
2.0922 2.0436 2.0212 2.0105
5.2228e–6 5.8915e–7 7.0171e–8 8.5554e–9
3.1481 3.0696 3.0360
ε
e1
r
e1 = u¯ ε (xs ) − u(xs ), e2 = 2u¯ ε/2 (xs ) − u¯ ε (xs ) − u(xs ),
1 8 e3 = u¯ ε (xs ) − 2u¯ ε/2 (xs ) + u¯ ε/4 (xs ) − u(xs ). 3 3 The results are shown in Table 1. We see clearly that they show a solid support to the expansion (2.8). 3.2. Integral across a shock; Corollary 2 The goal of this experiment is to numerically support our Corollary 2. In our test problem, the shock is formed at t = 1 where x = π2 , with shock speed 0. At t = π2 , after shock formation, we consider the integral 3π
I¯ε =
4
uε (x) dx π 4
1704
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Table 2 Errors and rates of convergence for I , the integral across a shock ε
e1
r
e2
r
e3
r
1/20 1/40 1/80 1/160 1/320 1/640
4.6644e–6 2.2077e–6 1.0774e–6 5.3253e–7 2.6477e–7 1.3201e–7
1.0791 1.0350 1.0166 1.0081 1.0040
2.4893e–7 5.2903e–8 1.2353e–8 2.9886e–9 7.3298e–10
2.2341 2.0982 2.0470 2.0275
1.2439e–8 1.1634e–9 1.3286e–10 1.8903e–11
2.3130 2.5746 2.7848
across the shock for the viscous equation, and compare it to the same integral for the invicid equation 3π
4 I=
u(x) dx. π 4
Again, we compute several I¯ε with the same choices of ε as in our previous test. Assume the expansion (2.39) is valid, we can again design a Richardson extrapolation technique, and the errors at each level are e1 = I¯ε − I, e2 = 2I¯ε/2 − I¯ε − I,
1 8 e3 = I¯ε − 2I¯ε/2 + I¯ε/4 − I. 3 3 The results are shown in Table 2. Again, the rates are as predicted by Corollary 2. 3.3. Recovering shock location with high order accuracy Based on Theorem 1 and Corollary 2, we can design an algorithm to recover the location of a shock with high order accuracy as follows. • We first roughly determine the shock location xs . This can be achieved by two tale-telling signs: (1) a local maximum in the size of the gradient |uεx (τ, ·)|; (2) the blow-up of the gradient as ε → 0, i.e. |uεx (τ, x)| = O(ε −1 ). • We construct an interval around x0 , containing the shock. In our experiment, we choose the interval [a, b] where a = x0 − |4ε log(ε)| and b = x0 + |4ε log(ε)|. • For the integral I across the shock, provided the expansion (2.39), we can extrapolate to get high order accuracy: Iv = 1/3I¯ε − 2I¯ε/2 + 8/3I¯ε/4 . • At the points a, b, we reconstruct two fifth order polynomials Pa (x), Pb (x) from extrapolated high order approximate solution v = 1/3uε − 2uε/2 + 8/3uε/4 .
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1705
Table 3 Errors and rate of convergence of recovered shock location at t = π/2 ε
eu
1/20 1/40 1/80 1/160
1.5214e–3 5.6003e–4 3.7014e–5 4.9121e–7
r
ev
r
1.4215 3.9198 6.2386
1.6173e–2 4.1072e–4 1.0218e–4 5.3349e–7
1.9644 2.0356 7.5516
Table 4 Errors and rate of convergence at t = 1, in the smooth region inside a centered rarefaction fan ε
e1
r
e2
r
e3
r
1/20 1/40 1/80 1/160 1/320 1/640
9.3135e–6 5.4262e–6 2.9793e–6 1.5711e–6 8.0841e–7 4.1028e–7
0.7793 0.8649 0.9232 0.9586 0.9784
1.5388e–6 5.3252e–7 1.6287e–7 4.5712e–8 1.2152e–8
1.5309 1.7090 1.8331 1.9113
1.9708e–7 3.9660e–8 6.6573e–9 9.6601e–10
2.3130 2.5746 2.7848
• We use the polynomials Pa (x), Pb (x) as high order approximation to the invicid Burgers’ equation at the left and the right of the shock, respectively. The approximated shock location x¯s can then be computed by solving the equation x¯s
b Pa (x) dx +
Pb (x) dx = Iv . x¯s
a
The errors and rates of convergence are given in Table 3, where ev = x¯s − xs is the error and xs = π2 is the exact shock location. We see clearly that high order rate of convergence is obtained. 3.4. In the smooth region inside a rarefaction wave This numerical experiment studies if an expansion of the form (2.8) is also valid in the smooth region in the middle of a rarefaction wave. We note that a theoretical proof of such an expansion is not provided in this paper. The difficulty lies in the fact that the backward characteristics for a smooth region in the rarefaction fan will join at a point where the initial discontinuity lies. The technique for proving Theorem 1 cannot be applied in this case. This is a topic of further research for the authors. Here, however, we present a preliminary numerical study which indicates that such an expansion is also valid. We choose our initial data with a discontinuity of an upward jump, u(x, 0) =
−1, x < 0.5, 1, x 0.5.
We consider the point at t = 1 and x = 0, which lies in the middle of the rarefaction fan. The experiment is carried out in a similar way as in our first test in Section 3.1. The errors and rates of convergence are presented in Table 4, which clearly indicates such an expansion (2.8) is valid.
1706
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Appendix A We prove here a technical lemma, which provide the basic ingredients in the proof of Theorem 1. Lemma 1. Consider the linear parabolic Cauchy problem z(0, x) = z¯ (x) zt + λ(t, x)z x = φ(t, x) + zxx ,
(A.1)
with λ uniformly bounded. Assume that there exist constants 0 < η 1, η > 1/2η and a smooth curve t → α(t) with the following properties z¯ (x) = 0, φ(t, x) = 0, λ(t, x) α(t) ˙ − η,
x α(0), x α(t), x ∈ α(t), α(t) + η .
(A.2) (A.3) (A.4)
Then the solution satisfies z(t, x) Ce−0 η ,
x α(t) + η .
(A.5)
Here the constants C, 0 > 0 depend on η, λL∞ , φL∞ , and on the local Lipschitz constant of λ in a neighborhood of (t, x), but not on η . Proof. By possibly replacing x by the shifted space variable y = x − α(t), it is not restrictive to . assume that α(t) = 0. For some constants η, λˆ > 0 we thus have λ(t, x) −η, x ∈ [0, η ], λ(t, x) λ, ˆ x ∈ R.
(A.6) (A.7)
The proof will be given in several steps. 1. Consider the fundamental solution Γ (t, x, τ, y) of the linear homogeneous equation zt + λ(t, x)z x = zxx (A.8) with initial data at time τ consisting of a unit mass concentrated at the point y. The first step in the proof is to show that, for any y 0, the total mass that creeps inside the half line x > η /2 is small. Indeed, Γ = Zx , where Z provides the solution to 0 if x < y, Zt + λ(t, x)Zx = Zxx , Z(τ, x) = (A.9) 1 if x > y. Using (A.6) we can construct a lower solution to (A.9) by setting
Z0 (x) =
⎧ 0 ⎪ ⎪ ⎨ 1 − e−ηx ⎪ P (x) ⎪ ⎩ P (η )
if x 0, if x ∈ [0, η /2], if x ∈ [η /2, η ], if x η .
(A.10)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1707
Here P (·) is the quadratic polynomial obtained as the second order Taylor expansion of the function f (x) = 1 − e−ηx at the point x = η /2. Moreover, . η 1 η = + η 2 η is the point where P attains its maximum. Defining the small constant
η2 e−ηη /2 . P , κ = − (η /2) = P 1 − e−ηη /2
(A.11)
it is straightforward to check that the function . Z∗ (t, x) = e−κ(t−τ ) Z0 (x) is a lower solution to (A.9). In particular, for every t > τ we obtain the upper bound ∞
Γ (t, x, τ, y) dx 1 − Z∗ (t, η /2) = 1 − e−κ(t−τ ) 1 − e−ηη /2
η /2
< κ(t − τ ) + e−ηη /2 Cη e−ηη /2 (1 + t − τ )
(A.12)
for a suitable constant Cη . In the case y 0, a better upper bound is obtained simply observing that, by (A.7), the solution to 0 if x < y, ˆ Wt + λWx = Wxx , W (τ, x) = (A.13) 1 if x > y. provides an upper solution to Z in (A.9). In particular, for t > τ and any y ∈ R one has ∞
∞
η /2
e−x /4(t−τ ) dx √ 2 π(t − τ ) 2
Γ (t, x, τ, y) dx η /2−λˆ (t−τ )−y
∞ = (η /2−y)(t−τ )−1/2 −λˆ (t−τ )1/2
e−x /4 √ dx. 2 π 2
Recall that, for s > 0, ∞ s
e−x /4 √ dx 2 π 2
∞ s
e−sx/4 2 2 √ dx = √ e−s /4 . πs 2 π
ˆ − τ ), the above estimates yield When y −2λ(t
(A.14)
1708
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
∞
∞
2
Γ (t, x, τ, y) dx √ −y/2 (t−τ )
η /2
e−x /4 √ dx 2 π
4 2 2 ˆ ˆ e−y /16(t−τ ) √ eλy/8 C eλy/9 −y √ π ( 2 t−τ ) ˆ −2π λy
(A.15)
for a suitable constant C . 2. The solution of the non-homogeneous problem (A.1) admits the integral representation 0 z(t, x) =
t 0 Γ (t, x, 0, y)¯z(y) dy +
−∞
(A.16)
Γ (t, x, τ, y)φ(τ, y) dy dτ. 0 −∞
Therefore, for every t > 0 using the bounds (A.14) and (A.15) we obtain ∞
z(t, x) dx ¯zL∞ ·
−η −2λˆ t
0 +
η /2
−η −2λˆ t
t
∞
Γ (t, x, 0, y) dx dy −∞
η /2 −η −2 λˆ (t−τ )
0
+ φL∞ ·
+
Γ (t, x, 0, y) dx dy dt −∞
−η −2λˆ (t−τ )
0
ˆ · Cη e ¯zL∞ · (η + 2λt)
−ηη /2
∞ η /2
−η −2λˆ t
(1 + t) +
−λˆ y/9
Ce
dy
−∞
t + φL∞ ·
ˆ − τ ) · Cη e−ηη /2 (1 + t − τ ) η + 2λ(t
0 −η −2 λˆ (t−τ )
+
−λˆ y/9
Ce
dy dτ.
(A.17)
−∞
As long as t (η )2 , from (A.17) we deduce ∞
z(t, x) dx C1 ¯zL∞ + φL∞ · (1 + η )3 e−1 η
η /2
C0 ¯zL∞ + φL∞ e−0 η for suitable constants 1 > 0 > 0 and C1 , C0 large enough.
(A.18)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1709
3. The estimate (A.18) provides an integral bound on the total mass which is carried by diffusion beyond the line x = η /2. In this last step, we transform this integral estimate (A.18) into a pointwise estimate, valid on a region where the drift coefficient λ = λ(t, x) is Lipschitz con. tinuous. Fix a point (t ∗ , x ∗ ) with t ∗ > 0, x ∗ η . Call t0 = max{t ∗ − 1, 0}, and assume that λ is Lipschitz continuous with constant L on the rectangle . ˆ x ∗ + 2λ]. ˆ Q = [t0 , t ∗ ] × [x ∗ − 2λ, . Set λ∗ = λ(t ∗ , x ∗ ) and consider the Green kernel √ ∗ 2 . G∗ (t, x) = e−(x−λ t) /4t /2 πt. Writing (A.1) in the equivalent form zt + λ(t, x) − λ∗ z x + λ∗ zx = φ(t, x) + zxx ,
z(0, x) = z¯ (x),
(A.19)
we can represent its solution as t
∗
z(t) = G (t − t0 ) ∗ z(t0 ) +
∗
t
G (t − s) ∗ φ(s) ds − t0
G∗x (t − s) ∗ λ(s) − λ∗ z(s) ds.
t0
(A.20) Hence η /2 ∞ z(t ∗ , x ∗ ) = + G∗ (t ∗ − t0 , x ∗ − y)z(t0 , y) dy −∞
η /2
t ∗ 0 +
G∗ (t ∗ − s, x ∗ − y)φ(s, y) dy ds
t0 −∞
t η /2 ∞ + G∗x (t ∗ − s, x ∗ − y) λ(s, y) − λ∗ z(s, y) dy ds. −
t0
−∞
(A.21)
η /2
Because of (A.18), it is clear that all of the above integrals become exponentially small as η → ∞. Indeed, the portion containing the singularity can also be estimated, observing that 1 t · G∗x (t, x) √ , 8πe This yields
C |x| · G∗x (t, x) √ . t
1710
t ∗
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711 x ∗ +λ∗(s−t ∗ )+1
G∗x (t ∗ − s, x ∗ − y) λ(s, y) − λ∗ z(s, y) dy ds
t0 x ∗ +λ∗ (s−t ∗ )−1
t ∗
max
t0
t ∗ t0
|x−x ∗ −λ∗ (t−t ∗ )|1
∗ ∗ G (t − s, x ∗ − y) λ(s, y) − λ∗ · x
x ∗ +λ∗(s−t ∗ )+1
z(s, y) dy dt
x ∗ +λ∗ (s−t ∗ )−1
∞ 1 C z(s, y) dy ds. · L· √ +√ ∗ (t − s) 8πe η /2
This establishes (A.5), with a possibly smaller constant 0 > 0.
2
References [1] S. Bianchini, A. Bressan, Vanishing viscosity solutions to nonlinear hyperbolic systems, Ann. of Math. 161 (2005) 223–342. [2] A. Bressan, Hyperbolic Systems of Conservation Laws. The One-Dimensional Cauchy Problem, Oxford Univ. Press, 2000. [3] A. Bressan, A. Marson, Error bounds for a deterministic version of the Glimm scheme, Arch. Ration. Mech. Anal. 142 (1998) 155–176. [4] A. Bressan, T. Yang, On the convergence rate of vanishing viscosity approximations, Comm. Pure Appl. Math. 57 (2004) 1075–1109. [5] B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Lect. Notes Comput. Sci. Eng., vol. 11, Springer-Verlag, 2000. [6] B. Cockburn, C. Johnson, C.W. Shu, E. Tadmor, in: A. Quarteroni (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, in: Lecture Notes in Math., vol. 1697, Springer-Verlag, 1998. [7] C. Dafermos, Polygonal approximations of solutions of the initial value problem for a conservation law, J. Math. Anal. Appl. 38 (1972) 33–41. [8] C. Dafermos, Generalized characteristics and the structure of solutions of hyperbolic conservation laws, Indiana Univ. Math. J. 26 (6) (1977) 1097–1119. [9] C. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, Springer-Verlag, Berlin, 1999. [10] R. DiPerna, Global existence of solutions to nonlinear hyperbolic systems of conservation laws, J. Differential Equations 20 (1976) 187–212. [11] W.H. Fleming, P.E. Souganidis, Asymptotic series and the method of vanishing viscosity, Indiana Univ. Math. J. 35 (1986) 425–447. [12] J. Glimm, Solutions in the large for nonlinear hyperbolic systems of equations, Comm. Pure Appl. Math. 18 (1965) 697–715. [13] J. Glimm, X. Li, Y. Liu, Conservative front tracking in one space dimension, in: Fluid Flow and Transport in Porous Media: Mathematical and Numerical Treatment, South Hadley, MA, 2001, in: Contemp. Math., vol. 295, Amer. Math. Soc., Providence, RI, 2002, pp. 253–264. [14] J. Glimm, X. Li, Y. Liu, Z. Xu, N. Zhao, Conservative front tracking with improved accuracy, SIAM J. Numer. Anal. 41 (2003) 1926–1947. [15] H. Holden, N.H. Risebro, Front Tracking for Hyperbolic Conservation Laws, Springer-Verlag, New York, 2002. [16] G. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [17] S. Kruzhkov, First order quasilinear equations with several space variables, Math. USSR Sb. 10 (1970) 217–243. [18] P. Lax, Hyperbolic systems of conservation laws II, Comm. Pure Appl. Math. 10 (1957) 537–566. [19] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, 1990. [20] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Univ. Press, 2002. [21] O. Oleinik, Discontinuous solutions of nonlinear differential equations, Amer. Math. Soc. Transl. Ser. 2 26 (1957) 95–172. [22] D. Serre, Systems of Conservation Laws I, II, Cambridge Univ. Press, 2000.
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1711
[23] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in: A. Quarteroni (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, in: Lecture Notes in Math., vol. 1697, Springer-Verlag, 1998, pp. 325–432. [24] J. Smoller, Shock Waves and Reaction–Diffusion Equations, Springer-Verlag, New York, 1983. [25] A. Szpiro, P. Dupuis, Second order numerical methods for first order Hamilton–Jacobi equations, SIAM J. Numer. Anal. 40 (3) (2002) 1136–1183. [26] E. Tadmor, Approximate solutions of nonlinear conservation laws, in: A. Quarteroni (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, 1997 CIME Course, Cetraro, Italy, June 1997, in: Lecture Motes in Math., vol. 1697, Springer-Verlag, 1998, pp. 1–149.
Vanishing viscosity approximation to hyperbolic conservation laws ✩ Wen Shen ∗ , Zhengfu Xu Department of Mathematics, Penn State University, University Park, PA 16802, USA Received 7 December 2006 Available online 8 February 2008
Abstract We study high order convergence of vanishing viscosity approximation to scalar hyperbolic conservation laws in one space dimension. We prove that, under suitable assumptions, in the region where the solution is smooth, the viscous solution admits an expansion in powers of the viscosity parameter ε. This allows an extrapolation procedure that yields high order approximation to the non-viscous limit as ε → 0. Furthermore, an integral across a shock also admits a power expansion of ε, which allows us to construct high order approximation to the location of the shock. Numerical experiments are presented to justify our theoretical findings. Published by Elsevier Inc. Keywords: Hyperbolic conservation; Vanishing viscosity; High order convergence; Numerical simulation
1. Introduction Consider a system of equations of conservation laws with viscosity in one space dimension ut + f (u)x = εuxx ,
u(0, x) = u(x). ¯
(1.1)
We want to study qualitatively how the solution of this system approximates the one without ✩
This work is supported in part by NSF grant No. 0619587.
* Corresponding author.
E-mail address: [email protected] (W. Shen). 0022-0396/$ – see front matter Published by Elsevier Inc. doi:10.1016/j.jde.2008.01.005
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1693
viscosity ut + f (u)x = 0,
u(0, x) = u(x), ¯
(1.2)
as the viscosity parameter ε tends to zero. It is well known that, even for a smooth initial data, the solution of (1.2) can develop shocks in finite time. This accounts for most of the difficulties in the theoretical and numerical studies of the solutions of (1.2). For a short list of some classical literatures on hyperbolic conservation laws and the recent progresses, see e.g. [2,9,10,18,22,24]. For general solutions with small total variation, various approximating algorithms are known to converge to the unique entropy weak solution. Namely: the Glimm scheme [12], wave-front tracking [2,7,15], and vanishing viscosity approximations [1]. In all cases, the rate of convergence is rather slow. For a fixed time τ > 0, as the step sizes x, t of the grid approach zero at the same rate, the analysis in [3] proved that uGlimm (τ ) − uexact (τ )L1 = 0. √ x→0 x |ln x| lim
A similar bound was established in [4] for viscous approximations. Namely, let A(u) = Df (u) be the Jacobian matrix of the flux function f , and uε be the solution to the parabolic equation uεt + A uε uεx = εuεxx
(1.3)
with same initial data (1.2). Then, for a fixed time τ > 0 one has ε u (τ ) − u(τ )
L1
= O(1) ·
√ ε|ln ε|.
(1.4)
Higher order numerical approximations to (1.2) has a large literature. Many methods have been proposed, see [5,6,19,20,23,26] and reference therein. Most of high order schemes achieve high order convergence rates for smooth solutions. However, in the presence of shocks, the convergence is often degraded to fractional order. According to Glimm’s theorem, solutions to hyperbolic conservation laws are functions with bounded variation. Their set of discontinuities may well be everywhere dense. In this generality, it is unlikely that a high order approximation methods will ever be found. However, in most of the relevant applications, the solutions are much better behaved, with a finite number of shock curves and a finite set of points where new shocks form or two shocks interact. Throughout the paper, we assume the solution is piecewise smooth. In order to achieve high convergence rates also in the presence of shocks, one possibility is to use high order front tracing algorithms, developed by Glimm and his co-authors [13,14]. In a neighborhood of each shock line t → x(t), the system of quasi-linear equations ut + A(u)ux = 0
(1.5)
must be supplemented by the Rankine–Hugoniot conditions. For a shock connecting the left and right states u− , u+ , respectively, these take the form f (u+ ) − f (u− ) = x(u ˙ + − u− ).
(1.6)
1694
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Assuming that the solution is smooth to the right and to the left of the shock, one can use a Taylor expansion of (1.6) and derive a high order numerical algorithm. For a solution containing various shocks, it is essential that all shock curves are traced with high order. In the end, one has to solve a family of N ordinary differential equations, one for each shock curve, coupled with the quasi-linear hyperbolic system (1.5) in the various domains bounded by the shock curves, with appropriate boundary conditions derived from (1.6). Despite its complexities, the feasibility of this approach has been demonstrated in [13]. In the present paper we pursue a different approach. At the present stage, our study is only a theoretical one, and we consider only the scalar equations in one space dimension. The idea is to construct viscous approximate solutions uε of (1.1), with different viscosity coefficients ε1 > ε2 > · · · > 0. A high order approximation to the non-viscous solution u of (1.2) is then sought by means of Richardson’s extrapolation technique. More precisely, fix a time τ > 0. Assuming that the limit solution u(τ, ·) is piecewise smooth, with shocks located at x1 < · · · < xN . On regions where u is smooth, under suitable assumptions, we prove that the viscous solution uε has an expansion of the form uε (x) = u(x) + v1 (x)ε + v2 (x)ε 2 + · · · .
(1.7)
Knowing the values of uε (x) for different values of ε yields a high order approximation to u(x). Furthermore, the location of the j th shocks xj can be recovered by computing a scalar functional of the form bj Ij (u) =
ej · u(τ, x) dx,
(1.8)
aj
where [aj , bj ] is an interval containing one shock xj in its interior. We prove that Ij also satisfies an power expansion of the form I (ε) = I0 + εI1 + ε 2 I2 + · · · + ε k Ik + O ε k+1 .
(1.9)
Combining this with an expansion for the values of u(τ, ·) to the right and to the left of the shock, we can construct a high order approximation of the location of the shock xj . In the literature, similar high order convergence results were obtained in [11], for vanishing viscosity approximations to Hamilton–Jacobi equations, restricted to open sets where the limit solution is smooth. A high order numerical method based on [11] is presented in [25]. Our paper provides a counterpart for those two papers. In Section 2 we present our main results, and give the proof. Numerical experiments supporting our theoretical findings are included in Section 3, which shows that high order accuracy is obtained. In Appendix A we give the proof to a technical lemma, which provides an essential ingredient in the proof of our main theorem. 2. The main result We consider here a scalar conservation law ut + f (u)x = 0,
x ∈ R, t ∈ [0, T ],
(2.1)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1695
and a bounded initial data u(0, x) = u(x). ¯
(2.2)
We assume that the flux function f is smooth and strictly convex, so that fuu (u) κ > 0 for all u ∈ R.
(2.3)
In addition, we consider the viscous Cauchy problem with the same initial data uεt + f uε x = εuεxx ,
uε (0, x) = u(x). ¯
(2.4)
It is well known [17] that under these hypotheses the Cauchy problem (2.1)–(2.2) admits a unique entropy solution u = u(t, x). Moreover, the assumption (2.3) implies the Oleinik inequality [21] u(t, y) − u(t, x)
y −x κτ
∀x < y, t > 0.
(2.5)
In the following, the right and left limits of u at a point (t, x) are written . u+ (t, x) = lim u(t, y), y→x+
. u− (t, x) = lim u(t, y). y→x−
On regions where the limit solution is smooth, the main result on high order convergence of viscous approximations is the following. Theorem 1. Assume that the solution u of (2.1)–(2.2) is uniformly Lipschitz continuous and infinitely many times continuously differentiable restricted to a domain of the form Ω = (t, x); t ∈ [0, T ], a(t) < x < b(t) , where a(·), b(·) are Lipschitz continuous curves whose derivatives satisfy a(t) ˙ f u+ t, a(t) ,
˙ f u− t, b(t) , b(t)
(2.6)
for almost every t ∈ [0, T ]. For any δ > 0, consider the subdomain Ωδ = (t, x); t ∈ [0, T − δ], a(t) + δ x b(t) − δ .
(2.7)
Then, for any k 1, the family of viscous solutions uε admits an expansion at ε = 0: uε (t, x) = u(t, x) + εu1 (t, x) + · · · + ε k uk (t, x) + Rk,ε (t, x),
(2.8)
where u1 , . . . , uk are smooth functions and the remainder satisfies lim ε −k · Rk,ε (t, x) = 0
ε→0+
uniformly as (t, x) ∈ Ωδ .
(2.9)
1696
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Fig. 1. Illustration of α(t) and β(t).
Remark 1. A typical situation where (2.6) holds is when Ω is the region bounded by a maximal backward characteristic a(·) and a minimal backward characteristic b(·), as defined in [8]. Proof of Theorem 1. The result will be proved in several steps. 1. Let δ > 0 be given. We begin by constructing smooth curves t → α(t) and t → β(t) such that α(t), α(t) + η ⊂ a(t), a(t) + δ ,
β(t) − η , β(t) ⊂ b(t) − δ, b(t) ,
(2.10)
and
f u(t, x) α(t) ˙ − 2η, x ∈ α(t), α(t) + η , ˙ + 2η, x ∈ β(t) − η , β(t) , f u(t, x) β(t)
(2.11)
for some constants η, η > 0. To construct the curve α(·) we proceed as follows. Consider the backward Cauchy problem x(t) ˙ = min a(t), ˙ f u t, x(t) − 4η ,
x(T ) = a(T ) + η.
Since the right-hand side of the ODE is Lipschitz continuous, for every small η > 0 we obtain a unique solution x η (·). Clearly, x η (t) > a(t). Moreover, as η → 0, we have x η (t) → a(t) uniformly for all t ∈ [0, T ]. We can thus choose η > 0 small enough so that x η (t) ∈ a(t), a(t) + δ/2 ,
t ∈ [0, T ].
By a mollification, we now replace the Lipschitz function x η (·) with a smooth function t → α(t) > a(t). Using (2.6) and the continuity of u, we can achieve the bound α(t) ˙ f u t, α(t) − 3η,
t ∈ [0, T ].
Still by continuity, we can choose η > 0 such that the first inequality in (2.11) holds. The construction of the curve β(·) is similar. See Fig. 1 for an illustration.
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1697
2. Next, recalling that a(t) < α(t) < β(t) < b(t) for all t ∈ [0, T ] and that u is smooth on the domain Ω, we construct a smooth function u˜ : [0, T − δ/2] × R → R such that u(t, ˜ x) = u(t, x),
x ∈ α(t), β(t) ,
(2.12)
and u(t, ˜ x) = 0,
|x| M,
(2.13)
for all t ∈ [0, T − δ/2] and some constant M large enough. Defining the smooth function . ˜ x) x , φ(t, x) = u˜ t (t, x) + f u(t, it is trivial to check that u˜ provides a solution to the balance law ˜ x = φ. u˜ t + f (u)
(2.14)
Notice that, by construction, we have φ(t, x) = 0 whenever x ∈ [a(t) + η, b(t) − η], and also when |x| M. 3. We now consider the viscous approximations u˜ ε , defined as the solutions to the viscous Cauchy problems with source term u˜ ε + f u˜ ε x = φ + ε u˜ xx ,
u˜ ε (0, x) = u(0, ˜ x).
(2.15)
Since the limit solution u˜ is globally smooth, it is well known that the functions u˜ ε admit an asymptotic expansion in terms of powers of ε: k,ε . u˜ ε = u0 + εu1 + ε 2 u2 + · · · + ε k uk + R
(2.16)
We recall here the basic construction. The functions uj can be inductively determined by requiring that u0 + εu1 + · · · + ε j uj t + f u0 + εu1 + · · · + ε j uj x = φ + ε u0 + εu1 + · · · + ε j uj xx + O ε j +1
(2.17)
with initial data ˜ x), u0 (0, x) = u(0,
uj (0, x) = 0,
j 1.
In particular, we have . ˜ u0 = u, while u1 is found by solving the linear equation ˜ 1 )x + f (u) ˜ u˜ x u1 = u˜ xx . (u1 )t + f (u)(u
(2.18)
1698
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
In general, uj satisfies a linear equation of the form ˜ j )x = A(t, x)uj + B(t, x), (uj )t + f (u)(u
(2.19)
where the coefficients A, B are smooth and depend on the previous functions u0 , . . . , uj −1 . The linear equation (2.19) can of course be solved by the method of characteristics. Recalling (2.11), we conclude that the restriction of all functions uj to the region . Ω α,β = (t, x); t ∈ [0, T − δ/2], x ∈ α(t), β(t)
(2.20)
is independent of the choice of the smooth function u. ˜ k,ε in (2.16), we introduce the functions To estimate the remainder term R . u˜ ε,k = u0 + εu1 + ε 2 u2 + · · · + ε k uk . k,ε satisfies a nonlinear equation of the form, with w as the For each k 1, the remainder R unknown, (u˜ k,ε + w)t + f (u˜ k,ε + w)(u˜ k,ε + w)x = φ + ε(u˜ k,ε + w)xx
(2.21)
with zero initial data. Recalling (2.17), one can write (2.21) in the equivalent form wt + f (u˜ k,ε )wx + Φk,ε (w)w · (u˜ k,ε + w)x = εwxx + ψk,ε ,
(2.22)
where Φk,ε is the smooth function defined by . Φk,ε (w) =
1
f θ u˜ k,ε + (1 − θ )w dθ.
0
Moreover, by the smoothness of the functions u0 , . . . , uk , for every ν 0 we have lim sup ε −k−1 ψk,ε C ν < ∞.
(2.23)
ε→0+
Therefore, for every k, ν 1 there exists a constant Ck,ν such that k,ε C ν Ck,ν ε k+1 R
(2.24)
for every ε > 0 sufficiently small. 4. In the remainder of the proof we compare the solution uε of (2.4) with the solution u˜ ε . Uniformly on the domain Ω α,β defined in (2.20), for any k 1 we claim that |u˜ ε − uε | = 0. ε→0 εk lim
(2.25)
Together, (2.24) and (2.25) will establish the theorem. The estimate (2.25) will be established by a homotopy method. For θ ∈ [0, 1], define w ε,θ as the solution to the following Cauchy problem
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1699
wt + f (w)x = θ φ + εwxx ,
(2.26)
w(0, x) = θ u(0, ˜ x) + (1 − θ )u(0, x).
(2.27)
Observe that w ε,0 = uε while w ε,1 = u˜ ε . To prove the estimate (2.25), it suffices to prove an estimate on the partial derivative ∂ ε,θ w (t, x) Ck ε k+1 ∂θ . uniformly on Ω α,β . Toward a proof of (2.28), we observe that the function z = a solution to the linear, non-homogeneous Cauchy problem
∂ ε,θ ∂θ w
(2.28) provides
zt + f w ε,θ z x = φ + εzxx ,
(2.29)
z(0, x) = u(0, ˜ x) − u(0, x).
(2.30)
Before working out details, we explain here the heart of the matter. According to (2.29)–(2.30), the initial data vanishes in a neighborhood of the interval [α(0), β(0)]. Moreover, the source term φ vanishes on a whole neighborhood of Ω α,β . If ε = 0, Eq. (2.29) reduces to a first order linear equation, and the method of characteristics immediately yields z = 0 on a neighborhood of the domain Ω α,β . On the other hand, for ε > 0, the small diffusion, which possibly takes mass inside the domain, is overwhelmed by convection, which transports mass with a speed that points strictly outward across the boundaries of Ω α,β . The total mass which is found inside the domain is thus o(ε k ), for every k 1. It is convenient here to rescale time and space coordinates. In terms of the new coordinates (t , x ) = (t/ε, x/ε), the linear Cauchy problem (2.19)–(2.30) takes the form zt + f w ε,θ z x = εφ + zx x
(2.31)
z(0, x ) = u(0, ˜ εx ) − u(0, εx ).
(2.32)
with
Notice that in this case the solution is uniformly Lipschitz continuous on the domain . Ωα,β = (t, x); t ∈ [0, T /ε], x ∈ α(t)/ε, β(t)/ε .
(2.33)
Moreover, for all ε > 0 sufficiently small, the drift λ = f (w ε,θ ) ≈ f (u) points strictly outward, on the regions (t , x ); t ∈ [0, T /ε], α(εt ) εx α(εt ) + η , (t , x ); t ∈ [0, T /ε], β(εt ) − η εx β(εt ) . The key fact is that, as ε → 0, the width of these regions grows as η /ε. By linearity, the solution of the Cauchy problem (2.31)–(2.32) is the sum of two solutions: one with initial data and source term supported in the region where x < α(t), and one with initial data and source term supported
1700
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
in the region where x > β(t). We can thus apply Lemma 1 (see Appendix A) to each of these solutions, obtaining the bound
z(t , x ) 2Ce−0 η /ε , x ∈ α(εt ) + η , β(εt ) − η . (2.34) ε ε Returning to the original variables, for x ∈ [α(t) + η , β(t) − η ] using (2.34) we conclude ε u (t, x) − u˜ ε (t, x)
1 ∂ ε,θ w (t, x) dθ 2Ce−0 η /ε . ∂θ 0
Since 0 , η > 0, this implies (2.25), completing the proof of Theorem 1.
2
Next corollary gives an estimate to the remainder term in the expansion (2.8). Corollary 1. In the same setting as Theorem 1, in the expansion (2.8) the remainder term satisfies the estimate lim ε −k · Rk,ε C ν (Ωδ ) = 0
ε→0+
(2.35)
for every k, ν 1. Proof. In view of (2.14), it suffices to show that the difference u˜ ε − uε converges to zero in C ν (Ωδ ) faster than any power of ε. In the rescaled variables (t , x ) = (t/ε, x/ε), the functions uε , u˜ ε satisfy the uniformly parabolic equations uεt + f uε x = uεx x , u˜ εt + f u˜ ε x = εθ φ + u˜ εx x . Moreover, their initial data are smooth restricted to the intervals [a(0)/ε, b(0)/ε]. By standard α,β parabolic estimates, all their derivatives remain uniformly bounded on the rescaled domains Ωε at (2.33), namely m n ε m n ε ∂ ∂ u (t , x ) Cm,n , ∂ ∂ u (t , x ) Cm,n , (t , x ) ∈ Ω α,β . ε t x t x In the original variables t, x, for every ν 0 this yields the estimates ε ε u ν+1 α,β Cν ε −ν−1 , u˜ ν+1 α,β Cν ε −ν−1 C (Ω ) C (Ω )
(2.36)
for a suitable constant Cν . The bounds (2.36) are far from optimal, but suffice for our purposes. Indeed, let integers ν, 1 be given. Using (2.9) with k = ν + + 1 one obtains ε u − u˜ ε 0 α,β = Ck ε ν++2 . (2.37) C (Ω ) By interpolation, (2.36) and (2.37) yield ε u − u˜ ε ν α,β = O ε . C (Ω ) Since ν, are arbitrary, the corollary is proved.
2
(2.38)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1701
Next corollary shows that an integral across a shock also admits a similar expansion. Corollary 2. Fix a time τ > 0 and assume that u is smooth on a neighborhood of two backward ¯ Then the integral characteristics passing through the points (τ, a) ¯ and (τ, b). . I (ε) =
b¯ uε (τ, x) dx a¯
admits an expansion of the form I (ε) = I0 + εI1 + ε 2 I2 + · · · + ε k Ik + O ε k+1 .
(2.39)
Proof. Call t → a(t), t → b(t) the two characteristic curves. Since Eq. (1.3) is in conservation form, we can write b(0) τ ε I (ε) = f u t, a(t) − f uε t, b(t) dt u(0, x) dx + 0
a(0)
τ +
ε ux t, b(t) − ux t, a(t) dt.
(2.40)
0
Even if the limit solution u contains shocks inside the region bounded by the characteristics a(·) and b(·), applying Corollary 1 with ν = 1 we obtain that the last integrands in (2.40) admit an asymptotic expansion in powers of ε. This concludes the proof. 2 The previous results are obtained assuming that all viscous solutions have the same initial data (2.2). One can show that the same conclusions hold provided that the initial data uε (0, ·) converge to u¯ in a suitable way. Corollary 3. Assume that, as ε → 0+, the initial data uε (0, ·) remain uniformly bounded in L∞ and converge to u¯ in L1loc . Moreover, assume that an expansion of the form ¯ + εu1 (0, x) + · · · + ε k uk (0, x) + O ε k+1 uε (0, x) = u(x)
(2.41)
holds as ε → 0+, in neighborhood of the point y. Then the conclusion of Theorem 1 remains valid. If expansions of the form (2.7) hold in a neighborhood of the two points a = a − τf (u(τ, a)) and b = b − τf (u(τ, b)), then the conclusion of Corollary 1 remains valid. The proof is the same as in Theorem 1, except that the initial data for z now is not ≡ 0 but satisfies an estimate z(0, x) C0 ε k . The conclusion does not change.
1702
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
3. Numerical experiments In this section we present some numerical experiments to verify our theoretical findings in Section 2. Assume some quantity u(ε) admits a power expansion in ε such that u(ε) = u¯ + εu1 + ε 2 u2 + ε 3 u3 + · · · . Knowing the values of u for different values of ε allows us to obtain high order approximation to the zero limit u¯ by Richardson extrapolation technique. Say if we computed u1 = u(ε) and u2 = u(ε/2), then 2u2 − u1 is a second order approximation to u, ¯ 2u2 − u1 = u¯ + ε 2 u˜ 2 + · · · . By choosing a sequence of ε, we can generate a triangular table. If the rates of convergence are as expected, it is an indication that the expansion is valid. In our numerical experiments we consider the viscous Burgers’ equation ut +
u2 2
= εuxx
(3.1)
x
with periodic initial condition u(x, 0) = 1 − sin(2πx). For any positive ε, solutions of (3.1) are smooth for all time. Theoretically viscous shocks will form at t = 1. Several experiments are carried out. 3.1. Extrapolation in smooth region; Theorem 1 In this computation we want to verify Theorem 1, i.e., an expansion as in (2.8) is valid in the smooth region that satisfies the assumptions in Theorem 1. Our numerical approximation is computed on the interval x ∈ [0, 1], and we impose corresponding periodic boundary conditions at x = 0 and x = 1. We consider three cases: at t = 0.5 (before shock formation), at t = 1 (right at shock formation) and at t = π2 (after shock formation). To compute the solutions to the viscous Burgers’ equation, we use high order numerical approximation: we use a fifth order finite difference WENO reconstruction in space direction and third order TVD Runge–Kutta in time discretization. We refer to [16] for the detail of high order finite difference WENO schemes. Grid size is chosen very small to give an accurate approximation. We compute solutions to viscous Burgers’ equation with a sequence of ε, chosen in a way such that εi+1 = 12 εi . Provided that the expansion (2.8) is valid, one can design a Richardson extrapolation technique to obtain higher order approximation to the zero limit. Let u¯ ε be the (numerical) solutions to the viscous Burgers’ equation, and u the (exact) solution of the invicid Burgers’ equation. At a point x = xs where the invicid solution is smooth, the errors at different levels of extrapolations are computed as
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1703
Table 1 Errors of extrapolation in smooth region e2
r
e3
r
t = 1/2, before viscous shock is formed 1/20 3.9922e–2 1/40 2.0330e–2 0.9735 1/80 1.0255e–2 0.9873 1/160 5.1498e–3 0.9937 1/320 2.5804e–3 0.9972 1/640 1.2915e–3 0.9985
7.3812e–4 1.8040e–4 4.4387e–5 1.0994e–5 2.7349e–6
2.0395 2.0000 2.0324 2.0065
5.4983e–6 9.5282e–7 1.3672e–7 1.8176e–8
2.5300 2.7999 2.9109
t = 1, the time when viscous shock is formed 1/20 2.6790e–2 1/40 1.3030e–2 1.0398 1/80 6.4292e–3 1.0191 1/160 3.1937e–3 1.0093 1/320 1.5917e–3 1.0046 1/640 7.9459e–4 1.0023
7.2933e–4 1.7187e–4 4.1733e–5 1.0283e–5 2.5525e–6
2.0852 2.0421 2.0208 2.0103
1.3944e–5 1.6477e–6 1.9936e–7 2.4479e–8
3.0811 3.0470 3.0257
t = π/2, after viscous shock is formed 1/20 9.5957e–3 1/40 4.6714e–3 1.0385 1/80 2.3060e–3 1.0184 1/160 1.1458e–3 1.0090 1/320 5.7115e–4 1.0044 1/640 2.8513e–4 1.0022
2.5282e–4 5.9290e–5 1.4380e–5 3.5425e–6 8.7922e–7
2.0922 2.0436 2.0212 2.0105
5.2228e–6 5.8915e–7 7.0171e–8 8.5554e–9
3.1481 3.0696 3.0360
ε
e1
r
e1 = u¯ ε (xs ) − u(xs ), e2 = 2u¯ ε/2 (xs ) − u¯ ε (xs ) − u(xs ),
1 8 e3 = u¯ ε (xs ) − 2u¯ ε/2 (xs ) + u¯ ε/4 (xs ) − u(xs ). 3 3 The results are shown in Table 1. We see clearly that they show a solid support to the expansion (2.8). 3.2. Integral across a shock; Corollary 2 The goal of this experiment is to numerically support our Corollary 2. In our test problem, the shock is formed at t = 1 where x = π2 , with shock speed 0. At t = π2 , after shock formation, we consider the integral 3π
I¯ε =
4
uε (x) dx π 4
1704
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Table 2 Errors and rates of convergence for I , the integral across a shock ε
e1
r
e2
r
e3
r
1/20 1/40 1/80 1/160 1/320 1/640
4.6644e–6 2.2077e–6 1.0774e–6 5.3253e–7 2.6477e–7 1.3201e–7
1.0791 1.0350 1.0166 1.0081 1.0040
2.4893e–7 5.2903e–8 1.2353e–8 2.9886e–9 7.3298e–10
2.2341 2.0982 2.0470 2.0275
1.2439e–8 1.1634e–9 1.3286e–10 1.8903e–11
2.3130 2.5746 2.7848
across the shock for the viscous equation, and compare it to the same integral for the invicid equation 3π
4 I=
u(x) dx. π 4
Again, we compute several I¯ε with the same choices of ε as in our previous test. Assume the expansion (2.39) is valid, we can again design a Richardson extrapolation technique, and the errors at each level are e1 = I¯ε − I, e2 = 2I¯ε/2 − I¯ε − I,
1 8 e3 = I¯ε − 2I¯ε/2 + I¯ε/4 − I. 3 3 The results are shown in Table 2. Again, the rates are as predicted by Corollary 2. 3.3. Recovering shock location with high order accuracy Based on Theorem 1 and Corollary 2, we can design an algorithm to recover the location of a shock with high order accuracy as follows. • We first roughly determine the shock location xs . This can be achieved by two tale-telling signs: (1) a local maximum in the size of the gradient |uεx (τ, ·)|; (2) the blow-up of the gradient as ε → 0, i.e. |uεx (τ, x)| = O(ε −1 ). • We construct an interval around x0 , containing the shock. In our experiment, we choose the interval [a, b] where a = x0 − |4ε log(ε)| and b = x0 + |4ε log(ε)|. • For the integral I across the shock, provided the expansion (2.39), we can extrapolate to get high order accuracy: Iv = 1/3I¯ε − 2I¯ε/2 + 8/3I¯ε/4 . • At the points a, b, we reconstruct two fifth order polynomials Pa (x), Pb (x) from extrapolated high order approximate solution v = 1/3uε − 2uε/2 + 8/3uε/4 .
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1705
Table 3 Errors and rate of convergence of recovered shock location at t = π/2 ε
eu
1/20 1/40 1/80 1/160
1.5214e–3 5.6003e–4 3.7014e–5 4.9121e–7
r
ev
r
1.4215 3.9198 6.2386
1.6173e–2 4.1072e–4 1.0218e–4 5.3349e–7
1.9644 2.0356 7.5516
Table 4 Errors and rate of convergence at t = 1, in the smooth region inside a centered rarefaction fan ε
e1
r
e2
r
e3
r
1/20 1/40 1/80 1/160 1/320 1/640
9.3135e–6 5.4262e–6 2.9793e–6 1.5711e–6 8.0841e–7 4.1028e–7
0.7793 0.8649 0.9232 0.9586 0.9784
1.5388e–6 5.3252e–7 1.6287e–7 4.5712e–8 1.2152e–8
1.5309 1.7090 1.8331 1.9113
1.9708e–7 3.9660e–8 6.6573e–9 9.6601e–10
2.3130 2.5746 2.7848
• We use the polynomials Pa (x), Pb (x) as high order approximation to the invicid Burgers’ equation at the left and the right of the shock, respectively. The approximated shock location x¯s can then be computed by solving the equation x¯s
b Pa (x) dx +
Pb (x) dx = Iv . x¯s
a
The errors and rates of convergence are given in Table 3, where ev = x¯s − xs is the error and xs = π2 is the exact shock location. We see clearly that high order rate of convergence is obtained. 3.4. In the smooth region inside a rarefaction wave This numerical experiment studies if an expansion of the form (2.8) is also valid in the smooth region in the middle of a rarefaction wave. We note that a theoretical proof of such an expansion is not provided in this paper. The difficulty lies in the fact that the backward characteristics for a smooth region in the rarefaction fan will join at a point where the initial discontinuity lies. The technique for proving Theorem 1 cannot be applied in this case. This is a topic of further research for the authors. Here, however, we present a preliminary numerical study which indicates that such an expansion is also valid. We choose our initial data with a discontinuity of an upward jump, u(x, 0) =
−1, x < 0.5, 1, x 0.5.
We consider the point at t = 1 and x = 0, which lies in the middle of the rarefaction fan. The experiment is carried out in a similar way as in our first test in Section 3.1. The errors and rates of convergence are presented in Table 4, which clearly indicates such an expansion (2.8) is valid.
1706
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
Appendix A We prove here a technical lemma, which provide the basic ingredients in the proof of Theorem 1. Lemma 1. Consider the linear parabolic Cauchy problem z(0, x) = z¯ (x) zt + λ(t, x)z x = φ(t, x) + zxx ,
(A.1)
with λ uniformly bounded. Assume that there exist constants 0 < η 1, η > 1/2η and a smooth curve t → α(t) with the following properties z¯ (x) = 0, φ(t, x) = 0, λ(t, x) α(t) ˙ − η,
x α(0), x α(t), x ∈ α(t), α(t) + η .
(A.2) (A.3) (A.4)
Then the solution satisfies z(t, x) Ce−0 η ,
x α(t) + η .
(A.5)
Here the constants C, 0 > 0 depend on η, λL∞ , φL∞ , and on the local Lipschitz constant of λ in a neighborhood of (t, x), but not on η . Proof. By possibly replacing x by the shifted space variable y = x − α(t), it is not restrictive to . assume that α(t) = 0. For some constants η, λˆ > 0 we thus have λ(t, x) −η, x ∈ [0, η ], λ(t, x) λ, ˆ x ∈ R.
(A.6) (A.7)
The proof will be given in several steps. 1. Consider the fundamental solution Γ (t, x, τ, y) of the linear homogeneous equation zt + λ(t, x)z x = zxx (A.8) with initial data at time τ consisting of a unit mass concentrated at the point y. The first step in the proof is to show that, for any y 0, the total mass that creeps inside the half line x > η /2 is small. Indeed, Γ = Zx , where Z provides the solution to 0 if x < y, Zt + λ(t, x)Zx = Zxx , Z(τ, x) = (A.9) 1 if x > y. Using (A.6) we can construct a lower solution to (A.9) by setting
Z0 (x) =
⎧ 0 ⎪ ⎪ ⎨ 1 − e−ηx ⎪ P (x) ⎪ ⎩ P (η )
if x 0, if x ∈ [0, η /2], if x ∈ [η /2, η ], if x η .
(A.10)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1707
Here P (·) is the quadratic polynomial obtained as the second order Taylor expansion of the function f (x) = 1 − e−ηx at the point x = η /2. Moreover, . η 1 η = + η 2 η is the point where P attains its maximum. Defining the small constant
η2 e−ηη /2 . P , κ = − (η /2) = P 1 − e−ηη /2
(A.11)
it is straightforward to check that the function . Z∗ (t, x) = e−κ(t−τ ) Z0 (x) is a lower solution to (A.9). In particular, for every t > τ we obtain the upper bound ∞
Γ (t, x, τ, y) dx 1 − Z∗ (t, η /2) = 1 − e−κ(t−τ ) 1 − e−ηη /2
η /2
< κ(t − τ ) + e−ηη /2 Cη e−ηη /2 (1 + t − τ )
(A.12)
for a suitable constant Cη . In the case y 0, a better upper bound is obtained simply observing that, by (A.7), the solution to 0 if x < y, ˆ Wt + λWx = Wxx , W (τ, x) = (A.13) 1 if x > y. provides an upper solution to Z in (A.9). In particular, for t > τ and any y ∈ R one has ∞
∞
η /2
e−x /4(t−τ ) dx √ 2 π(t − τ ) 2
Γ (t, x, τ, y) dx η /2−λˆ (t−τ )−y
∞ = (η /2−y)(t−τ )−1/2 −λˆ (t−τ )1/2
e−x /4 √ dx. 2 π 2
Recall that, for s > 0, ∞ s
e−x /4 √ dx 2 π 2
∞ s
e−sx/4 2 2 √ dx = √ e−s /4 . πs 2 π
ˆ − τ ), the above estimates yield When y −2λ(t
(A.14)
1708
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
∞
∞
2
Γ (t, x, τ, y) dx √ −y/2 (t−τ )
η /2
e−x /4 √ dx 2 π
4 2 2 ˆ ˆ e−y /16(t−τ ) √ eλy/8 C eλy/9 −y √ π ( 2 t−τ ) ˆ −2π λy
(A.15)
for a suitable constant C . 2. The solution of the non-homogeneous problem (A.1) admits the integral representation 0 z(t, x) =
t 0 Γ (t, x, 0, y)¯z(y) dy +
−∞
(A.16)
Γ (t, x, τ, y)φ(τ, y) dy dτ. 0 −∞
Therefore, for every t > 0 using the bounds (A.14) and (A.15) we obtain ∞
z(t, x) dx ¯zL∞ ·
−η −2λˆ t
0 +
η /2
−η −2λˆ t
t
∞
Γ (t, x, 0, y) dx dy −∞
η /2 −η −2 λˆ (t−τ )
0
+ φL∞ ·
+
Γ (t, x, 0, y) dx dy dt −∞
−η −2λˆ (t−τ )
0
ˆ · Cη e ¯zL∞ · (η + 2λt)
−ηη /2
∞ η /2
−η −2λˆ t
(1 + t) +
−λˆ y/9
Ce
dy
−∞
t + φL∞ ·
ˆ − τ ) · Cη e−ηη /2 (1 + t − τ ) η + 2λ(t
0 −η −2 λˆ (t−τ )
+
−λˆ y/9
Ce
dy dτ.
(A.17)
−∞
As long as t (η )2 , from (A.17) we deduce ∞
z(t, x) dx C1 ¯zL∞ + φL∞ · (1 + η )3 e−1 η
η /2
C0 ¯zL∞ + φL∞ e−0 η for suitable constants 1 > 0 > 0 and C1 , C0 large enough.
(A.18)
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1709
3. The estimate (A.18) provides an integral bound on the total mass which is carried by diffusion beyond the line x = η /2. In this last step, we transform this integral estimate (A.18) into a pointwise estimate, valid on a region where the drift coefficient λ = λ(t, x) is Lipschitz con. tinuous. Fix a point (t ∗ , x ∗ ) with t ∗ > 0, x ∗ η . Call t0 = max{t ∗ − 1, 0}, and assume that λ is Lipschitz continuous with constant L on the rectangle . ˆ x ∗ + 2λ]. ˆ Q = [t0 , t ∗ ] × [x ∗ − 2λ, . Set λ∗ = λ(t ∗ , x ∗ ) and consider the Green kernel √ ∗ 2 . G∗ (t, x) = e−(x−λ t) /4t /2 πt. Writing (A.1) in the equivalent form zt + λ(t, x) − λ∗ z x + λ∗ zx = φ(t, x) + zxx ,
z(0, x) = z¯ (x),
(A.19)
we can represent its solution as t
∗
z(t) = G (t − t0 ) ∗ z(t0 ) +
∗
t
G (t − s) ∗ φ(s) ds − t0
G∗x (t − s) ∗ λ(s) − λ∗ z(s) ds.
t0
(A.20) Hence η /2 ∞ z(t ∗ , x ∗ ) = + G∗ (t ∗ − t0 , x ∗ − y)z(t0 , y) dy −∞
η /2
t ∗ 0 +
G∗ (t ∗ − s, x ∗ − y)φ(s, y) dy ds
t0 −∞
t η /2 ∞ + G∗x (t ∗ − s, x ∗ − y) λ(s, y) − λ∗ z(s, y) dy ds. −
t0
−∞
(A.21)
η /2
Because of (A.18), it is clear that all of the above integrals become exponentially small as η → ∞. Indeed, the portion containing the singularity can also be estimated, observing that 1 t · G∗x (t, x) √ , 8πe This yields
C |x| · G∗x (t, x) √ . t
1710
t ∗
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711 x ∗ +λ∗(s−t ∗ )+1
G∗x (t ∗ − s, x ∗ − y) λ(s, y) − λ∗ z(s, y) dy ds
t0 x ∗ +λ∗ (s−t ∗ )−1
t ∗
max
t0
t ∗ t0
|x−x ∗ −λ∗ (t−t ∗ )|1
∗ ∗ G (t − s, x ∗ − y) λ(s, y) − λ∗ · x
x ∗ +λ∗(s−t ∗ )+1
z(s, y) dy dt
x ∗ +λ∗ (s−t ∗ )−1
∞ 1 C z(s, y) dy ds. · L· √ +√ ∗ (t − s) 8πe η /2
This establishes (A.5), with a possibly smaller constant 0 > 0.
2
References [1] S. Bianchini, A. Bressan, Vanishing viscosity solutions to nonlinear hyperbolic systems, Ann. of Math. 161 (2005) 223–342. [2] A. Bressan, Hyperbolic Systems of Conservation Laws. The One-Dimensional Cauchy Problem, Oxford Univ. Press, 2000. [3] A. Bressan, A. Marson, Error bounds for a deterministic version of the Glimm scheme, Arch. Ration. Mech. Anal. 142 (1998) 155–176. [4] A. Bressan, T. Yang, On the convergence rate of vanishing viscosity approximations, Comm. Pure Appl. Math. 57 (2004) 1075–1109. [5] B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Lect. Notes Comput. Sci. Eng., vol. 11, Springer-Verlag, 2000. [6] B. Cockburn, C. Johnson, C.W. Shu, E. Tadmor, in: A. Quarteroni (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, in: Lecture Notes in Math., vol. 1697, Springer-Verlag, 1998. [7] C. Dafermos, Polygonal approximations of solutions of the initial value problem for a conservation law, J. Math. Anal. Appl. 38 (1972) 33–41. [8] C. Dafermos, Generalized characteristics and the structure of solutions of hyperbolic conservation laws, Indiana Univ. Math. J. 26 (6) (1977) 1097–1119. [9] C. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, Springer-Verlag, Berlin, 1999. [10] R. DiPerna, Global existence of solutions to nonlinear hyperbolic systems of conservation laws, J. Differential Equations 20 (1976) 187–212. [11] W.H. Fleming, P.E. Souganidis, Asymptotic series and the method of vanishing viscosity, Indiana Univ. Math. J. 35 (1986) 425–447. [12] J. Glimm, Solutions in the large for nonlinear hyperbolic systems of equations, Comm. Pure Appl. Math. 18 (1965) 697–715. [13] J. Glimm, X. Li, Y. Liu, Conservative front tracking in one space dimension, in: Fluid Flow and Transport in Porous Media: Mathematical and Numerical Treatment, South Hadley, MA, 2001, in: Contemp. Math., vol. 295, Amer. Math. Soc., Providence, RI, 2002, pp. 253–264. [14] J. Glimm, X. Li, Y. Liu, Z. Xu, N. Zhao, Conservative front tracking with improved accuracy, SIAM J. Numer. Anal. 41 (2003) 1926–1947. [15] H. Holden, N.H. Risebro, Front Tracking for Hyperbolic Conservation Laws, Springer-Verlag, New York, 2002. [16] G. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [17] S. Kruzhkov, First order quasilinear equations with several space variables, Math. USSR Sb. 10 (1970) 217–243. [18] P. Lax, Hyperbolic systems of conservation laws II, Comm. Pure Appl. Math. 10 (1957) 537–566. [19] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, 1990. [20] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Univ. Press, 2002. [21] O. Oleinik, Discontinuous solutions of nonlinear differential equations, Amer. Math. Soc. Transl. Ser. 2 26 (1957) 95–172. [22] D. Serre, Systems of Conservation Laws I, II, Cambridge Univ. Press, 2000.
W. Shen, Z. Xu / J. Differential Equations 244 (2008) 1692–1711
1711
[23] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in: A. Quarteroni (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, in: Lecture Notes in Math., vol. 1697, Springer-Verlag, 1998, pp. 325–432. [24] J. Smoller, Shock Waves and Reaction–Diffusion Equations, Springer-Verlag, New York, 1983. [25] A. Szpiro, P. Dupuis, Second order numerical methods for first order Hamilton–Jacobi equations, SIAM J. Numer. Anal. 40 (3) (2002) 1136–1183. [26] E. Tadmor, Approximate solutions of nonlinear conservation laws, in: A. Quarteroni (Ed.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, 1997 CIME Course, Cetraro, Italy, June 1997, in: Lecture Motes in Math., vol. 1697, Springer-Verlag, 1998, pp. 1–149.