## Single Step Methods

A significant source of error when using difference-equation methods to generate approximations to the solution of initial values problems is incurred due to the round off error associated with finite digit arithmetic or inconsistencies in the initial conditions. A difference-equation method may be defined as *stable* if small changes or perturbations in the initial conditions or difference-equation produce correspondingly small changes in the subsequent approximations. The following definitions are useful in analysing the stability of difference-equation methods.

A one step difference-equation method with local truncation error \(\tau_i(h)\) at the \(i^{th}\) step is said to be *consistent* with the differential equation it approximates if:

It should be noted that this definition is a local definition and assumes that the approximation \(\omega_{i-1}\) and the exact solution \(y(t_{i-1})\) are identical.

A multistep difference-equation method with local truncation error \(\tau_i(h)\) at the \(i^{th}\) step is said to be *consistent* with the differential equation it approximates if:

where \(\alpha_i\) denotes the starting values supplied to the multistep method and \(y(t_i)\) denotes the exact value of the solution to the differential equation. The above definition implies that a multistep method will not be consistent unless the one step method generating the starting values is also consistent.

A one step or multistep difference-equation method is said to be *convergent* with respect to the differential equation it approximates if:

where \(y(t_i)\) denotes the exact value of the solution to the differential equation and \(\omega_i\) denotes the approximation obtained from the difference-equation method at the \(i^{th}\) step.

The stability of a one step method is determined by a Lipschitz condition. Suppose the initial value problem:

$$ \frac{dy}{dt} = f(t,y), \quad a \leq t \leq b, \quad y(a) = y_0 \\ \label{eq:stb4} $$is approximated by a one step difference method of the form:

\begin{equation} \begin{aligned} \omega_0 & = \alpha \\ \omega_{i+1} & = \omega_i + hT^n(t_i,\omega_i), \quad i=0,\dots,N-1 \\ \end{aligned} \label{eq:stb5} \end{equation}Suppose also that \(T^n(t_i,\omega_i)\) is continuous and satisfies a Lipschitz condition in the variable \(\omega\) on the set \(D=\{(t,\omega) | a \leq t \leq b,-\infty \lt \omega \lt \infty\}\) then the method is stable.

## Multistep Methods

Determining the stability of a multistep method is more involved. Associated with a multistep difference-equation method given by:

$$ \omega_0 = \alpha_0 \quad \omega_1 = \alpha_1, \dots, \omega_{m-1} = \alpha_{m-1} \\ \omega_{i+1} = a_{m-1}\omega_i + a_{m-2}\omega_{i-1} + \dots + a_0\omega_{i+1-m} + hF(t_i,h,\omega_{i+1},\omega_i,\dots,\omega_{i+1-m}) \\ \label{eq:stb6} $$is a polynomial termed the characteristic polynomial of the method, given by:

$$ P(\lambda) = \lambda^m - a_{m-1}\lambda^{m-1} - a_{m-2}\lambda^{m-2} - \dots - a_1\lambda - a_0 \\ \label{eq:stb7} $$The stability of a multistep method with respect to round off error is dictated by the magnitudes of the zeros of the characteristic polynomial. Consider applying the multistep method to the initial value problem:

$$ \frac{dy}{dt} = 0, \quad a \leq t \leq b, \quad y(a) = \alpha, \quad \alpha \ne 0 \\ \label{eq:stb8} $$This problem has the exact solution \(y(t) = \alpha\). By examining the equations of the multistep difference-equation method, the method will, in theory, produce the exact solution \(\omega_n = \alpha\) for all \(n\) assuming that no round off error occurs during any step. Since \(f(t,y) = 0\), it is assumed that \(F(t_i,h,\omega_{i+1},\omega_i,\dots,\omega_{i+1-m})=0\). The multistep difference equation method therefore becomes:

$$ \omega_{i+1} = a_{m-1}\omega_i + a_{m-2}\omega_{i-1} + \dots + a_0\omega_{i+1-m} \\ \label{eq:stb9} $$Since \(\lambda\) is one of the zeros of the characteristic polynomial associated with the multistep difference-equation. Then \(\omega_n = \lambda^n\) for each \(n\) is a solution to the characteristic polynomial since:

$$ \lambda^{i+1} - a_{m-1}\lambda^{i} - a_{m-2}\lambda^{i-1} - \dots - a_0\lambda^{i+1-m} = \lambda^{i+1-m}\left[\lambda^m - a_{m-1}\lambda^{m-1} - \dots - a_0 \right] = 0\\ \label{eq:stb10} $$If \(\lambda_1 \dots \lambda_m\) are distinct roots of the characteristic polynomial, then the general solution to the difference-equation may be written as:

$$ \omega_n = \sum_{i=1}^{m} c_i {\lambda_i}^n \\ \label{eq:stb11} $$for some unique collection of constants \(c_1, \dots, c_m\). If the roots are not distinct, the then above summation must be replaced with an appropriate expression but the result does not affect the following analysis.

Since the exact solution to the initial value problem is \(y(t) = \alpha\), the choice \(\omega_n = \alpha\) for all \(n\) is a solution to the associated difference-equation:

$$ 0 = \alpha - \alpha a_{m-1} - \alpha a_{m-2} - \dots - \alpha a_{0} = \alpha \left[1 - a_{m-1} - a_{m-2} - \dots - a_0 \right] \\ \label{eq:stb12} $$Therefore, if \(a_{m-1} = 1\) then \(a_{m-2}, \dots, a_0 = 0\) then characteristic equation may be written as:

