7.1 Definitions and examples

Consider a well-posed linear initial-value problem (see Definition 3)
ut(t,x) = P (t,x,∂∕∂x )u(t,x), x ∈ ℝn, t ≥ 0, (7.1 ) n u(0,x ) = f(x), x ∈ ℝ . (7.2 )

Definition 11. An approximation-discretization to the Cauchy problem (7.1View Equation, 7.2View Equation) is numerically stable if there is some discrete norm in space ∥ ⋅ ∥d and constants Kd, αd such that the corresponding approximation v satisfies

αdt ∥v (t,⋅)∥d ≤ Kde ∥f ∥d, (7.3 )
for high enough resolutions, smooth initial data f, and t ≥ 0.

Note:

From here on {x ,t } j k denotes some discretization of space and time. This includes both finite difference and spectral collocation methods, which are the ones discussed in Sections 8 and 9, respectively. In addition, we use the shorthand notation

vkj := v(tk,xj). (7.8 )

In order to gain some intuition into the general problem of numerical stability we start with some examples of simple, low-order approximations for a test problem. Consider uniform grids both in space and time

tk = kΔt, xj = jΔx, k = 0,1,2,..., j = 0,1,2,...N, (7.9 )
and the advection equation,
ut = aux, x ∈ [0,2π ], t ≥ 0, (7.10 )
on a periodic domain with 2π = N Δx, and smooth periodic initial data. Then the solution u can be represented by a Fourier series:
1 ∑ u(t,x ) = √---- ˆu(t,ω )eiωx, (7.11 ) 2π ω∈ℤ
where
∫ 2π ˆu (t,ω ) = √1--- e− iωxu (t,x)dx, (7.12 ) 2π 0
and the stability of the following schemes can be analyzed in Fourier space.

Example 33. The one-sided Euler scheme.
Eq. (7.10View Equation) is discretized with a one-sided FD approximation for the spatial derivative and evolved in time with the forward Euler scheme,

k+1 k k k vj---−-vj-= avj+1-−-vj. (7.13 ) Δt Δx
In Fourier space the approximation becomes
ˆvk+1(ω) = ˆq(ω )ˆvk(ω) = [ˆq(ω)]k+1 ˆv0(ω), (7.14 )
where
( iωΔx ) ˆq(ω ) = 1 + a λ e − 1 (7.15 )
is called the amplification factor and
λ = Δt-- (7.16 ) Δx
the Courant–Friedrich–Levy (CFL) factor.

Using Parseval’s identity, we find

∑ ∥v(tk,⋅)∥2 = |ˆq(ω)|2k|ˆv0(ω )|2, (7.17 ) ω∈ℤ
and therefore, we see that the inequality (7.3View Equation) can only hold for all k if
|ˆq(ω)| ≤ 1 for all ω ∈ ℤ. (7.18 )
For a > 0, this is the case if and only if the CFL factor satisfies
1 0 < λ ≤ --, (7.19 ) a
and in this case the well-posedness estimate (7.3View Equation) holds with Kd = 1 and αd = 0. The upper bound in condition (7.19View Equation) for this example is known as the CFL limit, and (7.18View Equation) as the von Neumann condition. If a = 0, ˆq(ω ) = 1, while for a < 0 the scheme is unconditionally unstable even though the underlying continuum problem is well posed.

Next we consider a scheme very similar to the previous one, but which turns out to be unconditionally unstable for a ⁄= 0, regardless of the direction of propagation.

Example 34. A centered Euler scheme.
Consider first the semi-discrete approximation to Eq. (7.10View Equation),

d- vj+1-−-vj−1 dtvj = a 2Δx ; (7.20 )
it is easy to check that it is stable for all values of Δx. Next discretize time through an Euler scheme, leading to
vk+j1 − vkj vkj+1 − vkj−1 ----------= a-----------. (7.21 ) Δt 2 Δx
The solution again has the form given by Eq. (7.14View Equation), now with
|ˆq(ω)| = |1 + iaλ sin(ωΔx )| ≥ 1. (7.22 )
At fixed time tk, the norm of the solution to the fully-discrete approximation (7.21View Equation) for arbitrary small initial data with ωΔx ∕∈ π ℤ grows without bound as the timestep decreases.

The semi-discrete centered approximation (7.20View Equation) and the fully-discrete centered Euler scheme (7.21View Equation) constitute the simplest example of an approximation, which is not fully-discrete stable, even though its semi-discrete version is. This is related to the fact that the Euler time integration is not locally stable, as discussed in Section 7.3.2.

The previous two examples were one-step methods, where vk+1 can be computed in terms of vk. The following is an example of a two-step method.

Example 35. Leap-frog.
A way to stabilize the centered Euler scheme is by approximating the time derivative by a centered difference instead of a forward, one-sided operator:

