Flux Splitting for stiff equations: A notion on stability

Flux Splitting for stiff equations: A notion on stability

Jochen Schütz J. Schütz and S. Noelle Institut für Geometrie und Praktische Mathematik, RWTH Aachen University
Templergraben 55, 52062 Aachen
Tel.: +49 241 80 97677
22email: {schuetz,noelle}@igpm.rwth-aachen.de
   Sebastian Noelle J. Schütz and S. Noelle Institut für Geometrie und Praktische Mathematik, RWTH Aachen University
Templergraben 55, 52062 Aachen
Tel.: +49 241 80 97677
22email: {schuetz,noelle}@igpm.rwth-aachen.de
Received: date / Accepted: date

For low Mach number flows, there is a strong recent interest in the development and analysis of IMEX (implicit/explicit) schemes, which rely on a splitting of the convective flux into stiff and nonstiff parts. A key ingredient of the analysis is the so-called Asymptotic Preserving (AP) property, which guarantees uniform consistency and stability as the Mach number goes to zero. While many authors have focussed on asymptotic consistency, we study asymptotic stability in this paper: does an IMEX scheme allow for a CFL number which is independent of the Mach number? We derive a stability criterion for a general linear hyperbolic system. In the decisive eigenvalue analysis, the advective term, the upwind diffusion and a quadratic term stemming from the truncation in time all interact in a subtle way. As an application, we show that a new class of splittings based on characteristic decomposition, for which the commutator vanishes, avoids the deterioration of the time step which has sometimes been observed in the literature.

Keywords: IMEX Finite Volume, Asymptotic Preserving, Flux Splitting, Modified Equation, Stability Analysis

AMS subject classification: 35L65, 76M45, 65M08

1 Introduction, Underlying Equations And Flux Splitting

In recent years there has been a renewed interest in the computation of singularly perturbed differential equations. These equations arise, e.g., in the simulation of low-speed fluid flows. Here one is interested in computing waves with vastly different speeds. The goal is to resolve slow waves accurately and efficiently with a large time step, while approximating the fast waves in a stable way, using the same time step.

There is vast literature on the computation of low-speed viscous and inviscid fluid flows. Arguably the first contribution within this field is Chorin’s algorithm Chorin (), who proposes to solve the incompressible Navier-Stokes equations using a projection method. Similar methods have also been used in, e.g., ColellaPao (). A different approach to reduce the stiffness occurring at low Mach numbers is to introduce preconditioning, i.e., to multiply the temporal derivative with a suitable matrix, see the pioneering work by Turkel TurkelPreconditioning (), and, built on this result, the works by Guillard et al. Guillard2 (); Guillard1 (); Guillard3 (). In Guillard1 (), the author identifies the main problem in approximating low-speed flows: Roughly speaking, the variation of pressure is of second order in the Mach number . However, ’traditional’ Godunov schemes tend to produce pressure variations that are of first order in the Mach number, thus for , there can be no uniform convergence. For an extensive additional analysis, in particular with respect to suitable initial conditions, we refer to Dellacherie Dellacherie (). We do not intend to give a fully exhaustive overview on this topic. For an overview on the treatment of low-speed flows, we refer to ChoiMerkle () and the references therein; a more recent survey was given in surveyAsymptotic ().

A class of algorithms that has found particular attention are the Asymptotic Preserving schemes introduced by Jin Jin99 (), built on work with Pareschi and Toscani Jin1998 (). For an excellent review article, consult Jin2012 (); we refer to ArNo12 (); CoDeKu12 (); DegLoNaNe12 (); DeTa (); HaJiLi12 (); ArNoLuMu12 (); JsAP13 () for various applications of this method in different contexts.

Many algorithms, especially those used within the Asymptotic Preserving schemes rely on identifying stiff and nonstiff parts of the underlying equation. This point is generally considered crucial, and the hope is that a well-chosen splitting guarantees a good behavior of the algorithm. A splitting is usually obtained by physical reasoning, see, e.g., the fundamental work by Klein Kl95 ().

