% -*-latex-*- Put EMACS into LaTeX-mode
% Verbal description for system ImplicitRC (ImplicitRC_desc.tex)
% Generated by MTT on Wednesday June 24 09:50:17 BST 1998.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This report describes the \emph{implicit} integration methods
available in MTT. They are introduced to provide simulation
of systems within the following context:
\begin{enumerate}
\item The system may be stiff with a mixture of slow and fast
(possibly due to approximating algebraic loops) subsystems.
\item The fast parts of the response are of no interest
\item A fixed sample interval is required -- possibly for real-time
simulation or control
\item The system is nonlinear.
\item The solution of nonlinear algebraic equations is to be avaided.
\end{enumerate}
The following sections consider the linear and nonlinear versions
respectively. The ideas are based on a standard textbook
\footnote{Press et al: \emph{Numerical Recipes in C}, 2nd edition,
1992. Cambridge, Section 16.6}.
\subsection{Implicit integration - the linear case}
\label{sec:linear}
Consider the \emph{linear} system:
\begin{equation}
\label{eq:sys}
\dot x = A x + B u
\end{equation}
For the purposes of simulation, it can be discretised (with sample
interval $\Delta t$) in at least two
ways:
\begin{enumerate}
\item $ \dot x \approx \frac{x_{i+1} - x_{i}}{\Delta t}$
\item $ \dot x \approx \frac{x_{i} - x_{i-1}}{\Delta t}$
\end{enumerate}
The former is gives rise to the \emph{forward} Euler or \emph{explicit}
integration scheme:
\begin{equation}
x_{i+1} = x_{i} + \Delta t \left [ A x_{i} + B_{i} u \right ]
\end{equation}
and the latter gives rise to the \emph{backward} Euler or \emph{implicit}
integration scheme:
\begin{equation}
x_{i} = x_{i-1} + \Delta t \left [ A x_{i} + B_{i} u \right ]
\end{equation}
which must be rewritten as:
\begin{equation}
x_{i} = \left [ I - \Delta t A \right ]^{-1} x_{i-1} + \Delta t B_{i} u
\end{equation}
for the purposes of implementation.
The explicit method gives simple implementation whereas the implicit
method requires matrix inversion. However, the explicit method is only
stable if:
\begin{equation}
\Delta t < \frac{2}{| \lambda |}
\end{equation}
where $\lambda$ is the \emph{largest} eigenvalue of $A$. If this
largest eigenvalue is real so $\lambda = \frac{1}{\tau}$ where $\tau$
is the \emph{smallest} system time constant:
\begin{equation}
\Delta t < 2 \tau
\end{equation}
If the system is stiff, that is it contains at least one small time
constant relative to the dominant time constants, Euler integration is
not feasible due to the very small sample interval $\Delta t$
required.
In contrast, the implicit method is stable.
\subsubsection{Example}
The acausal bond graph of system \textbf{ImplicitRC} is
displayed in Figure \Ref{ImplicitRC_abg} and its label
file is listed in Section \Ref{sec:ImplicitRC_lbl}.
The subsystems are listed in Section \Ref{sec:ImplicitRC_sub}.
The system represents two simple RC circuits in series with
differential equations as given in Section \Ref{sec:ImplicitRC_ode.tex} and
transfer function as given in Section \Ref{sec:ImplicitRC_tf.tex}.
For the purposes of this example the two time constants are $1$ and
$\epsilon=10^{-3}$ -- this is a stiff system. All of the simulations
use a sample interval of $\Delta t = 0.1$ ang the input is a unit
step. Section \Ref{sec:ImplicitRC_sro} shows the exact (computed from
the matrix exponential) solution, and Section \Ref{sec:ImplicitRC_odeso}
shows the solution by implicit integration.
The explicit solution is not shown, but was found to be unstable for
$\Delta t > 0.002$ as predicted.
\subsection{Implicit integration - the nonlinear case}}
\label{sec:nonlinear}
Consider the \emph{nonlinear} system:
\begin{equation}
\label{eq:sys}
\dot x = f(x,u)
\end{equation}
and suppose it can be linearised about any state and input to give:
\begin{equation}
A(x,u) = \frac{\partial f(x,u)}{\partial x}
\end{equation}
The corresponding \emph{implicit} scheme is:
\begin{equation}
x_{i} = x_{i-1} + \Delta t f(x_{i},u_{i})
\end{equation}
This is not easy to solve in general due to the set of non-linear
equations that need to be solved. To avoid this, consider a further
approximation:
\begin{equation}
f(x_{i},u_{i}) \approx f(x_{i-1},u_{i}) + A(x_{i-1},u_i) ( x_{i} - x_{i-1} )
\end{equation}
This then gives the \emph{semi-implicit} scheme
\begin{equation}
x_{i} = x_{i-1} + \Delta t \left [ f(x_{i-1},u_{i}) + A(x_{i-1},u_i) (
x_{i} - x_{i-1} ) \right ]
\end{equation}
which can be rewritten as:
\begin{equation}\label{eq:implicit}
x_{i} = \left [ I - \Delta t A(x_{i-1},u_i) \right ]^{-1}
\Delta t \left [ f(x_{i-1},u_{i}) - A(x_{i-1},u_i) x_{i-1} \right ]
\end{equation}
Because of the approximations invoved, Equation \ref{eq:implicit} is
not guarenteed to be stable. Nevertheless, it should do a much better
job than the corresponding \emph{explicit} method for reasonably
smooth systems.
This method is chosen by setting
\begin{verbatim}
METHOD='Implicit'
\end{verbatim}
in the MTT simpar.txt file.
A further approximation arises by setting $A(x_{i-1},u_i) =
A(x_{0},u_0)$ ie computing it one only at the beginning of the
simulation.
This method is chosen by setting
\begin{verbatim}
METHOD='ImplicitL'
\end{verbatim}
in the MTT simpar.txt file.
Both methods make use of the \textbf{smx} ``state-matrix with state $x$''
representation of MTT which is generated symbolically from the system
bond graph.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End: