5. Matrix Stability Analysis

A general one–step algorithm can be represented in the form

∑_{n=1}^{N−1} A_{mn}U_{n}^{j+1} = ∑_{n=1}^{N−1} B_{mn}U_{n}^{j} + c_{m}^{j}  for 1 ⩽ m ⩽ N − 1, 1 ⩽ j ⩽ J .  (5.1)

A and B are (N − 1) × (N − 1) matrices and we write the truncation errors R_{n}^{j} and the solution errors z_{n}^{j} ≡ u_{n}^{j} − U_{n}^{j} in the form of (N − 1)–dimensional vectors, so that \underline{Z}^{j} = (z_{1}^{j}, z_{2}^{j}, …, z_{N−1}^{j}) and ε\underline{d}^{j} = (R_{1}^{j}, R_{2}^{j}, …, R_{N−1}^{j}) where ε is small. Typically ε = O(k^{2} + kh^{2}). Then the error–vector obeys the equation

A\underline{Z}^{j+1} = B\underline{Z}^{j} + ε\underline{d}^{j} .

Assuming the matrix inverse A^{−1} exists (else the entire method collapses as we will be unable to find U_{n}^{j+1} from U_{n}^{j}) we can write

\underline{Z}^{j+1} = (A^{−1}B) \underline{Z}^{j} + εA^{−1}\underline{d}^{j} = P \underline{Z}^{j} + ε\underline{f}^{j} say.

The initial error, \underline{Z}^{0} = 0. Applying the above relation iteratively leads to

\underline{Z}^{j} = ε [(P)^{j−1}\underline{f}^{0} + (P)^{j−2}\underline{f}^{1} + … + P \underline{f}^{j−2} + \underline{f}^{j−1}] .  (5.2)

The matrix P describes the propagation of errors from one time–step to the next. Now there are j terms on the RHS of (5.2), and jε ⩽ Jε = O(kT) typically. For the method (5.1) to be stable, we must have \underline{Z}^{J} → 0 as J → ∞ and k → 0. We will be able to show this to be the case, provided repeated pre–multiplication by the matrix P does not lead to an increase in the errors. We therefore want to consider the vector (P)^{j}\underline{x} where \underline{x} is any vector and j is large.

Suppose the matrix P has eigenvectors \underline{e}_{m} and eigenvalues λ_{m}. Then

P \underline{e}_{m} = λ_{m}\underline{e}_{m} ,  (P)^{2}\underline{e}_{m} = λ_{m}^{2}\underline{e}_{m}  and (P)^{j}\underline{e}_{m} = (λ_{m})^{j}\underline{e}_{m} .

Thus, if one of the eigenvalues has modulus greater than one it is possible for |(P)^{j}\underline{x}| to increase without limit. If this happens the method (5.1) will be unstable (although see the subtlety overleaf.) Conversely, suppose that |λ_{m}| ⩽ 1 for all m. Then it is easy to show that when the eigenvectors of P are complete so that we may expand \underline{x} = ∑ a_{m}\underline{e}_{m}, then

|(P)^{j}\underline{x}| = | ∑_{m=1}^{N−1} a_{m}(P)^{j}\underline{e}_{m}| = |∑ a_{m}(λ_{m})^{j}\underline{e}_{m}| ⩽ ∑ |a_{m}||λ_{m}|^{j} ⩽ ∑ |a_{m}| .

This is bounded and so the error vector \underline{Z}^{j} in (5.2) remains small. Indeed, if all the eigenvalues are strictly less than one in modulus, then the propagated errors decrease as j increases. The solution error is then dominated by the local truncation error.

We have shown that the method (5.1) is stable provided max |λ_{m}| ⩽ 1 and the eigenvectors are complete. In fact this result is true for incomplete sets of eigenvectors also, but is harder to prove. The maximum value of |λ_{m}| is often called the spectral radius of P, and sometimes written as ρ(P). When calculating the eigenvalues, it is usually best to consider the equation |B − λA| = 0 rather than calculating P.

A subtlety: Bounded growth. Even if the spectral radius is a little greater than one, the method may still be stable, provided the errors grow in a bounded fashion. Consider a time interval 0 < t < T = Jk. The errors over this period are bounded if for all m and some constant G,

|λ_{m}|^{J} ⩽ G ⟺ |λ_{m}| ⩽ G^{1/J} = e^{(ln G)k/T} ≈ 1 + k \frac{ln G}{T} + O(k^{2}) as k → 0 .

The stability condition allowing for bounded growth is

