[Click here for a PDF of this post with nicer formatting]

In class today was a highlight of stability methods for linear multistep methods. To motivate the methods used, it is helpful to take a step back and review stability concepts for LDE systems.

By way of example, consider a second order LDE homogeneous system defined by

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:20}
\frac{d^2 x}{dt^2} + 3 \frac{dx}{dt} + 2 = 0.
\end{equation}

Such a system can be solved by assuming an exponential solution, say

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:40}
x(t) = e^{s t}.
\end{equation}

Substitution gives

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:60}
e^{st} \lr{ s^2 + 3 s + 2 } = 0,
\end{equation}

The polynomial part of this equation, the characteristic equation has roots \( s = -2, -1 \).

The general solution of \ref{eqn:stabilityLDEandDiscreteTime:20} is formed by a superposition of solutions for each value of \(s\)

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:80}
x(t) = a e^{-2 t} + b e^{-t}.
\end{equation}

Independent of any selection of the superposition constants \( a, b \), this function will not blow up as \( t \rightarrow \infty \).

This “stability” is due to the fact that both of the characteristic equation roots lie in the left hand Argand plane.

Now consider a discretized form of this LDE. This will have the form

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:100}
\begin{aligned}
0 &=
\inv{\lr{\Delta t}^2}
\lr{ x_{n+2} – 2 x_{n-1} + x_n } + \frac{3}{\Delta t} \lr{ x_{n+1} – x_n } + 2
x_n \\
&=
x_{n+2} \lr{
\inv{\lr{\Delta t}^2}
}
+
x_{n+1} \lr{
\frac{3}{\Delta t}
-\frac{2}{\lr{\Delta t}^2}
}
+
x_{n} \lr{
\frac{1}{\lr{\Delta t}^2}
-\frac{3}{\Delta t}
+ 2
},
\end{aligned}
\end{equation}

or

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:220}
0
=
x_{n+2}
+
x_{n+1} \lr{
3 \Delta t – 2
}
+
x_{n} \lr{
1 – 3 \Delta t + 2 \lr{ \Delta t}^2
}.
\end{equation}

Note that after discretization, each subsequent index corresponds to a time shift. Also observe that the coefficients of this discretized equation are dependent on the discretization interval size \( \Delta t \). If the specifics of those coefficients are ignored, a general form with the following structure can be observed

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:120}
0 =
x_{n+2} \gamma_0
+
x_{n+1} \gamma_1
+
x_{n} \gamma_2.
\end{equation}

It turns out that, much like the LDE solution by characteristic polynomial, it is possible to attack this problem by assuming a solution of the form

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:140}
x_n = C z^n.
\end{equation}

A time shift index change \( x_n \rightarrow x_{n+1} \) results in a power adjustment in this assumed solution. This substitution applied to \ref{eqn:stabilityLDEandDiscreteTime:120} yields

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:160}
0 =
C z^n
\lr{
z^{2} \gamma_0
+
z \gamma_1
+
1 \gamma_2
},
\end{equation}

Suppose that this polynomial has roots \( z \in \{z_1, z_2\} \). A superposition, such as

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:180}
x_n = a z_1^n + b z_2^n,
\end{equation}

will also be a solution since insertion of this into the RHS of \ref{eqn:stabilityLDEandDiscreteTime:120} yields

\begin{equation}\label{eqn:stabilityLDEandDiscreteTime:200}
a z_1^n
\lr{
z_1^{2} \gamma_0
+
z_1 \gamma_1
+
\gamma_2
}
+
b
z_2^n
\lr{
z_2^{2} \gamma_0
+
z_2 \gamma_1
+
\gamma_2
}
=
a z_1^n \times 0
+b z_2^n \times 0.
\end{equation}

The zero equality follows since \( z_1, z_2 \) are both roots of the characteristic equation for this discretized LDE.
In the discrete \( z \) domain stability requires that the roots satisfy the bound \( \Abs{z} < 1 \), a different stability criteria than in the continuous domain. In fact, there is no a-priori guarantee that stability in the continuous domain will imply stability in the discretized domain. Let's plot those z-domain roots for this example LDE, using \( \Delta t \in \{ 1/2, 1, 2 \} \). The respective characteristic polynomials are \begin{equation}\label{eqn:stabilityLDEandDiscreteTime:260} 0 = z^2 - \inv{2} z = z \lr{ z - \inv{2} } \end{equation} \begin{equation}\label{eqn:stabilityLDEandDiscreteTime:240} 0 = z^2 + z = z\lr{ z + 1 } \end{equation} \begin{equation}\label{eqn:stabilityLDEandDiscreteTime:280} 0 = z^2 + 4 z + 3 = (z + 3)(z + 1). \end{equation} These have respective roots \begin{equation}\label{eqn:stabilityLDEandDiscreteTime:300} z = 0, \inv{2} \end{equation} \begin{equation}\label{eqn:stabilityLDEandDiscreteTime:320} z = 0, -1 \end{equation} \begin{equation}\label{eqn:stabilityLDEandDiscreteTime:340} z = -1, -3 \end{equation} Only the first discretization of these three yields stable solutions in the z domain, although it appears that \( \Delta t = 1 \) is right on the boundary.