$$ \lambda^{i+1} - \lambda^i = \lambda^i(\lambda - 1) = 0\\ \label{eq:stb13} $$and therefore \(\lambda = 1\) is one of the zeros of the characteristic polynomial.

If all calculations were exact, then the constants \(c_1, \dots, c_m\) of the general solution to the difference-equation would therefore be zero with the exception of \(c_1 = \alpha\). In practice, \(c_2, \dots, c_m\) would not equal zero due to round off error. The round off error will grow exponentially unless \(|\lambda_i| \leq 1\) for each of the roots \(\lambda_2, \dots, \lambda_m\). This leads to the definition of the *root condition* which may be stated as:

If \(|\lambda_i| \leq 1\) for each \(i=1 \dots m\) and all roots with absolute value equal to one are simple roots, then the difference method is said to satisfy the root condition. In addition:

- Methods which satisfy the root condition and have \(\lambda = 1\) as the only root of the characteristic equation of magnitude equal to one are called strongly stable.
- Methods that satisfy the root condition and have more than one distinct root with magnitude equal to one are called weakly stable.
- Methods that do not satisfy the root condition are called unstable.

## Stability Analysis Due to Small Perturbations in Initial Values

The stability characteristics for the initial value problem discussed in this section determine the stability characteristics for initial value problems when \(f(t,y) \ne 0\). This is because the solution to the initial value problem discussed is embedded in the solution to any initial value problem.

At every time step, a consistent method incurs a small error. A numerical method is said to be stable if the error propagates between successive time steps without growing too readily. Consider applying a numerical method to the perturbed problem:

$$ \frac{dy}{dt} = f(t,y), \quad a \leq t \leq b, \quad y(a) = \alpha + \epsilon \\ \label{eq:stb14} $$Using a Taylor series expansion about the point \((t_0,y_0)\) the function \(f(t,y)\) can be approximated by:

$$ \frac{dy}{dt} = f(t,y(t)) \approx f(t_0,y_0) + (t-t_0) \left.\frac{\partial f}{\partial t}\right|_{(t_0,y_0)} + (y(t)-y_0)\left.\frac{\partial f}{\partial y}\right|_{(t_0,y_0)} \\ \label{eq:stb15} $$Defining \(\overline{y}(t)=y(t)-y_0\) and noting that \(\frac{d\overline{y}}{dt} = \frac{dy}{dy}\) since \(y_0\) is constant, then:

\begin{equation} \begin{aligned} \frac{d \overline{y}}{dt} & = \overline{y}(t)\left.\frac{\partial f}{\partial y}\right|_{(t_0,y_0)} + f(t_0,y_0) + (t-t_0)\left.\frac{\partial f}{\partial t}\right|_{(t_0,y_0)}\\ \frac{d \overline{y}}{dt} & = \lambda \overline{y}(t) + g(t) \\ \end{aligned} \label{eq:stb16} \end{equation}where

$$ \lambda = \left.\frac{\partial f}{\partial y}\right|_{(t_0,y_0)}\\ \label{eq:stb17} $$ $$ g(t) = f(t_0,y_0) + (t-t_0)\left.\frac{\partial f}{\partial t}\right|_{(t_0,y_0)} \\ \label{eq:stb18} $$for small \(t-t_0\). Stability analysis may therefore be restricted to equations of the form \(\frac{dy}{dt} = \lambda y(t) + g(t)\).

Let \(\epsilon_0 \dots \epsilon_{m-1}\) be small perturbations in the initial values supplied to a multistep method, and let:

\begin{equation} \begin{aligned} \omega_0 &= \alpha_0 \dots \omega_{m-1} = \alpha_{m-1} \\ \omega_{i+1} &= \sum_{k=0}^{m-1} a_{m-1-k} \omega_{i-k} + h \sum_{k=-1}^{m-1} b_{m-1-k} f(t_{i-k}, \omega_{i-k}) \\ \end{aligned} \label{eq:stb19} \end{equation}and

\begin{equation} \begin{aligned} \upsilon_0 &= \alpha_0 + \epsilon_0 \dots \upsilon_{m-1} = \alpha_{m-1} + \epsilon_{m-1} \\ \upsilon_{i+1} &= \sum_{k=0}^{m-1} a_{m-1-k} \upsilon_{i-k} + h \sum_{k=-1}^{m-1} b_{m-1-k} f(t_{i-k}, \upsilon_{i-k}) \\ \end{aligned} \label{eq:stb20} \end{equation}where \(f(t,y(t)) = \lambda y(t) + g(t)\). The difference due to the perturbation is:

\begin{equation} \begin{aligned} \upsilon_{i+1} - \omega_{i+1} &= \sum_{k=0}^{m-1} a_{m-1-k} (\upsilon_{i-k} - \omega_{i-k}) + h \sum_{k=-1}^{m-1} b_{m-1-k} \left[(\lambda \upsilon_{i-k} + g(t_{i-k})) - (\lambda \omega_{i-k} + g(t_{i-k})) \right] \\ \upsilon_{i+1} - \omega_{i+1} &= \sum_{k=0}^{m-1} a_{m-1-k} (\upsilon_{i-k} - \omega_{i-k}) + h \sum_{k=-1}^{m-1} b_{m-1-k} \lambda (\upsilon_{i-k} - \omega_{i-k}) \\ \end{aligned} \label{eq:stb21} \end{equation}and therefore \(g(t)\) does not affect the stability analysis. Stability analysis may therefore be restricted to considering the homogeneous problem \(\frac{dy}{dt} = \lambda y(t)\).