|λ_{m}| ⩽ 1 + O(k) .

Example: We consider the θ-method (4.2) for the equation u_{t} = u_{xx} for which

A = Tp [1 + 2rθ, −rθ, −rθ] and B = Tp [1 − 2r(1 − θ), (1 − θ)r, (1 − θ)r] ,

where we define the tri–diagonal Toeplitz matrix

Tp [a, b, c] = ( \begin{matrix}a & b & 0 & ... & 0 \\ c & a & b & ... & 0 \\ 0 & c & a & ... & ... \\ ... & ... & ... & ... & b \\ 0 & 0 & ... & c & a \end{matrix} ) .

We may express A = I − rθC and B = I + r(1 − θ)C where

I = Tp [1, 0, 0]  and C = Tp [−2, 1, 1] .  (5.3)

Then

B − λA = (1 − λ)I + [r(1 − θ) + λrθ]C and so
|B − λA| = 0 ⟺ |C − σI| = 0  where σ = \frac{λ − 1}{r(1 − θ) + λrθ} .  (5.4)

We have thus expressed the eigenvalues λ_{m} of P in terms of the eigenvalues σ_{m} of C. Regrouping (5.4) we have

λ_{m} = \frac{1 + r(1 − θ)σ_{m}}{1 − rθσ_{m}} .

Now it can be shown that the eigenvalues of Tp [a, b, c] are a + 2\sqrt{bc} cos(mπ/N) for 1 ⩽ m ⩽ N − 1. The eigenvalues of C are therefore

σ_{m} = −4 sin^{2} (\frac{mπh}{2})  for m = 1, …, N − 1  and h = \frac{1}{N} .

Substituting these values into (5.4) and using the condition |λ_{m}| ⩽ 1, leads to a stability condition identical to (4.9), which we obtained before from the exact solution, namely

|\frac{1 − 4r(1 − θ) sin^{2}(\frac{1}{2}mπh)}{1 + 4rθ sin^{2}(\frac{1}{2}mπh)}| ⩽ 1 or 2r(1 − 2θ) sin^{2} (\frac{mπh}{2}) ⩽ 1 .

In particular, we note that the implicit methods with θ ⩾ \frac{1}{2} are unconditionally stable.

Two-step Methods

Using a similar notation, the general two-step method can be written

A\underline{Z}^{j+1} = B\underline{Z}^{j} + C\underline{Z}^{j−1} + \underline{d}^{j}  (5.5)

where \underline{Z}^{j} and \underline{d}^{j} are (N − 1)–vectors and A, B, and C are (N − 1) × (N − 1) matrices. We introduce a (2N − 2)–vector \underline{W}^{j}, where

\underline{W}^{j} = ( \begin{matrix}\underline{Z} & j \\ − − − \\ \underline{Z} & j−1 \end{matrix} )  (5.6)

Then (5.5) reduces to \underline{W}^{j+1} = P \underline{W}^{j} + \underline{q}, where

P = ( \begin{matrix}A^{−1}B & | & A^{−1}C \\ − − − & − + − & − − − \\ I & | & O \end{matrix} )  and \underline{q} = ( \begin{matrix}A^{−1}\underline{d} \\ − − − \\ O \end{matrix} )  (5.7)

As before, we want ρ(P) ≡ |max e-value of P | ⩽ 1. The eigenvalues and eigenvectors are given by the equation

( \begin{matrix}A^{−1}B − λI & | & A^{−1}C \\ − − − − − & − + − & − − − \\ I & | & −λI \end{matrix} ) ( \begin{matrix}\underline{Z} & 1 \\ −− \\ \underline{Z} & 2 \end{matrix} ) = ( \begin{matrix}O \\ −− \\ O \end{matrix} ) .

Now the bottom half of these equations clearly imply \underline{Z}_{1} = λ\underline{Z}_{2}. Thus the top half require

(λA^{−1}B − λ^{2}I + A^{−1}C) \underline{Z}_{2} = 0 .

The requisite eigenvalue relation is therefore

|λ^{2}A − λB − C| = 0 .  (5.8)

The stability condition is that all (2N − 2) solutions of (5.8) must be such that |λ| ⩽ 1.

The matrix method is a general, powerful and precise method of determining stability of a finite difference scheme. Unlike the Fourier method, it takes into account the boundary conditions. In general, the eigenvalues may be harder to find than in our example, but it is often possible to determine whether or not they all lie in the circle |λ| ⩽ 1 without calculating them explicitly.