In general, this method provides only a subclass of numerically-stable approximations. However, it is a very practical one, since spatial and time stability are analyzed separately and stable semi-discrete approximations and appropriate time integrators can then be combined at will, leading to modularity in implementations.

Consider the approximation

for the initial value problem (7.1, 7.2). The scheme is semi-discrete stable if the solution to Eqs. (7.76, 7.77) satisfies the estimate (7.3).In the time-independent case, the solution to (7.76, 7.77) is

and stability holds if and only if there are constants and such that The von Neumann condition now states that there exists a constant , independent of spatial resolution (i.e., the size of the matrix ), such that the eigenvalues of satisfy This is the discrete-in-space version of the Petrovskii condition; see Lemma 2. As already pointed out, it is not always a sufficient condition for stability, unless can be uniformly diagonalized. Also, if the lower-order terms are dropped from the analysis then with independent of , and in order for (7.80) to hold for all (in particular small ), which is a stronger version of the semi-discrete von Neumann condition.Semi-discrete stability also follows if is semi-bounded, that is, there is a constant independent of resolution such that (cf. Eq. (3.25) in Theorem 1)

In that case, the semi-discrete approximation (7.76, 7.77) is numerically stable, as follows immediately from the following energy estimate arguments,For a large class of problems, which can be shown to be well posed using the energy estimate, one can construct semi-bounded operators by satisfying the discrete counterpart of the properties of the differential operator in Eq. (7.1) that were used to show well-posedness. This leads to the construction of spatial differential approximations satisfying the summation by parts property, discussed in Sections 8.3 and 9.4.

Now we consider explicit time integration for systems of the form (7.76, 7.77) with time-independent coefficients. That is, if there are points in space we consider the system of ordinary differential equations (ODEs)

where is an matrix.In the previous Section 7.3.1 we derived necessary conditions for semi-discrete stability of such systems. Namely, the von Neumann one in its weak (7.80) and strong (7.82) forms. Below we shall derive necessary conditions for fully-discrete stability for a large class of time integration methods, including Runge–Kutta ones. Upon time discretization, stability analyses of (7.85) require the introduction of the notion of the region of absolute stability of ODE solvers. Part of the subtlety in the stability analysis of fully-discrete systems is that the size of the system of ODEs is not fixed; instead, it depends on the spatial resolution. However, the obtained necessary conditions for fully-discrete stability will also turn out to be sufficient when combined with additional assumptions. We will also discuss sufficient conditions for fully-discrete stability using the energy method.

Suppose now that the system of ODEs (7.85) is evolved in time using a one-step explicit scheme,

and recall (cf. Eq. (7.49)) that a necessary condition for the stability of the fully-discrete system (7.87) under the assumption (7.48) that depends only on the CFL factor is that its eigenvalues satisfy for all spatial resolutions . Next, assume that the ODE solver is such that and notice that if is an eigenvalue of , then is an eigenvalue of , and (7.88) implies that a necessary condition for fully-discrete stability isThe necessary condition (7.92) can then be restated as:

Lemma 6 (Fully-discrete von Neumann condition for the method of lines.). Consider the semi-discrete system (7.85) and a one-step explicit time discretization (7.87) satisfying the assumptions (7.89). Then, a necessary condition for fully-discrete stability is that the spectrum of the scaled spatial approximation is contained in the region of absolute stability of the ODE solver ,

for all spatial resolutions .In the absence of lower-order terms and under the already assumed conditions (7.89) the strong von Neumann condition (7.82) then implies that must overlap the half complex plane . In particular, this is guaranteed by locally-stable schemes, defined as follows.

Definition 14. An ODE solver is said to be locally stable if its region of absolute stability contains an open half disc for some such that

As usual, the von Neumann condition is not sufficient for numerical stability and we now discuss an example, drawn from [268], showing that the particular version of Lemma 6 is not either.

Example 40. Consider the following advection problem with boundaries:

The corresponding semi-discrete scheme can be written in the form

for (notice that in the following expression the boundary point is excluded and not evolved) with the banded matrix followed by integration in time through the Euler method,Since is triangular, its eigenvalues are the elements of the diagonal; namely, , i.e., there is a single, degenerate eigenvalue .

The region of local stability of the Euler method is

which is a closed disk of radius in the complex plane centered at ; see Figure 2. The von Neumann condition then is or On the other hand, if the initial data has compact support, the numerical solution at early times is that of the periodic problem discussed in Example 33, for which the Fourier analysis gave a stability condition of

Theorem 11 (Kreiss–Wu [268]). Assume that

- A consistent semi-discrete approximation to a constant-coefficient, first-order IBVP is stable in the generalized sense (see Definition 12 and the following remarks).
- The resulting system of ODEs is integrated with a locally-stable method of the form (7.89), with stability radius .
- If is such that and then there is no such that , , and .

Then the fully-discrete system is numerically stable, also in the generalized sense, under the CFL condition

Remarks

- Condition (iii) can be shown to hold for any consistent approximation, if is sufficiently small [268].
- Explicit, one-step Runge–Kutta (RK) methods, which will be discussed in Section 7.5, are in
particular of the form (7.89) when applied to linear, time-independent problems. In fact, consider an
arbitrary, consistent, one-step, explicit ODE solver (7.87) of the form given in Eq. (7.89),
the integer is referred to as the number of stages.
Since the exact solution to Eq. (7.85) is and, in particular,

Eq. (7.106) must agree with the first terms of the Taylor expansion of , where is the order of the global^{28}truncation error of the ODE solver, defined through Therefore, we must have We then see that a scheme of order needs at least stages, , and The above expression in particular shows that when (i.e., when the second sum on the right-hand side is zero) the scheme is unique, with coefficients given by Eq. (7.109). In particular, for , such a scheme corresponds to the Euler one discussed in the example above; see Eq. (7.99). - As we will discuss in Section 7.5, in the nonlinear case it is possible to choose RK methods with if and only if .
- When , first and second-order RK are not locally stable, while third and fourth order are. The fifth-order Dormand–Prince scheme (also introduced in Section 7.5) is also locally stable.

Using the energy method, fully-discrete stability can be shown (resulting in a more restrictive CFL limit) for third-order Runge–Kutta integration and arbitrary dissipative operators [282, 409]:

Theorem 12 (Levermore). Suppose is dissipative, that is Eq. (7.83) with holds. Then, the third-order Runge–Kutta approximation , where

for the semi-discrete problem (7.85) is strongly stable, under the CFL timestep restriction .Notice that the restriction is not so severe, since one can always achieve it by replacing with . A generalization of Theorem 12 to higher-order Runge–Kutta methods does not seem to be known.

Living Rev. Relativity 15, (2012), 9
http://www.livingreviews.org/lrr-2012-9 |
This work is licensed under a Creative Commons License. E-mail us: |