Efficient numerical stability analysis of detonation waves in ZND
aa r X i v : . [ m a t h . NA ] N ov EFFICIENT NUMERICAL STABILITY ANALYSIS OF DETONATIONWAVES IN ZND
JEFFREY HUMPHERYS AND KEVIN ZUMBRUN
Abstract.
As described in the classic works of Lee–Stewart and Short–Stewart, the numericalevaluation of linear stability of planar detonation waves is a computationally intensive problemof considerable interest in applications. Reexamining this problem from a modern numericalEvans function point of view, we derive a new algorithm for their stability analysis, related toa much older method of Erpenbeck, that, while equally simple and easy to implement as thestandard method introduced by Lee–Stewart, appears to be potentially faster and more stable. Introduction
As described for example in [20, 21, 22, 34, 44, 45], the numerical stability analysis of detona-tion wave solutions of the Zeldovich–von Neumann–D¨oring (ZND), or reactive Euler equations,is a rich and computationally challenging problem. Planar detonation waves can often changestability as physical parameters are varied, undergoing interesting bifurcations to pulsating,spinning, and cellular solutions [12, 23, 2, 32, 35, 29, 47, 48, 49]. This motivates the numericalstudy of their stability, originated by Erpenbeck in [20, 21], both for its interest in its own rightand as a benchmark for more general time-evolution codes [12, 45].Due both to the number of physical parameters (four for a polytropic gas ) and the difficulty ofindividual computations, this problem has proven to be numerically intensive. In their classical1990 paper [34], in which they introduced the algorithm that has become the modern-daystandard, computing accurately for the first time the stability boundaries for one-dimensionaldetonations, Lee and Stewart conclude (p. 131 of the reference): “Finally, we point out thatthough our scheme is direct and easy to implement, complete investigation of the various regionsof parameter space is computationally intensive. Any equivalent or more efficient numericalmethod should be considered a valuable contribution and such approaches are needed to furtherexplore the parameter regimes of instability.”Despite these comments, the basic algorithm introduced by Lee-Stewart (or perhaps vari-ants thereof) as described in the 2006 survey [45] appears still to be the current state of theart. Of course, computational power has increased tremendously in the interim, making once-prohibitive computations now accessible. Nonetheless, it seems of interest to explore moreefficient algorithms if they can be found. Date : November 1, 2018.This work was supported in part by the National Science Foundation award numbers DMS-0607721 andDMS-0300487, and National Science Fountation CAREER award DMS-0847074. Thanks to Mark Williams forstimulating discussions regarding the numerical literature on stability of ZND detonations. Gas constant Γ = γ −
1, heat release coefficient q , activation energy E A , and detonation amplitude [20, 34, 55]. In particular, the computations of [34] were carried out in 1990 on a Cray X-MP/48 super-computer, with several hours required to produce individual figures. (For example, Fig. 9 of[34] tracking the top 6 unstable eigenvalues of detonations of a polytropic gas with gas constant γ = 1 . Hence, the challenge is trans-posed from the level of the national lab to the level of individual users, and from feasibility topractical ease of use. However, the impetus is no less real to reduce computation time fromhours to the minutes required for interactive numerical explorations, and such improvementwould undoubtedly lead to further advances in our understanding of detonation phenomena.Meanwhile, in parallel development, there has been considerable activity, centered aroundthe
Evans function [1, 39, 25], in the numerical evaluation of stability of viscous shock wavesand other traveling front or pulse and boundary layer solutions arising in a variety of equations[14, 15, 16, 13, 28, 4, 26, 27, 5, 7, 8, 18, 9], some of which problems- see, e.g., [27, 5, 7, 18]exhibit complexity rivalling that of detonations. The authors and collaborators have developeda general model-independent method and set of numerical principles for the treatment of suchproblems [28, 54], encoded in the MATLAB-based platform STABLAB [6], which performsextremely well on all of the above-described applications.At the same time, there has been a successful push to place detonation stability in a commonframework with stability of shock waves [50, 35, 36, 29, 49, 52, 55]. In particular, in [50, 29,49, 52, 55], the determination of stability of both viscous (reactive Navier–Stokes) and inviscid(reactive Euler or ZND) detonations has been reduced to the computation of an Evans functiondefined exactly as in the viscous shock and other cases described above. Thus, it is a naturalstep to study ZND stability within this common framework, using the general tools of [28, 54].In this paper, we do exactly that, proposing a new algorithm for the numerical determinationof stability of ZND detonations derived from the point of view of [28, 54].
Surprisingly, thoughboth are shooting methods, this is quite different from the Lee-Stewart algorithm currently instandard use , shooting from x = −∞ to x = 0 rather than from x = 0 to x = −∞ as in[34]; indeed, it is more closely related to the original algorithm of Erpenbeck [21]. The preciserelations between the various methods are described in Section 4.The advantage of shooting from −∞ to 0 is that we seek generalized eigenfunctions decayingexponentially at −∞ . Thus, in the forward direction ( −∞ → → −∞ ), the desired solution decays exponentially while error modesare exponentially amplified , a numerically undesirable situation (“numerical pitfall 1” of [54]).For this reason, we expect that our algorithm should be faster and better conditioned than theLee-Stewart algorithm currently in use. However, there are other aspects that cloud the issue,in particular the singular perturbation structure that arises in the high-activation energy or“square-wave” limit in which instabilities are often studied [20, 22, 23, 17, 2]. For this reason, A Cray X-MP/48 cost roughly $15-20M dollars in the mid-1980’s, having 2 processors with a 105 MHz clockspeed and a theoretical peak performance about 200 MFLOPS per processor or 400 MFLOPS total. A 2010 Mac Pro 8-core (2 quad-core Xeon processors) for example is a $4-5k system with a 2.5GHz clockspeed and a theoretical peak performance around 10 GFLOPS per core or 80 GFLOPS total. Hence, it hasroughly 200 times the processing power at a five thousandth the price (not even adjusting for inflation).
FFICIENT NUMERICAL STABILITY ANALYSIS 3 careful comparison of methods in physically relevant regimes is an important step before makingconclusions.In the present paper, we introduce the algorithm, and give some supporting numerical exper-iments for a simple model equation indicating the advantages of our approach. Followup workin [10, 11] indicates that, also in physically realistic settings, the algorithm performs favorablycompared to the current standard. Specifically, the standard adaptive-mesh version of the al-gorithm described here appears to outperform the fixed-mesh algorithm described in [34, 45] by2-3 orders of magnitude. Much of this improvement appears to be due to the difference betweenfixed and adaptive mesh. However, even compared to an adaptive-mesh version of the methodof Lee-Stewart, our algorithm appears to be 1-10 times faster, depending on the parameterregime: at the least, it is equivalent, and in some situations substantially more efficient.
Plan of the paper.
In Section 2, we review the ZND equations and detonation structure. InSection 3, we give a simple derivation of the Evans/Lopatinski function condition for detonationstability from a general point of view following [50, 29]. For clarity, we specialize in most of thediscussion to the single-species, ideal gas case with Arrhenius ignition dynamics, working in thesame framework as in [34]. The general case is discussed briefly in Remark 5.1. In Section 4,we determine the relation between the derived Evans/Lopatinski condition the related stabilitydeterminants of Erpenbeck [20] and Lee-Stewart [34]. In Section 5, we describe a proposednumerical implementation within the standard STABLAB package developed by the authorsand collaborators. Finally, in Section 6, we present numerical experiments for a simple modelindicating the advantages of integrating in the forward direction and factoring out expecteddecay at −∞ as prescribed in [28, 54].2. ZND detonations
The model.
In Eulerian coordinates the Zeldovich–von Neumann–D¨oring (ZND) equa-tions of reacting gas dynamics in one space dimension may be written as(2.1) ρ t + ( ρu ) x = 0( ρu ) t + ( ρu + p ) x = 0( ρ ˜ E ) t + (( ρ ˜ E + p ) u ) x = 0( ρY ) t + ( ρuY ) x = − ρϕ ( T ) KY, where ρ , u , p , ˜ E , T ∈ R represent density, velocity, pressure, total energy, and temperature,and Y = ( Y , . . . , Y r ) ∈ R r the mass fractions of reactants. Here, ˜ E = u / e is the non-reacting gas-dynamical energy E = u / e modified by chemical potential according to˜ e = e + qY, where e is the specific internal energy of the gas and qY is the specific chemical energy. Thematrix K ∈ R r × r and vector q ∈ R × r measure the rates of reaction and the heat released inreaction, respectively, and ϕ is an “ignition function” that is positive for T above some ignitiontemperature T i and zero for T ≤ T i , serving to “turn on” the reaction. The matrix − K isassumed to be stable, i.e., to have spectrum of strictly negative real part, so that reaction in a Alternatively, the equations may be written in terms of progress variables λ j = 1 − Y j [22, 34, 36]. JEFFREY HUMPHERYS AND KEVIN ZUMBRUN quiescent flow indeed proceeds to the completely burned state Y = 0. In the simplest case of asingle-species, exothermic reaction, Y ∈ R is a scalar, and K and q are positive constants.The system is closed by specifying equations of state (i.e., thermodynamic relations) p = p ( ρ, e, Y ) and T = T ( ρ, e, Y ) and the ignition function. Standard assumptions (in particular,the ones made in [34], etc.) are the ideal gas laws(2.2) p ( ρ, e ) = Γ ρe, T ( e ) = e/C v , where Γ, C v > ϕ ( T ) = exp (cid:18) − E A RT (cid:19) β ( T ) , where E A is the activation energy, R = γC V is the gas constant, and β is an artificial smoothcutoff function with the property that β ≡ T ≥ T i and β ≡ T ≤ T i . Under usualassumptions, the specific form of the function β plays no role in the analysis; see Remark 2.2. Remark 2.1.
More realistic rate laws r ( ρ, T, Y ) may be considered in place of the linear law r = − ρϕ ( T ) KY with little additional difficulty [34] ; however, we lose the explicit form of thereaction profile (5.3) computed in Section 5.1. In the single-species case, these are equivalent. Alternative formulation.
Subtracting q times the fourth equation of (2.1) from thethird equation, we obtain the alternative formulation(2.4) ρ t + ( ρu ) x = 0( ρu ) t + ( ρu + p ) x = 0( ρE ) t + (( ρE + p ) u ) x = ρqϕ ( T ) KY ( ρY ) t + ( ρuY ) x = − ρϕ ( T ) KY in terms of the usual gas-dynamical variables ρ , u , E . We alternate between the two formula-tions as convenient for the analysis.2.3. Detonation waves.
For temperatures T ≤ T i below igition level, equations (2.1) evi-dently reduces to the usual Euler equations of nonreactive gas dynamics, with the reactants Y convected passively by the velocity field u . In particular, so long as T ( ρ ± , e ± , Y ) ≤ T i , theysupport as traveling-wave solutions ordinary gas-dynamical shock waves( ρ, u, E, Y )( x − st ) = ( ( ρ + , u + , E + , Y ) x − st > ρ − , u − , E − , Y ) x − st ≤ s [ ρ ] = [ ρu ] , s [ ρu ] = [ ρu + p ] , s [ ρE ] = [( ρE + p ) u ] , [ Y ] = 0 , or, equivalently,(2.6) s [ ρ ] = [ ρu ] , s [ ρu ] = [ ρu + p ] , s [ ρ ˜ E ] = [( ρ ˜ E + p ) u ] , [ Y ] = 0 , The latter, standard modification circumvents the “cold-boundary difficulty” that the unburned state Y ≡ β ≡
1, and so steady traveling detonation waves do not exist.Though not mentioned, this assumption is also made implicitly in [34], etc.
FFICIENT NUMERICAL STABILITY ANALYSIS 5 where for an arbitrary function h ( ρ, u, ˜ E, Y ), [ h ] := h ( ρ + , u + , ˜ E + , Y + ) − h ( ρ − , u − , ˜ E − , Y − ) de-notes jump across the discontinuity. This also holds if there is no reactant, Y = (0 , . . . , Y = (0 , . . . , T + ≤ T i but T − ≥ T i , with u ± < s (alternatively, T + ≥ T i and T − ≤ T i , with u ± > s ), then there appears a different type of traveling-wavesolution known as a strong detonation , given by ( z = x − st )(2.7) ( ρ, u, E, Y )( z ) = ( ( ρ + , u + , E + , Y ) z > ρ, ¯ u, ¯ E, ¯ Y )( z ) z ≤ , where ¯ Y ( z ) satisfies the smooth traveling-profile ODE(2.8) (¯ ρ (¯ u − s ) ¯ Y ) ′ = − ¯ ρϕ ( ¯ T ) K ¯ Y on ( −∞ , Y (0) = Y , decaying to the completely burned state (0 , . . . , z → −∞ , with (¯ ρ, ¯ u, ¯ E ) = (¯ ρ, ¯ u, ¯ E )( ¯ Y ) determined through the generalized Rankine–Hugoniot relations(2.9) s ¯ ρ − ¯ ρ ¯ u = ( sρ − ρu ) ± s ¯ ρ ¯ u − (¯ ρ ¯ u + ¯ p ) = ( sρu − ( ρu + p )) ± s ¯ ρ ¯˜ E − (¯ ρ ¯˜ E + ¯ p )¯ u = ( sρ ˜ E − ( ρ ˜ E + p ) u ) ± obtained by integrating the remaining traveling-profile equations(2.10) ( s ¯ ρ − ¯ ρ ¯ u ) ′ = 0( s ¯ ρ ¯ u − (¯ ρ ¯ u + ¯ p )) ′ = 0( s ¯ ρ ¯˜ E − (¯ ρ ¯˜ E + ¯ p )¯ u ) ′ = 0from 0 to z (where z <
0) and recalling the Rankine–Hugoniot conditions (2.6) satisfied acrossthe jump at z = 0.That is, strong detonations moving to the right with respect to fluid velocity u (i.e., u < s ,where s is the speed of the detonation) have the structure of an initiating gas-dynamical shockcalled the Neumann shock, which rapidly compresses the gas, raising temperature to the pointof ignition, followed by a reaction zone (the profile (¯ ρ, ¯ u, ¯ E, ¯ Y )) resolving to the final burnedstate. This characteristic “detonation spike” in temperature and pressure profiles agrees wellwith observed features in laboratory experiments.Substituting into (2.8) the first relation in (2.9) and introducing the constant m := ( ρ ( s − u )) ± , we obtain the simplified reaction equation(2.11) Y ′ = m − ρϕ ( T ) KY that we will actually use to solve for the profile. Further simplifying (2.9), we obtain(2.12) s ¯ ρ − ¯ ρ ¯ u = ( sρ − ρu ) ± s ¯ ρ ¯ u − (¯ ρ ¯ u + ¯ p ) = ( sρu − ( ρu + p )) ± s ¯ ρ ¯ E − (¯ ρ ¯ E + ¯ p )¯ u + mq ¯ Y = ( sρE − ( ρE + p ) u + mqY ) ± . An application of the Implicit Function Theorem reveals that (2.9) (as, likewise, the originalODE (2.10)) may be solved for (¯ ρ, ¯ u, ¯ E ) in terms of ¯ Y so long as the gas-dynamical state (¯ ρ, ¯ u, ¯ E ) JEFFREY HUMPHERYS AND KEVIN ZUMBRUN remains noncharacteristic with respect to speed s , or, equivalently, the Rankine–Hugoniot re-lation (2.12) remains full rank in (¯ ρ, ¯ u, ¯ E ). For typical reactions and equations of state, inparticular ideal gas dynamics with single exothermic reaction, this condition holds for all so-lutions of (2.12) with ¯ Y j ≥
0, except for special limiting values of s for which the asymptoticstate ¯ Y = 0 is characteristic, or “sonic”; see, e.g., [35]. These limiting, characteristic wavesare called Chapman–Jouget detonations, and have a special place in the theory. The usual,noncharacteristic type are called overdriven detonations.For our present purposes, the main import of characteristicity is that the eigenvalue equationbecomes singular at x → −∞ in the coordinates we use here, complicating the discussion. Forsimplicity, we restrict hereafter to the overdriven case. The Chapman–Jouget case may betreated similarly using ideas of [34]; see Remark 5.1. Remark 2.2.
For the modified Arrhenius ignition function (2.3) , a standard assumption isthat ¯ T ≥ T i , all x ≤ , T + ≤ T i , so that β ≡ for x ≤ and β ≡ for x ≥ . Under thisassumption, the specific form of the cutoff β plays no role in the analysis. Linearized stability analysis: the Evans–Lopatinski determinant
We now carry out a linearized interface analysis, loosely following [29]. Setting V :=( ρ, u, e ) T , write (2.4) in abstract form as(3.1) F ( W ) t + F ( W ) x = R ( W ) ,W , F j , R ∈ R r , where(3.2) W := (cid:18) VY (cid:19) , F j := (cid:18) f j ( W ) Y g j ( V ) (cid:19) , R := (cid:18) QKY ψ ( W ) − KY ψ ( W ) (cid:19) , (3.3) f := ρρuρ ( e + u / , f := ρuρu + p ( ρ, e, Y )( ρ ( e + u /
2) + p ( ρ, e, Y )) u ,g := ρ, g = ρu, Q := (cid:18) · · · q · · · q r (cid:19) , ψ := ρφ ( T ( ρ, e, Y )) . with V , f j ∈ R , Y ∈ R r , g j , ψ ∈ R , Q ∈ R × r . Remark 3.1.
A minor departure from [19, 20, 50, 29] is to admit the possible dependence ofpressure and temperature on chemical makeup of the gas ( Y ) , an important feature in realisticmodeling of reactive flow. To investigate solutions in the vicinity of a discontinuous detonation profile, we postulateexistence of a single shock discontinuity at location X ( t ), and reduce to a fixed-boundaryproblem by the change of variables x → x − X ( t ). In the new coordinates, the problembecomes(3.4) F ( W ) t + ( F ( W ) − X ′ ( t ) F ( W )) x = R ( W ) , x = 0 , See also the related [50, 35, 36], and the original treatments in [19, 20, 34], etc.
FFICIENT NUMERICAL STABILITY ANALYSIS 7 with jump condition(3.5) X ′ ( t )[ F ( W )] − [ F ( W )] = 0 , [ h ( x, t )] := h (0 + , t ) − h (0 − , t ) as usual denoting jump across the discontinuity at x = 0.3.1. Linearization.
Without loss of generality, suppose for simplicity that the backgroundprofile ¯ W is a steady detonation, i.e., s = 0, hence ( ¯ W , ¯ X ) = ( ¯ W ,
0) is also a steady solution of(3.4)–(3.5). Linearizing (3.4)–(3.5) about the solution ( ¯
W , linearized equations (3.6) A ( W t − X ′ ( t ) ¯ W ′ ( x )) + ( A W ) x = CW, (3.7) X ′ ( t )[ F ( ¯ W )] − [ A W ] = 0 , x = 0 , (3.8) A j := ( ∂/∂W ) F j , C := ( ∂/∂W ) R. Reduction to homogeneous form.
As pointed out in [29], it is convenient for thestability analysis to eliminate the front from the interior equation (3.6). Therefore, we reversethe original transformation to linear order by the change of dependent variables(3.9) W → W − X ( t ) ¯ W ′ ( x ) , following the calculation W ( x − X ( t ) , t )) − W ( x, t ) ∼ X ( t ) W x ( x, t ) ∼ X ( t ) ¯ W ′ ( x ) . approximating to linear order the original, nonlinear transformation. Substituting (3.9) in(3.6)–(3.7), and noting that x -differentiation of the steady profile equation F ( ¯ W ) x = R ( ¯ W )gives(3.10) ( A ( ¯ W ) ¯ W ′ ( x )) x = C ( ¯ W ) ¯ W ′ ( x ) , we obtain modified, homogeneous interior equations(3.11) A W t + ( A W ) x = CW agreeing with those that would be obtained by a naive calculation without consideration of thefront, together with the modified jump condition(3.12) X ′ ( t )[ F ( ¯ W )] − X ( t )[ A ¯ W ′ ( x )] − [ A W ] = 0correctly accounting for front dynamics.The reduction to homogeneous interior equations puts the linearized problem in a standardlinear boundary-value-problem format for which stability may be investigated in straightforwardfashion by the construction of an Evans/Lopatinski determinant . Besides simplifying consider-ably Erpenbeck’s original derivation of his equivalent stability function [20], the homogeneousformat makes possible the application of standard numerical Evans function techniques for itsevaluation. This useful reduction was first carried out, in slightly different form, in [29]. Thetransformation (3.9) is of general use in interface problems, comprising the “good unknown” ofAlinhac [3]. A similar discussion in the simpler context of shock waves may be found in [24];however, in this case, ¯ W ′ ( x ) ≡
0, and so the transformation (3.9) does not make itself evident,nor do front dynamics modify (3.12).
JEFFREY HUMPHERYS AND KEVIN ZUMBRUN
The stability determinant.
Seeking normal mode solutions W ( x, t ) = e λt W ( x ), X ( t ) = e λt X , W bounded, of the linearized equations (3.11)–(3.12), we are led to the generalizedeigenvalue equations ( A W ) ′ = ( − λA + C ) W, x = 0 ,X ( λ [ F ( ¯ W )] − [ A ¯ W ′ ( x )]) − [ A W ] = 0 , where “ ′ ” denotes d/dx , or, setting Z := A W , to(3.13) Z ′ = GZ, x = 0 , (3.14) X ( λ [ F ( ¯ W )] − [ A ¯ W ′ ( x )]) − [ Z ] = 0 , with(3.15) G := ( − λA + C )( A ) − . Here, we are implicitly using the following elementary observation.
Lemma 3.2. A ( ¯ W ( x )) is invertible for all x such that ∂f /∂V is invertible (i.e. V is non-characteristic as a gas-dynamical state with Y held fixed).Proof. Similarly as in the discussion of existence of steady profiles, we may by subtracting Y times the first row of A from the block Y -row, reduce A to block upper-triangular form, withdiagonal blocks ∂f /∂V and g ( V, Y ) I r × r with g ( V, Y ) = ρu = 0. (cid:3) Remark 3.3.
As discussed in Section 2.1, this assumption is essentially necessary already forexistence of a steady profile. In particular, it is satisfied for the usual ideal gas equation of state.
We require also the following fundamental properties.
Lemma 3.4 ([19, 20, 29]) . On ℜ eλ > , the limiting (3 + r ) × (3 + r ) coefficient matrices G ± := lim z →±∞ G ( z ) have unstable subspaces of fixed rank: full rank r for G + and rank r for G − . Moreover, these subspaces have continuous limits as ℜ eλ → .Proof. Straightforward calculation using the fact that G ± are block upper-triangular in ( V, Y );see, e.g., [19, 20, 50, 29] in the case that f , g depend only on V . (cid:3) Corollary 3.5 ([50, 29]) . On ℜ eλ > , the only bounded solution of (3.13) for x > is the triv-ial solution W ≡ . For x < , the bounded solutions consist of an ( r + 2) -dimensional subspace Span { Z +1 , . . . , Z + r +2 } ( λ, x ) of exponentially decaying solutions, analytic in λ and tangent as x → −∞ to the subspace of exponentially decaying solutions of the limiting, constant-coefficientequations Z ′ = G − Z ; moreover, this subspace has a continuous limit as ℜ eλ → .Proof. The first observation is immediate, using the fact that G is constant for x >
0. Thesecond follows from asymptotic ODE theory, using the “gap” or “conjugation” lemmas of[25, 30], [38] together with the fact that G decays exponentially to its end state as x → −∞ .See [29, 52, 55] for details. (cid:3) FFICIENT NUMERICAL STABILITY ANALYSIS 9
Definition 3.6.
We define the Evans–Lopatinski determinant (3.16) D ( λ ) := det (cid:0) Z − ( λ, , · · · , Z − r +2 ( λ, , λ [ F ( ¯ W )] − [ A ¯ W ′ ( x )] (cid:1) = det (cid:0) Z − ( λ, , · · · , Z − r +2 ( λ, , λ [ F ( ¯ W )] + A ¯ W ′ (0 − ) (cid:1) , where Z − j ( λ, x ) are as in Corollary 3.5. The function D is exactly the stability function derived in a different form by Erpenbeck [20];see Section 4.2 below. The formulation (3.16) is of the standard form arising in the simplercontext of (nonreactive) shock stability [37, 19]. Evidently (by (3.14) combined with Corollary3.5), λ is a generalized eigenvalue/normal mode for ℜ eλ ≥ D ( λ ) = 0. Remark 3.7.
As noted in [52, 55] , consideration of the traveling-wave equation F ( W ) ′ = AW ′ = R ( W ) yields the simpler formula (3.17) D ( λ ) = det (cid:0) Z − ( λ, , · · · , Z − r +2 ( λ, , λ [ F ( ¯ W )] + R ( ¯ W (0 − )) (cid:1) . Dual formulation.
The ( n + r ) × ( n + r ) determinant (3.17) may be expressed moresuccinctly in dual form(3.18) D ( λ ) = ˜ Z − ( λ, · ( λ [ F ( ¯ W )] + R ( ¯ W (0 − )) , where ˜ Z − ( λ, x ) is the cross product Z − ∧ · · · ∧ Z − r +2 ( λ, x ) defined by˜ Z − · x = det (cid:0) Z − , · · · , Z − r +2 , x (cid:1) . The vector ˜ Z − may alternatively be characterized directly as the unique up to constant factorbounded solution on x ≤ Z ′ = − G ∗ ˜ Z, which, as x → −∞ is both exponentially decaying and tangent to the corresponding expo-nentially decaying one-dimensional subspace of bounded solutions of the limiting constant-coefficient equations ˜ Z ′ = − G ∗− ˜ Z . It may be specified analytically in λ by the additionalrequirement(3.20) ˜Π( ˜ Z − ) conj ( − M ) = ℓ ( λ ) ,M >
0, where ℓ is an analytically chosen left eigenvector of G − ( λ ) associated with the uniqueeigenvalue g − ( λ ) of negative real part and ˜Π the associated eigenprojection. Here, and else-where, conj denotes complex conjugate. By (3.20) together with the tangency property, ˜ Z − iswell-approximated at x = − M , for M > Z − ( − M ) = ℓ conj ( λ ) . This reduces the approximate evaluation of D ( · ) to the straightforward and extremely well-conditioned numerical problem of integrating a single exponentially growing (in forward direc-tion) mode from x = − M to x = 0. The stability of the computation derives from the fact thaterrors lying in other, exponentially decaying modes, are exponentially damped [54]. Alternate initialization.
Alternatively, following [14, 15, 16, 13], ˜ Z − may be specified byboundary conditions at −∞ , via(3.22) lim x →−∞ e g conj − x ˜ Z − ( x ) = ℓ conj ( λ ) , whence (3.21) becomes(3.23) ˜ Z − ( − M ) = e − g conj − M ℓ conj ( λ ) . This is the method that we prescribe here. It has the advantage of removing the dependenceof ˜ Z − on the artificial parameter M , allowing the flexible choice of M in different parameterregimes, as dictated by numerical considerations, while preserving analyticity. However, inpractice, there is usually not much difference between (3.21) and (3.23). In particular, if, as in[34], one is not interested in analyticity, then one may vary M freely in (3.21) as well.4. Relations to other methods
Relation to the method of Lee and Stewart.
Denoting by Z the solution on x ≤ Z (0) := λ [ F ( ¯ W )] + R ( W (0 − )),we have by standard duality properties that(4.1) ˜ Z − · Z ( λ, x ) ≡ D ( λ )is independent of x ≤
0, or ( ˜ Z − · Z ) ′ ( λ, x ) ≡
0. Taking x = − M and recalling (3.21), we arriveat the alternative Evans–Lopatinski approximation(4.2) D ( λ ) ∼ ℓ conj ( λ ) · Z ( λ, − M )used by Lee and Stewart [34], where ℓ conj · Z ( − M ) = 0 is their “nonradiative condition”enforcing boundedness of Z . The solution of Z from x = 0 to x = − M , on the other hand, isnumerically comparatively ill-conditioned in the vicinity of roots of D ( · ), since Z in this regimeis approximately exponentially decaying in the backward direction while errors are exponentiallygrowing. The version (3.18) is therefore much preferable from the numerical point of view, atleast when used (as here, and in [34]) as a shooting method.4.2.
Relation to the method of Erpenbeck.
Erpenbeck [21] computes ˜ Z − in much thesame way as we do here. However, in place of the homogeneous duality relation (4.1), he usesthe “inhomogeous Abel relation”(4.3) ( ˜ Z − · ˆ Z ) ′ ( λ, x )) = ˜ Z · λ ¯ W ′ ( x ) , valid for the solution ˆ Z of the inhomogeneous equation Z ′ = GZ + λ ¯ W ′ ( x ) with initial dataˆ Z (0) := λ [ ¯ W ] deriving from the unmodified equations (3.6)–(3.7), together with ¯ W ′ ( −∞ ) = 0,to evaluate D ( λ ) = Z −∞ ˜ Z − ( y ) · λ ¯ W ′ ( y ) dy + ˜ Z − (0) · λ [ ¯ W ] . Though it is mathematically equivalent to the homogeneous scheme described above, this hasthe disadvantage that it is difficult to implement adaptive control on truncation error simul-taneously for the ODE and quadrature steps. Indeed, the method is in general a bit morecumbersome to implement and understand than either of the previous two described methods.As a one-time cost, the latter is a rather minor point. However, the implications of the for-mer for performance appear to be significant. Our experience in similar Evans function-type More precisely, they solve the inhomogeneous equations Z ′ = GZ + λ ¯ W ′ ( x ) with initial data ˆ Z (0) := λ [ ¯ W ],and compute ℓ conj ( λ ) · Z ( λ, − M ) ∼ ℓ conj ( λ ) · ˆ Z ( λ, − M ), which is numerically equivalent. Here we are usingˆ Z − Z = ¯ W ′ ( x ) → x → −∞ . FFICIENT NUMERICAL STABILITY ANALYSIS 11 shooting computations [16, 28, 4, 26] of spectra of asymptotically constant-coefficient opera-tors is that a fixed-step scheme can be orders of magnitude slower than a comparable adaptivescheme; see [54] for a general discussion of performance of numerical Evans/Lopatinski solvers.Moreover, even in the solution of ˜ Z alone, the use of an adaptive solver without factoring outexpected decay is much less effective in our experience (“numerical pitfall 3” of [54]).4.3. Expression as boundary-value solver.
We mention in passing an alternative “localEvans function” formulation in the spirit of [34], suggested by Sandstede [41] as a generalmethod for numerical Evans function investigations using collocation/continuation rather thanshooting. By the analysis of the previous subsections, we may recast the eigenvalue equation(3.13)–(3.14) as in [34] as an overdetermined two-point boundary-value problem Z ′ = GZ with r + 4 boundary conditions(4.4) Z (0) := λ [ F ( ¯ W )] + R ( ¯ W (0 − )) , lim x →−∞ ℓ conj · Z ( x ) = 0 . Relaxing at random one of the r + 3 conditions at x = 0, say the requirement on the j thcoordinate, we generically obtain a well-posed boundary-value problem with the correct number r + 3 of boundary conditions; one of the coordinates will always suffice. More, the projectiveboundary-condition at x = −∞ is numerically “correct”, making this problem extremely well-conditioned for solution by collocation/continuation methods (see, e.g., [40]). Defining Z ( λ, x )to be the solution of this relaxed problem, we may then define a local, analytic Evans function˜ D ( λ ) := e j · ( Z ( λ, − ( λ [ F ( ¯ W )] + R ( ¯ W (0 − ))))that is numerically well-conditioned and vanishes if and only if λ is an eigenvalue. This gives asecond way to convert (4.4) into a numerically well-conditioned problem, though the speed andsimplicity of shooting is lost in this approach, along with global analyticity useful for windingnumber calculations. We shall not investigate this method here, but note that it could be usefulin extreme conditions such as the ultra-high activation energy limit [17].5. Numerical implementation
We now describe in detail the numerical algorithm proposed to compute (3.18), following thegeneral approach set out in [16, 28, 53, 54].5.1.
Computing the profile.
In Evans function computations, a delicate aspect is often thecomputation of the background nonlinear profile. We sidestep this issue by the explicit so-lution technique used in [19, 20, 34], modified slightly to accomodate the multi-species case(specifically, the simplified uniform ignition one considered here).Introducing the new variable y defined by(5.1) dy/dx = m − ρϕ ( T ) , y (0) = 0 , where m := ( ρ ( s − u )) ± , we reduce the reaction equation (2.11) to(5.2) dY /dy = KY, Y (0) = Y , obtaining an explicit solution(5.3) Y ( y ) = e Ky Y
02 JEFFREY HUMPHERYS AND KEVIN ZUMBRUN from which the full profile can be recovered through (2.12), either by explicit calculation, ascarried out for ideal gas dynamics in Appendix B, or, more generally, by Newton iteration.
Remark 5.1.
In the single-species case, (5.2) reduces to the change of coordinates x → y :=log Y used in [34] ; general, nonlinear rate laws, or Chapman–Jouget waves, may be accomodatedby a change of variables x → Y r for appropriate r , as discussed in [34] . Computing the stability determinant.
The linearized stability analysis can then becarried out in the variable y defined in (5.1), using the instantaneous change of variables formula(5.4) dx/dy = m/ (¯ ρϕ ( ¯ T )) , y ≤ . Remark 5.2.
Since the righthand side of (5.4) is uniformly positive and bounded, the variables x and y are equivalent in the sense that Cx ≤ y ≤ x/C for x , y ≤ , for some C > . Specifically, we solve from y = − M to y = 0 the ODE ( d/dy ) ˜ Z = − m/ (¯ ρϕ ( ¯ T )) G ∗ ( y ) ˜ Z, withinitial condition ˜ Z − ( − M ) = e − g conj − ( λ ) M ℓ conj ( λ ), M > ℓ ( λ )and limiting eigenvalue g = ( u − + c − ) − are as computed in eqs. (A.2) and (A.3) of AppendixA, the coefficient G ( λ, ¯ V , ¯ Y ) is as described in eqs. (3.15), (3.8), and (3.2)–(3.3), and the profile( ¯ V , ¯ Y )( y ) is as computed in Appendix B. As prescribed in (3.18), we may then compute thestability determinant D ( λ ) = ˜ Z − ( λ, · ( λ [ F ( ¯ W )] + R ( ¯ W ′ (0 − ))) . More precisely, we may solve the numerically more advantageous equations(5.5) ( d/dy ) ˆ Z = − ( m/ ¯ ρϕ ( ¯ T ))( G ( y ) + g − ( λ ) I ) ∗ ˆ Z, with initial conditions ˆ Z ( y ) := e − ( mg conj − / ¯ ρϕ ( ¯ T ) y ˜ Z ( y ), and compute D ( λ ) = ˆ Z − ( λ, · ( λ [ F ( ¯ W )] + R ( ¯ W (0 − ))) . This may readily be computed with good results by an adaptive solver such as the standardRK45; see [16, 28, 54] for further discussion.5.3.
Determination of stability: winding number vs. stability curves.
With an Evanssolver in hand, stability may be checked either by winding number computations as in [21, 4, 26],or by root-following methods based on the Implicit Function Theorem, as in [34]. In the firstmethod, a large semicircle S centered at the origin and lying in ℜ λ ≥ D , andthe number of zeros of D (unstable normal modes) lying within S computed using the principleof the argument, making use of the underlying analyticity of D . Unstable modes lying outside S may be excluded by a separate, asymptotic, argument based on high-frequency behavior of D [14, 15, 26]; for implementations in the context of ZND, see [55, 33] (analytical) or [10, 11](numerical). In the second method, individual roots are followed, avoiding the need to computearound a contour, but typically requiring an extra Newton iteration with each change in modelparameters; see, for example, [34, 44]. Both are by now completely standard.6. A simple model problem
We conclude by an examination of efficiency within the context of a simple but illustrativemodel problem. Consider the ODE(6.1) y ′ = A ( x, λ ) y, A ( x, λ ) = λ (cid:18) c e x − (cid:19) FFICIENT NUMERICAL STABILITY ANALYSIS 13 defined on −∞ < x ≤ x ∈ R , λ ∈ C , y ∈ C , with boundary conditions y ∼ e λx/ (1 , T as x → −∞ and y (0) = (1 , T , modeling a variable-coefficient eigenvalue problem of the formarising in ZND, where the coefficient c = 0 encodes rapidity of exponential decay. As forZND, the coefficient matrix is exponentially asymptotically constant as x → −∞ , with sizegrowing linearly in λ , and has a unique decaying mode as x → −∞ for all ℜ λ >
0, extendingcontinuously to ℜ λ = 0. Thus, we may expect somewhat similar behavior, at least away fromthe high-activation energy “square-wave” regime.In this context, our proposed algorithm consists of factoring out the expected decay e λx/ from the solution to obtain a “neutral” equation(6.2) ˆ y ′ = ˆ A ( x, λ )ˆ y, ˆ A ( x, λ ) = λ (cid:18) c e x − (cid:19) ,y := ˆ ye λx/ , then solving (6.2) from x = − M to x = 0 and checking whether ˆ y (0) lies parallelto (1 , T . For reasonable values of c , a computational domain of M = 5 is sufficient. Themethod of Lee-Stewart, consists roughly of integrating the original equation (6.1) from x = 0to x = − M ; the method of Erpenbeck consists roughly of integrating (6.1) from x = − M to x = 0 without first factoring out expected exponential decay. For comparison, we consideredalso a worst-case scenario with maximum amplification of error modes, integrating (6.2) from x = 0 to x = −∞ .We computed all with the adaptive-mesh RK45 algorithm (ode45) supported in MATLAB, with error tolerance set at the standard level 10 − used for Evans computations [26, 5, 7,8], measuring efficiency by the number of mesh points/function calls required to completethe computation. Extreme cases are λ real- the “best” case, with a spectral gap betweenexponentially growing and exponentially decaying modes at −∞ - and λ imaginary- the “worst”case from our standpoint, with neither spectral gap nor exponential decay. From the standpointof the Lee-Stewart method, the best and worst cases would appear to be reversed.The results, displayed in Tables 1 and 2 for a typical value c = 10, indicate that the proposednew algorithm performs 1-5 times faster than (adaptive versions of) either the Erpenbeck orLee-Stewart methods, depending on the value of λ , with particular improvement as | λ | becomeslarge. It should be noted, moreover, that this is only a comparison of speed (number of meshpoints) for the various methods to produce output with fixed truncation error. If we consideralso accuracy , i.e., convergence error, then the results could be expected to be more dramatic,since both Lee-Stewart and Erpenbeck methods are numerically less well-posed than the forward“neutral” algorithm that we propose. Appendix A. Calculation of ℓ In this appendix, we show how to calculate for general equations of state the initializingvector ℓ ( λ ) used in (3.21), the unique stable left eigenvector of the limiting coefficient matrix(A.1) G − = ( − λA − + C − )( A − ) − = − λf V − ( f V − ) − ( λf V − ( f V − ) − f Y − + QKψ − )( g − ) − − λg − − Kψ − )( g − ) − , In practice, faster than corresponding fixed-mesh methods [54, 10, 11]. mesh pointsforward integration backward integration λ c = 10 100 1000 c = 10 100 10001.0+ 0i 19 14 12 26 24 194.0+ 0i 43 29 19 94 92 8816.0+ 0i 107 76 51 363 361 35764.0+ 0i 261 191 138 1438 1436 1432256.0+ 0i 657 519 427 3177 3186 31920.4+ 0i 14 12 11 17 14 110.4+ 1i 17 13 12 30 27 180.4+ 4i 43 29 19 100 97 730.4+16i 111 77 51 385 382 2960.4+64i 317 224 177 1528 1523 11850.4+256i 1088 870 827 6104 6086 4738 Table 1.
Runs for Eq. (6.2). Forward corresponds to our proposed method,with expected decay factored out. Backward is a worst-case scenario not corre-sponding to any of the methods considered.mesh pointsforward integration backward integration λ c = 10 100 1000 c = 10 100 10001.0+ 0i 23 19 15 19 17 154.0+ 0i 61 58 56 52 50 4916.0+ 0i 181 181 181 186 184 18364.0+ 0i 719 719 719 723 721 721256.0+ 0i 2868 2868 2868 2873 2871 28700.4+ 0i 16 13 12 17 13 120.4+ 1i 20 17 15 20 17 150.4+ 4i 55 52 50 54 52 500.4+16i 196 194 193 197 195 1930.4+64i 765 765 765 775 771 7650.4+256i 3055 3055 3055 3084 3074 3055 Table 2.
Runs for Eq. (6.1). Forward corresponds to Erpenbeck method,backward to Lee–Stewart method.where for a general function h ( V, Y ), we use h − to denote h ( V − , Y − ). Here, we have stronglyused Y − = 0 to obtain the simple upper block-triangular form.By the upper block-triangular form of G − , and the fact that the lower right-hand blockhas spectrum of positive real part for ℜ eλ > g > g < − K is assumed to have spectrum of negative real part), we find that ℓ T mustbe of form ( ℓ TV , ℓ TY ), where ℓ V is the unique unstable eigenvector, associated with eigenvalue α , FFICIENT NUMERICAL STABILITY ANALYSIS 15 of the purely gas-dynamical matrix f V − ( f V − ) − , and(A.2) ℓ TY = ℓ TV ( λf V − ( f V − ) − f Y − + QKψ − )( λ ( g − − αg − ) I + Kψ − ) − . To determine α , ℓ V , and thereby ℓ Y , we observe that f V − ( f V − ) − , is related by similaritytranform M → ( f V − ) − M f V − to the inverse ( f V − ) − f V − of the hyperbolic convection matrix( f V ) − f V of the nonreactive Euler equations V t + ( f V ) − f V V x = 0 written in nonconservativeform in V coordinates with Y ≡
0. Thus, α − is an eigenvalue of ( f V ) − f V , i.e., a hyperboliccharacteristic speed of the non-reactive Euler equations, and ℓ TV = ˜ ℓ T ( f V − ) − , where ˜ ℓ is theassociated left characteristic direction (eigenmode).Noting that α − , as the unique positive characteristic at state V = V − , must be the largestcharacteristic speed, we have by standard formulae [46, 42, 43, 35, 51] or direct calculation(A.3) ℓ TV = ˜ ℓ T ( f V − ) − = ( p ρ − cu + ρ − p e ( u / − e ) , c − ρ − p e u, ρ − p e ) , determining ℓ T ( λ ) = ( ℓ TV , ℓ TY )( λ ) through (A.2). Note that ℓ V is independent of λ . For Y -independent equations of state, (A.2) simplifies considerably, to ℓ TY = ℓ TV QKψ − ( λ ( g − − αg − ) I + Kψ − ) − . Remark A.1.
Noting that the e -component ρ − p e of ℓ V does not vanish in (A.3) , we mayalternatively rescale by ρ/p e to obtain an analytic choice of form ℓ T = ( ∗ , ∗ , , ∗ ) convenient fornumerical solution. A.1.
Alternative, numerical computation.
Alternatively, an analytic choice of ℓ may bedetermined numerically by solution of Kato’s ODE [31] as described in [16, 28, 53, 54]. For ℜ λ bounded from zero, this involves finding numerically at each λ -value the unique stable left andright eigenvectors of G − and computing the associated eigenprojection for use in the Kato ODEas in the general problem-independent method of [16, 28, 53, 54]. At or near ℜ λ = 0, however,this method must be modified, since the stable eigenvector becomes neutral at ℜ λ = 0. Asimple resolution is to notice that, there, the eigenvalues of G − consist of a single eigenvaluewith strictly positive real part, which may be discarded, and three eigenvalues of form g j = α j λ ,where α j (see above) are hyperbolic characteristic speeds for the non-reactive Euler equations,of which the one for which g j /λ = α j < ℓ . Appendix B. Ideal gas profile
In this appendix, we explicitly solve (2.12) for the case of an ideal gas. Restricting to asteady shock, s = 0, and using the ideal gas law (2.2), we may rewrite (2.12) as(B.1) ¯ ρ ¯ u = ρ ± u ± := − m ¯ u + Γ ¯ e ¯ u = u ± + Γ e ± u ± := b ¯ u e + q ¯ Y = u ± e ± + qY ± := c. Combining the second two equations and simplifying gives (Γ+2)¯ u − b ¯ u +2Γ( c − q ¯ Y ) = 0 . Solving using the quadratic formula, we obtain(B.2) ¯ ρ = − m ¯ u , ¯ e = b ¯ u − ¯ u Γ , ¯ u = Γ + 1Γ + 2 b ± s(cid:18) Γ + 1Γ + 2 (cid:19) b + 2Γ( q ¯ Y − c )Γ + 2 , where we have chosen the negative solution branch for ¯ u in accordance with the fact that [ u ] > ρ ] <
0, for a right-moving gas-dynamical shock, so that ¯ u (0 − ) < u + . (Recallthat ¯ u (0 − ) and u + are the two branches of the square root for Y = 1, corresponding to thesolutions of the Rankine–Hugoniot conditions for a nonreacting gas-dynamical shock.) With(5.3), (B.2) gives an explicit expression for the profile as a function of variable y .For a given Neumann shock, there is a one-parameter family of possible endstates ( ρ, u, e ) − determined by the value of q , the maximum value of q corresponding to a Chapman–Jougetwave, for which the argument of the square root vanishes for y = 0. References [1] J. Alexander, R. Gardner and C.K.R.T. Jones. A topological invariant arising in the analysis oftraveling waves.
J. Reine Angew. Math. (1990) 167–212.[2] G. Abouseif and T.Y. Toong, Theory of unstable one-dimensional detonations,
Combust. Flame
Comm. Partial Differential Equations , 14(2):173–230, 1989.[4] B. Barker, J. Humpherys, , K. Rudd, and K. Zumbrun, Stability of viscous shocks in isentropic gasdynamics,
Comm. Math. Phys.
281 (2008), no. 1, 231–249.[5] B. Barker, J. Humpherys, and K. Zumbrun. Stability of isentropic parallel mhd shock layers.
J.Differential Equations , 249(9):2175–2213, 2010.[6] B. Barker, J. Humpherys, and K. Zumbrun. STABLAB: A MATLAB-based numerical library forEvans function computation. Available at: http://impact.byu.edu/stablab/.[7] B. Barker, O. Lafitte, and K. Zumbrun. Existence and stability of viscous shock profiles for 2-Disentropic MHD with infinite electrical resistivity.
Acta Math. Sci. Ser. B Engl. Ed. , 30(2):447–498,2010.[8] B. Barker, M. Lewicka, and K. Zumbrun. Existence and stability of viscoelastic shock profiles, 2010.To appear, Arch. Ration. Mech. Anal.[9] B. Barker, S. Shaw, S. Yarahmadian, and K. Zumbrun, Existence and stability of steady states of areaction convection diffusion equation modeling microtubule formation, to appear, J. Math. Biology.[10] B. Barker and K. Zumbrun, Numerical stability of ZND detonations for Majda’s model, in preparation(2010).[11] B. Barker and K. Zumbrun, Numerical stability of ZND detonations, in preparation.[12] A. Bourlioux, A. Majda, and V. Roytburd, Theoretical and numerical structure for unstable one-dimensional detonations.
SIAM J. Appl. Math.
51 (1991) 303–343.[13] T. J. Bridges, G. Derks, and G. Gottwald. Stability and instability of solitary waves of the fifth-orderKdV equation: a numerical framework.
Phys. D , 172(1-4):190–216, 2002.[14] L. Q. Brin.
Numerical testing of the stability of viscous shock waves.
PhD thesis, Indiana University,Bloomington, 1998.[15] L. Q. Brin. Numerical testing of the stability of viscous shock waves.
Math. Comp. , 70(235):1071–1088,2001.[16] L. Q. Brin and K. Zumbrun. Analytically varying eigenvectors and the stability of viscous shockwaves.
Mat. Contemp. , 22:19–32, 2002. Seventh Workshop on Partial Differential Equations, Part I(Rio de Janeiro, 2001).
FFICIENT NUMERICAL STABILITY ANALYSIS 17 [17] J. Buckmaster and J. Neves, One-dimensional detonation stability: the spectrum for infinite activationenergy.
Phys. Fluids 31 (1988) no. 12, 3572–3576.[18] N. Costanzino, J. Humpherys, T. Nguyen, and K. Zumbrun, Spectral stability of noncharacteristicisentropic Navier-Stokes boundary layers.
Arch. Ration. Mech. Anal.
192 (2009), no. 3, 537–587.[19] J. J. Erpenbeck. Stability of step shocks.
Phys. Fluids , 5:1181–1187, 1962.[20] J. J. Erpenbeck. Stability of steady-state equilibrium detonations.
Physics of Fluids , 5(5):604–614,1962.[21] J. J. Erpenbeck. Stability of idealized one-reaction detonations,
Phys. Fluids , 7 (1964).[22] W. Fickett and W. C. Davis.
Detonation: Theory and Experiment . Dover Publications, 2000.[23] Fickett and Wood, Flow calculations for pulsating one-dimensional detonations.
Phys. Fluids
Shock induced transitions and phase structures in general media , volume 52 of
IMA Vol.Math. Appl. , pages 93–110. Springer, New York, 1993.[25] R. A. Gardner and K. Zumbrun. The gap lemma and geometric criteria for instability of viscous shockprofiles.
Comm. Pure Appl. Math. , 51(7):797–855, 1998.[26] J. Humpherys, O. Lafitte, and K. Zumbrun. Stability of isentropic Navier-Stokes shocks in the high-Mach number limit.
Comm. Math. Phys. , 293(1):1–36, 2010.[27] J. Humpherys, G. Lyng, and K. Zumbrun. Spectral stability of ideal-gas shock layers.
Arch. Ration.Mech. Anal. , 194(3):1029–1079, 2009.[28] J. Humpherys and K. Zumbrun. An efficient shooting algorithm for Evans function calculations inlarge systems.
Phys. D , 220(2):116–126, 2006.[29] H. K. Jenssen, G. Lyng, and M. Williams. Equivalence of low-frequency stability conditions formultidimensional detonations in three models of combustion.
Indiana Univ. Math. J. , 54(1):1–64,2005.[30] T. Kapitula and B. Sandstede. Stability of bright solitary-wave solutions to perturbed nonlinearSchr¨odinger equations.
Phys. D , 124(1-3):58–103, 1998.[31] T. Kato, Perturbation theory for linear operators. Springer–Verlag, Berlin Heidelberg (1985).[32] A.R. Kasimov and D.S. Stewart, Spinning instability of gaseous detonations.
J. Fluid Mech.
J. Fluid Mech. , 216:103–132, 1990.[35] G. Lyng and K. Zumbrun. One-dimensional stability of viscous strong detonation waves.
Arch. Ration.Mech. Anal. , 173(2):213–277, 2004.[36] G. Lyng and K. Zumbrun. A stability index for detonation waves in Majda’s model for reacting flow.
Phys. D , 194(1-2):1–29, 2004.[37] A. Majda. The stability of multidimensional shock fronts.
Mem. Amer. Math. Soc. , 41(275):iv+95,1983.[38] C. Mascia and K. Zumbrun. Stability of large-amplitude viscous shock profiles of hyperbolic-parabolicsystems.
Arch. Ration. Mech. Anal. , 172(1):93–131, 2004.[39] 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.[40] B. Sandstede. Stability of travelling waves. In
Handbook of dynamical systems, Vol. 2 , pages 983–1055.North-Holland, Amsterdam, 2002.[41] B. Sandstede. private communication, 1999.[42] D. Serre.
Systems of conservation laws. 1 . Cambridge University Press, Cambridge, 1999. Hyperbol-icity, entropies, shock waves, Translated from the 1996 French original by I. N. Sneddon.[43] D. Serre.
Systems of conservation laws. 2 . Cambridge University Press, Cambridge, 2000. Geometricstructures, oscillations, and initial-boundary value problems, Translated from the 1996 French originalby I. N. Sneddon. [44] M. Short and D. S. Stewart. The multi-dimensional stability of weak-heat-release detonations.
J. FluidMech. , 382:109–135, 1999.[45] D. S. Stewart and A. R. Kasimov, On the State of Detonation Stability Theory and Its Applicationto Propulsion,
Journal of Propulsion and Power , 22:6, 1230-1244, 2006.[46] J. Smoller.
Shock waves and reaction-diffusion equations . Springer-Verlag, New York, second edition,1994.[47] B. Texier and K. Zumbrun. Galloping instability of viscous shock waves.
Phys. D , 237(10-12):1553–1601, 2008.[48] B. Texier and K. Zumbrun. Hopf bifurcation of viscous shock waves in compressible gas dynamicsand MHD.
Arch. Ration. Mech. Anal. , 190(1):107–140, 2008.[49] B. Texier and K. Zumbrun. Transition to instability of viscous detonation waves is generically asso-ciated with Hopf bifurcation to time-periodic galloping solutions, to appear,
Comm. Math. Physics [50] K. Zumbrun. Multidimensional stability of planar viscous shock waves. In
Advances in the theory ofshock waves , volume 47 of
Progr. Nonlinear Differential Equations Appl. , pages 307–516. Birkh¨auserBoston, Boston, MA, 2001.[51] K. Zumbrun. Stability of large-amplitude shock waves of compressible Navier-Stokes equations. In
Handbook of mathematical fluid dynamics. Vol. III , pages 311–533. North-Holland, Amsterdam, 2004.With an appendix by Helge Kristian Jenssen and Gregory Lyng.[52] K. Zumbrun. Stability of viscous detonations in the ZND limit. to appear, Arch. Rational Mech. Anal.[53] K. Zumbrun. A local greedy algorithm and higher order extensions for global numerical continuationof analytically varying subspaces. to appear, Quart. Appl. Math.[54] K. Zumbrun. Numerical error analysis for evans function computations: a numerical gap lemma,centered-coordinate methods, and the unreasonable effectiveness of continuous orthogonalization.preprint, 2009.[55] K. Zumbrun. High-frequency asymptotics and stability of ZND detonations in the high-overdrive andsmall-heat release limits. preprint, 2010.
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 ::