3.1 Absolute Stability

(1)

3 Stability and Stiffness

3.1 Absolute Stability

Zero-stability is the property of a numerical method for solving IVPs to generate bounded approximations on a fixed interval as h → 0.

In contrast, the term absolute stability is concerned with the behavior of such a scheme for long integration times with fixed stepsize h > 0.

Qualitatively: given an IVP and a (convergent) numerical method, how small must we choose the stepsize in order for the numerical

approximation to have at least the correct “qualitative bahavior“?

Stiff ODEs TU Bergakademie Freiberg 127

Example 1: The IVP

y

⁰

(t) = − sin t, y(0) = 1 (1) has the solution y(t) = cos t. The local discretization error of the explicit Euler method for (1) at a point t is

T (t) = h

2 y

⁰⁰

(t) + O(h

²

) = − h

2 cos t + O(h

²

).

Our standard bound for the global discretization error on an interval [0, t

_end

] yields, for t

_end

= 2,

| e

n

| ≤ t

_end

max

0≤t≤T

| T (t) | ≈ h max

0≤t≤tend

| cos t | = h.

To achieve an error | e | ≤ 10

⁻³

when integrating up to t

_end

= 2, a stepsize of h = 10

⁻³

should suffice, and for N = 2000 we obtain

y

N

= −4.166014e − 01, |y(2) − y

N

| = 4.547667e − 04.

Example 2:

Consider the modification of (1)

y

⁰

(t) = − sin t

+λ(y−cost),

y(0) = 1, (2) with the same solution y(t) = cos t. For λ = − 10, we obtain in this case

y

N

= −4.170721e − 01, |y(2) − y

N

| = 1.611611e − 05.

Attention:

The theoretical bound on the error from pg. 128 is no longer valid!

(2)

Example 3:

Consider the previous example with λ = −2100. In this case one obtains y

N

= − 1.452516e + 76, | y(2) − y

N

| = 1.452516e + 76.

Different stepsizes for this IVP

N h error at t = 2

2000 0.001 1.4e + 76 2050 0.000976 5.8e + 35 2090 0.000957 1.0e + 02 2095 0.000955 4.9e − 03 2100 0.00095 3.2e − 07 2500 0.0008 7.9e − 08 5000 0.0004 4.0e − 08

3.1.1 Runge-Kutta Methods

Solutions

y(t)

of

y⁰

= Ay with constant A ∈

Rⁿ^×ⁿ

converge to

0

(as t → ∞) if Re λ < 0 for all eigenvalues λ of A.

We seek conditions under which a RKM generates approximate solutions with the same asymptotic behavior.

Therefore we apply the considered OSM/RKM to Dahlquist’s test equation

y

⁰

= λy

and analyze the behavior of the numerical approximation.

In general: any OSM applied to the test equation yields approximations y

n+1

= R(ˆ h)y

n

with a stability function R of ˆ h = λh depending on the method.

Therefore, for h fixed, we have lim

n→∞

y

n

= 0 for all y

0

if, and only if,

| R(ˆ h) | < 1.

Definition 3.1.

We define the region of absolute stability R

A

of a OSM by R

A

:= ˆ h ∈

C

: |R(ˆ h)| < 1 .

A OSM is called absolutely stable, if R

A

⊇ { ˆ h : Re (ˆ h) < 0 }.

(3)

For the implicit Euler method, R(ˆ h) = 1

1 − ˆ h so that R

A

=

n

ˆ h : | 1 − ˆ h | > 1

o

For the trapezoidal rule,

R(ˆ h) = 1 +

¹₂

ˆ h

1 −

¹₂

ˆ h so that R

A

= { ˆ h : Re ˆ h < 0 } . Both methods are therefore absolutely stable.

Draw a picture of both stability regions to support this statement.

Applying an m-stage RKM, given by A, b, c from the Butcher tableau, to the test equation y

⁰

= λy yields the linear system

k

j

= λ y

n

+ h

Xm

l=1

α

j,l