Having obtained a splitting into stiff and nonstiff parts, the nonstiff part is then treated explicitly, and the stiff one implicitly. This procedure naturally leads to so-called IMEX schemes as introduced in Ascher1997 (). We refer to Bo07 (); RuBosc12 () for an interesting discussion on the quality of these schemes in the asymptotic limit.

As far as the authors can see, a fully nonlinear asymptotic stability analysis for the non-isentropic Euler equations is still out of reach. In this work, we attempt to reveal an important structural stability property of flux splittings via the considerably simpler modified equation analysis WaHy74 () for a prototype, linear system of conservation laws.

More specifically we derive the modified parabolic system of equations of second order and investigate under what conditions its solutions are bounded in the -norm. For simple problems, one can investigate these conditions analytically. This approach is closely related to the classical von-Neumann stability analysis (see, e.g., Lax1961 (); Richtmyer1957 ()). Strang str64 () showed that, under some assumptions, it is enough to consider only linearized problems, so the approach used in this work is actually more general than it seems at first sight.

Throughout the paper, we consider the linear hyperbolic system of conservation laws


with a constant matrix that has distinct eigenvalues and a full set of corresponding eigenvectors. For simplicity, we set and consider smooth periodic solutions with initial data . Furthermore, we assume that the matrix is a function of a parameter such that with , some eigenvalues of diverge towards infinity.

The motivation to consider (1) stems from considering linearized versions of classical systems of conservation laws


e.g., the (non-dimensionalized) Euler equations at low Mach number , with characteristic quantities density, momentum and total energy, , and


Its linearization around a state yields a matrix with eigenvalues

where denotes the speed of sound. is the ratio of specific heats, frequently taken to be for air. Note that two eigenvalues tend to infinity as .

To highlight the difficulties posed by eigenvalues of multiple scales we briefly discuss a standard, explicit finite volume scheme for (2)

with consistent numerical flux . From the stability conditions by Courant, Friedrichs and Lewy Courant1928 (), it is known that explicit schemes are only stable under a condition, which is typically given by


In the limit as , one mainly wants to resolve the advective wave traveling with speed . Given the restrictive condition (4), this would imply that one needs steps to advect a signal across a single grid cell. For small , this is prohibitively inefficient, and for many schemes also prohibitively dissipative. However, using


as advective condition would result in an unstable scheme.

One potential remedy is to use fully implicit or mixed implicit / explicit (IMEX) methods. The latter class of methods requires a splitting of the flux into components with ’slow’ and ’fast’ waves. More precisely, in the context of (1), one seeks matrices and , such that


with the following conditions posed on and :

Definition 1

The splitting (6) is called admissible, if

  • both and induce a hyperbolic system, i.e., they have real eigenvalues and a complete set of eigenvectors.

  • the eigenvalues of are bounded independently of .

In this work, we give a recipe for identifying stable classes of flux splittings for the linear hyperbolic system (1). We use the well-known modified equation analysis as a tool for (heuristically) investigating (linear) -stability.

The paper is outlined as follows: In Section 2, we introduce the lowest-order IMEX scheme, while in Section 3, we investigate this scheme using the modified equation approach. In Section 4, we introduce so-called characteristic splittings, which are a basic ingredient for a uniformly stable scheme. Based on those sections, we show our main result in Section 5, Theorem 5.1: Characteristic splittings are stable in the sense as explained in Section 3 with a time step size independently of . Section 6 shows an example of a scheme that is only stable with a time step size decreasing with , i.e., the allowable such that the scheme is stable behaves as , which, for small is obviously not the desired effect. In Section 7, we apply our analysis to the linearized Euler equations and show some numerical results that substantiate the theory developed in this paper. Finally, Section 8 offers conclusions and possible future work.

2 IMEX Discretization

Based on a splitting as given in (6), we introduce a straightforward first-order IMEX discretization for (1) based on nonstiff and stiff numerical fluxes and . We assume that the temporal domain is subdivided as

with constant spacing ; and that we have a subdivision

