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

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

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

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

Substitution gives

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

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$$

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

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

\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}

or

\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
}.

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

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

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

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

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

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

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

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

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

\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.

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 $$\label{eqn:stabilityLDEandDiscreteTime:260} 0 = z^2 - \inv{2} z = z \lr{ z - \inv{2} }$$ $$\label{eqn:stabilityLDEandDiscreteTime:240} 0 = z^2 + z = z\lr{ z + 1 }$$ $$\label{eqn:stabilityLDEandDiscreteTime:280} 0 = z^2 + 4 z + 3 = (z + 3)(z + 1).$$ These have respective roots $$\label{eqn:stabilityLDEandDiscreteTime:300} z = 0, \inv{2}$$ $$\label{eqn:stabilityLDEandDiscreteTime:320} z = 0, -1$$ $$\label{eqn:stabilityLDEandDiscreteTime:340} z = -1, -3$$ 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.