multiphysicsL16.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{multiphysicsL16}
%\renewcommand{\dirname}{notes/ece1254/}
%\newcommand{\keywords}{ECE1254H}
%\input{../peeter_prologue_print2.tex}
%
%\beginArtNoToc
%\generatetitle{ECE1254H Modeling of Multiphysics Systems.  Lecture 16: LMS systems and stability.  Taught by Prof.\ Piero Triverio}
%\chapter{LMS systems and stability}
\label{chap:multiphysicsL16}
%\section{Disclaimer}
%
%Peeter's lecture notes from class.  These may be incoherent and rough.
%
\section{Residual for LMS methods.}
\index{residual}
\index{linear multistep methods}
%
\paragraph{Mostly on slides:} 12_ODS.pdf
%
Residual is illustrated in \cref{fig:lecture16:lecture16Fig1}, assuming that the iterative method was accurate until \( t_{n} \)
%
\imageFigure{../figures/ece1254-multiphysics/lecture16Fig1}{Residual illustrated.}{fig:lecture16:lecture16Fig1}{0.3}
%
\paragraph{Summary}
%
\begin{itemize}
\item [FE]: \( R_{n+1} \sim \lr{ \Delta t}^2 \).   This is of order \( p = 1 \).
\item [BE]: \( R_{n+1} \sim \lr{ \Delta t}^2 \).   This is of order \( p = 1 \).
\item [TR]: \( R_{n+1} \sim \lr{ \Delta t}^3 \).   This is of order \( p = 2 \).
\item [BESTE]: \( R_{n+1} \sim \lr{ \Delta t}^4 \).   This is of order \( p = 3 \).
\end{itemize}
%
\section{Global error estimate.}
\index{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 \).
%
\section{Stability.}
\index{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 \cref{fig:lecture16:lecture16Fig2}, \( \Real(\lambda) < 0 \) is assumed to ensure stability.
%
\imageFigure{../figures/ece1254-multiphysics/lecture16Fig2}{Stable system.}{fig:lecture16:lecture16Fig2}{0.2}
%
Linear stability theory can be thought of as asking the question: ``Is the solution of \cref{eqn:multiphysicsL16:40} computed by my LMS method also stable?''
%
Application of \cref{eqn:multiphysicsL16:20} to \cref{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
%
\begin{itemize}
\item discrete time system
\item FIR filter
\end{itemize}
%
The numerical solution \( x_n \) will be stable if \cref{eqn:multiphysicsL16:120} is stable.
%
A characteristic equation associated with \cref{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 \cref{fig:lecture16:lecture16Fig3}
%
\imageFigure{../figures/ece1254-multiphysics/lecture16Fig3}{Stability.}{fig:lecture16:lecture16Fig3}{0.3}
%
Observe that the \( \gamma's \) are dependent on \( \Delta t \).
%
Some of the details associated with this switch to discrete time systems have been assumed.  A refresher, by example, of those ideas can be found in \cref{chap:stabilityLDEandDiscreteTime}.
%FIXME: There's a lot of handwaving here that could use more strict justification.  Check if the text covers this in more detail.
%
\index{difference!forward Euler}
\makeexample{Forward Euler stability.}{example:multiphysicsL16:160}{
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{subequations}
\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}
\end{subequations}
%
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}
\boxedEquation{eqn:multiphysicsL16:260}{
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 \cref{fig:lecture16:lecture16Fig4}.
%
\imageFigure{../figures/ece1254-multiphysics/lecture16Fig4}{Stability region of FE.}{fig:lecture16:lecture16Fig4}{0.2}
%
All poles of my system must be inside the stability region in order to get stable \( \gamma \).
}
%
%\EndNoBibArticle