Stability of viscous shocks in isentropic gas dynamics
aa r X i v : . [ m a t h . A P ] M a y STABILITY OF VISCOUS SHOCKS IN ISENTROPICGAS DYNAMICS
BLAKE BARKER, JEFFREY HUMPHERYS,KEITH RUDD, AND KEVIN ZUMBRUN
Abstract.
In this paper, we examine the stability problem for viscousshock solutions of the isentropic compressible Navier–Stokes equations,or p -system with real viscosity. We first revisit the work of Matsumuraand Nishihara, extending the known parameter regime for which small-amplitude viscous shocks are provably spectrally stable by an optimizedversion of their original argument. Next, using a novel spectral energyestimate, we show that there are no purely real unstable eigenvalues forany shock strength.By related estimates, we show that unstable eigenvalues are con-fined to a bounded region independent of shock strength. Then throughan extensive numerical Evans function study, we show that there is nounstable spectrum in the entire right-half plane, thus demonstrating nu-merically that large-amplitude shocks are spectrally stable up to Machnumber M ≈ ≤ γ ≤
3. This strongly suggests that shocksare stable independent of amplitude and the adiabatic constant γ . Wecomplete our study by showing that finite-difference simulations of per-turbed large-amplitude shocks converge to a translate of the originalshock wave, as expected. Introduction
Consider the isentropic compressible Navier-Stokes equations in one spa-tial dimension expressed in Lagrangian coordinates, also known as the p -system with real viscosity: v t − u x = 0 ,u t + p ( v ) x = (cid:16) u x v (cid:17) x , (1)where, physically, v is the specific volume, u is the velocity, and p ( v ) is thepressure law. We assume that p ( v ) is adiabatic, that is, a γ -law gas satisfying p ( v ) = a v − γ for some constants a > γ ≥
1. In the thermodynamicalrarified gas approximation, γ > γ = ( n + 2) /n , where n is the number of internal degrees of freedom of anindividual particle [3]: n = 3 ( γ = 1 . ... ) for monatomic, n = 5 ( γ = 1 . Date : Last Updated: March 23, 2007.This work was supported in part by the National Science Foundation award numbersDMS-0607721 and DMS-0300487. for diatomic gas. More generally, γ is usually taken within 1 ≤ γ ≤ Conclusions and results of numerical investigations:
We carry outour numerical experiments far into the hypersonic shock regime, exploring upthrough Mach number M ≈ ≤ γ ≤
3. Particular attention is given
TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 3 to the monatomic and diatomic cases, γ = 5 / γ = 7 /
4, respectively. Inall cases, our results are consistent with spectral stability (hence also linearand nonlinear stability [23, 22, 31]). This strongly suggests that viscousshock profiles in an isentropic are spectrally stable independent of bothamplitude and γ . Our bounds on the size of unstable eigenvalues, which areindependent of shock strength, may be viewed as a first step in establishingsuch a result analytically. Extensions and open questions.
The present study is not a numericalproof. However, as discussed in [6], it could be converted to numerical proofby the implementation of interval arithmetic and a posteriori error estimatesfor numerical solution of ODE. This would be an interesting direction forfuture investigation. A crucial step in carrying out numerical proof by in-terval arithmetic is by analytical estimates special to the problem at handto sufficiently reduce the computational domain to make the computationsfeasible in realistic time. This we have carried out by the surprisingly strongestimates of Section 5 and Appendices B–C.A second very interesting direction would be to establish stability in thelarge-amplitude limit by a singular-perturbation analysis, an avenue we in-tend to follow in future work. Together, these two projects would give acomplete, rigorous proof of stability for arbitrary-amplitude viscous shocksolutions of (1) on the physical range γ ∈ [1 , Background
By a viscous shock profile of (1), we mean a traveling wave solution v ( x, t ) = ˆ v ( x − st ) ,u ( x, t ) = ˆ u ( x − st ) , moving with speed s and having asymptotically constant end-states ( v ± , u ± ).As an alternative, we can translate x → x − st , and consider instead sta-tionary solutions of v t − sv x − u x = 0 ,u t − su x + ( a v − γ ) x = (cid:16) u x v (cid:17) x . (2)Under the rescaling ( x, t, v, u ) → ( − εsx, εs t, v/ε, − u/ ( εs )), where ε is cho-sen so that 0 < v + < v − = 1, our system takes the form v t + v x − u x = 0 ,u t + u x + ( av − γ ) x = (cid:16) u x v (cid:17) x , (3)where a = a ε − γ − s − . Thus, the shock profiles of (3) satisfy the ordinarydifferential equation v ′ − u ′ = 0 ,u ′ + ( av − γ ) ′ = (cid:18) u ′ v (cid:19) ′ , BARKER, HUMPHERYS, RUDD, AND ZUMBRUN subject to the boundary conditions ( v ( ±∞ ) , u ( ±∞ )) = ( v ± , u ± ). This sim-plifies to v ′ + ( av − γ ) ′ = (cid:18) v ′ v (cid:19) ′ . By integrating from −∞ to x , we get the profile equation(4) v ′ = v ( v − a ( v − γ − , where a is found by setting x = + ∞ , thus yielding the Rankine-Hugoniotcondition(5) a = − v + − v − γ + − v γ + − v + − v γ + . Evidently, a → γ − in the weak shock limit v + →
1, while a ∼ v γ + in thestrong shock limit v + → Remark 2.1.
Since the profile equation (4) is first order scalar, it has amonotone solution. Since v + < v − , we have that ˆ v x < for all x ∈ R . By linearizing (3) about the profile (ˆ v, ˆ u ), we get the eigenvalue problem λv + v ′ − u ′ = 0 ,λu + u ′ − (cid:18) h (ˆ v )ˆ v γ +1 v (cid:19) ′ = (cid:18) u ′ ˆ v (cid:19) ′ , (6)where(7) h (ˆ v ) = − ˆ v γ +1 + a ( γ −
1) + ( a + 1)ˆ v γ We say that a shock profile of (1) is spectrally stable if the linearized system(6) has no spectrum in the closed deleted right half-plane P = {ℜ e ( λ ) ≥ } \ { } , that is, there are no growth or oscillatory modes for (6). We remark that atraveling wave profile always has a zero-eigenvalue associated with its trans-lational invariance. This generally negates the possibility of good uniformbounds in energy estimates, and so we employ the standard technique (see[16, 32]) of transforming into integrated coordinates. This goes as follows:Suppose that ( v, u ) is an eigenfunction of (6) with a non-zero eigenvalue λ . Then ˜ u ( x ) = Z x −∞ u ( z ) dz, ˜ v ( x ) = Z x −∞ v ( z ) dz, and their derivatives decay exponentially as x → ∞ . Thus, by substitutingand then integrating, (˜ u, ˜ v ) satisfies (suppressing the tilde) λv + v ′ − u ′ = 0 , (8a) λu + u ′ − h (ˆ v )ˆ v γ +1 v ′ = u ′′ ˆ v . (8b)This new eigenvalue problem differs spectrally from (6) only at λ = 0, hencespectral stability of (6) is implied by spectral stability of (8). Moreover, TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 5 since (8) has no eigenvalue at λ = 0, one can expect to have a better chanceof developing a successful spectral energy method to prove stability. Wedemonstrate this in the following section.3. Stability of small-amplitude shocks
In this section we further the work in [24] by extending slightly the knownparameter regime for which small-amplitude viscous shocks are provablyspectrally stable. We also show that this method cannot be extended anyfurther for larger amplitudes.
Theorem 3.1 ([24]) . Viscous shocks of (1) are spectrally stable whenever (9) v γ +1+ aγ ! + 2( γ − v γ +1+ aγ ! − ( γ − ≥ In particular, as v + → (hence aγ → ), the left-hand side of (9) ap-proaches γ and so the inequality is satisfied. Therefore, small-amplitudeviscous shocks of (1) are spectrally stable.Proof. We note that h (ˆ v ) >
0. By multiplying (8b) by both the conjugate¯ u and ˆ v γ +1 /h (ˆ v ) and integrating along x from ∞ to −∞ , we have Z R λu ¯ u ˆ v γ +1 h (ˆ v ) + Z R u ′ ¯ u ˆ v γ +1 h (ˆ v ) − Z R v ′ ¯ u = Z R u ′′ ¯ u ˆ v γ h (ˆ v ) . Integrating the last three terms by parts and appropriately using (8a) tosubstitute for u ′ in the third term gives us Z R λ | u | ˆ v γ +1 h (ˆ v ) + Z R u ′ ¯ u ˆ v γ +1 h (ˆ v ) + Z R v ( λv + v ′ ) + Z R ˆ v γ | u ′ | h (ˆ v ) = − Z R (cid:18) ˆ v γ h (ˆ v ) (cid:19) ′ u ′ ¯ u. We take the real part and appropriately integrate by parts to get ℜ e ( λ ) Z R (cid:20) ˆ v γ +1 h (ˆ v ) | u | + | v | (cid:21) + Z R g (ˆ v ) | u | + Z R ˆ v γ h (ˆ v ) | u ′ | = 0 , where g (ˆ v ) = − "(cid:18) ˆ v γ +1 h (ˆ v ) (cid:19) ′ + (cid:18) ˆ v γ h (ˆ v ) (cid:19) ′′ . Thus, to prove stability it suffices to show that g (ˆ v ) ≥ v + , γh (ˆ v ) − ˆ vh ′ (ˆ v ) = aγ ( γ −
1) + ˆ v γ +1 and(10) ˆ v γ − ˆ v x = aγ − h (ˆ v ) . (11) BARKER, HUMPHERYS, RUDD, AND ZUMBRUN
Using (10) and (11), we abbreviate a few intermediate steps below: g (ˆ v ) = − ˆ v x (cid:20) ( γ + 1)ˆ v γ h (ˆ v ) − ˆ v γ +1 h ′ (ˆ v ) h (ˆ v ) + dd ˆ v (cid:20) γ ˆ v γ − h (ˆ v ) − ˆ v γ h ′ (ˆ v ) h (ˆ v ) ˆ v x (cid:21)(cid:21) = − ˆ v x (cid:20) ˆ v γ (( γ + 1) h (ˆ v ) − ˆ vh ′ (ˆ v )) h (ˆ v ) + dd ˆ v (cid:20) γh (ˆ v ) − ˆ vh ′ (ˆ v ) h (ˆ v ) ( aγ − h (ˆ v )) (cid:21)(cid:21) = − a ˆ v x ˆ v γ − h (ˆ v ) × (cid:2) γ ( γ + 1)ˆ v γ +2 − a + 1) γ ( γ − v γ +1 + ( a + 1) γ ( γ − v γ + aγ ( γ + 2)( γ − v − a ( a + 1) γ ( γ − (cid:3) = − a ˆ v x ˆ v γ − h (ˆ v ) [( γ + 1)ˆ v γ +2 + ˆ v γ ( γ −
1) (( γ + 1)ˆ v − ( a + 1) γ ) (12) + aγ ( γ − γ + 2)ˆ v − a ( a + 1) γ ( γ − ≥ − a ˆ v x ˆ v γ − h (ˆ v ) [( γ + 1)ˆ v γ +2 + aγ ( γ − γ + 2)ˆ v − a ( a + 1) γ ( γ − ≥ − γ a ˆ v x ( γ + 1)2 h (ˆ v ) v + v γ +1+ aγ ! + 2( γ − v γ +1+ aγ ! − ( γ − . (13)Thus from (9), we have spectral stability. (cid:3) (cid:3) We note that the hypothesis in (9) is not sharp. Indeed, one can showfrom (12), that a stronger condition could be given as( γ + 1)ˆ v γ +2 + ˆ v γ ( γ −
1) (( γ + 1)ˆ v − ( a + 1) γ ) (14) + aγ ( γ − γ + 2)ˆ v − a ( a + 1) γ ( γ − ≥ , which is sharp in the following sense: When this condition fails to be true,then g (ˆ v ) is no longer nonnegative, and thus the energy method fails. InFigure 1(a), we see that g (ˆ v ) dips on the left-hand side when this inequalityis compromised. We remark further that near v + , the left-hand side of (14)is monotone increasing in ˆ v . Thus, (14) holds if and only if( γ + 1) v γ +2+ + v γ + ( γ −
1) (( γ + 1) v + − ( a + 1) γ ) (15) + aγ ( γ − γ + 2) v + − a ( a + 1) γ ( γ − ≥ . We see from Figure 1(b) that (15) is only a marginal improvement over(9). However, since (15) is sharp, we cannot hope to prove large-amplitudespectral stability using this approach. Instead, we proceed by a combinedanalytical and numerical approach as in [6, 7, 8], first showing that unstableeigenvalues can occur only in a bounded set, then searching for eigenvaluesin this set by computing the Evans function numerically. Before doing so,
TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 7 −4 −0.1−0.08−0.06−0.04−0.0200.020.040.060.080.1 Stable Regime
M=2M=5M=10 (a) (b)
Figure 1.
In (a), we have a graph of g (ˆ v ) against ˆ v for v + = 1 × − and γ = 2 .
0. Note that g (ˆ v ) dips down belowzero on the left-hand size. Hence, the energy estimate willnot generalize beyond the small-amplitude regime. In (b) wegraph the stability boundaries (dark lines) given by (9) and(15), where γ and v + are the horizontal and vertical axes,respectively. We see that in our scaling (15) is only a modestimprovement over (9). The dotted lines correspond from topto bottom as the parameter regimes for Mach numbers 2, 5,and 10 (see the appendix to see how to determine the Machnumber). Hence, this energy estimate does not even hold forshocks at Mach 2 and γ > . No real unstable eigenvalues
In this section we use a novel spectral energy estimate to show that thereare no purely real unstable eigenvalues for any shock strength. We note thatthis result is stronger than that which could be given by the Evans functionstability index (sometimes called the orientation index), which only measuresthe parity of unstable real eigenvalues, see [4, 15, 31]. The fact that thisholds for all shock strengths is interesting because it is among the strongeststatements about large-amplitude spectral stability that has been proven todate.
Theorem 4.1.
Viscous shocks of (1) have no unstable real spectra.Proof.
We multiply (8b) by the conjugate ¯ v and integrate along x from ∞ to −∞ . This gives Z R λu ¯ v + Z R u ′ ¯ v − Z R h (ˆ v ) v ′ ¯ v ˆ v γ +1 = Z R u ′′ ¯ v ˆ v . BARKER, HUMPHERYS, RUDD, AND ZUMBRUN
Notice that on the real line, ¯ λ = λ . Using (8a) to substitute for λv in thefirst term and for u ′′ in the last term, we get Z R u (¯ u ′ − ¯ v ′ ) + Z R u ′ ¯ v − Z R h (ˆ v ) v ′ ¯ v ˆ v γ +1 = Z R ( λv ′ + v ′′ )¯ v ˆ v . Separating terms and simplifying gives Z R u ¯ u ′ + 2 Z R u ′ ¯ v − Z R h (ˆ v ) v ′ ¯ v ˆ v γ +1 = λ Z R v ′ ¯ v ˆ v + Z R v ′′ ¯ v ˆ v . We further simplify by substituting for u ′ in the second term and integratingthe last terms by parts to give, Z R u ¯ u ′ + 2 Z R ( λv + v ′ )¯ v − Z R h (ˆ v ) v ′ ¯ v ˆ v γ +1 = λ Z R v ′ ¯ v ˆ v + Z R ˆ v x ˆ v v ′ ¯ v − Z R | v ′ | ˆ v , which yields Z R u ¯ u ′ + 2 λ Z R | v | + 2 Z R v ′ v − aγ Z R v ′ ¯ v ˆ v γ +1 + Z R | v | ˆ v = λ Z R v ′ ¯ v ˆ v . By taking the real part, we arrive at λ Z R (cid:18) − ˆ v x ˆ v (cid:19) | v | − aγ ( γ + 1)2 Z R ˆ v x ˆ v γ +2 | v | + Z R | v ′ | ˆ v = 0 . This is a contradiction when λ ≥ (cid:3) (cid:3) We remark that the absence of positive real eigenvalues limits the admis-sible onset of instability to Hopf-like bifurcations where a pair of conjugateeigenvalues crosses the imaginary axis. In the following section, we givean upper bound on the spectral frequencies that are admissible. In otherwords, we show that if a pair of conjugate eigenvalues cross the imaginaryaxis, they must do so within these bounds.5.
High-frequency bounds
In this section, we prove high-frequency spectral bounds. This is an im-portant step because it gives a ceiling as to how far along both the imaginaryand real axes that one must explore for spectrum when doing Evans functioncomputations. Indeed if we want to check for roots of the Evans function inthe unstable half-plane, say using the argument principle, then we need onlycompute within these bounds. If no roots are found therein, then we havestrong numerical evidence that the given shock is spectrally stable. (Indeed,at the expense of further effort, such a calculation may be used as the basisof numerical proof, as described in [6, 7].) In this section we show that thehigh-frequency bounds are quite strong, only allowing unstable eigenvaluesto persist in a relatively small triangle adjoining the origin. Moreover, thesebounds are independent of the shock amplitude . TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 9
Lemma 5.1.
The following identity holds for ℜ eλ ≥ : ( ℜ e ( λ ) + |ℑ m ( λ ) | ) Z R ˆ v | u | − Z R ˆ v x | u | + Z R | u ′ | ≤ √ Z R h (ˆ v )ˆ v γ | v ′ || u | + Z R ˆ v | u ′ || u | . (16) Proof.
We multiply (8b) by ˆ v ¯ u and integrate along x . This yields λ Z R ˆ v | u | + Z R ˆ vu ′ ¯ u + Z R | u ′ | = Z R h (ˆ v )ˆ v γ v ′ ¯ u. We get (16) by taking the real and imaginary parts and adding them to-gether, and noting that |ℜ e ( z ) | + |ℑ m ( z ) | ≤ √ | z | . (cid:3) (cid:3) Lemma 5.2.
The following identity holds for ℜ eλ ≥ : (17) Z R | u ′ | = 2 ℜ e ( λ ) Z R | v | + ℜ e ( λ ) Z R | v ′ | ˆ v + 12 Z R (cid:20) h (ˆ v )ˆ v γ +1 + aγ ˆ v γ +1 (cid:21) | v ′ | Proof.
We multiply (8b) by ¯ v ′ and integrate along x . This yields λ Z R u ¯ v ′ + Z R u ′ ¯ v ′ − Z R h (ˆ v )ˆ v γ +1 | v ′ | = Z R v u ′′ ¯ v ′ = Z R v ( λv ′ + v ′′ )¯ v ′ . Using (8a) on the right-hand side, integrating by parts, and taking the realpart gives ℜ e (cid:20) λ Z R u ¯ v ′ + Z R u ′ ¯ v ′ (cid:21) = Z R (cid:20) h (ˆ v )ˆ v γ +1 + ˆ v x v (cid:21) | v ′ | + ℜ e ( λ ) Z R | v ′ | ˆ v . The right hand side can be rewritten as(18) ℜ e (cid:20) λ Z R u ¯ v ′ + Z R u ′ ¯ v ′ (cid:21) = 12 Z R (cid:20) h (ˆ v )ˆ v γ +1 + aγ ˆ v γ +1 (cid:21) | v ′ | + ℜ e ( λ ) Z R | v ′ | ˆ v . Now we manipulate the left-hand side. Note that λ Z R u ¯ v ′ + Z R u ′ ¯ v ′ = ( λ + ¯ λ ) Z R u ¯ v ′ − Z R u (¯ λ ¯ v ′ + ¯ v ′′ )= − ℜ e ( λ ) Z R u ′ ¯ v − Z R u ¯ u ′′ = − ℜ e ( λ ) Z R ( λv + v ′ )¯ v + Z R | u ′ | . Hence, by taking the real part we get ℜ e (cid:20) λ Z R u ¯ v ′ + Z R u ′ ¯ v ′ (cid:21) = Z R | u ′ | − ℜ e ( λ ) Z R | v | . This combines with (18) to give (17). (cid:3) (cid:3)
Remark 5.3.
Lemma 5.2 is a special case of the high-frequency boundsgiven in [22, 31] . We also note that (17) follows from a “Kawashima-type”estimate as described in [22, 31] . Lemma 5.4.
For h (ˆ v ) as in (7) , we have (19) sup ˆ v (cid:12)(cid:12)(cid:12)(cid:12) h (ˆ v )ˆ v γ (cid:12)(cid:12)(cid:12)(cid:12) = γ − v + − v γ + ≤ γ. Proof.
Defining(20) H (ˆ v ) := h (ˆ v )ˆ v − γ = − ˆ v + a ( γ − v − γ + ( a + 1) , we have H ′ (ˆ v ) = − − aγ ( γ − v − γ − < < v + ≤ ˆ v ≤ v − = 1, hencethe maximum of H on ˆ v ∈ [ v + , v − ] is achieved at ˆ v = v + . Substituting (5)into (20) and simplifying yields (19). (cid:3) We complete this section by proving our high-frequency bounds.
Theorem 5.5.
Any eigenvalue λ of (8) with nonnegative real part satisfies (21) ℜ e ( λ ) + |ℑ m ( λ ) | ≤ ( √ γ + 12 ) . Proof.
Using Young’s inequality twice on right-hand side of (16) togetherwith (19), we get( ℜ e ( λ )+ |ℑ m ( λ ) | ) Z R ˆ v | u | − Z R ˆ v x | u | + Z R | u ′ | ≤ √ Z R h (ˆ v )ˆ v γ | v ′ || u | + Z R ˆ v | u ′ || u |≤ θ Z R h (ˆ v )ˆ v γ +1 | v ′ | + ( √ θ Z R h (ˆ v )ˆ v γ ˆ v | u | + ǫ Z R ˆ v | u ′ | + 14 ǫ Z R ˆ v | u | < θ Z R h (ˆ v )ˆ v γ +1 | v ′ | + ǫ Z R | u ′ | + (cid:20) γ θ + 14 ǫ (cid:21) Z R ˆ v | u | . Assuming that 0 < ǫ < θ = (1 − ǫ ) /
2, this simplifies to( ℜ e ( λ ) + |ℑ m ( λ ) | ) Z R ˆ v | u | + (1 − ǫ ) Z R | u ′ | < − ǫ Z R h (ˆ v )ˆ v γ +1 | v ′ | + (cid:20) γ θ + 14 ǫ (cid:21) Z R ˆ v | u | . Applying (17) yields( ℜ e ( λ ) + |ℑ m ( λ ) | ) Z R ˆ v | u | < (cid:20) γ − ǫ + 14 ǫ (cid:21) Z R ˆ v | u | , or equivalently, ( ℜ e ( λ ) + |ℑ m ( λ ) | ) < (4 γ − ǫ − ǫ (1 − ǫ ) . Setting ǫ = 1 / (2 √ γ + 1) gives (21). (cid:3) (cid:3) In the following section, we do an extensive Evans function study to ex-plore numerically the rest of the right-half plane in an effort to locate thepresence of unstable complex eigenvalues.
TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 11 Evans function computation
In this section, we numerically compute the Evans function to locateany unstable eigenvalues, if they exist, in our system. The Evans function D ( λ ) is analytic to the right of the essential spectrum and is defined as theWronskian of decaying solutions of the eigenvalue equation for the linearizedoperator (8) (see [9, 10, 11, 12, 13, 1]). In a spirit similar to the characteristicpolynomial, we have that D ( λ ) = 0 if and only if λ is in the point spectrumof the linearized operator (8). While the Evans function is generally toocomplex to compute explicitly, it can readily be computed numerically, evenfor large systems [19].Since the Evans function is analytic in the region of interest, we can nu-merically compute the winding number in the right-half plane. This allowsus to systematically locate roots (and hence eigenvalues) within. As a re-sult, spectral stability can be determined, and in the case of instability, onecan produce bifurcation diagrams to illustrate and observe its onset. Thisapproach was first used by Evans and Feroe [13] and has been advancedfurther in various directions (see for example [26, 25, 2, 6, 7, 5, 19]).We begin by writing (8) as a first-order system W ′ = A ( x, λ ) W , where(22) A ( x, λ ) = λ
10 0 1 λ ˆ v λ ˆ v f (ˆ v ) − λ , W = uvv ′ , ′ = ddx , and f (ˆ v ) = ˆ v − ˆ v − γ h (ˆ v ), with h as in (7). Note that eigenvalues of (8)correspond to nontrivial solutions of W ( x ) for which the boundary condi-tions W ( ±∞ ) = 0 are satisfied. We remark that since ˆ v is asymptoticallyconstant in x , then so is A ( x, λ ). Thus at each end-state, we have theconstant-coefficient system(23) W ′ = A ± ( λ ) W. Hence solutions that satisfy the needed boundary condition must emergefrom the unstable manifold W − ( x ) at x = −∞ and the stable manifold W +2 ( x ) ∧ W +3 ( x ) at x = ∞ . In other words, eigenvalues of (8) correspondto the values of λ for which these two manifolds intersect, or more precisely,when D ( λ ) := det( W − W +2 W +3 ) | x =0 is zero.As an alternative, we consider the adjoint formulation of the Evans func-tion [26, 4], where instead of computing the 2-dimensional stable manifold,we find the single trajectory f W +1 ( x ) coming from the unstable manifold of(24) f W ′ = − f W A ( x, λ )at x = ∞ . Note that f W +1 ( x ) is biorthogonal to both W +2 ( x ) and W +3 ( x )since ( f W ( x ) W ( x )) ′ = 0 for all x and the initial data of f W is biorthogo-nal to that of W . Hence, the original manifolds intersect when f W +1 and W − are biorthogonal. Therefore, the adjoint Evans function takes the form D + ( λ ) := ( f W +1 W − ) | x =0 . L γ = 1 . γ = 1 . γ = 1 . γ = 1 .
88 1.23(-1) 1.16(-1) 1.08(-1) 1.04(-1)10 2.07(-2) 1.46(-2) 1.75(-2) 1.78(-2)12 2.00(-3) 1.40(-3) 9.85(-4) 7.20(-4)14 6.90(-4) 5.31(-4) 4.73(-4) 4.71(-4)
Table 1.
Relative errors in D ( λ ) are computed by tak-ing the maximum relative error for 60 contour points eval-uated along the semicircle in Figure 2(a). Samples weretaken for varying L and γ , leaving v + fixed at v + = 10 − (Mach M ≈ L = 8 , , , ,
16 and γ = 1 . , . , . , .
8. Relative errors were computed us-ing the next value of L as the baseline.To further improve the numerical efficiency and accuracy of the shootingscheme, we rescale W and f W to remove exponential growth/decay at infinity,and thus eliminate potential problems with stiffness. Specifically, we let W ( x ) = e µ − x V ( x ), where µ − is the growth rate of the unstable manifold at x = −∞ , and we solve instead V ′ ( x ) = ( A ( x, λ ) − µ − I ) V ( x ). We initialize V ( x ) at x = −∞ to be the eigenvector of A − ( λ ) corresponding to µ − .Similarly, it is straightforward to rescale and initialize f W ( x ) at x = ∞ .Numerically, we use a finite domain for x , replacing the end states x = ±∞ with x = ± L , for sufficiently large L . Some care needs to be taken,however, to make sure that we go out far enough to produce good results.In Appendices B and C, we explore the decay rates of the profile ˆ v and A ( x, λ ) and combine our analysis with numerical convergence experimentsto conclude that our domain [ − L, L ] is sufficiently large. Our experiments,described below, were primarily conducted using L = 12, with relative errorbounds ranging mostly between 10 − and 10 − . Larger choices of L wereneeded on the high end of the Mach scale, going up to L = 18 in some cases,to get the relative errors down to 10 − . In Table 1, we provide a sample ofrelative errors in D ( λ ) for large-amplitude shocks.We remark also that in order to produce analytically varying Evans func-tion output, the initial data. V ( − L ) and e V ( L ), must be chosen analytically.The method of Kato [20, pg. 99], also described in [8], does this well by re-placing the eigenvectors of (23) with analytically defined spectral projectors(see also [5, 18]).Before we can compute the Evans function, we first need to computethe traveling wave profile. We use both Matlab’s ode45 routine, whichis the adaptive fourth-order Runge-Kutta-Fehlberg method (RKF45), andMatlab’s bvp4c routine, which is an adaptive Lobatto quadrature scheme.Both methods work well and produce essentially equivalent profiles.Our experiments were carried out on the range( γ, M ) ∈ [1 , × [1 . , . TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 13 (a) (b)
Figure 2.
The graph of the contour and its mapping viathe Evans function. In (a), we have a contour, which is alarge semi-circle aligned on the imaginary axis on the right-half plane (horizontal axis real, vertical axis imaginary). In(b) we have the image of the contour mapped by the Evansfunction. Note that the winding number of this graph iszero. Hence, there are no unstable eigenvalues in the semi-circle. Together with the high-frequency bounds, this impliesspectral stability. Our computation was carried out for γ =5 / v + = 1 × − . This corresponds to a monatomicgas with Mach number M ≈ γ ∈ [1 , M ∈ [1 , . γ, M ) ∈ [1 , × [1 , M = 1 far into the hypersonic shock regime,and encompassing the entire physically relevant range of γ .For each value of γ , the Mach number M was varied on logarithmic scalewith regular mesh from M = 1 . M = 3000. We used 50 mesh pointsfor γ = 1 . . k , where k = 1 , , . . . ,
20. For the special values γ = 1 . .
666 (monatomic and diatomic cases), we did a refined study with 400mesh points.For each value of ( γ, M ), we computed the Evans function along semi-circular contours that contain the triangular region found in the previoussection via our high-frequency bounds; see Figure 2(a). The ODE calcula-tions for individual λ were carried out using Matlab’s ode45 routine, whichis the adaptive 4th-order Runge-Kutta-Fehlberg method (RKF45), and afterscaling out the exponential decay rate of the constant-coefficient solution atspatial infinity, as described in [6, 7, 8, 5, 19]. This method is known tohave excellent accuracy [5, 19]; in addition, the adaptive refinement givesautomatic error control. Typical runs involved roughly 300 mesh points,with error tolerance set to AbsTol = 1e-6 and
RelTol = 1e-8 . Values of λ were varied on a contour with 60 points. As a check on winding numberaccuracy, it was tested a posteriori that the change in argument of D foreach step was less than π/
25. Recall, by Rouch´e’s Theorem, that accuracyis preserved so long as the argument varies by less than π along each meshinterval.In all the cases that we examined, the winding number was zero. Thisindicates that the shocks we considered are spectrally stable, and in viewof [21, 23, 22], nonlinear stability follows. Moreover, in light of the largeMach numbers considered, this is highly suggestive of stability for all shockstrengths. 7. Numerical evolution of strong shocks
In this section, we use a standard finite-difference method to simulateperturbed large-amplitude viscous shocks, and show that they converge, asexpected, to a translate of the original profile.We do this with a nonlinear Crank-Nicholson scheme with a Newton solverto deal with the nonlinearity. This method provides second-order accuracyand will allow for larger time steps than a naive explicit scheme. By differ-entiating the viscosity term in (3), we have v t + v x − u x = 0 ,u t + u x − aγv − γ − v x = u xx v − u x v x v . By implementing the Crank-Nicholson averaging, we obtain v n +1 j − v nj ∆ t + 14∆ x ( v n +1 j +1 − v n +1 j − + v nj +1 − v nj − ) − x ( u n +1 j +1 − u n +1 j − + u nj +1 − u nj − ) = 0and u n +1 j − u nj ∆ t + 14∆ x ( u n +1 j +1 − u n +1 j − + u nj +1 − u nj − ) − a x γ ( v nj ) − γ − ( v n +1 j +1 − v n +1 j − + v nj +1 − v nj − ) − x ) v nj ( u n +1 j +1 − u n +1 j + u n +1 j − + u nj +1 − u nj + u nj − )+ 116(∆ x ) ( v nj ) ( u n +1 j +1 − u n +1 j − + u nj +1 − u nj − ) × ( v n +1 j +1 − v n +1 j − + v nj +1 − v nj − ) = 0 , where n and j are, respectively, the discretized temporal and spacial indices.To cope with the nonlinearities, we use the Newton solver (cid:18) U n +1 V n +1 (cid:19) k +1 = (cid:18) U n +1 V n +1 (cid:19) k − (cid:18) D u F D v FD u G D v G (cid:19) − k , (cid:18) F ( U n +1 , V n +1 ) G ( U n +1 , V n +1 ) (cid:19) k TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 15
10 15 20 25 30 35−1−0.500.511.52 10 15 20 25 30 35−1−0.500.511.5210 15 20 25 30 35−1−0.500.511.52 10 15 20 25 30 35−1−0.500.511.52
Figure 3.
Snapshots of the evolution of a perturbed vis-cous shock wave solution generated by our extended Crank-Nicholson scheme (top curve: v against x ; bottom curve: u against x , with time t fixed and increasing from fig-ure to figure). The parameters used were γ = 1 . v + = 9 × − . This corresponds to a diatomic gas withMach number M ≈ F ( U n +1 , V n +1 ) and G ( U n +1 , V n +1 ) are the above finite differenceschemes, and D u F , D v F , D u G , and D v G are their corresponding partialderivatives. Hence, we use the previous time step as our initial guess inNewton’s method and then iterate until convergence.In Figure 3, we see the evolution of a perturbed viscous shock. We useda perturbation with a positive mass so that we could observe the conver-gence to a translate of the original profile, thus numerically demonstratingnonlinear stability. Appendix A. Mach number for the p -system The Mach number is defined as M = u + − σc + , where u + is the downwind velocity, σ is the shock speed, and c + is the down-wind shock speed, all in Eulerian coordinates. By considering the conserva-tion of mass equation, we have ρ t + ( ρu ) x = 0. Hence, the jump conditionis given by σ [ ρ ] = [ ρu ], which implies, in the original scaling for (1), that σ = u + v − − u − v + v − − v + . Hence, M = u + − σc + = v + ( u − − u + ) c + ( v − − v + ) = v + [ u ] c + [ v ] = − s v + c + or M = (cid:18) u + − σc + (cid:19) = s v − p ′ ( v + ) = s v γa v − γ − = v γ +3+ γ s a . To express this in our scaling, which is given in (3), we need to swap thepluses and minuses. Noting that 0 < v + < v − = 1, we simplify to get M = 1 γv γ + − v γ + − v + = 1 γa . Recalling that a ∼ v γ + as v + →
0, we find that(25) v + ∼ ( γM ) − γ as M → ∞ . In particular,(26) | log v + | ∼ γ − (log γ + 2 log M ) . Appendix B. Profile bounds
Denote profile equation (4) as v ′ = H ( v, v + ) := v ( v − a ( v − γ − Lemma B.1.
For γ ≥ , ≤ x < , (27) 1 ≤ − x γ − x ≤ γ. Proof.
By convexity of f ( x ) = x γ , the difference quotient (27) is increasingin x , bounded above by f ′ (1) = γ , and below by its value at x = 0. (cid:3) Proposition B.2.
For γ ≥ and v + ≤ v ≤ , v + ≤ , (28) − γ ( v − v + ) ≤ H ( v, v + ) ≤ −
34 ( v − v + ) . For γ ≥ , v + ≤ γ , and ≤ v ≤ v − = 1 , (29) 12 ( v − v − ) ≤ H ( v, v + ) ≤ ( v − v − ) . TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 17
Proof.
By (5), H ( v, v + ) = v (cid:16) ( v − − ( v + − v − γ − v − γ + − (cid:17) = v (cid:16) ( v − v + ) + (cid:16) − v + − v γ + (cid:17)(cid:16)(cid:0) v + v (cid:1) γ − (cid:17)(cid:17) = ( v − v + ) (cid:16) v − (cid:16) − v + − v γ + (cid:17)(cid:16) − (cid:0) v + v (cid:1) γ − (cid:0) v + v (cid:1) (cid:17)(cid:17) . Applying (27) with x = v + v , we obtain (28) from v − γ ≤ H ( v, v + ) v − v + ≤ v − (1 − v + )Similarly, we may obtain (29) by the calculation H ( v, v + ) v − v (cid:16) − (cid:16) v + v (cid:17) γ (cid:16) − v + − v γ + (cid:17)(cid:16) − v γ − v (cid:17)(cid:17) ≥ v − γv + ≥ − . (cid:3) Corollary B.3.
For γ ≥ , < v + ≤ , and ˆ v (0) := v + + , the solution ˆ v of (4) satisfies | ˆ v ( x ) − v + | ≤ (cid:16) (cid:17) e − x x ≥ , (30a) | ˆ v ( x ) − v − | ≤ (cid:16) (cid:17) e x +122 x ≤ . (30b) Proof.
Observing that (ˆ v − v + )(0) = , we obtain (30a) by PropositionB.2 and the comparison principle for first-order scalar ODE. Likewise, (30b)follows by (29) together with the observation that, by convexity of H , | H | is bounded below by estimates obtained at ˆ v = v + + and ˆ v = of and , respectively, so that ˆ v traverses [ v + + , ] over an x -interval of length ≤ / / = 12. (cid:3) Remark B.4.
From Proposition B.2 and Corollary B.3, we obtain the re-markable fact that, in the scaling we have chosen, ˆ v decays up to first de-rivative to endstates v ± at a uniform exponential rate independent of shockstrength, despite the apparent singularity at v = v + . Appendix C. Initialization error
Lemma C.1.
For A as in (22) , | · | the Euclidean ( ℓ ) matrix operatornorm, | A ( x, λ ) − A + ( λ ) | ≤ (cid:16) | λ | + 1 + γ ( γ − v − (cid:17) e − x , x ≥ , (31a) | A ( x, λ ) − A − ( λ ) | ≤ (cid:16) | λ | + 1 + 2 γ ( γ − (cid:17) e x +122 , x ≤ . (31b) Proof.
By (22), | A ( x, λ ) − A ± ( λ ) | ≤ λ | ˆ v − v ± | + | f (ˆ v ) − f ( v ± ) | . As computedin the proof of Lemma 5.4, f ′ (ˆ v ) = − − aγ ( γ − v − γ − . Applying (27) tothe expression for a in (5), we obtain v γ + ≤ a ≤ γv γ + , so that | f ′ (ˆ v ) | ≤ γ ( γ − v − ≤ γ ( γ − v − , yielding (31a) by the Mean Value Theorem and (30a). Bound (31b) followssimilarly. (cid:3) Theorem C.2 (Simplified Gap Lemma [15, 6, 7]) . Let e V + and e µ + be a lefteigenvector and associated eigenvalue of − A + ( λ ) and suppose that (32) | e ( − A + − e µ + ) x | ≤ C e − ˆ ηx x ≤ , | ( A − A + )( x ) | ≤ C e − ηx x ≥ , with ≤ ˆ η < η . Then, there exists a solution W = e e µ + x e V ( x, λ ) of (24) with (33) | e V ( x, λ ) − e V + ( λ ) || e V + ( λ ) | ≤ C C e − ηx ( η − ˆ η )(1 − ǫ ) x ≥ L provided ( η − ˆ η ) − C C e − ηL ≤ ǫ . Similar estimates hold for solutions of (22) on x ≤ .Proof. Writing e V ′ = e V ( − A + − e µ + )+ e V ( − A + A + ) and imposing the limitingbehavior e V (+ ∞ , λ ) = e V + , we obtain by Duhamel’s Principle e V ( x ) = e V + − Z + ∞ x e V ( y ) e ( − A + − e µ + )( x − y ) ( A − A + ) dy, from which the result follows by a straightforward Contraction Mappingargument using (32). See [15, 6, 7] for details. (cid:3) Lemma C.3.
For λ satisfying (21) , γ ∈ [1 , , and v + ≤ , (32) (a) holdswith η = 1 / γ , ˆ η = 1 / γ , and C = 10 , and similarly for x ≤ .Proof. Using the inverse Laplace transform representation e ( − A + − e µ + ) x = 12 πi I Γ e zt ( z + A + + e µ + ) − dz, where Γ is a contour enclosing the eigenvalues of ( − A + − e µ + ) and distanceˆ η away, and estimating the resolvent norm | ( z + A + + e µ + ) − | by Kramer’srule, we obtain the stated crude bound, and similarly for x ≤ (cid:3) For error tolerance θ := 10 − k ( | log θ | ∼ k ), define L − ( θ ) := 2 (cid:16) | log 10 − | + | log(2 γ + 7 + 2 γ ( γ − | + | log θ | (cid:17) + 12 ,L + ( θ, v + ) := 43 (cid:16) | log 10 − | + | log(2 γ + 7 + γ ( γ − v − ) | + | log θ | (cid:17) . TABILITY OF VISCOUS SHOCKS IN ISENTROPIC GAS DYNAMICS 19
Corollary C.4.
For λ satisfying (21) , we have relative error bounds (34) | W − ( − L − , λ ) − V − e µ − x || V − e µ − x | , | f W +1 ( L + , λ ) − e V + e e µ + x || e V + e e µ + x | ≤ θ. Proof.
This follows by (33) with (31), Lemma C.3, and (21). (cid:3)
Remark C.5.
Combining (34) with (26) , we find that L + ∼
43 (2 log M + 4 + | log 10 − | + | log θ | ) suffices to obtain relative initialization error less than θ for γ ∈ [1 , and λ in the computational region (21) . For θ = 10 − , we thus obtain θ toleranceup to M = 3 , ∼ for L + ∼ . This is conservative, as we have madelittle effort to optimize bounds, but still within the realm of our experiments.Our numerical convergence studies indicate that L ± = 18 in fact suffices for − accuracy. References [1] J. Alexander, R. Gardner, and C. Jones. A topological invariant arising in the stabilityanalysis of travelling waves.
J. Reine Angew. Math. , 410:167–212, 1990.[2] J. C. Alexander and R. Sachs. Linear instability of solitary waves of a Boussinesq-typeequation: a computer assisted computation.
Nonlinear World , 2(4):471–507, 1995.[3] G. K. Batchelor.
An introduction to fluid dynamics . Cambridge Mathematical Library.Cambridge University Press, Cambridge, paperback edition, 1999.[4] S. Benzoni-Gavage, D. Serre, and K. Zumbrun. Alternate Evans functions and viscousshock waves.
SIAM J. Math. Anal. , 32(5):929–962 (electronic), 2001.[5] T. J. Bridges, G. Derks, and G. Gottwald. Stability and instability of solitary wavesof the fifth-order KdV equation: a numerical framework.
Phys. D , 172(1-4):190–216,2002.[6] L. Q. Brin.
Numerical testing of the stability of viscous shock waves.
PhD thesis,Indiana University, Bloomington, 1998.[7] L. Q. Brin. Numerical testing of the stability of viscous shock waves.
Math. Comp. ,70(235):1071–1088, 2001.[8] L. Q. Brin and K. Zumbrun. Analytically varying eigenvectors and the stability ofviscous shock waves.
Mat. Contemp. , 22:19–32, 2002. Seventh Workshop on PartialDifferential Equations, Part I (Rio de Janeiro, 2001).[9] J. W. Evans. Nerve axon equations. I. Linear approximations.
Indiana Univ. Math.J. , 21:877–885, 1971/72.[10] J. W. Evans. Nerve axon equations. II. Stability at rest.
Indiana Univ. Math. J. ,22:75–90, 1972/73.[11] J. W. Evans. Nerve axon equations. III. Stability of the nerve impulse.
Indiana Univ.Math. J. , 22:577–593, 1972/73.[12] J. W. Evans. Nerve axon equations. IV. The stable and the unstable impulse.
IndianaUniv. Math. J. , 24(12):1169–1190, 1974/75.[13] J. W. Evans and J. A. Feroe. Traveling waves of infinitely many pulses in nerveequations.
Math. Biosci. , 37:23–50, 1977.[14] R. Gardner and C. K. R. T. Jones. A stability index for steady state solutions ofboundary value problems for parabolic systems.
J. Differential Equations , 91(2):181–203, 1991.[15] R. A. Gardner and K. Zumbrun. The gap lemma and geometric criteria for instabilityof viscous shock profiles.
Comm. Pure Appl. Math. , 51(7):797–855, 1998. [16] J. Goodman. Nonlinear asymptotic stability of viscous shock profiles for conservationlaws.
Arch. Rational Mech. Anal. , 95(4):325–344, 1986.[17] P. Howard and K. Zumbrun. Pointwise estimates and stability for dispersive-diffusiveshock waves.
Arch. Ration. Mech. Anal. , 155(2):85–169, 2000.[18] J. Humpherys, B. Sandstede, and K. Zumbrun. Efficient computation of analyticbases in Evans function analysis of large systems.
Numer. Math. , 103(4):631–642,2006.[19] J. Humpherys and K. Zumbrun. An efficient shooting algorithm for evans functioncalculations in large systems.
Physica D , 220(2):116–126, 2006.[20] T. Kato.
Perturbation theory for linear operators . Classics in Mathematics. Springer-Verlag, Berlin, 1995. Reprint of the 1980 edition.[21] C. Mascia and K. Zumbrun. Pointwise Green function bounds for shock profiles ofsystems with real viscosity.
Arch. Ration. Mech. Anal. , 169(3):177–263, 2003.[22] C. Mascia and K. Zumbrun. Stability of large-amplitude viscous shock profiles ofhyperbolic-parabolic systems.
Arch. Ration. Mech. Anal. , 172(1):93–131, 2004.[23] C. Mascia and K. Zumbrun. Stability of small-amplitude shock profiles of symmetrichyperbolic-parabolic systems.
Comm. Pure Appl. Math. , 57(7):841–876, 2004.[24] A. Matsumura and K. Nishihara. On the stability of travelling wave solutions of aone-dimensional model system for compressible viscous gas.
Japan J. Appl. Math. ,2(1):17–25, 1985.[25] R. L. Pego, P. Smereka, and M. I. Weinstein. Oscillatory instability of traveling wavesfor a KdV-Burgers equation.
Phys. D , 67(1-3):45–65, 1993.[26] R. L. Pego and M. I. Weinstein. Eigenvalues, and instabilities of solitary waves.
Philos.Trans. Roy. Soc. London Ser. A , 340(1656):47–94, 1992.[27] J. Smoller.
Shock waves and reaction-diffusion equations . Springer-Verlag, New York,second edition, 1994.[28] B. Texier and K. Zumbrun. Relative Poincar´e-Hopf bifurcation and galloping insta-bility of traveling waves.
Methods Appl. Anal. , 12(4):349–380, 2005.[29] B. Texier and K. Zumbrun. Galloping instability of viscous shock waves. Preprint,2006.[30] B. Texier and K. Zumbrun. Hopf bifurcation of viscous shock waves in compressiblegas- and magnetohydrodynamics. Preprint, 2006.[31] K. Zumbrun. Multidimensional stability of planar viscous shock waves. Lecture notes,TMR summer school, 1999.[32] K. Zumbrun and P. Howard. Pointwise semigroup methods and stability of viscousshock waves.
Indiana Univ. Math. J. , 47(3):741–871, 1998.
Department of Mathematics, Brigham Young University, Provo, UT 84602
E-mail address : [email protected] Department of Mathematics, Indiana University, Bloomington, IN 47402
E-mail address ::