1 Introduction

We have seen that time dependent problems are approximated by a vector of values, ${\rm\bf U}^n$ at each time step and further the form of the iteration for these values was

$${\rm\bf B}_1{\rm\bf U}^{n+1}={\rm\bf B}_0{\rm\bf U}^n+{\rm\bf f}^n,$$ (96)

where ${\rm\bf B}_0,{\rm\bf B}_1$ are matrices. If ${\rm\bf B}_1={\rm\bf I}$ we called the method explicit, otherwise we referred to an implicit method. If the PDE involved two or more space dimensions, it can still be dealt with under this framework, elements of vector ${\rm\bf U}$ have just a more complicated mapping onto the physical mesh. Thus in principle we have iterations of the form

$${\rm\bf U}^{n+1}={\rm\bf B}_1^{-1}{\rm\bf B}_0{\rm\bf U}^n+{\rm\bf B}_1^{-1}{\rm\bf f}^n,$$ (97)

to calculate the approximation time step by time step. Thus we need first to understand some things about iterated sequences.

Suppose we have the sequence

$$2y_{n+2}-5y_{n+1}+2y_n=0,$$ (98)

with $y_0=a$, $y_1=b$. We look for a solution $y_n\sim\lambda^n$ so that $2\lambda^2-5\lambda+2=0,.$ The roots are $\lambda=2, {1\over 2}$ and the general solution is $y_n=A2^{-n}+B2^n$. Using the initial conditions

$$y_n={{4a-2b}\over 3}2^{-n}+{{2b-a}\over 3}2^n.$$ (99)

Now suppose we choose the inital values $a,b$ such that $2b-a=0$, then $y_n=a2^{-n}\to 0,\ {\rm as}\ n\to\infty,$ and it might appear that we have a sensible computation to determine $y_n$. However suppose, because of finite machine accuracy, that $y_0=a\pm\epsilon$, $y_1=b\pm\epsilon$, then

$$y_n=(a\pm {2\over 3}\epsilon)2^{-n}\pm\epsilon 2^n,$$ (100)

and for example, if $\epsilon =10^{-6}$ (typical single precison accuracy) then $\epsilon 2^n\sim 1$ after only $n\sim 20$ iterations! Hence we could not carry out any practical computations using this iteration without the calulations being swamped by the growing solution even though the `solution' should have no component of this mode.

Consider now a matrix example,

$${\rm\bf U}^{n+1}={\rm\bf A}{\rm\bf U}^n,$$ (101)

where

$${\rm\bf A}=\left(\begin{array}{r r} 2&1\\ 1&1\\ \end{array}\right).$$ (102)

We look for a solution where $U^n\sim\lambda{\rm\bf z}$ so that ${\rm\bf A}{\rm\bf z}=\lambda{\rm\bf z}$, and $\lambda$ is an eigenvalue of the iteration matrix ${\rm\bf A}$. In this case the eigenvalues are $\lambda_1=(2-\sqrt{5})/2<1$, $\lambda_2=(3+\sqrt{5})/2>1,$ and if we denote the eigenvectors by ${\rm\bf z}_1,{\rm\bf z}_2$, the general solution is

$${\rm\bf U}^n=A\lambda_1^n{\rm\bf z}_1+B\lambda_2^n{\rm\bf z}_2,$$ (103)

where $A$ and $B$ are constants. Again, if we choose the inital data so that $B$ is zero, machine rounding errors will still produce a growing component which will swamp the calculation very quickly.

Thus one difficulty with calculating iterated sequences is that growing modes lead to rapid breakdown of a calculation and this `instability' may not be a part of the original physical system, but introduced as an artefact of the discretisation of the continuous system. Often stability is a more challenging practical problem than accuracy.