7.2 The von Neumann condition

Consider a discretization for a linear system with variable, time-independent coefficients such that
k+1 k v = Qv , (7.40 ) v0 = f, (7.41 )
where k v denotes the gridfunction k k v = {vj : j = 0,1,...,N } and Q is called the amplification matrix. We assume that Q is also time-independent. Then
k k v = Q f (7.42 )
and the approximation (7.40View Equation, 7.41View Equation) is stable if and only if there are constants Kd and αd such that
∥Qk ∥d ≤ Kde αdtk (7.43 )
for all k = 0,1,2,... and high enough resolutions.

In practice, condition (7.43View Equation) is not very manageable as a way of determining if a given scheme is stable since it involves computing the norm of the power of a matrix. A simpler condition based on the eigenvalues {q } i of Q as opposed to the norm of Qk is von Neumann’s one:

|qi| ≤ eαdΔt for all eigenvalues qi of Q and all Δt > 0. (7.44 )
This condition is necessary for numerical stability: if qi is an eigenvalue of Q, k qi is an eigenvalue of k Q and
|qki | ≤ ∥Qk ∥ ≤ Kdeαdtk = Kde αdkΔt. (7.45 )
That is,
|qi| ≤ K1 ∕keαdΔt, (7.46 ) d
which, in order to be valid for all k, implies Eq. (7.44View Equation).

As already mentioned, in order to analyze numerical stability, one can drop lower-order terms. Doing so typically leads to Q depending on Δt and Δx only through a quotient (the CFL factor) of the form (with p = 1 for hyperbolic equations)

λ = -Δt---, (7.47 ) (Δx )p
Q (Δt, Δx ) = Q (λ). (7.48 )
Then for Eq. (7.44View Equation) to hold for all Δt > 0 while keeping the CFL factor fixed (in particular, for small Δt > 0), the following condition has to be satisfied:
|qi| ≤ 1 for all eigenvalues qi of Q, (7.49 )
and one has a stronger version of the von Neumann condition, which is the one encountered in Example 33; see Eq. (7.18View Equation).

7.2.1 The periodic, scalar case

We return to the periodic scalar case, such as the schemes discussed in Examples 33, 34, 35, and 36 with some more generality. Suppose then, in addition to the linearity and time-independent assumptions of the continuum problem, that the initial data and discretization (7.40View Equation, 7.41View Equation) are periodic on the interval [0,2π]. Through a Fourier expansion we can write the grid function f = (f(x0),f (x1),...,f(xN )) corresponding to the initial data as

1 ∑ f = √---- fˆ(ω )eiω, (7.50 ) 2 π ω∈ℤ
where iω iωx0 iωx1 iωx e = (e ,e ,...,e N). The approximation becomes
k --1--∑ ˆ k iω v = √2-π- f(ω)Q e . (7.51 ) ω∈ℤ
Assuming that Q is diagonal in the basis eiω, such that
iω iω Qe = ˆq(ω)e , (7.52 )
as is many times the case, we obtain, using Parseval’s identity,
( ) ∑ 1∕2 ∥vk∥ = |fˆ(ω)|2|ˆq(ω )k|2 . (7.53 ) ω∈ℤ
If
|qˆ(ω )| ≤ eαdΔt for all eigenvalues ω ∈ ℤ and Δt > 0, (7.54 )
for some constant αd then
( ) ∑ 1∕2 ∥vk ∥ ≤ eαdkΔt |fˆ(ω)|2 = eαdkΔt∥f∥ = eαdtk∥f ∥ (7.55 ) ω
and stability follows. Conversely, if the scheme is stable and (7.52View Equation) holds, (7.54View Equation) has to be satisfied. Take
f = eiω, (7.56 )
then
|ˆqk(ω )|∥f ∥ = ∥vk∥ ≤ K eαdtk∥f∥, (7.57 ) d
or
1∕k α Δt |ˆq(ω )| ≤ Kd e d (7.58 )
for arbitrary k, which implies (7.54View Equation). Therefore, provided the condition (7.52View Equation) holds, stability is equivalent to the requirement (7.54View Equation) on the eigenvalues of Q.

7.2.2 The general, linear, time-independent case

However, as mentioned, the von Neumann condition is not sufficient for stability, neither in its original form (7.44View Equation) nor in its strong one (7.49View Equation), unless, for example, Q can be uniformly diagonalized. This means that there exists a matrix T such that

