7.3 The method of lines

A convenient approach both from an implementation point of view as well as for analyzing numerical stability or constructing numerically-stable schemes is to decouple spatial and time discretizations. That is, one first analyzes stability under some spatial approximation assuming time to be continuous (semi-discrete stability) and then finds conditions for time integrators to preserve stability in the fully-discrete case.

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.

7.3.1 Semi-discrete stability

Consider the approximation

vt(t) = Lv, t > 0 (7.76 ) v(0) = f (7.77 )
for the initial value problem (7.1View Equation, 7.2View Equation). The scheme is semi-discrete stable if the solution to Eqs. (7.76View Equation, 7.77View Equation) satisfies the estimate (7.3View Equation).

In the time-independent case, the solution to (7.76View Equation, 7.77View Equation) is

Lt v(t) = e f (7.78 )
and stability holds if and only if there are constants Kd and αd such that
∥eLt∥d ≤ Kde αdt for all t ≥ 0. (7.79 )
The von Neumann condition now states that there exists a constant αd, independent of spatial resolution (i.e., the size of the matrix L), such that the eigenvalues ℓi of L satisfy
max Re (ℓi) ≤ αd. (7.80 ) i
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 L can be uniformly diagonalized. Also, if the lower-order terms are dropped from the analysis then
L = --1---&tidle;L (7.81 ) (Δx )p
with &tidle; L independent of Δx, and in order for (7.80View Equation) to hold for all Δx (in particular small Δx),
max Re(ℓi) ≤ 0, (7.82 ) i
which is a stronger version of the semi-discrete von Neumann condition.

Semi-discrete stability also follows if L is semi-bounded, that is, there is a constant α d independent of resolution such that (cf. Eq. (3.25View Equation) in Theorem 1)

2 ⟨v, Lv ⟩d + ⟨Lv, v⟩d ≤ 2αd ∥v∥d for all v. (7.83 )
In that case, the semi-discrete approximation (7.76View Equation, 7.77View Equation) is numerically stable, as follows immediately from the following energy estimate arguments,
d 2 d 2 --∥v ∥d = --⟨v,v ⟩d = ⟨Lv, v⟩d + ⟨v, Lv ⟩d ≤ 2αd∥v ∥d. (7.84 ) dt dt

For a large class of problems, which can be shown to be well posed using the energy estimate, one can construct semi-bounded operators L by satisfying the discrete counterpart of the properties of the differential operator P in Eq. (7.1View Equation) 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.

7.3.2 Fully-discrete stability

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

vt = Lv, (7.85 )
where L is an N × N 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.80View Equation) and strong (7.82View Equation) 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.85View Equation) 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 N 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.

Necessary conditions.
Recall the von Neumann condition for the semi-discrete system (7.85View Equation): if ℓi is an eigenvalue of L, a necessary condition for semi-discrete stability is [cf. Eq. (7.80View Equation)]

maxi Re (ℓi) ≤ αd. (7.86 )
for some αd independent of N.

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

k+1 k v = Qv , (7.87 )
and recall (cf. Eq. (7.49View Equation)) that a necessary condition for the stability of the fully-discrete system (7.87View Equation) under the assumption (7.48View Equation) that Q depends only on the CFL factor is that its eigenvalues {qi} satisfy
maxi |qi| ≤ 1 (7.88 )
for all spatial resolutions N. Next, assume that the ODE solver is such that
28 28 Q = R (ΔtL ), where R is a polynomial in ΔtL, default As we will discuss later, this includes explicit, one-step Runge–Kutta methods. (7.89)
and notice that if ℓi is an eigenvalue of L, then
ri := R (Δt ℓi) (7.90 )
is an eigenvalue of Q,
{r } ⊂ {q }, (7.91 ) i i
and (7.88View Equation) implies that a necessary condition for fully-discrete stability is
max |ri| ≤ 1. (7.92 ) i

Definition 13. The region of absolute stability of a one-step time integration scheme R is defined as

𝒮ℛ := {z ∈ ℂ : |R (z)| ≤ 1}. (7.93 )

The necessary condition (7.92View Equation) can then be restated as:

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

