## Disclaimer

Peeter’s lecture notes from class. These may be incoherent and rough.

## Residual for LMS methods

### Mostly on slides:

12_ODS.pdf

Residual is illustrated in fig. 1, assuming that the iterative method was accurate until $$t_{n}$$

### Summary

• [FE]: $$R_{n+1} \sim \lr{ \Delta t}^2$$. This is of order $$p = 1$$.
• [BE]: $$R_{n+1} \sim \lr{ \Delta t}^2$$. This is of order $$p = 1$$.
• [TR]: $$R_{n+1} \sim \lr{ \Delta t}^3$$. This is of order $$p = 2$$.
• [BESTE]: $$R_{n+1} \sim \lr{ \Delta t}^4$$. This is of order $$p = 3$$.

## Global error estimate

Suppose $$t \in [0, 1] s$$, with $$N = 1/{\Delta t}$$ intervals. For a method with local error of order $$R_{n+1} \sim \lr{ \Delta t}^2$$ the global error is approximately $$N R_{n+1} \sim \Delta t$$.

## Stability

Recall that a linear multistep method (LMS) was a system of the form

\begin{equation}\label{eqn:multiphysicsL16:20}
\sum_{j=-1}^{k-1} \alpha_j x_{n-j} = \Delta t \sum_{j=-1}^{k-1} \beta_j f( x_{n-j}, t_{n-j} )
\end{equation}

Consider a one dimensional test problem

\begin{equation}\label{eqn:multiphysicsL16:40}
\dot{x}(t) = \lambda x(t)
\end{equation}

where as in fig. 2, $$\Re(\lambda) < 0$$ is assumed to ensure stability.

Linear stability theory can be thought of as asking the question: “Is the solution of \ref{eqn:multiphysicsL16:40} computed by my LMS method also stable?”

Application of \ref{eqn:multiphysicsL16:20} to \ref{eqn:multiphysicsL16:40} gives

\begin{equation}\label{eqn:multiphysicsL16:60}
\sum_{j=-1}^{k-1} \alpha_j x_{n-j} = \Delta t \sum_{j=-1}^{k-1} \beta_j \lambda x_{n-j},
\end{equation}

or
\begin{equation}\label{eqn:multiphysicsL16:80}
\sum_{j=-1}^{k-1} \lr{ \alpha_j – \Delta \beta_j \lambda }
x_{n-j} = 0.
\end{equation}

With

\begin{equation}\label{eqn:multiphysicsL16:100}
\gamma_j = \alpha_j – \Delta \beta_j \lambda,
\end{equation}

this expands to
\begin{equation}\label{eqn:multiphysicsL16:120}
\gamma_{-1} x_{n+1}
+
\gamma_{0} x_{n}
+
\gamma_{1} x_{n-1}
+
\cdots
+
\gamma_{k-1} x_{n-k} .
\end{equation}

This can be seen as a

• discrete time system
• FIR filter

The numerical solution $$x_n$$ will be stable if \ref{eqn:multiphysicsL16:120} is stable.

A characteristic equation associated with \ref{eqn:multiphysicsL16:120} can be defined as

\begin{equation}\label{eqn:multiphysicsL16:140}
\gamma_{-1} z^k
+
\gamma_{0} z^{k-1}
+
\gamma_{1} z^{k-2}
+
\cdots
+
\gamma_{k-1} = 0.
\end{equation}

This is a polynomial with roots $$z_n$$ (poles). This is stable if the poles satisfy $$\Abs{z_n} < 1$$, as illustrated in fig. 3

Observe that the $$\gamma’s$$ are dependent on $$\Delta t$$.

FIXME: There’s a lot of handwaving here that could use more strict justification. Check if the text covers this in more detail.

### Example: Forward Euler stability

For $$k = 1$$ step.

\begin{equation}\label{eqn:multiphysicsL16:180}
x_{n+1} – x_n = \Delta t f( x_n, t_n ),
\end{equation}

the coefficients are $$\alpha_{-1} = 1, \alpha_0 = -1, \beta_{-1} = 0, \beta_0 =1$$. For the simple function above

\begin{equation}\label{eqn:multiphysicsL16:200}
\gamma_{-1} = \alpha_{-1} – \Delta t \lambda \beta_{-1} = 1
\end{equation}
\begin{equation}\label{eqn:multiphysicsL16:220}
\gamma_{0} = \alpha_{0} – \Delta t \lambda \beta_{0} = -1 – \Delta t \lambda.
\end{equation}

The stability polynomial is

\begin{equation}\label{eqn:multiphysicsL16:240}
1 z + \lr{ -1 – \Delta t \lambda} = 0,
\end{equation}

or

\begin{equation}\label{eqn:multiphysicsL16:260}
\boxed{
z = 1 + \delta t \lambda.
}
\end{equation}

This is the root, or pole.

For stability we must have

\begin{equation}\label{eqn:multiphysicsL16:280}
\Abs{ 1 + \Delta t \lambda } < 1,
\end{equation}

or
\begin{equation}\label{eqn:multiphysicsL16:300}
\Abs{ \lambda – \lr{ -\inv{\Delta t} } } < \inv{\Delta t},
\end{equation}

This inequality is illustrated roughly in fig. 4.

All poles of my system must be inside the stability region in order to get stable $$\gamma$$.