Λ = T −1QT = diag(q0,...,qN ) (7.59 )
is diagonal and the condition number of T with respect to the same norm,
−1 κd(T ) := ∥T ∥d∥T ∥d (7.60 )
is bounded
κ (T ) ≤ K (7.61 ) d d
for some constant Kd independent of resolution (an example is that of Q being normal, QQ ∗ = Q ∗Q). In that case
vk = T ΛkT − 1f (7.62 )
and
k k αdkΔt αdtk ∥v ∥d ≤ κ(T )maxi |qi|∥f ∥d ≤ Kde ∥f ∥d = Kde ∥f∥d. (7.63 )

Next, we discuss two examples where the von Neumann condition is satisfied but the resulting scheme is unconditionally unstable. The first one is for a well-posed underlying continuum problem and the second one for an ill-posed one.

Example 38. An unstable discretization, which satisfies the von Neumann condition for a trivially–well-posed problem [228Jump To The Next Citation Point].

Consider the following system on a periodic domain with periodic initial data

( ) u1 ut = 0, u = u , (7.64 ) 2
discretized as
vk+1 − vk ( 0 1) ----------= − Δx D20vk (7.65 ) Δt 0 0
with D0 given by Eq. (7.32View Equation). The Fourier transform of the amplification matrix and its k-th power are
( 1 λsin2(ωΔx )) ( 1 kλ sin2 (ω Δx )) Qˆ = , Qˆk = . (7.66 ) 0 1 0 1
The von Neumann condition is satisfied, since the eigenvalues are 1. However, the discretization is unstable for any value of λ > 0. For the unit vector T e = (0, 1), for instance, we have
∘ -------------------- |ˆQke | = 1 + (kλ)2sin4(ω Δx ), (7.67 )
which grows without bound as k is increased for sin (ωΔx ) ⁄= 0.

The von-Neumann condition is clearly not sufficient for stability in this example because the amplification matrix not only cannot be uniformly diagonalized, but it cannot be diagonalized at all because of the Jordan block structure in (7.65View Equation).

Example 39. Ill-posed problems are unconditionally unstable, even if they satisfy the von Neumann condition. The following example is drawn from [107Jump To The Next Citation Point].

Consider the periodic Cauchy problem

u = Au , (7.68 ) t x
where u = (u1,u2)T, A is a 2 × 2 constant matrix, and the following discretization. The right-hand side of the equation is approximated by a second-order centered derivative plus higher (third) order numerical dissipation (see Section 8.5)
3 2 2 Aux → AD0v − 𝜖I(Δx ) D+D − v, (7.69 )
where I is the 2 × 2 identity matrix, 𝜖 ≥ 0 an arbitrary parameter regulating the strength of the numerical dissipation and D+, D − are first-order forward and backward approximations of d∕dx,
D+vj := vj+1-−-vj, D − vj := vj-−-vj−1. (7.70 ) Δx Δx
The resulting system of ordinary differential equations is marched in time (method of lines, discussed in Section 7.3) through an explicit method: the iterated Crank–Nicholson (ICN) one with an arbitrary but fixed number of iterations p (see Example 37).

If the matrix A is diagonalizable, as in the scalar case of Example 37, the resulting discretization is numerically stable for λ ≤ 2∕a and p = 2,3, 6,7,10,11,..., even without dissipation. On the other hand, if the system (7.68View Equation) is weakly hyperbolic, as when the principal part has a Jordan block,

( ) a 1 A = 0 a , (7.71 )
one can expect on general grounds that any discretization will be unconditionally unstable. As an illustration, this was explicitly shown in [107] for the above scheme and variations of it. In Fourier space the amplification matrix and its k-th power take the form
( c b) (ck kck−1b ) Qˆ = , ˆQk = k , (7.72 ) 0 c 0 c
with coefficients c,b depending on {a, λ,ωΔx, 𝜖} such that for an arbitrary small initial perturbation at just one gridpoint,
0 T 0 T v0 = (0,2π𝜖) , vj = (0,0) otherwise, (7.73 )
the solution satisfies
k 2 5∕4 0 2 ∥v ∥d ≥ Ck ∥v ∥d for some constant C, (7.74 )
and is therefore unstable regardless of the value of λ and 𝜖. On the other hand, the von Neumann condition |a| ≤ 1 is satisfied if and only if
0 ≤ 𝜖λ ≤ 1∕8. (7.75 )
Notice that, as expected, the addition of numerical dissipation cannot stabilize the scheme independent of its amount. Furthermore, adding dissipation with a strength parameter 𝜖 > 1∕(8λ) violates the von Neumann condition (7.75View Equation) and the growth rate of the numerical instability worsens.
  Go to previous page Go up Go to next page