also with constant spacing and cell midpoints . As is customary, we denote a numerical approximation to by . Furthermore, the vector is denoted by . Now we can introduce a (classical) first-order IMEX scheme:

Definition 2

A sequence is a solution to an IMEX discretization, given that



Here, nonstiff and stiff numerical fluxes are defined by


with (positive) numerical viscosities and .

Remark 1

and are given in the so-called viscosity form of a numerical flux, see, e.g., GR1 (). More generally, one can also consider matrices for and instead of scalars. This plays a role in preconditioned schemes; here, it is omitted for the sake of simplicity.

For fixed , consistency analysis of the scheme is well-known, see, e.g., Crouzeix1980 (); GR1 (); Kroener (). However, we have to consider both and as small parameters. The crucial point is that we restrict our analysis to cases where the magnitude of and its first and second derivatives are independent of . (Especially, no derivative behaves as or worse.) This assumption is reasonable, as only those solutions allow for an asymptotic limit as . Similar assumptions have been made in KlMa81 ().

Lemma 1

Let denote the vector with solution to (1) whose derivatives can be bounded independently of ; and let . Furthermore, let . Then, the local truncation error is of order , i.e., there holds for all and for all ,


Apply a Taylor expansion to (7) and note that all the derivatives of can be bounded independently of . ∎

Remark 2

Note that the condition that the derivatives of can be bounded independently of is indeed a condition on the initial datum . Not every initial data gives rise to such a solution. More precisely, this means that as , there is indeed a solution to (1) that does not blow up. In the context of the low-Froude or low-Mach number limit, such a choice of initial conditions is often called well-prepared initial data. For a discussion of the influence of initial conditions on the limit solution, we refer to the work by Klainerman and Majda KlMa81 (). However, also for examples from ordinary differential equations, there are analogues, see, e.g., Bo07 (); HaiWan ().

3 Modified Equation Analysis

In this section, we derive the modified equation WaHy74 () corresponding to (7). As we consider a periodic setting, we can solve the resulting parabolic system explicitly using Fourier series. Using Plancherel’s theorem, we investigate the stability of the modified equation. This yields a practical criterion for the stability of the IMEX scheme.

We start by deriving the modified equation corresponding to (7).

Theorem 3.1

Let be a smooth solution of


Furthermore, we consider vectors and . Then, for fixed and , the IMEX scheme (7) is a second order accurate discretization of (10), i.e.


It is well-known that the modified equation for a first-order discretization is a parabolic equation, i.e., we expect to fulfill


for a (yet unknown) viscosity matrix that is in class . Using (11), one can conclude that


Note that this holds due to , which we will from now on exploit frequently. To simplify the presentation, we slightly abuse our notation, and write for . Using (13) at position ,





From (12),

while . Therefore,


Now we plug (14), (15) and (16) into (7) to obtain, always at position ,

This is if and only if fulfills (11) with


Note again that we have repeatedly used the assumption . This proves the claim. ∎

In the sequel, we show how (10) can be used to determine a necessary condition under what condition the IMEX scheme (7) is stable. We begin by deriving an exact solution to (11) using a Fourier ansatz. Note that is a matrix.

Lemma 2

Let be given by


Furthermore, let be a solution to


Then, admits a representation


with fulfilling the system of ordinary differential equations




and initial conditions


The proof exploits direct computations and starts with assuming that the representation (20) is correct. Thus, plugging (20) into (19), one obtains

Exploiting the linear independence of for different , one obtains (21). ∎

Remark 3
  • Every periodic smooth function can be written as in (18).

  • For future reference, we call the frequency matrices of the modified equation (10).

The following corollary is a direct consequence from the theory of ordinary differential equations, and Plancherel’s theorem.

Corollary 1

We consider the setting as in Lemma 2. Then,


holds for a positive constant if

for all eigenvalues of with .

Remark 4

One might argue that Corollary 1 is not needed in the sense that for every matrix with positive definite, there holds (23). However, this is not a necessary condition. Consider, e.g., the pair of matrices and . Obviously, is not positive definite (note that for, e.g., ), however, the eigenvalues of have negative real part, and consequently, the complete system (19) is stable. (A tedious computation reveals that the eigenvalues of are .)