, k

l

!

(j = 1, . . . , m),

which can be reformulated as

(I

m

− ˆ hA)

k

= λy

ne.

with ˆ h = λh and

e

= [1, . . . , 1]

^T

. Hence we get the approximation sequence

y

n+1

=

h

1 + ˆ h

b^T

(I

m

− ˆ hA)

⁻¹ei

y

n

and the stability function

R(ˆ h) = 1 + ˆ h

b^T

(I

m

− ˆ hA)

⁻¹e.

Applying the Neumann series we can rewrite the formula for sufficiently small ˆ h as

R(ˆ h) = 1 +

X∞ j=1

h ˆ

^jb^T

A

^j⁻¹e.

If the method has order p, then (cf. pg. 94) R(ˆ h) =

Xp

j=0

1 j! ˆ h

^j

+

X∞ j=p+1

h ˆ

^jb^T

A

^j−1e.

For explicit RKM, A is a strictly lower triangular matrix, s. t. A

^j

= O for j ≥ m. Therefore the series breaks down after m + 1 terms and

R(ˆ h) = 1 +

Xm

j=1

ˆ h

^jb^T

A

^j⁻¹e

=

Xp

j=0

1 j! ˆ h

^j

+

Xm

j=p+1

h ˆ

^jb^T

A

^j⁻¹e.

(4)

Explicit m-stage RKMs of order m (1 ≤ m ≤ 4) all have the stability function

R(ˆ h) =

Xm

j=0

1 j! ˆ h

^j

.

Therefore R

A

is independent of the coefficients for these methods.

Moreover, no explicit RKM has an unbounded region of absolute stability, as R is a polynomial in this case.

Conclusion:

Explicit RKM are never absolutely stable! In practice they can sometimes require unrealistic small stepsizes.

−4 −3 −2 −1 0 1 2

−3

−2

−1 0 1 2 3

Verfahren dritter Ordnung von Heun

−4 −3 −2 −1 0 1 2

−3

−2

−1 0 1 2 3

klassisches Rung−Kutta−Verfahren

Stability regions for Heun’s third order method and the classical Runge-Kutta method.

3.1.2 Linear Multistep Methods

When a LMM is applied to the test equation y

⁰

= λy the resulting numerical approximation satisfies the difference equation

Xr

j=0

α

j

y

n+j

= h

Xr

j=0

β

j

λy

n+j

or

_r

X

j=0

(α

j

− ˆ hβ

j

)y

n+j

= 0, ˆ h = hλ.

Its solutions depend on the roots ζ

1

, . . . , ζ

r

of its characteristic polynomial (cf. pg. 118 ff.)

π(ζ) = π(ζ; ˆ h) :=

Xr

j=0

(α

j

− hβ ˆ

j

)ζ

^j

= ρ(ζ) − ˆ hσ(ζ),

which is called the stability polynomial of the associated LMM.

(5)

A LMM is therefore absolutely stable for a given value of ˆ h if the zeros of the stability polynomials for this value of h ˆ satisfy the conditions for bounded solutions (cf. pg. 122).

Examples for stability polynomials:

Method π(ζ; ˆ h)

Euler (explicit) ζ + (1 + ˆ h) Euler (implicit) (1 − ˆ h)ζ − 1 trapezoidal rule (1 − ˆ h/2)ζ − (1 + ˆ h/2)ζ midpoint rule ζ

²

− 2ˆ hζ − 1 There are no explicit absolutely stable LMMs.

Implicit LMM, which are absolutely convergent, have order of convergence at most 2.

ï20 ï10 0 10 20 30 40

ï25 ï20 ï15 ï10 ï5 0 5 10 15 20 25

Re

Im

BDF1 BDF2 BDF3 BDF4 BDF5 BDF6

BDF methods: The region R

A

of absolute stability of the BDF methods (r = 1, . . . , 6) is the exterior of the corresponding closed curve.

3.2 Stiff Differential Equations

