from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

6.8.Zero-stability of multistep methods¶

For one-step methods such as Runge–Kutta, Theorem 6.2.1 guarantees that the method converges and that the global error is of the same order as the local truncation error. For multistep methods, however, a new wrinkle is introduced.

Example 6.8.1 (Instability)

It is straightforward to check that the two-step method LIAF, defined by

u_{i+1} = -4u_i + 5u_{i-1} + h(4f_i + 2f_{i-1}),

(6.8.1)

is third-order accurate. Let’s apply it to the ridiculously simple IVP $u'=u$ , $u(0)=1$ , whose solution is $e^t$ .

We’ll measure the error at the time $t=1$ .

du_dt = lambda t, u: u
u_exact = exp
a, b = (0.0, 1.0)

def LIAF(du_dt, tspan, u0, n):
    a, b = tspan
    h = (b - a) / n
    t = linspace(a, b, n+1)
    u = np.tile(np.array(u0), (n+1, 1))
    u[1] = u_exact(t[1])    # use an exact starting value
    f = copy(u)
    f[0] = du_dt(t[0], u[0])
    for i in range(n):
        f[i] = du_dt(t[i], u[i])
        u[i + 1] = -4 * u[i] + 5 * u[i-1] + h * (4 * f[i] + 2 * f[i-1])

    return t, u.T

n = [5, 10, 20, 40, 60]
results = PrettyTable(["n", "error"])
for j in range(5):
    t, u = LIAF(du_dt, [a, b], [1.0], n[j])
    err = abs(u_exact(b) - u[0, -1])
    results.add_row([n[j], err])
print(results)

+----+------------------------+
| n  |         error          |
+----+------------------------+
| 5  |    293.502838171541    |
| 10 |   452905.58324991807   |
| 20 |   2196838145485.221    |
| 40 | 1.0437865522698013e+26 |
| 60 | 6.628042953187168e+39  |
+----+------------------------+

There is no convergence in sight! A graph of the last numerical attempt yields a clue:

semilogy(t, abs(u[0]), "-o")
xlabel("$t$"), ylabel("$|u|$")
title(("LIAF solution"));

It’s clear that the solution is growing exponentially in time.

The source of the exponential growth in Example 6.8.1 is not hard to identify. Recall that we can rewrite (6.8.1) as $\rho(\mathcal{Z})u_{i-1}=h \sigma(\mathcal{Z})u_{i-1}$ using the forward shift operator $\mathcal{Z}$ :

(\mathcal{Z}^2 + 4\mathcal{Z} - 5) u_{i-1} = h(4\mathcal{Z} + 2)f_{i-1}.

(6.8.2)

(See (6.6.7), using $k=2$ here.) Next, suppose that $h$ is negligible in (6.8.2). Then the numerical solution of LIAF is roughly defined by

(\mathcal{Z}^2 + 4\mathcal{Z} - 5) u_{i-1} = 0.

(6.8.3)

The graph in Example 6.8.1 strongly suggests that for small $h$ , $|u_i|\approx c \alpha^i$ for some $\alpha>1$ as $m$ gets large. So we are motivated to try defining

u_i = c z^i

(6.8.4)

for all $i$ and see if we can prove that it is an exact solution. The beauty of this choice is that for all $i$ ,

\mathcal{Z} u_i = u_{i+1} = z u_i.

(6.8.5)

Hence, (6.8.3) becomes

z^2 + 4z - 5 = 0.

(6.8.6)

Therefore, as $h\to 0$ , the two roots of $z^2+4z-5$ will each correspond to an approximate solution in the form (6.8.4) of the LIAF method. These roots are $z=1$ and $z=-5$ , and the growth curve at the end of Example 6.8.1 is approximately $|(-5)^i|$ .

6.8.1Zero-stability¶

Here is the crucial property that LIAF lacks.

Without zero-stability, any truncation or roundoff error will get exponentially amplified and eventually overwhelm convergence to the exact solution.

The following theorem concisely summarizes when we can expect zero-stability.

Proof

(Partial proof, omitting repeated roots.) As explained above, the values produced by the numerical method approach solutions of the difference equation $\rho(\mathcal{Z})u_{i-k+1}=0$ . We consider only the case where the roots $r_1,\ldots,r_k$ of $\rho(z)$ are simple. Then $u_i=(r_j)^i$ is a solution of $\rho(\mathcal{Z})u_i=0$ for each $j=1,\ldots,k$ . By linearity,

u_i = c_1 (r_1)^i + c_2 (r_2)^i + \cdots + c_k (r_k)^i

(6.8.7)

is a solution for any values of $c_1,\ldots,c_k$ . These constants are determined uniquely by the starting values $u_0,\ldots,u_{k-1}$ (we omit the proof). Now, if all the roots satisfy $|r_j|\le 1$ , then

|u_i| \le \sum_{j=1}^k |c_j| |r_j|^i \le \sum_{j=1}^k |c_j|,

(6.8.8)

independently of $h$ and $i$ . This proves zero-stability. Conversely, if some $|r_j|>1$ , then $|u_i|$ cannot be bounded above by a constant independent of $i$ . Since $b=t_i$ , $i\to\infty$ at $t=b$ as $h\to 0$ , so zero-stability cannot hold.

A multiple root of $\rho$ introduces a modification of (6.8.4) that is considered in Exercise 6.8.4.

6.8.2Dahlquist theorems¶

It turns out that lacking zero-stability is the only thing that can go wrong for a consistent multistep method.

The Dahlquist equivalence theorem is one of the most important and celebrated in the history of numerical analysis. It can be proved more precisely that a zero-stable, consistent method is convergent in the same sense as Theorem 6.2.1, with the error between numerical and exact solutions being of the same order as the local truncation error, for a wide class of problems.

You may have noticed that the Adams and BD formulas use only about half of the available data from the past $k$ steps, i.e., they have many possible coefficients set to zero. For instance, a $k$ -step AB method uses only the $f_j$ -values and has order $k$ . The order could be made higher by also using $u_j$ -values, like the LIAF method does for $k=2$ . Also like the LIAF method, however, such attempts are doomed by instability.

The lesson of Theorem 6.8.3 is that accuracy is not the only important feature, and trying to optimize for it leads to failure. New lessons on the same theme appear in Absolute stability.

Zero-stability of multistep methods

6.8.Zero-stability of multistep methods¶

6.8.1Zero-stability¶

6.8.2Dahlquist theorems¶

6.8.3Exercises¶