As already pointed out in the introduction, we have the following important remark concerning commutative matrices:

Remark 5

The real part of the eigenvalues of is not affected by the terms stemming from if matrices and can be simultaneously diagonalized. This is the motivation for introducing so called characteristic splittings in the following section.

4 Characteristic Splitting

In this section, we introduce a new class of splittings that, with our analysis to be presented, turns out to be uniformly stable in without any additional stabilization terms. The splitting relies on a characteristic decomposition of the matrix , i.e., can be decomposed into


for an invertible and . The idea of the characteristic splitting is to split the matrix into stiff and nonstiff parts. We make this more precise in the following definition:

Definition 3

Let be decomposed as in (24), and let be split into

in such a way that and are diagonal matrices and define an admissible splitting of in the sense of Definition 1. (Consequently, the entries of can be bounded independently of .) Subsequently, the characteristic splitting is defined by


Obviously, the splitting of is admissible in the sense of Definition 1.

Let us make the following remark about characteristic splittings:

Remark 6
  • As the system (1) is hyperbolic for all , one can always obtain an admissible characteristic splitting by choosing as , and . Obviously, the eigenvalues of both and are real, and those of are trivially independent of . Such a decomposition would lead to a fully explicit scheme for , which is desirable as for nonstiff equations, those schemes usually are less diffusive than implicit ones.

  • Note that still depends on , and so, even if is independent of , generally is not (but its eigenvalues are).

  • As one can see from Section 5 and especially Section 6, a crucial part in our analysis is the fact that and commute, i.e., . In this way, the ’bad’ modes that can potentially destroy uniform stability are ruled out.

In the sequel, we apply this concept to a prototype matrix , given by


Its eigenvalues are

and for simplicity, we consider to be positive, i.e., is the largest eigenvalue.

In order to be fully explicit for , we use a characteristic splitting with

Consequently, we can derive matrices and via (25) as

The focus of this paper is on uniform stability as , where the fast wave speeds tend to infinity. As outlined in the introduction, the goal is to overcome the inefficiency of a fully explicit scheme due to condition (4), or the instability due to condition (5). In the following (cf. Theorem 5.1), we derive upper bounds on the nonstiff number that assure stability (in a sense to be made more precise) of IMEX scheme (7) for a characteristic splitting.

5 Stability Of Characteristic Flux Splittings

Now, we combine Theorem 3.1 and Corollary 1 to obtain a necessary criterion under what circumstances the IMEX scheme (7) is stable. We start with the general case, and subsequently consider the prototype equation.

5.1 General case

We consider the characteristic splitting (25) in the light of Corollary 1. For a generic splitting with commuting matrices and , the frequency matrix (see (22) and (17)) can be written as


Note that , since constant (in space) solutions of the modified equations are constant in time also. Therefore we need to analyze only the case . As we rely on a characteristic splitting, all the matrices occurring in (29) can be written as for some diagonal matrix . Thus, it is easy to see that the real part of the eigenvalues of is given by

where and are eigenvalues to and , respectively. Claiming that is negative leads to


This leads to a time step restriction that only depends on the explicit part. We summarize this in the following lemma:

Lemma 3

Let be restricted by


Then, (30) and thus Corollary 1 hold.

This is a good result, as can be bounded independently of , and therefore, (31) is also independent of . We summarize this in the following theorem.

Theorem 5.1

The characteristic splitting as introduced in (25) is such that holds for all , eigenvalue to , with a restriction on the time step size that is independent of .

In the sequel, we consider a prototype system in more detail to obtain quantitative information.

5.2 Characteristic Splitting of prototype equation

In this section, we consider the prototype matrix from Section 4, as it allows for easy and explicit computations. The non-dimensionalized advective number corresponding to is denoted by


Note that using , (31) reads