Stiffness is a property of differential equations with strong implications for their practical solution using numerical methods.

Unfortunately, a precise mathematical definition of stiffness covering all occurrences of this phenomenon has not been found.

We therefore resort to a series of examples of stiff differential equations.

(6)

Example 1. The two IVPs

y⁰

=

−2 1 1 −2

y

+

2 sin t 2(cos t − sin t)

,

y(0) =

2

3 , (IVP

1

)

y⁰

=

− 2 1 998 − 999

y

+

2 sin t 999(cos t − sin t)

,

y(0) =

2

3 , (IVP

2

) possess the same solution

y(t) = exp(−t)

2

2 + sin t

cos t

.

0 2 4 6 8 10

−1

−0.5 0 0.5 1 1.5 2 2.5 3

y₁ y₂

We solve both on t ∈ [0, 10] using

Matlab’s ODE solverode45¹

, specifying a relative tolerance of 0.01.

0 2 4 6 8 10

−1

−0.5 0 0.5 1 1.5 2 2.5 3

System 1, ode45: 64 steps.

y₁ y₂

0 2 4 6 8 10

−1

−0.5 0 0.5 1 1.5 2 2.5 3

System 2, ode45: 12044 steps.

y₁ y₂

1based on the embedded Runge-Kutta pair of Dormand and Prince of order 4 and 5 for adaptive stepsize selection

However, with this method solving (IVP

2

) is roughly 200 times more costly than for (IVP

1

):

IVP h

min

h

max

h

_∅

# steps

(IVP

1

) 1.19e − 1 2.03e − 1 1.56e − 1 64 (IVP

2

) 6.03e − 4 2.88e − 3 8.30e − 3

12044

Applying instead

Matlab’s ODE solverode23tb²

, again with a relative error tolerance of 0.01, one obtains

IVP h

min

h

max

h

_∅

# steps

(IVP

1

) 1.72e − 1 5.81e − 1 4.00e − 1 25 (IVP

2

) 2.13e − 1 7.04e − 1 4.55e − 1 22 This time both problems are solved without difficulty.

2based on an implicit 2-stage RKM with the trapezoidal rule and the BDF2 method as its first and second stages, respectively

(7)

0 2 4 6 8 10

−1

−0.5 0 0.5 1 1.5 2 2.5 3

System 1, ode23tb: 25 steps.

y₁ y2

0 2 4 6 8 10

−1

−0.5 0 0.5 1 1.5 2 2.5 3

System 2, ode23tb: 22 steps.

y₁ y2

Results of

Matlab’s solverode23tb

for both systems.

The reason for the different behavior of the two IVPs becomes apparent when we look at the general solutions of each associated ODE system.

For (IVP

1

) this is

y(t) =

κ

1

exp( − t) 1

1 + κ

2

exp( − 3t) 1

− 1

+ sin t

cos t

,

whereas for (IVP

2

)

y(t) =

κ

1

exp( − t)

1 1

+ κ

2

exp( − 1000t)

1 −998

+ sin t

cos t

. Both are inhomogeneous linear systems with constant coefficients, and the eigenvalues of the coefficient matrix are

λ

1

= − 1, λ

2

= − 3 for (IVP

1

), λ

1

= − 1, λ

2

= − 1000 for (IVP

2

).

In both cases imposing the initial condition leads to κ

1

= 2 and κ

2

= 0 (!).

In both ODE solutions we have steady-state and (decaying) transient components

[sin t, cos t]

^T

and exp(λ

1

t)

a1

+ exp(λ

2

t)

a2

. Fastest decaying transient (as t → ∞) determined by

λ

2

= − 3 for (IVP

1

) and λ

2

= − 1000 for (IVP

2

).

Even though neither is present in the solutions (due to κ

2

= 0), this rate determines the minimal required stepsize.

Changing the initial data to

y(0) = [0,

1]

^T

in (IVP

2

) results in

κ

1

= κ

2

= 0, i.e., both transient components are then completely

invisible in the resulting exact solution