σ(ΔtL ) ⊂ 𝒮 ℛ, (7.94 )
for all spatial resolutions N.

In the absence of lower-order terms and under the already assumed conditions (7.89View Equation) the strong von Neumann condition (7.82View Equation) then implies that 𝒮 ℛ must overlap the half complex plane {z ∈ ℂ : Re (z ) ≤ 0 }. In particular, this is guaranteed by locally-stable schemes, defined as follows.

Definition 14. An ODE solver R is said to be locally stable if its region of absolute stability 𝒮 ℛ contains an open half disc D− (r) := {μ ∈ ℂ : |μ| < r,Re(μ) < 0} for some r > 0 such that

D − (r) ⊂ 𝒮ℛ. (7.95 )

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

Example 40. Consider the following advection problem with boundaries:

ut = − ux, x ∈ [0,1 ], t ≥ 0, u(t,0) = 0, x = 0, t ≥ 0, u(0,x) = f (x), x ∈ [0,1 ], t = 0,
where f is smooth and compactly supported on [0,1 ], and the one-sided forward Euler discretization, cf. Example 33, with injection of the boundary condition (see Section 10):
k+1 k k k vj---−-vj-= − vj-−-vj−1, for j = 1,2,...N, Δt Δx vk0 = 0, for j = 0.

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

-d dt v = Lv (7.96 )
for (notice that in the following expression the boundary point is excluded and not evolved)
T v = (v1,v2,...,vN) , (7.97 )
with L the banded matrix
( ) − 1 0 0 ⋅⋅⋅ 0 | 1 − 1 0 ⋅⋅⋅ 0 | 1 || || L = ----| 0 1 − 1 ⋅⋅⋅ 0 | , (7.98 ) Δx |( ... ... ... ... 0 |) 0 0 ⋅⋅⋅ 1 − 1
followed by integration in time through the Euler method,
k+1 k v = (1 + ΔtL )v . (7.99 )

Since L is triangular, its eigenvalues are the elements of the diagonal; namely, { ℓi} = {− 1∕Δx }, i.e., there is a single, degenerate eigenvalue q = − 1∕Δx.

The region of local stability of the Euler method is

SE = {z ∈ ℂ : |1 + z | ≤ 1}, (7.100 )
which is a closed disk of radius 1 in the complex plane centered at z0 = − 1; see Figure 2View Image. The von Neumann condition then is
|| || |1 + qΔt | = |1 −-Δt-|≤ 1 (7.101 ) | Δx |
Δt--≤ 2. (7.102 ) Δx
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
Δt-- Δx ≤ 1. (7.103 )

Sufficient conditions.
Under additional assumptions, fully-discrete stability does follow from semi-discrete stability if the time integration is locally stable:

Theorem 11 (Kreiss–Wu [268Jump To The Next Citation Point]). Assume that

  1. 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).
  2. The resulting system of ODEs is integrated with a locally-stable method of the form (7.89View Equation), with stability radius r > 0.
  3. If α ∈ ℝ is such that |α| < r and
    R (iα ) = eiϕ, ϕ ∈ ℝ, (7.104)
    then there is no β ∈ ℝ such that |β | < r, iϕ R (iβ) = e, and β ⁄= α.

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

Δt∥L ∥d ≤ λ for any λ < r. (7.105 )


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 L [282, 409]:

Theorem 12 (Levermore). Suppose L is dissipative, that is Eq. (7.83View Equation) with αd = 0 holds. Then, the third-order Runge–Kutta approximation v(tk+1) = R3(ΔtL )v (tk), where

z2- z3- R3(z) = 1 + z + 2 + 6 , (7.111 )
for the semi-discrete problem (7.85View Equation) is strongly stable,
∥R3 (ΔtL )∥d ≤ 1 (7.112 )
under the CFL timestep restriction Δt ∥L ∥d ≤ 1.

Notice that the restriction α = 0 d is not so severe, since one can always achieve it by replacing L with L − αdI. A generalization of Theorem 12 to higher-order Runge–Kutta methods does not seem to be known.

  Go to previous page Go up Go to next page