In the sequel, we give stronger bounds on the allowable time step size. Given , we consider the frequency matrix for the characteristic splitting introduced in (25). One potential advantage of the characteristic splitting is that one can compute the eigenvalues explicitly, as all the matrices commute:

Lemma 4

The real part of the eigenvalues of are given by


With the notation introduced in (24) and (25), see also (29), there holds

where the vector is given by

Thus, one can conclude that the eigenvalues of are given in . Sorting the eigenvalues conveniently, starting with the one in the middle, one can conclude that their Real parts are given by formulae (33a) and (33b).∎

The problem under consideration has two asymptotics, namely the one associated to , and the other one associated to (which automatically includes ). We immediately obtain the following condition for the negativity of the first eigenvalue:

Lemma 5

under the condition

Remark 7

Condition (34) is a condition for the nonstiff number from (32). It depends on both and . The coefficient encodes the upwind viscosity of the explicit numerical flux (8), and is usually chosen as , the largest eigenvalue of the nonstiff matrix. There is more freedom to choose the viscosity coefficient of the implicit numerical flux (9), and limiting choices are either (the largest eigenvalue of the nonstiff matrix) or . In both cases, , which gives the sufficient stability condition

This is independent of .

Now we discuss . Obviously, for fixed and , , which directly yields the following Proposition:

Proposition 1

Let be fixed. Then, there exists an , such that for all and all , .

However, this is not the full asymptotics. We therefore change the point of view: Given a fixed , for which is ? The following lemma provides the crucial estimate:

Lemma 6

We define


Now, consider two cases as follows:

  1. Let . Then, holds unconditionally.

  2. Let . Then, holds for


For to be negative, it suffices to show that

We substitute and obtain


(36) is trivially fulfilled, if the right-hand side is not positive, i.e., if

This proves the first claim that for all , both are negative.

Now let . In this case, the right-hand side of (36) is positive, and therefore one has a restriction on . One can compute


Note that for any , there holds . Consequently, (37) is fulfilled if

This proves the lemma. ∎

With Lemmas 5 and 6 we obtain the following

Proposition 2

Recall the definition of from (34). For , and if




From the previous considerations, we can conclude:

Corollary 2

We choose . Then, there holds for all , eigenvalue to , if


Consider expression (38) and plug in the definition of both from (39) and from (34). We consider case first, and obtain

Ergo, is sufficient for (38), which implies (40). Similarly, for , one can easily show that is fulfilled given that (40) holds.∎

Figure 1: Plot of function from (35).
Remark 8
  • The function , see (35), has been plotted in Figure 1. Consider the stability constraint (38), together with the definition of in (39). From Figure 1, one can tell that in (39) is met for a sufficiently small . This then again implies that for ’small’ , one has stability that only depends on the convective number , see (34), as for this case. This is an impressive result in the sense that the only stability restriction for small depends on the slow waves.

  • In Figure 2, we plotted the numerically determined maximum allowable values such that the real parts of the eigenvalues of are negative. For the particular computations, we choose , , , or and determine the maximum , such that and One can infer from this figure that the stability restriction one has to impose on the ratio can be made independent of . Note that if there is no known explicit formula for the eigenvalues of , one could still investigate linearized splittings of, e.g., the Euler equations by simply computing all up to a given value of , getting a first glimpse of possible (in)stabilities in the splitting.

Figure 2: Determined number versus a-priori estimates. Left: , Right: .

6 On The Non-Uniform Stability Of Some Splittings

In this section, we consider a splitting that does not induce a scheme that is uniformly stable. This means that there is no bound on , independently of , such that holds for all and all .

In particular, we consider the (non-characteristic) splitting of matrix from (26) into

The eigenvalues of the splitting matrices are

Obviously, this only yields a hyperbolic splitting in the sense of Definition 1 if we restrict ourselves to . (We are only interested in the case , so this can be done without loss of generality, as it could also be circumvented by a reparametrization of .) For , the splitting is admissible.

The frequency matrix in this case has eigenvalues which can only be computed via an extremely tedious calculation, or with the aid of Maple. Expanded in terms of the asymptotic sequence , these eigenvalues are given by