y(t) = [sin

t, cos t]

^T

.

Even then, however,

ode45

would take excessively small stepsizes.

(8)

For

ode45, the region

R

A

of absolute stability satisfies R

A

∩

R

∼ [ − 3, 0] (stability interval)

To ensure hλ ∈ R

A

for both eigenvalues of (IVP

1

) it is sufficient to require h < 1. The observed average stepsize of h

_∅

∼ 0.156 is a consequence of the requested accuracy.

By contrast, for (IVP

2

) ensuring λh ∈ R

A

for all λ necessitates a stepsize of h < 0.003. In this case the observed h

_∅

∼ 0.008 results from the stability requirement.

Finally,

ode23t

is an absolutely stable method, and therefore hλ ∈ R

A

for any h > 0 for both (IVP

1

) and (IVP

2

) when this method is used.

Definition 3.2.

For a system of ODEs

y⁰

= Ay +

b(t)

the quantity

λ

max

∈Λ(A)

|Re λ | / min

λ∈Λ(A)

|Re λ | is called its stiffness ratio.

Statement 1:

A linear system of ODEs with constant coefficients is called stiff, when all its eigenvalues have negative real part and its stiffness ratio is large.

To see that Statement 1 is not a useful definition, consider yet another IVP with the same solution as (IVP

1

) and (IVP

2

):

y⁰=

−2 1

−1.999 0.999

y+

2 sint

−0.999(cost−sint)

, y(0) = 2

3

. (IVP3)

The eigenvalues of its coefficient matrix are λ

1

= − 1 and λ

2

= − 0.001, giving the same stiffness ratio of 1000 as for (IVP

2

).

However,

ode45

happily solves (IVP

3

) with reasonable stepsizes :

IVP h

min

h

max

h

_∅

# steps

(IVP

3

) 1.19e − 1 2.50e − 1 2.27e − 1 44 (ode45)

(IVP

3

) 1.25e − 1 4.92e − 1 4.55e − 1 22 (ode23tb)

The stiffness ratio alone can therefore not always determine whether or

not a system is stiff.

(9)

This might suggest modifying Statement 1 to Statement 2:

A system is stiff if max |Re λ | is large, e.g., if max |Re λ | 1.

However, the change of variables t 7→ 0.001 t transforms (IVP

2

) to an equivalent problem with identical stiffness properties, but for which max |Re λ | = 1.

Consequently one settles for pragmatic “definitions“ of stiffness which include

A system is stiff if stability requirements – and not accuracy requirements – determine the stepsize.

A system is stiff, if certain components decay much faster than others.

Stiff equations are problems for which explicit methods don’t work.

[Hairer and Wanner, Vol. II, p. 2]

Example 2. The ODE

y

⁰

(t) = λ[y(t) − g(t)] + g

⁰

(t) (∗) has the general solution

y(t) = γ exp(λt) + g(t), γ ∈

R

,

which, for Re λ < 0, consists of a decaying transient γ exp(λt) and a steady-state component g(t).

Setting λ = − 10, g(t) = arctan t and initial condition y(0) = 0, we approximate the exact solution y(t) = arctan t for t ∈ [0, 5] using

explicit Euler (R

A

= { ˆ h : | ˆ h + 1| < 1}) and implicit Euler (R

A

=

C

\ { ˆ h : | ˆ h − 1 | > 1 }).

−5 0 5

−3

−2

−1 0 1 2 3

Allgemeine Loesung von (*)

General solution

(10)

−1 0 1 2

h=0.3 (explizites Euler−Verfahren)

−1 0 1 2

0 1 2 3 4 5

−1 0 1 2

Explicit Euler,h= 0.1,0.2,0.3

−1 0 1 2

h=0.3 (implizites Euler−Verfahren)

−1 0 1 2

h=0.5 (implizites Euler−Verfahren)

0 1 2 3 4 5

−1 0 1 2

h=1 (implizites Euler−Verfahren)

Implicit Euler,h= 0.3,0.5,1