( ) vkj+1= vkj−1 + aλ vkj+1 − vkj− 1 . (7.23 )
Enlarging the system by introducing
( k ) wk := vj (7.24 ) j vkj−1
it can be cast into the one-step method
( ) ˆwk+1 = Qˆ(ω)wˆk = ˆQ (ω)k+1 ˆw0, with ˆQ (ω) = 2iaλ sin(ω Δx ) 1 . (7.25 ) 1 0
By a similar procedure, a general multi-step method can always be reduced to a one-step one. Therefore, in the stability results below we can assume without loss of generality that the schemes are one-step.

In the above example the amplification matrix ˆQ (ω ) can be diagonalized through a transformation that is uniformly bounded:

( ) ( ) μ+ 0 −1 k μk+ 0 −1 ˆQ(ω ) = ˆT(ω ) 0 μ Tˆ (ω), ˆQ (ω ) = ˆT(ω ) 0 μk Tˆ (ω ), (7.26 ) − −
with μ ± = z ± (1 + z2)1∕2, z := iaλ sin(ω Δx ), and
( μ μ ) ˆT(ω ) = + − . (7.27 ) 1 1
The eigenvalues μ ± are of unit modulus, |μ± | = 1. In addition, the norms of Tˆ(ω) and its inverse are
ˆ ∘ ---------- ˆ− 1 -----1------ |T (ω)| = 2 (1 + |z|), |T (ω )| = ∘ ----------. (7.28 ) 2(1 − |z|)
Therefore, the condition number of Tˆ(ω) can be bounded for all ω:
( )1∕2 ( )1 ∕2 |Tˆ(ω)| ⋅ |ˆT −1(ω)| = 1-+-|z| ≤ 1-+-|a|λ- < ∞, (7.29 ) 1 − |z| 1 − |a|λ
provided that
-1- λ < |a|, (7.30 )
and it follows that the Leap-frog scheme is stable under the condition (7.30View Equation).

The previous examples were explicit methods, where the solution vkj+1 (or wk+j1) can be explicitly computed from the one at the previous timestep, without inverting any matrices.

Example 36. Crank–Nicholson.
Approximating Eq. (7.10View Equation) by

( Δt ) ( Δt ) 1 − a---D0 vkj+1 = 1 + a---D0 vkj, (7.31 ) 2 2
with
-1--- D0vj := 2Δx (vj+1 − vj− 1) , (7.32 )
defines an implicit method. Fourier transform leads to
[ λ ] [ λ ] 1 − ia--sin(ωx ) ˆvkj+1 = 1 + ia -sin(ωx ) ˆvkj. (7.33 ) 2 2
The expressions inside the square brackets on both sides are different from zero and have equal magnitude. As a consequence, the amplification factor in this case satisfies
|qˆ(ω )| = 1 for all ω ∈ ℤ and λ > 0, (7.34 )
and the scheme is unconditionally stable at the expense of having to invert a matrix to advance the solution in time.

Example 37. Iterated Crank–Nicholson.
Approximating the Crank–Nicholson scheme through an iterative scheme with a fixed number of iterations is usually referred to as the Iterated Crank–Nicholson (ICN) method. For Eq. (7.10View Equation) it proceeds as follows [414Jump To The Next Citation Point]:

The resulting discretization is numerically stable for λ ≤ 2∕a and p = 2,3,6,7,10, 11,... iterations, and unconditionally unstable otherwise. In the limit p → ∞ the ICN scheme becomes the implicit, unconditionally-stable Crank–Nicholson scheme of the previous example. For any fixed number of iterations, though, the method is explicit and stability is contingent on the CFL condition λ ≤ 2 ∕a. The method is unconditionally unstable for p = 4,5,8,9,12, 13,... because the limit of the amplification factor approaching one in absolute value [cf. Eq. (7.34View Equation)] as p increases is not monotonic. See [414] for details and [380] for a similar analysis for “theta” schemes.

Similar definitions to the one of Definition 11 are introduced for the IBVP. For simplicity we explicitly discuss the semi-discrete case. In analogy with the definition of a strongly–well-posed IBVP (Definition 9) one has

Definition 12. A semi-discrete approximation to the linearized version of the IBVP (5.1View Equation, 5.2View Equation, 5.3View Equation) is numerically stable if there are discrete norms ∥ ⋅ ∥d at Σ and ∥ ⋅ ∥∂,d at ∂ Σ and constants Kd = Kd (T) and 𝜀d = 𝜀d(T) ≥ 0 such that for high-enough resolution the corresponding approximation v satisfies

∫ t ⌊ ∫t ⌋ 2 2 2⌈ 2 ( 2 2 ) ⌉ ∥v(t,⋅)∥d + 𝜀d ∥v(s,⋅)∥∂,dds ≤ Kd ∥f ∥d + ∥F (s,⋅)∥d + ∥g(s,⋅)∥ ∂,d ds , (7.39 ) 0 0
for all t ∈ [0,T ]. If the constant 𝜀 d can be chosen strictly positive, the problem is called strongly stable.

In addition, the semi-discrete version of Definitions 6 and 7 lead to the concepts of strong stability in the generalized sense and boundary stability, respectively, which we do not write down explicitly here. The definitions for the fully-discrete case are similar, with time integrals such as those in Eq. (7.39View Equation) replaced by discrete sums.


  Go to previous page Go up Go to next page