\chapter{Analysis of Modeling Equations}
\begin{flushright}
\parbox{5.5in}{Folk wisdom says that if you leave a cold alone, it lasts seven days; if you treat it, it will go away in a week. Nonetheless, people have sworn by such ``cures'' as this: hang your hat on the bedpost, drink gin until you see two hats, then pull the covers up until the siege is over.~\cite[page~52]{rd91}}
\end{flushright}

Most of dynamical systems analyzed analyzed by us in this book are modeled by ordinary differential equations. A main use of a mathematical model is to predict the system transient behavior. Unlike linear systems, where closed-form solutions can be written in terms of the system's eigenvalues\index{eigenvalue} and eigenvectors\index{eigenvector}, finding analytical, exact, solutions to nonlinear differential equations\index{differential equation!nonlinear} can be very difficult or impossible. However, we can approximate the solutions of differential equations with difference equations\index{difference equation} whose solutions can be obtained easily using a computer. In this chapter, after discussing methods for solving state equations modeling linear systems, we present a number of methods for solving a more general class of differential equations. The first class of methods allows us to graphically determine solutions to second-order differential equations. Then, we discuss numerical techniques for solving differential equations.
%After that, we present two methods of linear approximation of nonlinear systems.

Linear systems can be modeled using convolutions and in the case of lumped systems, we can use state-space equations. Analytical computation of convolutions can be quite involved. An efficient method of obtaining a solution to a model of a lumped system is to convert the transfer function  into a state-space model and then solve the state equation. We begin with a time-invariant case and then proceed with the time-varying case.


%\section*{Solving the State Equations}\addcontentsline{toc}{section}{Solving the State Equations}
\section{State Equation Solution of Uncontrolled System}
%\addcontentsline{toc}{subsection}{Solution of Uncontrolled System}
Given a dynamical system's time-invariant model of the form
\begin{equation}
\dot{\b x}(t)=\b A \b x(t)
\label{eq:xdobaxeq}
\end{equation}
subject to an initial condition
\begin{equation}
\b x(0)=\b x_0,
\label{eq:inicoax}
\end{equation}
where $\b A \in \mb{\Bb R}^{n \times n}$. The solution to~(\ref{eq:xdobaxeq}) subject to~(\ref{eq:inicoax}) is
\begin{equation}
\fbox{$\displaystyle \b x(t)=e^{\b A t}\b x_0 $}
\label{eqn:sol_uncon_sys_ti}
\end{equation}
The above is a generalization of the situation when $\b A$ is a scalar. The exponential matrix\index{exponential matrix} $e^{\b A t}$ is defined as
\[
e^{\b At}=\b I_n + t\b A + \frac{t^2}{2!}\b A^2+\frac{t^3}{3!}\b A^3+\cdots
\]
To verify that~\eqref{eqn:sol_uncon_sys_ti} solves~\eqref{eq:xdobaxeq} subject to~\eqref{eq:inicoax}, we use the fact that
\[
\frac{d}{dt}e^{\b At} = \b Ae^{\b At}= e^{\b At} \b A.
\]

Suppose that the initial condition $\b x(0)=\b x_0$ is replaced with a more general condition
\[
\b x(t_0)=\b x_0,
\]
where $t_0$ is a given initial time instant. Then, the solution to $\dot{\b x}(t)=\b A \b x(t)$ subject to $\b x(t_0)=\b x_0$ is
\[
\fbox{$\displaystyle\b x(t)=e^{\b A\left(t-t_0\right)}\b x_0 $}
\]
The matrix $e^{\b A\left(t-t_0\right)}$ is often written as
\[
e^{\b A\left(t-t_0\right)}=\b \Phi \left(t,t_0\right)
\]
and is called the {\em state transition matrix\/}\index{state transition matrix}, because it relates the state at any instant of time $t_0$ to the state at any other time $t$.


\section{Computing State Transition Matrix of LTI Systems}

For an overview of methods to calculate the state transition matrix\index{state transition matrix!of LTI system} for a linear, time-invariant~(LTI) system\index{LTI system}, we refer the reader to the classic paper on the subject by Moler and Van Loan~\cite{moler78}, where the numerical issues related to the computation of the matrix exponential are discussed. In what follows, we present a few simple methods to compute the matrix exponential.

\subsection{The Inverse Laplace Transform Method}
In the case when the matrix $\b A$ is of low dimensions, we can use the formula that results from applying the Laplace transform to the time-invariant system\index{time-invariant!system} $\dot{\b x}=\b A \b x,\quad \b x(0)=\b x_0$. The Laplace transform of the above matrix differential equation is
\[
s\b X(s)-\b x_0= \b A\b X(s).
\]
Hence
\[
\b X(s)=\left(s\b I_n-\b A\right)^{-1}\b x_0,
\]
and the inverse transform yields
\begin{equation}
\fbox{$\displaystyle e^{\b At}={\cal L}^{-1}\left\{\left(s\b I_n-\b A\right)^{-1}\right\} $}
\label{eq:inlatrofunco}
\end{equation}
where ${\cal L}^{-1}$ denotes the inverse Laplace transform\index{inverse Laplace transform}\index{Laplace transform} operator.
\begin{example}\label{exam:iLtr}\rm
We use the inverse Laplace method to compute $e^{\b A t}$ for
\[
\b A=\left[\begin{array}{cc}
0 & 3\\
2 & 1
\end{array}\right].
\]
We have
\begin{eqnarray}
e^{\b A t} &=& {\cal L}^{-1}\left(\left[ s\b I_2 - \b A  \right]^{-1}  \right)\nonumber\\
&=& {\cal L}^{-1}\left(\left[\begin{array}{cc}
s & -3\\
-2 & s-1
\end{array}\right]^{-1}\right)\nonumber\\
&=& {\cal L}^{-1}\left(\frac{1}{(s-3)(s+2)}\left[\begin{array}{cc}
s-1 & 3\\
2 & s
\end{array}\right]
\right)\nonumber\\
&=& \left[\begin{array}{cc}
0.4 e^{3t} + 0.6 e^{-2t} & 0.6 \left( e^{3t} - e^{-2t} \right)\\
0.4 \left( e^{3t} - e^{-2t} \right) & 0.6 e^{3t} + 0.4 e^{-2t}
\end{array}\right].
\label{eq:iLmex}
\end{eqnarray}
\end{example}
\markex

\subsection{Diagonalization Method}

A convenient way to calculate $e^{\b A t}$ is when the matrix $\b A$ has distinct eigenvalues\index{eigenvalues!distinct}\index{distinct eigenvalues}. In such a case the $n$-by-$n$ matrix $\b A$ has $n$ linearly independent eigenvectors as per the following theorem.
\begin{thm}\label{thm:dist_eigenval}
Suppose the characteristic equation,
\begin{eqnarray*}
\Delta (\lambda)&=&\det[\lambda \b I_n - \b A]\\
&=& \lambda^n + a_{n-1}\lambda^{n-1} + \cdots + a_1 \lambda + a_0\\
&=& 0,
\end{eqnarray*}
has $n$ distinct roots, $\lambda_1, \lambda_2, \ldots , \lambda_n$. Then, there exist $n$ linearly independent\index{linearly independent!vectors}\index{vectors!linearly independent} vectors $\b v_1, \b v_2, \ldots , \b v_n$ such that
\[
\b A\b v_i = \lambda_i \b v_i, \quad i=1,2,\ldots , n.
\]
\end{thm}
{\bf Proof}\quad See Exercise~\ref{ex:dist_eigenval_pf}.\\
\hspace*{5.5in} $\Box$

Let
\[
\b T=\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right].
\]
The matrix $\b T$ is invertible because the eigenvectors $\b v_i$, $i=1,2,\ldots, n$ are linearly independent. Then, because $\b T^{-1}\b T = \b I_n$, we have
\begin{eqnarray}
\b T^{-1}\b A\b T &=& \b T^{-1}\b A \left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right]\nonumber\\
&=& \b T^{-1} \left[\begin{array}{cccc}
\b A\b v_1 & \b A\b v_2 & \cdots & \b A\b v_n
\end{array}\right]\nonumber\\
&=& \b T^{-1} \left[\begin{array}{cccc}
\lambda_1\b v_1 & \lambda_2\b v_2 & \cdots & \lambda_n\b v_n
\end{array}\right]\nonumber\\
&=& \b T^{-1}\b T\left[\begin{array}{ccccc}
\lambda_1 & 0 & 0 & \cdots & 0\\
0 & \lambda_2 & 0 & \cdots & 0\\
\vdots & \vdots & {} & \ddots & \vdots\\
0 & 0 & 0 & \cdots & \lambda_n
\end{array}\right]\nonumber\\
&=& \b{\Lambda}.
\label{eq:diag_lambda}
\end{eqnarray}
We use the above result to compute the transition matrix $e^{\b A t}$. Recall that
\[
e^{\b At} = \b I_n + \b A t + \frac{\b A^2 t^2}{2!}  + \frac{\b A^3 t^3}{3!} + \cdots
\]
Note that
\[
\b{\Lambda}^k = \left[\begin{array}{ccccc}
\lambda_1^k & 0 & 0 & \cdots & 0\\
0 & \lambda_2^k & 0 & \cdots & 0\\
\vdots & \vdots & {} & \ddots & \vdots\\
0 & 0 & 0 & \cdots & \lambda_n^k
\end{array}\right],\quad k\ge 1
\]
Hence
\[
e^{\b{\Lambda} t} = \left[\begin{array}{ccccc}
e^{\lambda_1 t} & 0 & 0 & \cdots & 0\\
0 & e^{\lambda_2 t} & 0 & \cdots & 0\\
\vdots & \vdots & {} & \ddots & \vdots\\
0 & 0 & 0 & \cdots & e^{\lambda_n t}
\end{array}\right].
\]
From~\eqref{eq:diag_lambda}, we obtain
\[
\b A=\b T\b{\Lambda}\b T^{-1}.
\]
Therefore,
\begin{eqnarray*}
\b A^k &=& \left( \b T\b{\Lambda}\b T^{-1} \right)^k\\
&=&\left( \b T\b{\Lambda}\b T^{-1} \right)\left( \b T\b{\Lambda}\b T^{-1} \right)\times \cdots \times \left( \b T\b{\Lambda}\b T^{-1} \right)\\
&=& \b T \b{\Lambda}^k \b T^{-1},
\end{eqnarray*}
and hence
\begin{eqnarray}
e^{\b A t} &=& e^{\b T\b{\Lambda}\b T^{-1}}\nonumber\\
&=& \b I_n + \b T\b{\Lambda}\b T^{-1} t + \frac{\left( \b T\b{\Lambda}\b T^{-1} \right)^2 t^2}{2!} + \frac{\left( \b T\b{\Lambda}\b T^{-1} \right)^3 t^3}{3!} + \cdots \nonumber\\
&=& \b T\left( \b I_n + \b{\Lambda} t + \frac{\b{\Lambda}^2 t^2}{2!} + \frac{\b{\Lambda}^3 t^3}{3!}\right) \b T^{-1} \nonumber\\
&=& \b T e^{\b{\Lambda} t}\b T^{-1}.
\label{eq:diag_m_for}
\end{eqnarray}
\begin{example}\label{exam:diag_met}\rm
We use the diagonalization method to compute $e^{\b A t}$ for
\[
\b A=\left[\begin{array}{cc}
0 & 3\\
2 & 1
\end{array}\right].
\]
Note that this is the same matrix $\b A$ as in Example~\ref{exam:iLtr}. The eigenvalues of $\b A$ are distinct so the diagonalization method\index{diagonalization method} is applicable here. We first compute an eigenvector\index{eigenvector}, $\b v_1$, corresponding to the eigenvalue $\lambda_1=3$ by solving the eigenvalue-eigenvector\index{eigenvalue-eigenvector!equation} equation,
\[
\left(\b A- \lambda_1 \b I_2  \right)\b v_1 = \b 0.
\]
That is,
\[
\left(\b A- \lambda_1 \b I_2  \right)\b v_1 =\left[\begin{array}{cc}
-3 & 3\\
2 & -2
\end{array}\right]\left[\begin{array}{c}
v_{11}\\
v_{21}
\end{array}\right]= \left[\begin{array}{c}
0\\
0
\end{array}\right].
\]
We take
\[
\b v_1=\left[\begin{array}{c}
1\\
1
\end{array}\right].
\]
We next compute an eigenvector corresponding to the eigenvalue $\lambda_2=-2$ by solving the equation
\[
\left(\b A- \lambda_2 \b I_2  \right)\b v_2 =\left[\begin{array}{cc}
2 & 3\\
2 & 3
\end{array}\right]\left[\begin{array}{c}
v_{12}\\
v_{22}
\end{array}\right]= \left[\begin{array}{c}
0\\
0
\end{array}\right].
\]
We take
\[
\b v_2=\left[\begin{array}{c}
-3\\
2
\end{array}\right].
\]
We use the above eigenvectors to form the transformation matrix\index{transformation matrix}
\[
\b T=\left[\begin{array}{cc}
\b v_1 & \b v_2
\end{array}\right] = \left[\begin{array}{cc}
1 & -3\\
1 & 2
\end{array}\right].
\]
Its inverse is
\[
\b T^{-1}=\frac{1}{5}\left[\begin{array}{cc}
2 & 3\\
-1 & 1
\end{array}\right].
\]
We are now ready to compute $e^{\b A t}$ using~\eqref{eq:diag_m_for},
\begin{eqnarray*}
e^{\b At} &=& \b Te^{\Lambda t} \b T^{-1}\\
&=& \frac{1}{5}\left[\begin{array}{cc}
1 & -3\\
1 & 2
\end{array}\right]\left[\begin{array}{cc}
e^{3 t} & 0\\
0 & e^{-2 t}
\end{array}\right]\left[\begin{array}{cc}
2 & 3\\
-1 & 1
\end{array}\right]\\
&=&\frac{1}{5}\left[\begin{array}{cc}
2e^{3 t} + 3 e^{-2 t} & 3\left( e^{3 t}- e^{-2 t} \right)\\
2\left( e^{3 t}- e^{-2 t} \right) & 3e^{3 t} + 2 e^{-2 t}
\end{array}\right].
\end{eqnarray*}
The result is the same as when we used the inverse Laplace transform method.
\end{example}
\markex

The above example illustrates that the transformation $\b T$ producing the diagonal form of $\b A$ is a powerful computational tool. The existence of a diagonalizing $\b T$ is not always guaranteed. The inability to transform $\b A$ to diagonal form is connected with the existence of repeated eigenvalues\index{eigenvalues!repeated}\index{repeated eigenvalues} that we discuss next.

\section{Computing Functions of a Square Matrix}

In this section, we discuss methods of computing functions of a square matrix that is not necessarily diagonalizable. The methods discussed allow us to easily compute the matrix exponential and hence the state transition matrix\index{transition matrix} of LTI systems. A major tool used by the methods discussed is the so called Jordan canonical form\index{Jordan form} that we discuss next.

\subsection{The Jordan Canonical Form}

If a square $n$-by-$n$ matrix $\b A$ has repeated eigenvalues, then it may not be diagonalizable. Although some matrices possessing repeated eigenvalues can be transformed into diagonal form, in many cases it is only possible to obtain the so-called ``Jordan canonical form'' or just the ``Jordan from,''\index{Jordan canonical form|see{Jordan form}}\index{Jordan form} for short. To illustrate the point, we consider the case when the characteristic equation, $\Delta (\lambda)=(\lambda-\lambda_1)^n$, that is, the matrix $\b A$ has just one eigenvalue with multiplicity $n$. For simplicity, we analyze the case when $n=4$ and suppose that
\[
\mb{rank} \left[\begin{array}{c}
\b A - \lambda_1 \b I_4
\end{array}\right] = 3.
\]
That is, the nullity of the matrix $\left[\begin{array}{c}
\b A - \lambda_1 \b I_4
\end{array}\right]$ is one, which means that there is only one eigenvector corresponding to the eigenvalue $\lambda_1$. Thus such a matrix $\b A$ cannot be diagonalized. We need three more linearly independent vectors to represent $\b A$ in a new basis, where it has a particularly simple form that is as close as possible to being diagonal.

To proceed, we need a notion of a generalized eigenvector that we define next.
\begin{defn}\label{defn:gene_eigenvec}
A vector $\b v$ is called a generalized eigenvector\index{generalized eigenvector}\index{generalized eigenvector!of grade $k$}\index{eigenvector!generalized} of grade $k$ if
\[
\left(\b A - \lambda \b I_n \right)^k \b v = \b 0
\]
and
\[
\left( \b A - \lambda \b I_n \right)^{k-1} \b v \ne \b 0.
\]
\end{defn}
\markdefn

For $k=1$, we have
\[
\left( \b A-\lambda \b I_n \right) \b v = \b 0,
\]
which means that $\b v$ is an ordinary eigenvector corresponding to the eigenvalue $\lambda$.

Suppose now that $\b A\in \mb{\Bb R}^{4 \times 4}$ and $\b v$ is a generalized eigenvector of grade $n=4$ corresponding to the repeated eigenvalue $\lambda_1$ whose multiplicity is also $4$. Then by Definition~\ref{defn:gene_eigenvec},
\[
\left(\b A - \lambda_1 \b I_4 \right)^4 \b v = \b 0\quad\mb{but}\quad \left(\b A - \lambda_1 \b I_4 \right)^3 \b v \ne \b 0.
\]
Let
\[
\b v_1=\left(\b A - \lambda_1 \b I_4 \right)^3 \b v.
\]
Then,
\begin{eqnarray*}
\left(\b A - \lambda_1 \b I_4 \right) \b v_1 &=& \left(\b A - \lambda_1 \b I_4 \right)^4 \b v\\
&=& \b 0.
\end{eqnarray*}
Hence
\begin{equation}
\b A\b v_1 = \lambda_1 \b v_1.
\label{eq:first_in_chain}
\end{equation}
We also have to have
\[
\left(\b A - \lambda_1 \b I_4 \right)^2 \b v \ne \b 0,
\]
otherwise we would have $\left(\b A - \lambda_1 \b I_4 \right)^3 \b v=\b 0$ as well, which would mean that $\b v$ is not a generalized eigenvector of grade 4.

Let
\[
\b v_2=\left(\b A - \lambda_1 \b I_4 \right)^2 \b v.
\]
Note that
\begin{eqnarray*}
\left(\b A - \lambda_1 \b I_4 \right) \b v_2 &=& \left(\b A - \lambda_1 \b I_4 \right)^3 \b v\\
&=& \b v_1.
\end{eqnarray*}
Therefore,
\begin{equation}
\b A\b v_2 = \b v_1 + \lambda_1 \b v_2.
\label{eq:second_in_chain}
\end{equation}

We proceed in a similar way by letting
\[
\b v_3=\left( \b A - \lambda_1 \b I_4 \right) \b v.
\]
We observer that
\begin{eqnarray*}
\left(\b A - \lambda_1 \b I_4 \right) \b v_3 &=& \left(\b A - \lambda_1 \b I_4 \right)^2 \b v\\
&=& \b v_2.
\end{eqnarray*}
Hence,
\begin{equation}
\b A\b v_3 = \b v_2 + \lambda_1 \b v_3.
\label{eq:third_in_chain}
\end{equation}
Finally, we let
\[
\b v_4=\b v
\]
and we note that
\begin{eqnarray*}
\left(\b A - \lambda_1 \b I_4 \right) \b v_4 &=& \left(\b A - \lambda_1 \b I_4 \right) \b v\\
&=& \b v_3.
\end{eqnarray*}
Therefore,
\begin{equation}
\b A\b v_4 = \b v_3 + \lambda_1 \b v_4.
\label{eq:fourth_in_chain}
\end{equation}
Combining~\eqref{eq:first_in_chain}--\eqref{eq:fourth_in_chain}, we obtain
\begin{eqnarray}
\b A\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \b v_3 & \b v_4
\end{array}\right] &=& \left[\begin{array}{cccc}
\lambda_1\b v_1 & \b v_1 + \lambda_1\b v_2 & \b v_2 + \lambda_1\b v_3 & \b v_3 + \lambda_1\b v_4
\end{array}\right]\nonumber\\
&=& \left[\begin{array}{cccc}
\b v_1 & \b v_2 & \b v_3 & \b v_4
\end{array}\right]\left[\begin{array}{cccc}
\lambda_1 & 1 & 0 & 0\\
0 & \lambda_1 & 1 & 0\\
0 & 0 & \lambda_1 & 1\\
0 & 0 & 0 & \lambda_1
\end{array}\right].
\label{eq:pre_simtoJ}
\end{eqnarray}
Let
\[
\b T=\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \b v_3 & \b v_4
\end{array}\right].
\]
The matrix $\b T$ is invertible by construction, that is, the set of vectors $\left\{ \b v_1, \b v_2, \b v_3, \b v_4 \right\}$ forms a basis for $\mb{\Bb R}^4$. Therefore we can represent~\eqref{eq:pre_simtoJ} as
\begin{eqnarray}
\b T^{-1}\b A\b T &=& \left[\begin{array}{cccc}
\lambda_1 & 1 & 0 & 0\\
0 & \lambda_1 & 1 & 0\\
0 & 0 & \lambda_1 & 1\\
0 & 0 & 0 & \lambda_1
\end{array}\right]\nonumber\\
&=& \b J.
\label{eq:sim_to_J}
\end{eqnarray}
The representation  of $\b A$ in the basis $\left\{ \b v_1, \b v_2, \b v_3, \b v_4 \right\}$ has the form of a Jordan block of order 4.

\begin{example}\label{exam:Jord_formthre_by_three}\rm
Consider the matrix
\[
\b A=\left[\begin{array}{ccc}
5 & 2 & 4\\
0 & 5 & 3\\
0 & 0 & 5
\end{array}\right].
\]
The characteristic polynomial is $\Delta (\lambda)=(\lambda - 5)^3$ and
\[
\mb{rank} \left.\left[\begin{array}{c}
\b A - \lambda\b I_3
\end{array}\right]\right\arrowvert_{\lambda=5} = 2,
\]
that is, the nullity\index{nullity!of the characteristic matrix} of the characteristic matrix, $\left.\left[\begin{array}{c}
\b A - \lambda\b I_3
\end{array}\right]\right\arrowvert_{\lambda=5}$, is 1. We proceed to find $\b v_1$ such that
\begin{eqnarray*}
\left.\left[ \b A-\lambda \b I_3 \right]\right\arrowvert_{\lambda=5} \b v_1 &=&  \left[\begin{array}{ccc}
0 & 2 & 4\\
0 & 0 & 3\\
0 & 0 & 0
\end{array}\right]\left[\begin{array}{c}
v_{11}\\
v_{21}\\
v_{31}
\end{array}\right]\\
&=& \b 0.
\end{eqnarray*}
We take
\[
\b v_1=\left[\begin{array}{c}
1\\
0\\
0
\end{array}\right].
\]
Then we proceed to find $\b v_2$ such that
\[
\left.\left[ \b A-\lambda \b I_3 \right]\right\arrowvert_{\lambda=5} \b v_2=\b v_1,
\]
that is, we solve for $\b v_2$ satisfying
\[
\left[\begin{array}{ccc}
0 & 2 & 4\\
0 & 0 & 3\\
0 & 0 & 0
\end{array}\right]\left[\begin{array}{c}
v_{12}\\
v_{22}\\
v_{32}
\end{array}\right]=\left[\begin{array}{c}
1\\
0\\
0
\end{array}\right].
\]
We can take
\[
\b v_2=\left[\begin{array}{c}
0\\
1/2\\
0
\end{array}\right].
\]
Finally we find $\b v_3$ such that
\[
\left.\left[ \b A-\lambda \b I_3 \right]\right\arrowvert_{\lambda=5} \b v_3=\b v_2,
\]
that is, we solve for $\b v_3$ satisfying
\[
\left[\begin{array}{ccc}
0 & 2 & 4\\
0 & 0 & 3\\
0 & 0 & 0
\end{array}\right]\left[\begin{array}{c}
v_{13}\\
v_{23}\\
v_{33}
\end{array}\right]=\left[\begin{array}{c}
0\\
1/2\\
0
\end{array}\right].
\]
We take
\[
\b v_3=\left[\begin{array}{c}
0\\
-2/6\\
1/6
\end{array}\right].
\]
The transformation matrix $\b T$ has the form
\[
\b T=\left[\begin{array}{ccc}
\b v_1 & \b v_2 & \b v_3
\end{array}\right]=\left[\begin{array}{ccc}
1 & 0 & 0\\
0 & 1/2 & -2/6\\
0 & 0 & 1/6
\end{array}\right],
\]
and
\[
\b T^{-1}=\left[\begin{array}{ccc}
1 & 0 & 0\\
0 & 2 & 4\\
0 & 0 & 6
\end{array}\right].
\]
The matrix $\b A$ in the new basis has the form
\begin{eqnarray*}
\b T^{-1}\b A\b T &=& \left[\begin{array}{ccc}
1 & 0 & 0\\
0 & 2 & 4\\
0 & 0 & 6
\end{array}\right] \left[\begin{array}{ccc}
5 & 2 & 4\\
0 & 5 & 3\\
0 & 0 & 5
\end{array}\right]\left[\begin{array}{ccc}
1 & 0 & 0\\
0 & 1/2 & -2/6\\
0 & 0 & 1/6
\end{array}\right]\\
&=& \left[\begin{array}{ccc}
5 & 1 & 0\\
0 & 5 & 1\\
0 & 0 & 5
\end{array}\right],
\end{eqnarray*}
which is a third order Jordan block.
\end{example}
\markex


If the rank of the characteristic matrix\index{characteristic matrix} evaluated at  $\lambda_1$ is $n-2$, then the equation
\[
\left[ \b A-\lambda_1 \b I_n \right] \b v=\b 0,
\]
has two linearly independent solutions. Therefore, $\b A$ has two linearly independent eigenvectors. Let $\b u_1$ and $\b v_1$ be linearly independent eigenvectors associated with the eigenvalue $\lambda_i$. In order to transform $\b A$ into its Jordan form, we need to find $n-2$ generalized eigenvectors. In  this case we have two chains of generalized eigenvectors, $\{\b u_1, \b u_2,\ldots , \b u_k  \}$ and $\{\b v_1, \b u_2,\ldots , \b v_l \}$, where $k+l=n$. In the basis
\[
\{\b u_1, \b u_2,\ldots , \b u_k, \b v_1, \b u_2,\ldots , \b v_l \}
\]
the matrix $\b A$ will have a block diagonal form
\[
\b T^{-1}\b A\b T=\left[\begin{array}{cc}
\b J_1 & \b O\\
\b O & \b J_2
\end{array}\right],
\]
where $\b J_1$ is a Jordan block of order $k$ and $\b J_2$ is a Jordan block\index{Jordan block} of order $l$.
\begin{example}\label{exam:two_blks_Jord_form}\rm
Let
\[
\b A=\left[\begin{array}{ccc}
3 & 1 & -1\\
0 & 6 & -3\\
0 & 3 & 0
\end{array}\right].
\]
The characteristic polynomial\index{characteristic polynomial}
\[
\Delta (\lambda)=(\lambda-3)^3
\]
and
\begin{eqnarray*}
\mb{rank}\left.\left[ \b A-\lambda \b I_3 \right]\right\arrowvert_{\lambda=3} &=& \left[\begin{array}{ccc}
0 & 1 & -1\\
0 & 3 & -3\\
0 & 3 & -3
\end{array}\right]\\
&=& 1
\end{eqnarray*}
As two linearly independent solutions to
\[
\left.\left[ \b A-\lambda \b I_3 \right]\right\arrowvert_{\lambda=3}\b v=\b 0,
\]
we can take
\[
\b u_1=\left[\begin{array}{c}
1\\
3\\
3
\end{array}\right]\quad\mb{and}\quad \b v_1=\left[\begin{array}{c}
0\\
1\\
1
\end{array}\right].
\]
Now we have to find one generalized eigenvector\index{generalized eigenvector} that together with $\b u_1$ and $\b v_1$ would form a basis for $\mb{\Bb R}^3$ and enable us to construct a transformation that brings $\b A$ into its Jordan form. We try first $\b u_1$ to construct $\b u_2$ such that
\[
\left.\left[ \b A-\lambda \b I_3 \right]\right\arrowvert_{\lambda=3} \b u_2=\b u_1,
\]
that is, we solve for $\b u_2$ satisfying
\[
\left[\begin{array}{ccc}
0 & 1 & -1\\
0 & 3 & -3\\
0 & 3 & -3
\end{array}\right]\left[\begin{array}{c}
u_{12}\\
u_{22}\\
u_{32}
\end{array}\right]=\left[\begin{array}{c}
1\\
3\\
3
\end{array}\right].
\]
We can take
\[
\b u_2=\left[\begin{array}{c}
1\\
1\\
0
\end{array}\right].
\]
We form the transformation matrix
\[
\b T=\left[\begin{array}{ccc}
\b u_1 & \b u_2 & \b v_1
\end{array}\right]=\left[\begin{array}{ccc}
1 & 1 & 0\\
3 & 1 & 1\\
3 & 0 & 1
\end{array}\right].
\]
Its inverse is
\[
\b T^{-1}=\left[\begin{array}{ccc}
1 & -1 & 1\\
0 & 1 & -1\\
-3 & 3 & -2
\end{array}\right].
\]
It is then easy to check that
\[
\b T^{-1}\b A\b T=\left[\begin{array}{cc}
\b J_1 & \b O\\
\b O & \b J_2
\end{array}\right]=\left[\begin{array}{ccc}
3 & 1 & 0\\
0 & 3 & 0\\
0 & 0 & 3
\end{array}\right],
\]
where
\[
\b J_1=\left[\begin{array}{cc}
3 & 1\\
0 & 3
\end{array}\right]\quad\mb{and}\quad \b J_2=3.
\]
\end{example}
\markex



In general, the Jordan canonical form\index{Jordan canonical form} is a block\index{block!diagonal matrix} diagonal matrix\index{block diagonal matrix}, where each diagonal block is a Jordan block corresponding to a possibly different eigenvalue, that is,
\begin{eqnarray}
\b T^{-1}\b A\b T &=& \left[\begin{array}{ccccccccccc}
\lambda_1 & 1 & 0 & \cdots & 0 & {} & {} & {} & {} &{} &{} \\
0 & \lambda_1 & 1 & \cdots & 0 & {} & {} & {} & {} &{} &{} \\
0 & 0 & \lambda_1 & \cdots & 0 & {} & {} & {} & {} &{} &{} \\
\vdots &\vdots &\vdots & \ddots & \vdots &{} &{}& {}& {} & {} & {}\\
0 & 0 & 0 & \cdots & \lambda_1 & {} & {} & {} & {} &{} &{} \\
%{} & {} & {} & {} &{} &{} & {} & {} & {} & {} &{} \\
{} & {} & {} & {} &{} &\ddots & {} & {} & {} & {} &{} \\
{} & {} & {} & {} &{} &{} &\lambda_r & 1 & 0 & \cdots & 0\\
{} & {} & {} & {} &{} &{} &0 & \lambda_r & 1 & \cdots & 0\\
{} & {} & {} & {} &{} &{} &0 & 0 & \lambda_r & \cdots & 0\\
{} &{}& {}& {} & {} & {} &\vdots &\vdots &\vdots & \ddots & \vdots\\
{} & {} & {} & {} &{} &{} &0 & 0 & 0 & \cdots & \lambda_r
\end{array}\right]\nonumber\\
&=& \left[\begin{array}{cccc}
\b J_1 & \b O & \cdots &\b O\\
\b O & \b J_2 & {} &\b O\\
\vdots & {} & \ddots & \vdots\\
\b O & \b O& \cdots & \b J_q
\end{array}\right],
\label{eq:gen_Jord_can_form}
\end{eqnarray}
where each $r_i\times r_i$ Jordan block\index{Jordan block} has the form $\b J_i=\lambda_i$ if $r_i=1$, or, for $r_i>1$ the corresponding Jordan block has a non-diagonal\index{non-diagonal} structure. Note that
\[
n=r_1 + r_2 + \cdots + r_q
\]
and
\[
q \ge r,
\]
which means that there may be more than one Jordan blocks associated with one eigenvalue as was the case in Example~\ref{exam:two_blks_Jord_form}.
\begin{example}\label{exam:two_eig_Jord_form}\rm
To illustrate the general case depicted in~\eqref{eq:gen_Jord_can_form}, we consider, as in Chen~\cite[p.~60]{ctchen99}, a $5$-by-$5$ matrix $\b A$ with repeated eigenvalue $\lambda_1$ of multiplicity 4 and a simple eigenvalue $\lambda_2$, that is, the characteristic polynomial of $\b A$ is
\[
\Delta(\lambda)=(\lambda-\lambda_1)^4(\lambda-\lambda_2).
\]
Depending on the nullity\index{nullity} of the matrix $[\b A-\lambda_1\b I_5]$, the Jordan form of the matrix $\b A$ make take one of the following forms:
\begin{itemize}
\item If the nullity of $[\b A-\lambda_1\b I_5]$ is 1, then there is only one eigenvector\index{eigenvector}, say $\b u_1$ associated with $\lambda_1$ and we need to compute three generalized eigenvectors to form a basis $\{\b u_1, \b u_2, \b u_3, \b u_4, \b v_1\}$, where $\b v_1$ is an eigenvector associated with $\lambda_2$. The matrix $\b A$ in this basis will have the form:
    \[
    \hat{\b A}_1=\b T_1^{-1}\b A\b T_1 = \left[\begin{array}{cccccc}
\lambda_1 & 1 & 0 & 0 & \vdots & 0\\
0 & \lambda_1 & 1 & 0 &  \vdots & 0\\
0 & 0 & \lambda_1 & 1 &   \vdots & 0\\
0 & 0 & 0 &\lambda_1 &    \vdots & 0\\
\cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\
0 & 0 & 0 & 0 & \vdots & \lambda_2
\end{array}\right].
    \]
Note that we can re-order the basis vectors as $\{\b v_1, \b u_1, \b u_1, \b u_3, \b u_4\}$. In this basis the matrix $\b A$ representation has the form
\[
\hat{\b A}_2=\b T_2^{-1}\b A\b T_2 = \left[\begin{array}{cccccc}
\lambda_2 & \vdots & 0 & 0 & 0 & 0\\
\cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\
0 & \vdots & \lambda_1 & 1 & 0 &   0\\
0 & \vdots & 0 & \lambda_1 & 1  & 0 \\
0 & \vdots & 0 & 0 & \lambda_1 &  1\\
0 & \vdots & 0 & 0 & 0 & \lambda_1
\end{array}\right].
\]
    \item If the nullity of $[\b A-\lambda_1\b I_5]$ is 2, then there are two eigenvectors associated with $\lambda_1$, that is, there are two Jordan blocks associated with $\lambda_1$. The matrix $\b A$ may take one of the following forms:
\[
\hat{\b A}_3=\b T_3^{-1}\b A\b T_3 = \left[\begin{array}{ccccccc}
\lambda_1 & 1 & 0 & \vdots & 0 &{}&  0 \\
0 & \lambda_1 & 1 & \vdots & 0 &{} & 0 \\
0 & 0 & \lambda_1 & \vdots & 0  & {} &0 \\
\cdots & \cdots & \cdots & \cdots & \cdots & \cdots &{} \\
0 & 0 & 0 & \vdots & \lambda_1 &  \vdots &0\\
{} & {} & {} &  \cdots  & \cdots & \cdots & \cdots\\
0  & 0 & 0 & {}  &0& \vdots &\lambda_2
\end{array}\right]
\]
or
\[
\hat{\b A}_4=\b T_4^{-1}\b A\b T_4 = \left[\begin{array}{ccccccc}
\lambda_1 & 1 & \vdots & 0  & 0 &{}&  0 \\
0 & \lambda_1 & \vdots & 0  & 0 &{} & 0 \\
\cdots & \cdots & \cdots & \cdots & \cdots & \cdots &{} \\
0 & 0 & \vdots & \lambda_1 & 1 & \vdots & 0 \\
0 & 0 &  \vdots & 0 & \lambda_1 &  \vdots &0\\
{} & {} & \cdots &  \cdots  & \cdots & \cdots & \cdots\\
0  & 0 & {} & 0  &0& \vdots &\lambda_2
\end{array}\right].
\]
    \item If the nullity of $[\b A-\lambda_1\b I_5]$ is 3, then there are three eigenvectors associated with $\lambda_1$, that is, there are three Jordan blocks associated with $\lambda_1$. The matrix $\b A$ may take the following form:
        \[
\hat{\b A}_5=\b T_5^{-1}\b A\b T_5 = \left[\begin{array}{cccccccc}
\lambda_1 & 1 & \vdots & 0  & {} & 0 &{}&  0 \\
0 & \lambda_1 & \vdots & 0  & {} & 0 &{} & 0 \\
\cdots & \cdots & \cdots & \cdots & \cdots & {} &{} & {}\\
0 & 0  & \vdots & \lambda_1 &  \vdots & 0 &{} &0\\
{} & {} & \cdots & \cdots & \cdots & \cdots &\cdots \\
0 & 0 &  {} & 0&  \vdots & \lambda_1 &  \vdots  & 0\\
{} & {} & {} & {}  & \cdots & \cdots & \cdots&\cdots\\
0  & 0 & {}& 0 & {}  &0& \vdots &\lambda_2
\end{array}\right].
\]
\item Finally, if the nullity of $[\b A-\lambda_1\b I_5]$ is 4, then there are four eigenvectors associated with $\lambda_1$, that is, there are four Jordan blocks associated with $\lambda_1$. The Jordan form of the matrix matrix $\b A$ is a diagonal matrix,
    \[
  \hat{\b A}_6=\b T_6^{-1}\b A\b T_6 = \left[\begin{array}{cccccc}
\lambda_1 & 0 & 0 & 0 & 0 \\
0 & \lambda_1 & 0 & 0 &  0 \\
0 & 0 & \lambda_1 & 0 &   0 \\
0 & 0 & 0 &\lambda_1 &    0 \\
0 & 0 & 0 & 0 &  \lambda_2
\end{array}\right].
    \]
\end{itemize}
\end{example}
\markex

To summarize our discussion, we need the following definitions.
\begin{defn}
For a given eigenvalue $\lambda_i$ of an $n$-by-$n$ matrix $\b A$, its algebraic multiplicity\index{algebraic multiplicity}\index{algebraic!multiplicity} is the multiplicity of it in the characteristics polynomial $\Delta (\lambda)$, that is, the algebraic multiplicity of $\lambda_i$ is the sum of orders of all Jordan blocks associated with $\lambda_i$.
\end{defn}\markdefn
\begin{defn}
For a given eigenvalue $\lambda_i$ of an $n$-by-$n$ matrix $\b A$, its geometric multiplicity\index{geometric multiplicity}\index{geometric!multiplicity} is the nullity of the associated characteristic matrix\index{characteristic matrix} $\left[\begin{array}{c}
\b A - \lambda_i\b I_n
\end{array}\right]$, that is, the geometric multiplicity of $\lambda_i$ equals the number of Jordan blocks associated with it.
\end{defn}\markdefn

Then we have the following theorem:
\begin{thm}
A given $n$-by-$n$ matrix $\b A$ is diagonalizable if and only if geometric multiplicities of its eigenvalues are the same as their algebraic multiplicities, that is, if and only if $\b A$ has $n$ linearly independent\index{linearly independent} eigenvectors.
\end{thm}

We close our discussion on the Jordan canonical form with an observation by Owens~\cite[p.~34]{owens81} that in almost all situations of practical relevance, the matrix $\b A$ has distinct eigenvalues. So in a sense the cases described in this section can be considered ``pathological''\index{pathological} from the applications point of view. In the literature, see for example~\cite[p.~61]{ctchen99}, if $\b A$ cannot be diagonalized then it is called to be defective\index{defective}. Calculation of the Jordan canonical is numerically ill conditioned\index{ill conditioned!numerically}. One can try MATLAB's function \verb!jordan!. However, the matrix $\b A$ must be known exactly, which in MATLAB's language means that its elements must be integers or ratios of small integers. Finally, to diagonalize a square matrix\index{square matrix!diagonalize} using MATLAB, one can call the function \verb!eig!. If, however, $\b A$ is not diagonalizable,\index{not diagonalizable} then the function \verb!eig! will give incorrect solution.


\subsection{Computing Polynomial Functions of a Square Matrix}

If $\b A$ is an $n$-by-$n$ square matrix then its $k$-th power is
\[
\b A^k=\underbrace{\b A\b A \times \cdots \times \b A}_{\mb{$k$ times}},
\]
where
\[
\b A^0=\b I_n.
\]
Suppose now that we have a polynomial
\[
f(\lambda)=\lambda^n + \alpha_{n-1}\lambda^{n-1} +\cdots + \alpha_1 \lambda + \alpha_0.
\]
Then we define $f(\b A)$ as
\[
f(\b A)=\b A^n +\alpha_{n-1}\b A^{n-1} + \cdots + \alpha_1 \b A + \alpha_0 \b I_n.
\]
For example, if
\[
f(\lambda)=\lambda^2 + 2\lambda -3=(\lambda +3)(\lambda -1),
\]
then
\[
f(\b A)=\b A^2 + 2\b A - 3\b I_n= (\b A + 3\b I_n)(\b A-\b I_n).
\]
Now if the matrix $\b A$ is block diagonal, for example, such as
\[
\b A=\left[\begin{array}{ccc}
\b A_1 & \b O & \b O\\
\b O & \b A_2 & \b O\\
\b O & \b O & \b A_3
\end{array}\right],
\]
where $\b A_i$ are square matrices and $\b O$ denotes a matrix of zeros\index{matrix of zeros} of appropriate dimensions, then we have
\[
\b A^k=\left[\begin{array}{ccc}
\b A_1^k & \b O & \b O\\
\b O & \b A_2^k & \b O\\
\b O & \b O & \b A_3^k
\end{array}\right]
\]
and
\[
f(\b A)=\left[\begin{array}{ccc}
f(\b A_1) & \b O & \b O\\
\b O & f(\b A_2) & \b O\\
\b O & \b O & f(\b A_3)
\end{array}\right].
\]
Consider a similar matrix\index{similar matrix}\index{matrix!similar} $\hat{\b A}$ of the given matrix $\b A$, that is,
\[
\hat{\b A}=\b T^{-1}\b A \b T,
\]
where $\b T\in \mb{\Bb R}^{n\times n}$ is an invertible matrix. The above can be re-written as
\[
\b A=\b T\hat{\b A}\b T^{-1},
\]
and we have
\begin{eqnarray*}
\b A^k &=& \left( \b T\hat{\b A}\b T^{-1} \right)\left( \b T\hat{\b A}\b T^{-1} \right)\times \cdots \times \left( \b T\hat{\b A}\b T^{-1} \right)\\
&=& \b T\hat{\b A}^k\b T^{-1},
\end{eqnarray*}
or
\begin{eqnarray*}
\hat{\b A}^k &=& \left( \b T^{-1}\b A\b T \right)\left( \b T^{-1}\b A\b T \right)\times \cdots \times \left( \b T^{-1}\b A\b T \right)\\
&=& \b T^{-1}\b A^k\b T.
\end{eqnarray*}
Thus, we have
\begin{equation}
f\left(\b A\right)= \b Tf( \hat{\b A})\b T^{-1}
\label{eq:ahatAtinv}
\end{equation}
or
\begin{equation}
f(\hat{\b A})= \b T^{-1} f\left(\b A\right)\b T.
\label{eq:tinvaAt}
\end{equation}


Now let $\lambda_i$ be the $i$-th eigenvalue of $\b A$ and suppose that the algebraic multiplicity of $\lambda_i$ is $n_i$. Then the characteristic polynomial of $\b A$ can be written as,
\begin{equation}
\Delta(\lambda)=\det (\lambda\b I_n - \b A)= \prod_{i=1}^r \left(\lambda - \lambda_i \right)^{n_i}.
\label{eq:char_poly_of_A}
\end{equation}

We now discuss a very useful property of Jordan blocks that we utilize in the proof of the very important theorem of linear algebra, the theorem of Cayley and Hamilton\index{Cayley and Hamilton!theorem of}\index{Cayley and Hamilton theorem}. Suppose that we have a fourth order Jordan block,
\[
\b J=\left[\begin{array}{cccc}
\lambda_i & 1 & 0 & 0\\
0 & \lambda_i & 1 & 0\\
0 & 0 & \lambda_i & 1\\
0 & 0 & 0 & \lambda_i
\end{array}\right].
\]
Then we have
\[
(\b J-\lambda_i \b I)=\left[\begin{array}{cccc}
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & 0 & 0
\end{array}\right].
\]
Let us now compute successive powers of $(\b J-\lambda \b I)$. We obtain
\[
(\b J-\lambda_i \b I)^2=\left[\begin{array}{cccc}
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0
\end{array}\right],\quad
(\b J-\lambda_i \b I)^3=\left[\begin{array}{cccc}
0 & 0 & 0 & 1\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0
\end{array}\right],
\]
and
\[
(\b J-\lambda_i \b I)^k=\b O\quad\mb{for}\quad k\ge 4.
\]
A matrix with above property is called a {\it nilpotent} matrix\index{nilpotent matrix}.

In general, if we have a Jordan block of order $r_i$ then
\begin{equation}
(\b J-\lambda_i \b I)^k=\b O\quad\mb{for}\quad k\ge r_i.
\label{eq:nilpotentJ}
\end{equation}

We have all we need to state and prove the theorem of Cayley and Hamilton.
\begin{thm}\label{thm:CHthm}(Cayley-Hamilton theorem)\index{Cayley and Hamilton theorem}
Let
\[
\Delta (\lambda) = \det [\lambda \b I_n - \b A]
= \prod_{i=1}^r (\lambda - \lambda_i)^{n_i},
\]
be the characteristic polynomial of a matrix $\b A$. A square matrix satisfies its characteristic polynomial equation,
\[
\Delta(\b A) = \prod_{i=1}^r (\b A - \lambda_i\b I_n)^{n_i}= \b O.
\]
\end{thm}
{\bf Proof}\quad Suppose that $\hat{\b A}$ is the Jordan canonical form of $\b A$. Thus there exists a non-singular matrix $\b T$ such that
\[
\hat{\b A}=\b T^{-1}\b A\b T=\left[\begin{array}{cccc}
\b J_1 & \b O & \cdots &\b O\\
\b O & \b J_2 & {} &\b O\\
\vdots & {} & \ddots & \vdots\\
\b O & \b O& \cdots & \b J_q
\end{array}\right].
\]
Let us now consider the matrix polynomial $\Delta(\b A)$. We use~\eqref{eq:ahatAtinv} to obtain
\begin{eqnarray*}
\Delta(\b A) &=& \b T \Delta( \hat{\b A})\b T^{-1}\\
&=&\b T\prod_{i=1}^r \left(\hat{\b A}- \lambda_i\b I_n \right)^{n_i} \b T^{-1}\\
%&=& \b T\left(\hat{\b A}- \lambda_1\b I_n \right)^{n_1}\times \left(\hat{\b A}- \lambda_2\b I_n \right)^{n_2}\times \cdots \times \left(\hat{\b A}- \lambda_r\b I_n \right)^{n_r}\b T^{-1}\\
&=& \b T \prod_{i=1}^r\left[\begin{array}{cccc}
\b J_1-\lambda_i\b I_{r_1} & \b O & \cdots &\b O\\
\b O & \b J_2-\lambda_i\b I_{r_2} & {} &\b O\\
\vdots & {} & \ddots & \vdots\\
\b O & \b O& \cdots & \b J_q-\lambda_i\b I_{r_q}
\end{array}\right]^{n_i}\b T^{-1}\\
&=& \b T \prod_{i=1}^r\left[\begin{array}{cccc}
(\b J_1-\lambda_i\b I_{r_1})^{n_i} & \b O & \cdots &\b O\\
\b O & (\b J_2-\lambda_i\b I_{r_2})^{n_i} & {} &\b O\\
\vdots & {} & \ddots & \vdots\\
\b O & \b O& \cdots & (\b J_q-\lambda_i\b I_{r_q})^{n_i}
\end{array}\right]\b T^{-1}.
\end{eqnarray*}
The order of any Jordan block associated with $\lambda_i$ cannot exceed the corresponding algebraic multiplicity. Taking into account the nilpotent property given by~\eqref{eq:nilpotentJ}, we obtain
\begin{eqnarray}
\Delta(\b A) &=&\b T \left[\begin{array}{cccc}
\b O & \b O & \cdots &\b O\\
\b O & (\b J_2-\lambda_1\b I_{r_2})^{n_1} & {} &\b O\\
\vdots & {} & \ddots & \vdots\\
\b O & \b O& \cdots & (\b J_q-\lambda_1\b I_{r_q})^{n_1}
\end{array}\right]\times \nonumber\\
{}&{}&\quad \cdots \times \left[\begin{array}{cccc}
(\b J_1-\lambda_r\b I_{r_1})^{n_r} & \b O & \cdots &\b O\\
\b O & (\b J_2-\lambda_r\b I_{r_2})^{n_r} & {} &\b O\\
\vdots & {} & \ddots & \vdots\\
\b O & \b O& \cdots & \b O
\end{array}\right]\b T^{-1}\nonumber\\
&=& \b T \left[\begin{array}{cccc}
\b O & \b O & \cdots &\b O\\
\b O & \b O & {} &\b O\\
\vdots & {} & \ddots & \vdots\\
\b O & \b O& \cdots & \b O
\end{array}\right]\b T^{-1}\nonumber\\
&=&\b O,
\label{eq:C_H_thm_pf_end}
\end{eqnarray}
which completes the proof of the theorem.\\
\hspace*{5.5in} $\Box$

It follows from the Cayley-Hamilton theorem that
\[
\Delta(\b A)=\b A^n +\alpha_{n-1}\b A^{n-1} + \cdots + \alpha_1 \b A + \alpha_0 \b I_n=\b O,
\]
that is,
\begin{equation}
\b A^n=-\alpha_{n-1}\b A^{n-1} - \cdots - \alpha_1 \b A - \alpha_0 \b I_n.
\label{eq:Aton-alin_combo}
\end{equation}
Thus $\b A^n$ can be expressed as a linear combination\index{linear combination!of matrices} of the matrices
\[
\{\b I_n, \b A, \ldots, \b A^{n-1}  \}.
\]
If we multiply both sides of~\eqref{eq:Aton-alin_combo} by $\b A$, then we obtain
\begin{equation}
\b A^{n+1}=-\alpha_{n-1}\b A^{n} - \cdots - \alpha_1 \b A^2 - \alpha_0 \b A.
\label{eq:Atonplusone}
\end{equation}
Thus, $\b A^{n+1}$ can be expressed as a linear combination\index{linear combination} of
\[
\{\b A, \b A^2, \ldots, \b A^{n}  \}.
\]
Using~\eqref{eq:Aton-alin_combo}, we conclude that $\b A^{n+1}$ can be expressed as a linear combination of $\{\b I_n, \b A, \ldots, \b A^{n-1}  \}$. Using the same argument, we conclude that any $\b A^k$, $k\ge 0$, can be expressed as a linear combination of $\{\b I_n, \b A, \ldots, \b A^{n-1}  \}$. Therefore, any matrix polynomial $f(\b A)$ can also be expressed as
\[
f(\b A)=\beta_0\b I_n + \beta_1\b A +\cdots + \beta_{n-1}\b A^{n-1}
\]
for some $\beta_i$, $i=0,1,\ldots, n-1$. The above property is used to efficiently compute matrix polynomials of $\b A$, which can then be used to compute the matrix exponential\index{matrix exponential} $e^{\b At}$ as well as the state transition matrix for discrete-time\index{discrete-time!systems} systems. In what follows, we present such a method and illustrate it with examples. We follow the presentation of Chen~\cite[pp.~63--67]{ctchen99}.


Suppose we are given a polynomial $f(\lambda)$ and we are interested in calculating $f(\b A)$. We begin with dividing $f(\lambda)$ by $\Delta (\lambda)$ to obtain
\begin{equation}
f(\lambda)=q(\lambda)\Delta(\lambda) + h(\lambda),
\label{eq:didefbydelat}
\end{equation}
where $q(\lambda)$ is the quotient, $h(\lambda)$ is the remainder, and
\[
\deg_{\lambda} h(\lambda) < n.
\]
We evaluate $f(\b A)$ using~\eqref{eq:didefbydelat} and the Cayley-Hamilton theorem\index{Cayley-Hamilton theorem} to obtain
\begin{eqnarray*}
f(\b A)&=&q(\b A)\Delta(\b A) + h(\b A)\\
&=& q(\b A)\b O + h(\b A)\\
&=& h(\b A).
\end{eqnarray*}
Thus instead of evaluating the matrix polynomial $f(\b A)$ that may contain very high powers of $\b A$, we can find $f(\b A)$ by computing $h(\b A)$, where the highest power of $\b A$ is less than $n$. In addition, rather than performing the division~\eqref{eq:didefbydelat}\index{division}, we can find $h(\lambda)$ by assuming that it has the form
\[
h(\lambda)=\beta_0 + \beta_1 \lambda + \cdots + \beta_{n-1}\lambda^{n-1}
\]
and then solve $n$ equations for $\beta_i$, $i=0,1,\ldots,n-1$. Consider first the case when $\b A$ has distinct eigenvalues. To obtain equations involving the coefficients $\beta_i$, we note that
\[
f(\lambda_i)=q(\lambda_i)\Delta(\lambda_i) + h(\lambda_i) = h(\lambda_i),\quad i=1,2,\ldots, n.
\]
\begin{example}\label{ex:dist_eig_eqna_for_beta}\rm
For a given
\[
f(\lambda)=\lambda^{7},
\]
we compute
\[
f(\b A)=\b A^{7},
\]
where
\begin{equation}
\b A=\left[\begin{array}{cc}
0 & 2\\
1 & 1
\end{array}\right].
\label{eq:simpleA}
\end{equation}
The characteristic polynomial of $\b A$ is
\[
\Delta (\lambda)=\lambda^2-\lambda-2=(\lambda+1)(\lambda-2).
\]
Let
\[
h(\lambda)=\beta_0 + \beta_1\lambda.
\]
Then we have two equations:
\[
f(\lambda_i)=h(\lambda_i),\quad i=1,2,
\]
that is,
\[
f(-1)=(-1)^{7}=-1=\beta_0 -\beta_1
\]
and
\[
f(2)=2^{7}=128=\beta_0+2\beta_1 .
\]
We represent the above linear equations for $\beta_0$ and $\beta_1$ in the matrix format as
\[
\left[\begin{array}{cc}
1 & -1\\
1 & 2
\end{array}\right]\left[\begin{array}{c}
\beta_0\\
\beta_1
\end{array}\right]=\left[\begin{array}{c}
-1\\
128
\end{array}\right].
\]
Solving, we obtain $\beta_0=42$ and $\beta_1=43$. Hence
\[
f(\b A)=\b A^7=42\b I_2 + 43\b A=\left[\begin{array}{cc}
42 & 86\\
43 & 85
\end{array}\right].
\]
\end{example}
\markex
\begin{example}\label{ex:mat_exp_with_h}\rm
We can use the above method to evaluate $e^{\b At}$, where $\b A$ is given by~\eqref{eq:simpleA}. In other words, let $f(\lambda)=e^{\lambda t}$; we wish to find
\[
f(\b A)=e^{\b At}.
\]
Let $h(\lambda)=\beta_0 + \beta_1\lambda$. Then
\[
f(-1)=h(-1)
\]
yields
\[
e^{-t}=\beta_0-\beta_1,
\]
and
\[
f(2)=h(2)
\]
yields
\[
e^{2t}=\beta_0+2\beta_1.
\]
We thus have two equations for $\beta_0$ and $\beta_1$,
\[
\left[\begin{array}{cc}
1 & -1\\
1 & 2
\end{array}\right]\left[\begin{array}{c}
\beta_0\\
\beta_1
\end{array}\right]=\left[\begin{array}{c}
e^{-t}\\
e^{2t}
\end{array}\right].
\]
Solving the above gives
\[
\left[\begin{array}{c}
\beta_0\\
\beta_1
\end{array}\right]=\frac{1}{3}\left[\begin{array}{cc}
2 & 1\\
-1 & 1
\end{array}\right]\left[\begin{array}{c}
e^{-t}\\
e^{2t}
\end{array}\right]=\frac{1}{3}\left[\begin{array}{c}
2e^{-t}+e^{2t}\\
-e^{-t}+e^{2t}
\end{array}\right].
\]
Hence
\[
e^{\b At}=\beta_0\b I_2 + \beta_1\b A=\frac{1}{3}\left[\begin{array}{cc}
2e^{-t}+e^{2t} & -2e^{-t}+2e^{2t}\\
-e^{-t}+e^{2t} &  e^{-t}+2e^{2t}
\end{array}\right].
\]
\end{example}
\markex

When the characteristic polynomial has repeated zeros\index{repeated zeros}, we use the derivatives of $f(\lambda)$ and $h(\lambda)$ to obtain additional equations to determine the coefficients of $h(\lambda)$ as per the following theorem.
\begin{thm} Chen~\cite[p.~64]{ctchen99}\label{thm:handfon sepecofA}
We are given $f(\lambda)$ and an $n\times n$ matrix $\b A$ with characteristic polynomial
\[
\Delta (\lambda) = \prod_{i=1}^m (\lambda-\lambda_i)^{n_i},
\]
where $n=\sum_{i=1}^m n_i$. Let
\[
h(\lambda) = \beta_0 + \beta_1\lambda + \cdots + \beta_{n-1}\lambda^{n-1}.
\]
It is a polynomial of degree $n-1$ with $n$ unknown coefficients. These $n$ unknowns are to be solved from the following set of $n$ equations:
\[
f^{(l)}(\lambda_i)=h^{(l)}(\lambda_i)\quad\mb{for}\quad l=0,1,\ldots , n_i-1 \quad\mb{and}\quad i=1,2,\ldots , m
\]
where
\[
f^{(l)}(\lambda_i)=\left.\frac{d^l f(\lambda)}{d \lambda^l}\right.\arrowvert_{\lambda = \lambda_i}
\]
and $h^{(l)}(\lambda_i)$ is similarly defined. Then we have
\[
f(\b A)=h(\b A)
\]
and $h(\lambda)$ is said to be equal $f(\lambda)$ on the spectrum of $\b A$.
\end{thm}
\begin{example}\label{ex:mat_exp_sim_multi_ex}\rm
We compute $e^{\b A t}$ for
\[
\b A=\left[\begin{array}{ccc}
3 & 0 & 1\\
-1 & 1 & 2\\
-2 & 0 & 0
\end{array}\right].
\]
The characteristic polynomial is
\[
\Delta (\lambda)=(\lambda -1)^2(\lambda-2).
\]
Let $f(\lambda)=e^{\lambda t}$ and let $h(\lambda)=\beta_0 + \beta_1 \lambda + \beta_2 \lambda^2$. Then using Theorem~\ref{thm:handfon sepecofA}, we obtain
\begin{eqnarray*}
f(1)&=& h(1),\quad e^t = \beta_0 + \beta_1 + \beta_2,\\
f'(1)&=& h'(1), \quad te^t = \beta_1 + 2\beta_2,\\
f(2)&=& h(2),\quad e^{2t} = \beta_0 + 2\beta_1 + 4\beta_2.
\end{eqnarray*}
Thus we have a set of equations
\[
\left[\begin{array}{ccc}
1 & 1 & 1\\
0 & 1 & 2\\
1 & 2 & 4
\end{array}\right]\left[\begin{array}{c}
\beta_0\\
\beta_1\\
\beta_2
\end{array}\right]=\left[\begin{array}{c}
e^t\\
te^t\\
e^{2t}
\end{array}\right].
\]
Solving the above gives
\[
\left[\begin{array}{c}
\beta_0\\
\beta_1\\
\beta_2
\end{array}\right]=\left[\begin{array}{ccc}
0 & -2 & 1\\
2 & 3 & -2\\
-1 & -1 & 1
\end{array}\right]\left[\begin{array}{c}
e^t\\
te^t\\
e^{2t}
\end{array}\right]=\left[\begin{array}{c}
-2te^t + e^{2t}\\
2e^t + 3te^t -2 e^{2t}\\
-e^t - t e^t + e^{2t}
\end{array}\right].
\]
Hence
\begin{eqnarray*}
e^{\b A t} &=& h(\b A)\\
&=& \beta_0 \b I_3 + \beta_1\b A + \beta_2 \b A^2\\
&=& \left(-2te^t + e^{2t} \right)\left[\begin{array}{ccc}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{array}\right]+ \left(2e^t + 3te^t -2 e^{2t} \right)\left[\begin{array}{ccc}
3 & 0 & 1\\
-1 & 1 & 2\\
-2 & 0 & 0
\end{array}\right]\\
&{}& + \left(-e^t - t e^t + e^{2t} \right)\left[\begin{array}{ccc}
7 & 0 & 3\\
-8 & 1 & 1\\
-6 & 0 & -2
\end{array}\right]\\
&=& \left[\begin{array}{ccc}
-e^t+ 2e^{2t} & 0 & -e^t+ e^{2t}\\
6e^t+5te^t - 6e^{2t} & e^t & 3e^t+5te^t- 3e^{2t}\\
2e^t-2e^{2t} & 0 & 2e^t - e^{2t}
\end{array}\right].
\end{eqnarray*}
\end{example}
\markex
\begin{example}\label{ex:mat_exp_Jord_block}\rm
Suppose we are given the Jordan block\index{Jordan block} of order $n_i$,
\[
\b J = \left[\begin{array}{cccccc}
\lambda_i & 1 & 0 &  \cdots & 0 & 0\\
0 & \lambda_i & 1 &  \cdots & 0 & 0\\
0 & 0 & \lambda_i &  \cdots & 0 & 0\\
\vdots & \vdots & {} & \ddots & {}& \vdots\\
0 & 0 & 0 & \cdots & \lambda_i & 1\\
0 & 0 & 0 & \cdots & 0 & \lambda_i
\end{array}\right].
\]
The characteristic polynomial\index{characteristic polynomial} of $\b J$ is $\Delta (\lambda)=(\lambda - \lambda_i)^{n_i}$. For computational simplicity we select the $(n_i-1)$-st order polynomial $h(\lambda)$ of the form
\[
h(\lambda)=\beta_0 + \beta_1(\lambda-\lambda_i) + \beta_2(\lambda - \lambda_i)^2 + \cdots + \beta_{n_i - 1} (\lambda - \lambda_i)^{n_i-1}.
\]
The condition $f(\lambda)=h(\lambda)$ on the spectrum of $\b J$ yields
\[
\beta_0 = f(\lambda_i),\quad \beta_1 = f'(\lambda_i),\quad \beta_2 = \frac{f''(\lambda_i)}{2!},\ldots , \beta_{n_i - 1}= \frac{f^{(n_i-1)}(\lambda_i)}{(n_i-1)!}.
\]
Hence
\begin{eqnarray*}
f(\b J)&=&f(\lambda_i) \b I_{n_i} + \frac{f'(\lambda_i)}{1!}(\b J-\lambda_i\b I_{n_i}) + \frac{f''(\lambda_i)}{2!}(\b J-\lambda_i\b I_{n_i})^2 + \cdots\\
 {}&{}&\quad \cdots + \frac{f^{(n_i-1)}(\lambda_i)}{(n_i-1)!}(\b J-\lambda_i\b I_{n_i})^{n_i-1}.
\end{eqnarray*}
Using the nilpotent\index{nilpotent} property of the matrix $(\b J-\lambda_i\b I_{n_i})$, we obtain
\[
f(\b J)=\left[\begin{array}{ccccc}
f(\lambda_i) & f'(\lambda_i)/1! & f''(\lambda_i)/2! & \cdots &  f^{(n_i-1)}(\lambda_i)/(n_i-1)!\\
0 & f(\lambda_i) & f'(\lambda_i)/1! &\cdots & f^{(n_i-2)}(\lambda_i)/(n_i-2)! \\
0 & 0 & f(\lambda_i) &   \cdots & f^{(n_i-3)}(\lambda_i)/(n_i-3)! \\
%0 & 0 & 0 & f(\lambda_i) & \cdots  & 0\\
\vdots & \vdots & \vdots &  {} &  \vdots\\
0 & 0 & 0 & \ddots &  f'(\lambda_i)/1!\\
0 & 0 & 0 & \cdots  & f(\lambda_i)
\end{array}\right].
\]
In particular, if $f(\lambda)=e^{\lambda t}$, then we have
\[
e^{\b J t}=\left[\begin{array}{ccccc}
e^{\lambda_i t} & te^{\lambda_i t} & t^2e^{\lambda_i t}/2! & \cdots &  t^{n_i - 1}e^{\lambda_i t}/(n_i - 1)!\\
0 & e^{\lambda_i t} & te^{\lambda_i t} &\cdots & t^{n_i - 2}e^{\lambda_i t}/(n_i - 2)! \\
0 & 0 & e^{\lambda_i t} &   \cdots & t^{n_i - 3}e^{\lambda_i t}/(n_i - 3)! \\
%0 & 0 & 0 & f(\lambda_i) & \cdots  & 0\\
\vdots & \vdots & \vdots &  {} &  \vdots\\
0 & 0 & 0 & \ddots &  te^{\lambda_i t}\\
0 & 0 & 0 & \cdots  & e^{\lambda_i t}
\end{array}\right].
\]
\end{example}
\markex

\section{State Equation Solution of Controlled System}
%\addcontentsline{toc}{subsection}{Solution of Controlled System}

Consider a model of a controlled dynamical system\index{controlled dynamical system}\index{dynamical system!controlled},
\begin{equation}
\dot{\b x}(t)=\b A\b x(t) + \b B\b u(t)
\label{eq:cosymoli}
\end{equation}
subject to a given initial condition $\b x(0)=\b x_0$, where $\b B$ is an $n \times m$ matrix, and $m<n$ as is usually in practice. Premultiplying both sides of~(\ref{eq:cosymoli}) by the integrating factor, $e^{-\b At}$, yields
\[
e^{-\b At}\dot{\b x}(t)=e^{-\b At}\b A \b x(t) + e^{-\b At} \b B\b u(t).
\]
We then rearrange the above equation as
\begin{equation}
e^{-\b At}\dot{\b x}(t)-e^{-\b At}\b A \b x(t) = e^{-\b At} \b B\b u(t).
\label{eq:xdaxburear}
\end{equation}
Noting that
\[
\frac{d}{dt}\left( e^{-\b At}\b x(t) \right)=-\b Ae^{-\b At}\b x(t)+e^{-\b At}\dot{\b x}(t),
\]
we can represent~(\ref{eq:xdaxburear}) in the form
\[
\frac{d}{dt}\left(e^{-\b At}\b x(t)\right)=e^{-\b At}\b B\b u(t).
\]
Integrating both sides of the resulting equation gives
\[
e^{-\b At}\b x(t)-\b x(0)=\int_{0}^{t}e^{-\b A \tau}\b B\b u(\tau)d\tau.
\]
Hence
\[
\b x(t)=e^{\b At}\b x(0)+\int_{0}^{t}e^{\b A(t-\tau)}\b B\b u(\tau)d\tau.
\]
If the initial condition is $\b x(t_0)=\b x_0$, then integration from $t_0$ to $t$ yields
\begin{equation}
\fbox{$\displaystyle \b x(t)=e^{\b A(t-t_0)}\b x(t_0)+\int_{t_0}^{t}e^{\b A(t-\tau)}\b B\b u(\tau)d\tau $}
\label{eq:coreofcotrls}
\end{equation}

As in the case of the uncontrolled system, we can also use the Laplace transform technique to obtain the solution to~(\ref{eq:cosymoli}) subject to the initial condition $\b x(0)=\b x_0$. Indeed, taking the Laplace transform of~(\ref{eq:cosymoli}) and performing some manipulations yields
\[
\b X(s)=\left(sI_n-\b A\right)^{-1}\b x(0)+\left(sI_n-\b A\right)^{-1}\b B\b U(s).
\]
Hence
\[
\b x(t)={\cal L}^{-1}\left\{\left(sI_n-\b A\right)^{-1}\b x(0)\right\}+{\cal L}^{-1}\left\{\left(sI_n-\b A\right)^{-1}\b B\b U(s)\right\}.
\]
If $\b x_0=\b x(t_0)$, then
\[
\b x(t)={\cal L}^{-1}\left\{e^{-t_0s}\left(sI_n-\b A\right)^{-1}\b x(t_0)\right\}+{\cal L}^{-1}\left\{e^{-t_0s}\left(sI_n-\b A\right)^{-1}\b B\b U(s)\right\}.
\]


\section{From Continuous to Discrete Models}
Our objective here is to present a method for constructing linear discrete-time\index{discrete-time!linear} models from given linear continuous-time\index{continuous-time!linear} models.
%The obtained discrete models will be used to perform computations to generate control commands.

We use a sample-and-hold device that transforms a continuous signal, $f(t)$, into the staircase signal,
\[
f(kh),\quad kh\le t< (k+1)h,
\]
where $h$ is the sampling period\index{sampling period}. The sample and zero-order hold~(ZOH)\index{ZOH|see{zero-order hold}}\index{zero-order hold (ZOH)} operation is illustrated in Figure~\ref{fig:sample_hold}.
\begin{figure}[h]
\centering \epsfig{file=sample_zoh_2009.eps,width=12cm}
\caption{Sample and zero-order hold~(ZOH) element operating on a continuous function.}
\label{fig:sample_hold}
\end{figure}

Suppose that we are given a continuous-time model,
\begin{eqnarray}
\dot{\b x}(t) &=& \b A\b x(t) + \b B\b u(t),\quad \b x_0=\b x(t_0)\label{eq:state_eqn_forc2d}\\
\b y(t) &=& \b C\b x(t).
\label{eqn:out_eqn_for_c2d}
\end{eqnarray}
The solution to the state equation is
\begin{equation}
\b x(t)=e^{\b A(t-t_0)} \b x(t_0) +\int_{t_0}^t e^{\b A(t-\tau)} \b B \b u(\tau) d \tau.
\label{eq:contin_model}
\end{equation}
We assume that the input to the system is generated by a sample and zero-order device and has the form,
\[
\b u(t)=\b u(k), \quad kh\le t< (k+1)h.
\]
Let $t_0=kh$ and $t=(k+1)h$ and let us use shorthand notation,
\[
\b x[k] = \b x(kh)
\]
Then taking into account that $\b u(t)$ is constant on the interval $[kh, (k+1)h)$ we represent~(\ref{eq:contin_model}) as
\begin{eqnarray}
\b x[k+1] &=& e^{\b Ah}\b x[k]+\int_{kh}^{(k+1)h} e^{\b A(kh+h-\tau)}\b B\b u[k] d\tau\nonumber\\
&=& e^{\b Ah}\b x[k]+\int_{kh}^{(k+1)h} e^{\b A(kh+h-\tau)}\b B d\tau\, \b u[k].
\label{eq:pre_chang_lim}
\end{eqnarray}
Consider now the second term on the right-hand side of the above equation. Let
\[
\eta=kh+h-\tau.
\]
Then we can represent~(\ref{eq:pre_chang_lim}) as
\begin{eqnarray}
\b x[k+1] &=& e^{\b Ah}\b x[k]+\int_{kh}^{(k+1)h} e^{\b A(kh+h-\tau)}\b B d\tau\, \b u[k]\nonumber\\
&=& e^{\b Ah}\b x[k]+\int_{0}^{h} e^{\b A\eta}\b B d\eta\, \b u[k]\nonumber\\
&=& \b A_d \b x[k] +\b B_d\b u[k],
\label{eq:final_d_time}
\end{eqnarray}
where
\[
\b A_d=e^{\b Ah}\quad\mb{and}\quad \b B_d= \int_{0}^{h} e^{\b A\eta}\b B d\eta.
\]
The discrete output equation has the form
\[
\b y[k] = \b C\b x[k].
\]

We point out that the above discrete model\index{discrete model}, employing the zero-order hold, is exact\index{exact!discrete model} in the sense that the state, input, and output of this model are exactly the same at the sampling times as the state, input, and output of the underlying continuous-time model.


\begin{example}\label{ex:c2d}\rm
We construct an exact linear discrete state-space model for the continuous system model,
\begin{eqnarray*}
\dot{\b x}(t) &=& \left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\b x(k) + \left[\begin{array}{c}
2\\
1
\end{array}\right]u(t)\\
y(t) &=& \left[\begin{array}{cc}
1 & 0
\end{array}\right]\b x(t).
\end{eqnarray*}
The sampling period\index{sampling period} is $h=1$~sec.

We first evaluate $\b A_d=e^{\b Ah}$. Note that
$\b A^2=\b O$. Hence
\begin{eqnarray*}
\b A_d &=& \b I_2+\b Ah\\
&=& \left[\begin{array}{cc}
1 & 1\\
0 & 1
\end{array}\right].
\end{eqnarray*}
Next, we evaluate
\[
\b B_d=\int_0^h e^{\b A \tau}\b B d\tau.
\]
We obtain
\begin{eqnarray*}
\b B_d&=&\int_0^h e^{\b A \tau}\b B d\tau\\
&=&\int_0^h \left[\begin{array}{cc}
1 & \tau\\
0 & 1
\end{array}\right]\left[\begin{array}{c}
2\\
1
\end{array}\right]d\tau\\
&=& \int_0^1 \left[\begin{array}{c}
2+\tau\\
1
\end{array}\right]d\tau\\
&=& \left. \left[\begin{array}{c}
2\tau+\frac{1}{2}\tau^2\\
\tau
\end{array}\right] \right\arrowvert_0^1\\
&=& \left[\begin{array}{c}
2.5\\
1
\end{array}\right].
\end{eqnarray*}
The discrete model then has the form
\begin{eqnarray*}
\b x[k+1] &=& \left[\begin{array}{cc}
1 & 1\\
0 & 1
\end{array}\right]\b x[k] + \left[\begin{array}{c}
2.5\\
1
\end{array}\right]\\
y[k] &=& \left[\begin{array}{cc}
1 & 0
\end{array}\right]\b x[k].
\end{eqnarray*}
\end{example}

\section{State Plane Analysis}
The {\em state plane\/}\index{state plane}, which is a state space of two-dimensional\index{two-dimensional state space|see{state plane}} systems, is also called the {\em phase plane\/}\index{phase plane|see{state plane}}. Analysis in the state plane is applicable to linear and nonlinear systems modeled by second-order ordinary differential equations. The state-plane methods\index{state plane!methods} are graphical procedures for solving such equations. Using state-plane methods one can graphically determine the transient response of a second-order dynamical system. We now introduce relevant terms. Consider a class of second-order systems whose dynamics are modeled by the differential equation,
\begin{equation}
\ddot{x}=f(x,\dot{x}).
\label{eq:secoordernot}
\end{equation}
Define the state variables,
\[
x_1=x\quad\mb{and}\quad x_2=\dot{x}.
\]
Then, equation~(\ref{eq:secoordernot}) can be represented as
\begin{equation}
\left.\begin{array}{lll}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& f(x_1,x_2)
\end{array}\right\}
\label{eq:stspreofsoode}
\end{equation}
The plane with coordinates $x_1$ and $x_2$ is labeled as the state plane or phase plane. A solution of~(\ref{eq:secoordernot}) can be illustrated by plotting $x$ versus $t$, or by plotting $x_2$ versus $x_1$ using $t$ as a parameter. To each value of the state $\b x(t)=\left[\begin{array}{cc}
x_1(t) & x_2(t)
\end{array}\right]^T$ there corresponds a point in the $(x_1,x_2)$-plane called the {\em representative point\/}\index{representative point (RP)} (RP). As $t$ varies, the RP\index{RP|see{representative point}} describes a curve in the state plane called a {\em trajectory\/}\index{trajectory}. A family of trajectories is called a {\em phase portrait}\index{phase portrait}.

\subsection{Examples of Phase Portraits}
\begin{example}
\label{example_ph}\rm
Consider the spring-mass system\index{spring-mass system}\index{system!spring-mass} shown in Figure~\ref{fig_circ_ph}, where $x$ is the displacement of the mass $M$ from equilibrium\index{equilibrium}, and $k$ is the linear spring\index{linear spring} constant.
\begin{figure}
\centering
\includegraphics[scale=0.55]{circ_ph.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/circ_ph.bmp}
%\hskip1.4in\special{bmp:c:/mybook/circ_ph.bmp x=4in y=3in}
\caption{A spring-mass system of Example~\ref{example_ph}.}
\label{fig_circ_ph}
\end{figure}
The spring is assumed to be linear, that is, it obeys Hooke's law\index{Hooke's law}. Applying Newton's law\index{Newton's law}, we obtain the equation of motion of the mass,
\[
M\ddot{x}+kx=0.
\]
Let $\omega^2=k/M$. Then, the above equation can be represented in standard form,
\begin{equation}
\ddot{x}+\omega^2x=0.
\label{eq:stasecode}
\end{equation}
Observe that
\[
\ddot{x} = \dot{x}\frac{\ddot{x}}{\dot{x}}=\dot{x}\frac{\frac{d\dot{x}}{dt}}{\frac{dx}{dt}}= \dot{x}\frac{d\dot{x}}{dx}.
\]
Using the above, we represent equation~(\ref{eq:stasecode}) as
\[
\dot{x}\frac{d\dot{x}}{dx}+\omega^2x=0.
\]
Separating the variables yields
\[
\frac{\dot{x}}{\omega^2}d\dot{x}+xdx=0.
\]
Integrating the above gives
\[
\left( \frac{\dot{x}}{\omega}\right)^2+x^2=A^2,
\]
where $A$ is the integration constant\index{integration constant} determined by the initial conditions. We see that a phase portrait\index{phase portrait} of equation~(\ref{eq:stasecode}) consists of a family of circles centered at the origin of the $(x,\frac{\dot{x}}{\omega})$-plane as illustrated in Figure~\ref{fig_circ_tr}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{circ_tr.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/circ_tr.bmp}
%\hskip1.3in\special{bmp:c:/mybook/circ_tr.bmp x=4in y=3in}
\caption{A phase portrait for the spring-mass system of Example~\ref{example_ph}.}
\label{fig_circ_tr}
\end{figure}
Different circular trajectories, in this plane, correspond to different values of the constant $A$ representing initial conditions. Note that
\[
x=\omega\int\frac{\dot{x}}{\omega}dt,
\]
where $\omega > 0$. Hence, for positive values of $\dot{x}/\omega$ the displacement $x$ increases. This means that for the values $\dot{x}/\omega > 0$ the arrows on the trajectories indicate $x$ moving to the right. For negative values of $\dot{x}/\omega$, that is, for the values of $\dot{x}/\omega$ below the horizontal axis, the arrows on the trajectories indicate $x$ decreasing by going to the left. Therefore, as time increases the RP moves in a clockwise direction along the trajectory.
\end{example}\markex
\begin{example}
\label{example_sat}\rm
In this example, we consider a model of a satellite\index{satellite!model of} shown in Figure~\ref{fig_sat}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{satelli.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/satelli.bmp}
%\hskip1.1in\special{bmp:c:/mybook/satelli.bmp x=4in y=3in}
\caption{A model of a rigid satellite of Example~\ref{example_sat}.}
\label{fig_sat}
\end{figure}
We assume that the satellite is rigid, and in a frictionless environment\index{frictionless environment}. It can rotate about an axis perpendicular to the page as a result of torque $\tau$ applied to the satellite\index{satellite!thrusters} by firing the thrusters\index{thrusters!satellite}. The system input is the applied torque, and the system output is the attitude angle $\theta$. The satellite's moment of inertia is $I$. Applying the rotational version of Newton's second law\index{Newton's second law!rotational version}, we obtain the satellite's equation of motion,
\[
I\frac{d^2\theta}{dt^2}=\tau,
\]
which can be represented in the following equivalent form:
\[
\frac{d^2\theta}{dt^2}=\frac{1}{I}\tau.
\]
Let
\[
\frac{1}{I}\tau=u,
\]
then the equation describing the satellite's motion becomes
\begin{equation}
\frac{d^2\theta}{dt^2}=u.
\label{eq:doubintsy}
\end{equation}
We assume that when the thrusters fire, the thrust is constant, that is, $u=\pm U$. For simplicity, let $U=1$. Using the state variables $x_1=\theta$ and $x_2=\dot{\theta}$, we obtain the satellite's state-space model,
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2
\end{array}\right]=\left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\left[\begin{array}{c}
x_1\\
x_2
\end{array}\right]+\left[\begin{array}{c}
0\\
1
\end{array}\right]u,
\label{eq:stsmodoubint}
\end{equation}
where $u=\pm 1$. We first consider the case when $u=1$. With this control, the system equations are:
\[
\left.\begin{array}{lll}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& 1.
\end{array}\right\}
\]
Solving the second equation yields
\[
x_2=t+c_1,
\]
where $c_1$ is an integration constant. We substitute $x_2=t+c_1$ into $\dot{x}_1=x_2$ to get
\[
\dot{x}_1=t+c_1.
\]
Solving the above equation, we obtain
\begin{eqnarray*}
x_1 &=& \int(t+c_1)dt\\
&=& \frac{1}{2}(t+c_1)^2+c_2\\
&=& \frac{1}{2}x^2_2+c_2,
\end{eqnarray*}
where $c_2$ is an integration constant. In the state plane the above equation represents a family of parabolas\index{family of parabolas!phase portrait} open toward the positive $x_1$ axis. The representative point moves upward along these parabolas because $\dot{x}_2=1$, that is, $\dot{x}_2>0$. Typical system trajectories for this case are shown in Figure~\ref{fig:bbol1}.
\begin{figure}
\centering
\includegraphics[scale=0.55]{bbol1.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/bbol1.bmp}
%\hskip1in\special{bmp:c:/mybook/bbol1.bmp x=4in y=3in}
\caption{Trajectories of the system~(\ref{eq:stsmodoubint}) for $u=1$.}
\label{fig:bbol1}
\end{figure}

Similarly, when $u=-1$, we have
\[
\dot{x}_2=-1.
\]
Solving the above differential equation, yields
\[
x_2=-t+\tilde{c}_1,
\]
where $\tilde{c}_1$ is an integration constant. Substituting $x_2=-t+\tilde{c}_1$ into $\dot{x}_1=x_2$, and then solving the resulting equation yields
\[
x_1=-\frac{1}{2}x^2_2+\tilde{c}_2,
\]
where $\tilde{c}_2$ is an integration constant. In the state plane, this equation represents a family of parabolas open toward the negative $x_1$ axis. The representative point moves downward along these parabolas because $\dot{x}_2=-1$, that is, $\dot{x}<0$. Typical system trajectories for this case are shown in Figure~\ref{fig:bbol2}.
\begin{figure}
\centering
\includegraphics[scale=0.55]{bbol2.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/bbol2.bmp}
%\hskip1in\special{bmp:c:/mybook/bbol2.bmp x=4in y=3in}
\caption{Trajectories of the system~(\ref{eq:stsmodoubint}) for $u=-1$.}
\label{fig:bbol2}
\end{figure}
\end{example}\markex
In the above examples, systems' trajectories were constructed analytically. Analytical methods are useful for systems modeled by differential equations that can be easily solved. If the system of differential equations cannot be easily solved analytically, we can use graphical or numerical methods. We now describe a graphical method for solving second-order, or a system of two first-order, differential equations.

\subsection{The Method of Isoclines}
An {\em isocline\/}\index{isocline} is a curve along which the trajectories' slope is constant. In other words, an isocline is a locus\index{locus!set of points}, a set of points, of the trajectories' constant slope\index{constant slope}. To obtain equations that describe isoclines, we consider a system of two first-order differential equations of the form
\begin{equation}
\left.\begin{array}{lll}
\dot{x}_1 &=& \frac{dx_1}{dt}=f_1(x_1,x_2)\\
\dot{x}_2 &=& \frac{dx_2}{dt}=f_2(x_1,x_2).
\end{array}\right\}
\label{eq:systfirsto}
\end{equation}
Dividing the second of the above equations by the first one, we obtain
\begin{equation}
\frac{dx_2}{dx_1}=\frac{f_2(x_1,x_2)}{f_1(x_1,x_2)}.
\label{eq:ftbyf1afo}
\end{equation}
Thus, we eliminated the independent variable $t$ from the set of the first-order differential equations given by~(\ref{eq:systfirsto}). In equation~(\ref{eq:ftbyf1afo}), we consider $x_1$ as the independent variable, and $x_2$ as dependent. Let
\[
m=m(x_1,x_2)=\frac{dx_2}{dx_1}.
\]
Note that $m$ is just the slope of the tangent to the trajectory passing through the point $\left[\begin{array}{cc}
x_1 & x_2
\end{array}\right]^T$. The locus of the trajectory constant slope
\[
\frac{dx_2}{dx_1}=m(x_1,x_2)=\mb{constant},
\]
can be obtained from the equation
\[
f_2(x_1,x_2)=mf_1(x_1,x_2).
\]
The curve that satisfies the above equation is an isocline corresponding to the trajectories' slope $m$ because a trajectory crossing the isocline will have its slope equal to $m$. The idea of the method of isoclines is to construct several isoclines in the state plane. This then allows one to construct a field of local tangents\index{local tangents!field of} $m$. Then, the trajectory passing through any given point in the state plane is obtained by drawing a continuous curve following the directions of the field\index{directions of the field}.
\begin{example}
\label{example_phcirc}\rm
Consider the spring-mass system\index{spring-mass system} from Example~\ref{example_ph}. Defining the state variables $x_1=x$ and $x_2=\dot{x}$, we obtain the state equations
\[
\left.\begin{array}{lll}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& -\omega^2 x_1.
\end{array}\right\}
\]
Let, for simplicity, $\omega=1$. Then,
\[
\frac{dx_2}{dx_1} = m= -\frac{x_1}{x_2}.
\]
From the above we obtain the equation of isoclines
\[
x_2=-\frac{1}{m}x_1.
\]
This equation describes a family of straight lines through the origin. In Figure~\ref{fig_isocline}, we show some of the isoclines and a typical trajectory for the above spring-mass system.
\begin{figure}
\centering
%\includegraphics[scale=0.55]{isocline.eps}
\includegraphics[scale=0.5]{circ_tr_2.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/isocline.bmp}
%\hskip1.4in\special{bmp:c:/mybook/isocline.bmp x=4in y=3in}
\caption{Isoclines and a typical trajectory for the spring-mass system of Example~\ref{example_phcirc}.}
\label{fig_isocline}
\end{figure}
Note that in this example the trajectory slope is $m$, while the isocline slope is $-1/m$. This means that trajectories of the spring-mass system are perpendicular to the isoclines, that is, the trajectories are circles centered at the origin of the state plane.
\end{example}\markex
\begin{example}\rm
Given a second-order differential equation
\begin{equation}
\ddot{x}=-4\eta \dot{x}-4x,
\label{eq:secodietap}
\end{equation}
where $\eta$ is a parameter. We will find values of $\eta$ for which there are isoclines along which the trajectory slope and the isocline slope\index{isocline!slope} are equal. Such isoclines, if they exist, are called the {\em asymptotes}\index{asymptotes|see{isocline}}.

Let $x_1=x$ and $x_2=\dot{x}$. Then, the second-order differential equation given by~(\ref{eq:secodietap}) can be represented as a system of first-order differential equations,
\[
\left.\begin{array}{ccl}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& -4\eta x_2 - 4x_1
\end{array}\right\}
\]
The isocline equation is
\[
\frac{dx_2}{dx_1}=\frac{-4 \eta x_2 -4x_1}{x_2}=m,
\]
or equivalently
\begin{equation}
x_2=-\left(\frac{4}{m + 4 \eta} \right)x_1.
\label{eq:isoeqxtxo}
\end{equation}
A trajectory slope is equal to an isocline slope if and only if
\[
-\frac{4}{m+4 \eta}=m,
\]
that is,
\begin{equation}
m^2+4 \eta m+4=0.
\label{eq:secalwm}
\end{equation}
The discriminant\index{discriminant} of~(\ref{eq:secalwm}) is
\[
\Delta=4\sqrt{\eta ^2-1}.
\]
Hence, the condition for the existence of asymptotes is
\begin{equation}
|\eta| \ge 1.
\label{eq:conforasym}
\end{equation}
\end{example}\markex
\begin{example}\rm
Consider the system of differential equations,
\[
\left.\begin{array}{lll}
\dot{x}_1 &=& x_1+x_2+1\\
\dot{x}_2 &=& -x_2 + 2.
\end{array}\right\}
\]
\begin{itemize}
\item[(i)] Find the equation of the isocline corresponding to the trajectory slope $m=1$;
\item[(ii)] Does the given system have asymptotes? If yes, then write the equations of the asymptotes.
\end{itemize}

%\markdefn

The trajectory slope is
\[
\frac{d x_2}{d x_1}=m=\frac{-x_2+2}{x_1+x_2+1}.
\]
Hence, the equation of the isocline corresponding to the trajectory slope $m=1$ is obtained by solving
\[
\frac{-x_2+2}{x_1+x_2+1}=1.
\]
We obtain
\[
x_2=-\frac{1}{2}x_1+\frac{1}{2}.
\]
The trajectory slope is
\begin{equation}
m=\frac{-x_2+2}{x_1+x_2+1}.
\label{eq:eqformb}
\end{equation}
From~(\ref{eq:eqformb}), we obtain an equation of isoclines of the form
\[
x_2=-\frac{m}{m+1}x_1-\frac{m-2}{m+1}.
\]
An asymptote is the isocline\index{isocline!slope} for which its slope and the trajectory\index{trajectory!slope} slope are equal. Therefore, we can obtain the equations of asymptotes\index{asymptotes!equations} by solving
\[
m=-\frac{m}{m+1}.
\]
We obtain
\[
m^2+2m=m(m+2)=0.
\]
Hence, the equations of asymptotes are
\[
x_2=2\quad\mb{and}\quad x_2=-2x_1-4.
\]
\end{example}\markex
We can use a computer to draw phase portraits interactively\index{phase portrait!draw interactively}. Below is an example of MATLAB's script that was used to sketch a phase portrait of the van der Pol equation\index{van der Pol equation} shown in Figure~\ref{fig:php_pol}.
\begin{center}
\rule{3.5in}{0.01in}
\end{center}
\begin{verbatim}
t0=0;                      % initial time
tf=20;                     % final time
tspan=tf-t0;
x0=[-4 -4]';               % first initial condition
button=1;
p=4*[-1 0;1 0];            % preparation to draw internal axes
clf;plot(p(:,1),p(:,2))
hold on
plot(p(:,2),p(:,1))
axis(4*[-1 1 -1 1])
    while(button==1)

% Select a new initial condition and press the left button
% Press the right button if you wish to stop

      [t,x]=ode45(@my_xdot,tspan,x0); % invoke your own m-file
      plot(x(:,1),x(:,2))
      [x1,x2,button]=ginput(1);
      x0=[x1 x2]';
    end
\end{verbatim}
\begin{center}
\rule{3.5in}{0.01in}
\end{center}
\begin{figure}
\centering
\includegraphics[scale=0.55]{php_pol.eps}
%\vspace{3.5in}
%\centerbmp{3.97in}{3.5in}{c:/mybook/php_pol.bmp}
%\hskip1in\special{bmp:c:/mybook/php_pol.bmp x=3.97in y=3.5in}
\caption{A phase-plane portrait of the van der Pol equation for $\mu=1$.}
\label{fig:php_pol}
\end{figure}
The van der Pol equation has the form
\begin{equation}
\ddot{x}-\mu(1-x^2)\dot{x}+x=0,
\label{eq:vderpol}
\end{equation}
where $\mu$ is a parameter. We can equivalently represent the van der Pol equation as
\begin{equation}
\left.\begin{array}{ccl}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& \mu(1-x_1^2)x_2-x_1
\end{array}\right\}
\label{eq:eqrepvdpol}
\end{equation}
An inspection of the phase portrait of the van der Pol equation reveals the existence of a closed trajectory in the state plane. This solution is an example of a limit cycle\index{limit cycle!van der Pol equation}.
\begin{defn}
A limit cycle\index{limit cycle!definition of} is a closed trajectory in the state plane, such that no other closed trajectory can be found arbitrarily close to it.
\end{defn}\markdefn

Note that the differential equation modeling the spring-mass system does not have a limit cycle because the solutions are not isolated. Limit cycles can only appear in nonlinear systems; they do not exist in linear systems.
\begin{thm}[Bendixson]\index{Bendixson's theorem}
If $f_1$ and $f_2$ in~(\ref{eq:ftbyf1afo}) have continuous partial derivatives in a domain $\cal A$ bounded by a closed path $\b C$ and if $\partial f_1/\partial x_1+\partial f_2/\partial x_2$ does not vanish and is of constant sign in that domain, then~(\ref{eq:ftbyf1afo}) can have no limit cycle in that domain.
\end{thm}
{\bf Proof}\quad We prove the theorem by contradiction. Equation~(\ref{eq:ftbyf1afo}) is satisfied on any of its trajectory. Therefore, if~(\ref{eq:ftbyf1afo}) has a limit cycle then clearly
\begin{equation}
f_1dx_2-f_2dx_1=0.
\label{eq:rotzerdiv}
\end{equation}
We now apply Green's theorem\index{Green's theorem}, stated on page~\pageref{thm:green}, to the vector field
\[
\left[\begin{array}{cc}
-f_2 & f_1
\end{array}\right]^T
\]
to obtain
\begin{equation}
\int_{\b C}f_1\,dx_2-f_2dx_1=\int\!\!\!\int_{\cal A}\left( \frac{\partial f_1}{\partial x_1}+\frac{\partial f_2}{\partial x_2}\right)dx_1dx_2.
\label{eq:divthviag}
\end{equation}
By~(\ref{eq:rotzerdiv}), the right-hand side of~(\ref{eq:divthviag}) must vanish, which contradicts the hypothesis. Thus, if $\partial f_1/\partial x_1+\partial f_2/\partial x_2$ does not vanish and is of constant sign in a domain $\cal A$ bounded by a closed path $\b C$, then no limit cycle can exist in $\cal A$. The proof of the theorem is complete.\\
\hspace*{5.5in}$\Box$
\begin{example}\rm
We apply Bendixson's theorem to the van der Pol equation\index{van der Pol equation!Bendixson's theorem} represented by~(\ref{eq:eqrepvdpol}). We have
\[
f_1=x_2,\quad\mb{and}\quad f_2=\mu(1-x_1^2)x_2-x_1.
\]
We calculate $\partial f_1/\partial x_1+\partial f_2/\partial x_2$ to get
\[
\partial f_1/\partial x_1+\partial f_2/\partial x_2=\mu(1-x_1^2).
\]
It follows from the Bendixson theorem that there are no limit cycles in the region $|x_1|< 1$ because $\partial f_1/\partial x_1+\partial f_2/\partial x_2>0$ in that region.
\end{example}\markex
Bendixson's theorem gives us a necessary condition for a closed curve to be a limit cycle. It is useful for the purpose of establishing the nonexistence of a limit cycle. In general, it is very difficult to establish the existence of a limit cycle\index{limit cycle!existence of}. Sufficient conditions for the existence of a limit cycle are discussed, for example, by Hochstadt~\cite[Chapter 7]{hochstadt75}. A method for predicting limit cycles, called the describing function method\index{describing function method!limit cycle existence}, is presented in Section~\ref{sec:descfunme}.

\section{Numerical Techniques}\label{section:numeric_techniq}

In this section, we discuss numerical techniques\index{numerical techniques!solving differential equations} for solving differential equations. We begin our presentation of numerical techniques that allow us to determine approximate solutions\index{numerical techniques!approximate solutions} to the given system of differential equations with a method based on a Taylor series expansion\index{numerical techniques!Taylor series expansion}. Much of this section presentation is based on the Class Notes of the late Professor~V.~B.~Haas~\cite{haas80}.

\subsection{The Method of Taylor Series}
The basis for many numerical methods\index{numerical methods} for solving differential equations is the Taylor formula\index{Taylor formula!numerical methods}. Suppose we are given a state-space model of the form
\begin{equation}
\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t)), \quad \mb{\bm$x$}(t_0)=\mb{\bm$x$}_0.
\label{eq:taydiffe}
\end{equation}
The method, we discuss, depends on the fact that if $\mb{\bm$f$}(t,\mb{\bm$x$},\mb{\bm$u$})$ is differentiable with respect to its $n+m+1$ variables $x_1,\ldots, x_n,u_1,\ldots ,u_m, t$ in some $(n+m+1)$-dimensional region ${\cal D}$, and if $\mb{\bm$x$}(t)$ is a solution of~(\ref{eq:taydiffe}), where $(t_0,\mb{\bm$x$}_0, \b u(t_0))$ lies in ${\cal D}$, then we can expand $\mb{\bm$x$}(t)$ into a Taylor series\index{Taylor series} to obtain
\begin{eqnarray*}
\mb{\bm$x$}(t)&=&\mb{\bm$x$}_0+\dot{\mb{\bm$ x$}}(t_0)(t-t_0)+\ddot{\mb{\bm$ x$}}(t_0)\frac{(t-t_0)^2}{2}+\cdots \\
&=&\mb{\bm$x$}_0+\mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$ u$}(t_0))(t-t_0)\\
{}&{}&+\left[\frac{\partial \mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0))}{\partial t}+\sum_{i=1}^{n} \frac{\partial \mb{\bm$f$}(t_0,\mb{\bm$ x$}_0,\mb{\bm$u$}(t_0))}{\partial x_i} \dot{x}_i \right.\\
{}&{}&+ \left.\sum_{j=1}^{m} \frac{\partial\mb{\bm$ f$}(t_0,\mb{\bm$ x$}_0,\mb{\bm$ u$}(t_0))}{\partial u_j} \dot{u}_j \right] \frac{(t-t_0)^2}{2}+\cdots
\end{eqnarray*}
If we stop after the $q$-th term of the Taylor series solution\index{Taylor series solution}, then the remainder\index{Taylor series!remainder}\index{remainder!Taylor series} after $q$ terms of a Taylor series gives us an estimate of the error\index{Taylor series!estimate of the error}\index{estimate of the error!Taylor series}. We can expect the Taylor series solution to give a fairly accurate approximation of the solution of the state-space equations only in a small interval about $t_0$. Thus, if $\Delta t = t - t_0$ is some small time interval such that the Taylor series converges for $t \in \left[\begin{array}{cc}
t_0,&t_0 + \Delta t
\end{array}\right]$, we can use the expansion
\[
\mb{\bm$x$}(t)=\mb{\bm$x$}_0+\frac{d\mb{\bm$ x$}(t_0)}{dt}(t-t_0)+\frac{d^2\mb{\bm$ x$}(t_0)}{dt^2}\frac{(t-t_0)^2}{2}+\cdots
\]
to evaluate $\mb{\bm$ x$}(t_1)=\mb{\bm$ x$}(t_0+\Delta t)$. Then, we can start over again solving $\dot{\mb{\bm$ x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t))$ with the new initial condition $\mb{\bm$x$}(t_1)$.
We illustrate the method with simple examples.
\begin{example} \rm
Consider the first-order differential equation\index{first-order differential equation}
\[
\frac{dx}{dt}=t^2-x^2,\quad x(1)=1.
\]
Expanding $x=x(t)$ into a Taylor series about $x(1)=1$ yields
\[
x(t)=x(1)+\frac{dx(1)}{dt}(t-1)+\frac{d^2x(1)}{dt^2} \frac{(t-1)^2}{2}+\frac{d^3x(1)}{dt^3} \frac{(t-1)^3}{3!}+\frac{d^4x(1)}{dt^4} \frac{(t-1)^4}{4!}\cdots
\]
Because $\frac{dx}{dt}=t^2-x^2$, we obtain by successive differentiations
\begin{eqnarray*}
\frac{d^2x}{dt^2}&=&2t-2x\frac{dx}{dt}, \\
\frac{d^3x}{dt^3}&=&2-2\left( \frac{dx}{dt} \right) ^2-2x\frac{d^2x}{dt^2}, \\
\frac{d^4x}{dt^4}&=&-6\frac{dx}{dt}\frac{d^2x}{dt^2}-2x\frac{d^3x}{dt^3},\\
& \vdots &
\end{eqnarray*}
Evaluating the above at $t=1$, yields
\begin{eqnarray*}
\frac{dx(1)}{dt}&=&0,\\
\frac{d^2x(1)}{dt^2}&=&2,\\
\frac{d^3x(1)}{dt^3}&=&-2,\\
\frac{d^4x(1)}{dt^4}&=&4,\\
& \vdots &
\end{eqnarray*}
Hence,
\[
x(t)=1+(t-1)^2-\frac{(t-1)^3}{3}+\frac{(t-1)^4}{6}+\cdots .
\]
\end{example}
\markex
\begin{example} \rm
In this example, we consider the second-order differential equation\index{second-order differential equation}
\[
\frac{d^2x}{dt^2}+3\frac{dx}{dt}+6tx=0,
\]
subject to the initial conditions $x(0)=1$, and $\frac{dx(0)}{dt}=1$. Using $\frac{d^2x}{dt^2}=-3\frac{dx}{dt}-6tx$, we obtain by successive differentiation\index{successive differentiation}
\begin{eqnarray*}
\frac{d^3x}{dt^3}&=&-3\frac{d^2x}{dt^2}-6t\frac{dx}{dt}-6x,\\
\frac{d^4x}{dt^4}&=&-3\frac{d^3x}{dt^3}-6t\frac{d^2x}{dt^2}-12\frac{dx}{dt},\\
\frac{d^5x}{dt^5}&=&-3\frac{d^4x}{dt^4}-6t\frac{d^3x}{dt^3}-18\frac{d^2x}{dt^2},\\
& \vdots &
\end{eqnarray*}
Taking into account the initial conditions, we evaluate all coefficients at $t=0$ to get
\begin{eqnarray*}
\frac{d^2x(0)}{dt^2}&=&-3,\\
\frac{d^3x(0)}{dt^3}&=&3,\\
\frac{d^4x(0)}{dt^4}&=&-21,\\
\frac{d^5x(0)}{dt^5}&=&117,\\
& \vdots &
\end{eqnarray*}
Hence,
\begin{eqnarray*}
x(t)&=&x(0)+\frac{dx(0)}{dt}t+\frac{d^2x(0)}{dt^2} \frac{t^2}{2}+\frac{d^3x(0)}{dt^3} \frac{t^3}{3!}+\frac{d^4x(0)}{dt^4} \frac{t^4}{4!}+\frac{d^5x(0)}{dt^5} \frac{t^5}{5!}+\cdots\\
&=&1+t-\frac{3}{2}t^2+\frac{1}{2}t^3-\frac{7}{8}t^4+\frac{117}{120}t^5+\cdots .
\end{eqnarray*}
\end{example}
\markex
We now illustrate the method of Taylor series when applied to a system of first-order differential equations.
\begin{example} \rm
We consider the following system of two first-order equations
\[
\left.\begin{array}{lll}
\frac{dx_1}{dt}&=&f_1(t,x_1,x_2)\\
\frac{dx_2}{dt}&=&f_2(t,x_1,x_2)
\end{array}\right\}
\]
The initial conditions at $t=t_0$ are: $x_1(t_0)=x_{10}$, $x_2(t_0)=x_{20}$. Having the initial conditions, we can determine the values $\frac{dx_1(t_0)}{dt}$ and $\frac{dx_2(t_0)}{dt}$ from the original equations. We then differentiate $\frac{dx_1}{dt}$ and $\frac{dx_2}{dt}$ to obtain
\begin{eqnarray*}
\frac{d^2x_1}{dt^2}&=&\frac{d}{dt}f_1\left(t,x_1,x_2,\frac{dx_1}{dt},\frac{dx_2}{dt}\right)\\
\frac{d^2x_2}{dt^2}&=&\frac{d}{dt}f_2\left(t,x_1,x_2,\frac{dx_1}{dt},\frac{dx_2}{dt}\right).
\end{eqnarray*}
Now, $\frac{d^2x_1}{dt^2}$ and $\frac{d^2x_2}{dt^2}$ can be calculated by substituting known quantities on the right-hand sides of the above equations. By successive differentiation and substitution, we can determine $\frac{d^3x_1}{dt^3}$ and $\frac{d^3x_2}{dt^3}$ as well as higher derivatives at $t=t_0$. The solution to the given system is then computed as
\begin{eqnarray*}
x_1(t)&=&x_{10}+\frac{dx_1(t_0)}{dt}(t-t_0)+\frac{d^2x_1(t_0)}{dt^2}\frac{\left(t-t_0\right)^2}{2}+\frac{d^3x_1(t_0)}{dt^3}\frac{\left(t-t_0\right)^3}{3!}+\cdots \\
x_2(t)&=&x_{20}+\frac{dx_2(t_0)}{dt}(t-t_0)+\frac{d^2x_2(t_0)}{dt^2}\frac{\left(t-t_0\right)^2}{2}+\frac{d^3x_2(t_0)}{dt^3}\frac{\left(t-t_0\right)^3}{3!}+\cdots
\end{eqnarray*}
\end{example}
\markex
The Taylor series method can be programmed on a computer. However, this is not very efficient method from numerical point of view. In the following, we discuss other methods for numerical solution of the state equations that are related to the method of Taylor series\index{Taylor series}. We first present two simple numerical techniques known as the forward and backward Euler methods.

\subsection{Forward Euler Method}
Suppose that we wish to obtain the solution of the modeling equation
\[
\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t))
\]
in the time interval $\left[\begin{array}{cc}
t_0,&t_f
\end{array}\right]$, subject to the initial condition $\mb{\bm$x$}(t_0)=\mb{\bm$x$}_0$, and the input vector $\mb{\bm$u$}(t)$. We divide the interval $\left[\begin{array}{cc}
t_0,&t_f
\end{array}\right]$ into $N$ equal subintervals of width
\[
h=\frac{t_f-t_0}{N}.
\]
We call $h$ the {\em step length}\index{step length}. We set
\[
t_k=t_0+kh.
\]
To proceed further, we approximate the derivative\index{approximate the derivative!forward Euler} $\frac{d\mb{\bm$x$}(t)}{dt}$ at time $t_k$ by
\[
\frac{d\mb{\bm$x$}(t_k)}{dt}\approx \frac{\mb{\bm$x$}(t_{k+1})-\mb{\bm$x$}(t_k)}{h},
\]
and write
\[
\frac{\mb{\bm$x$}(t_{k+1})-\mb{\bm$x$}(t_k)}{h}=\mb{\bm$ f$}(t_k,\mb{\bm$x$}(t_k),\mb{\bm$u$}(t_k)).
\]
Hence, for $k=1,2,\ldots ,N$, we have
\begin{equation}
\fbox{$\displaystyle \mb{\bm$x$}(t_{k+1})=\mb{\bm$x$}(t_k)+h\mb{\bm$f$}(t_k,\mb{\bm$ x$}(t_k),\mb{\bm$u$}(t_k))$}
\label{eq:forward_Euler_meth}
\end{equation}
The above is known as the {\em forward Euler\/}\index{forward Euler} algorithm. We see that if $h$ is sufficiently small, we can approximately determine the state at time $t_1=t_0+h$ from the initial condition $\mb{\bm$x$}_0$ to get
\[
\mb{\bm$x$}(t_1)=\mb{\bm$x$}_0+h\mb{\bm$f$}(t_0,\mb{\bm$ x$}_0,\mb{\bm$u$}(t_0)).
\]
Once we have determined the approximate solution at time $t_1$, we can determine the approximate solution at time $t_2=t_1+h$, and so on.
The forward Euler method can also be arrived at when instead of considering the differential equation $\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t))$, we rather start with its equivalent integral representation
\begin{equation}
\mb{\bm$x$}(t)=\mb{\bm$x$}_0+\int_{t_0}^{t}\mb{\bm$f$}(\tau,\mb{\bm$x$}(\tau),\mb{\bm$u$}(\tau))d \tau.
\label{eq:eqintrepz}
\end{equation}
We then use the {\em rectangular rule\/}\index{forward Euler!rectangular rule}\index{rectangular rule!forward Euler} for numerical integration that can be stated as follows. The area under the curve $z=g(t)$ between $t=t_k$ and $t=t_k+h$ is approximately equal to the area of the rectangle $ABCD$ as shown in Figure~\ref{fig_rectint}, that is,
\begin{equation}
\int_{t_k}^{t_k+h}g(t)dt \approx hg(t_k).
\label{eq:recrulofint}
\end{equation}
\begin{figure}
\centering
\includegraphics[scale=0.55]{rectint.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/rectint.bmp}
%\hskip1.25in\special{bmp:c:/mybook/rectint.bmp x=4in y=3in}
\caption{Illustration of the rectangular rule of integration.}
\label{fig_rectint}
\end{figure}
Applying the rectangular rule of integration to~(\ref{eq:eqintrepz}), we obtain the forward Euler method.

\begin{example}\label{ex:handle_g_and_For_E}\rm  We combine basic MATLAB's handle graphics\index{MATLAB's handle graphics} commands with the forward Euler method to write a simple program that animates the motion\index{motion!animation} of a simple pendulum\index{simple pendulum!motion animation}. First, we present a few basic definitions from handle graphics. We start with a description of a handle\index{handle}---it is a floating-point\index{handle!floating-point number} number that MATLAB\index{MATLAB} assigns to every object in the figure window such as a line, tex, label, axes, figure, etc. Thus a handle can be viewed as an object's ID\index{handle!object's ID}\index{MATLAB's handle!object's ID}.

Each object has different properties. If we wish to change a property\index{property!handle}, we need to get hold of the object in question. We use object's handle to get hold of it and then change the object's property\index{property!change} using the following commands,
\begin{verbatim}
get(handle,'PropertyName')
set(handle,'PropertyName','PropertyValue')
\end{verbatim}
In our animation of the motion of a pendulum\index{simple pendulum!motion animation}, we use the rotation matrix\index{rotation matrix} that we discuss next. Suppose that we have a point mass on a circle of radius $r$ as shown in Figure~\ref{fig:Rota_anima}. The point mass is connected to the circle center; we call it hinge. The line segment, bar, connecting the hinge and the mass forms the angle $\alpha$. Therefore, the point mass has the coordinates
\[
\b x_{present}=\left[\begin{array}{c}
r\cos (\alpha)\\
r\sin (\alpha)
\end{array}\right].
\]
\begin{figure}
\centering
\includegraphics[scale=0.35]{Rota_anima.eps}
\caption{Rotating a point mass from $\alpha$ to $\alpha+\theta$.}
\label{fig:Rota_anima}
\end{figure}
Suppose that the bar was rotated by an angle $\theta$ from its initial displacement angle $\alpha$. The new angle of the bar, with the positive $x$-axis, is then $\alpha + \theta$. We use the following well-known trigonometric identities\index{trigonometric identities}:
\begin{eqnarray*}
\cos(\alpha + \theta) &=& \cos (\alpha)\cos (\theta) - \sin (\alpha) \sin (\theta)\\
\sin(\alpha + \theta) &=& \sin (\alpha)\cos (\theta) + \cos (\alpha) \sin (\theta),
\end{eqnarray*}
to obtain the mass new coordinates,
\begin{eqnarray*}
\b x_{next} &=& \left[\begin{array}{cc}
\cos(\theta) & -\sin(\theta)\\
\sin(\theta) & \cos(\theta)
\end{array}\right]\b x_{present}\\
&=& \left[\begin{array}{cc}
\cos(\theta) & -\sin(\theta)\\
\sin(\theta) & \cos(\theta)
\end{array}\right]\left[\begin{array}{c}
r\cos (\alpha)\\
r\sin (\alpha)
\end{array}\right].
\end{eqnarray*}
Equipped with the above background results, we are ready to animate the motion of the simple pendulum from Example~\ref{ex:simP_pend_Lagrange}. Let
\[
x_1=\theta\quad\mb{and}\quad x_2=\dot{\theta}.
\]
Then, we represent the second-order differential equation~\eqref{eq:sec_ord_mode_simp_pen} as a system of first-order differential equations,
\begin{eqnarray*}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& -5\sin(x_1).
\end{eqnarray*}
We approximate the above system of differential equations using Euler's forward\index{Euler's forward} method,
\begin{eqnarray*}
\frac{x_{1next}-x_{1present}}{h} &=& x_{2present}\\
\frac{x_{2next}-x_{2present}}{h} &=& -5\sin(x_{1present}),
\end{eqnarray*}
where $h$ is the step size. Combining the above with the handle graphics gives the following simple MATLAB program\index{MATLAB program!animation} that animates the simple pendulum motion. This code is based on the code of Example~3 in Pratap~\cite[p.~207]{pratap06}.
\begin{verbatim}
theta=input('Initial angular displacement of the pendulum in deg= ');
theta=theta*pi/180;           % convert from degrees to radians
tfinal=input('Final time tfinal = ');
data_init=[0 -2]';            % pendulum at rest
R=[cos(theta) -sin(theta);sin(theta) cos(theta)];    % rotation matrix
data=R*data_init;             % initial pendulum position
axis(2*[-1 1 -1 1])
axis('equal')
axis('off')
bar=line('xdata',[0 data(1)],'ydata',[0 data(2)]','linewidth',7,...
    'erase','xor');
mass=line('xdata',data(1),'ydata',data(2),'linestyle','.',...
    'markersize',50,'erase','xor');
hinge=line('xdata',0,'ydata',0,'linestyle','o','markersize',12,...
    'erase','xor');
h=0.0025;     % step-size for solving differential
t=0;           % initial time
thetadot=0;    % initial angular speed

% solve the differential eqns using the Euler's method
while(t<tfinal)
   t=t+h;
   theta = theta + thetadot*h;
   thetadot=thetadot - 5*sin(theta)*h;
   R=[cos(theta) -sin(theta);sin(theta) cos(theta)];
   datanew=R*data_init;

% change the property values of the bar and hinge objects
   set(bar,'xdata',[0 datanew(1)],'ydata',[0 datanew(2)]);
   set(mass,'xdata',datanew(1),'ydata',datanew(2));
   set(hinge,'xdata',0,'ydata',0)
   drawnow;
end
\end{verbatim}
\end{example}\markex

\subsection{Backward Euler Method}

The {\em backward Euler\/}\index{backward Euler} method differs from the forward Euler method in the way we approximate the derivative\index{approximate the derivative!backward Euler},
\[
\frac{d\mb{\bm$x$}(t_k)}{dt}\approx \frac{\mb{\bm$x$}(t_{k})-\mb{\bm$x$}(t_{k-1})}{h},
\]
This approximation gives the backward Euler formula,
\[
\fbox{$\displaystyle \mb{\bm$x$}(t_{k+1})=\mb{\bm$x$}(t_k)+h\mb{\bm$f$}(t_{k+1},\mb{\bm$ x$}(t_{k+1}),\mb{\bm$u$}(t_{k+1}))$}
\]
In the above formula, $\b x(t_{k+1})$ is a function of itself. For this reason, we say that the backward Euler method is an implicit integration algorithm. Implicit integration\index{implicit integration!backward Euler} algorithms require additional computation to solve for $\b x(t_{k+1})$. For further discussion of this issue, we refer to Parker and Chua~\cite{parker89}.

Euler's methods are rather elementary, and thus, not as accurate as some of the more sophisticated techniques that we are going to discuss next.

\subsection{Predictor-Corrector Method}
The predictor-corrector\index{predictor-corrector} method, also known as Heun's method\index{Heun's method|see{predictor-corrector}}, is based on the trapezoidal rule\index{trapezoidal rule!predictor-corrector} for numerical integration. The trapezoidal rule states that the area under the curve $z=g(t)$ between $t=t_k$ and $t=t_k+h$ is approximately equal to the area of the trapezoid $ABCD$ as depicted in Figure~\ref{fig_trapint}, that is,
\begin{equation}
\int_{t_k}^{t_k+h}g(t)dt\approx h \frac{g(t_k)+g(t_k+h)}{2}.
\label{eq:trapruofint}
\end{equation}
\begin{figure}
\centering
\includegraphics[scale=0.55]{trapint.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/trapint.bmp}
%\hskip1.25in\special{bmp:c:/mybook/trapint.bmp x=4in y=3in}
\caption{The trapezoidal rule of integration.}
\label{fig_trapint}
\end{figure}
Suppose we wish to find $\mb{\bm$x$}(t)$ at $t=t_{1} > t_0$. Applying the trapezoidal rule\index{trapezoidal rule!of integration} of integration to~(\ref{eq:taydiffe}), we obtain
\[
\mb{\bm$x$}(t_1)=\mb{\bm$x$}_0+\frac{h}{2}\left(\mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0)) + \mb{\bm$f$}(t_1,\mb{\bm$x$}(t_1),\mb{\bm$u$}(t_1)) \right).
\]
The above is an implicit algorithm and it may be viewed as a merging of the forward and backward Euler algorithms. If the function $\mb{\bm$f$}$ is nonlinear, then we will, in general, not be able to solve the above equation for $\mb{\bm$x$}(t_1)$ exactly. We attempt to obtain $\mb{\bm$x$}(t_1)$ iteratively. Let
\begin{eqnarray*}
h&=&t_{k+1}-t_{k},\\
\mb{\bm$x$}^{*}_1 &=& \mb{\bm$x$}_0+h\mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0)),\\
\mb{\bm$m$}_0 &=& \mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0)),\\
\mb{\bm$m$}_1 &=& \mb{\bm$f$}(t_1,\mb{\bm$x$}^{*}_1,\mb{\bm$u$}(t_1)).
\end{eqnarray*}
A first approximation to $\mb{\bm$x$}(t_1)$ is
\[
\fbox{$\displaystyle \mb{\bm$x$}^{*}_1 = \mb{\bm$x$}_0 + h\mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0))$}
\]
and is called a {\em predictor\/}\index{predictor}. The next approximation is
\[
\fbox{$\displaystyle \mb{\bm$x$}(t_1)=\mb{\bm$x$}_0+\frac{h}{2}\left(\mb{\bm$m$}_0 + \mb{\bm$m$}_1 \right)$}
\]
and is called a {\em corrector}\index{corrector}.


Suppose now that we wish to find $\b x(t)$ at $t=t_f>t_0$. We first divide the interval $\left[\begin{array}{cc}
t_0,&t_{f}
\end{array} \right]$ into $N$ equal subintervals using evenly spaced subdivision points to obtain
\[
t_0 < t_1 < t_2 < \ldots < t_{N}=t_{f}.
\]
Once we have found $\mb{\bm$x$}(t_1)$, we can apply the same procedure to find $\mb{\bm$x$}(t_2)$, and so on.
\begin{example}
\label{ex:pred_cor}
\rm
We use the predictor-corrector\index{predictor-corrector} method to find $x(1)$, where the step length $h=0.5$, for the following differential equation
\begin{equation}
\frac{dx}{dt}=tx^2-x,\quad x(0)=1.
\label{eq:predcorrex}
\end{equation}
The calculations are summarized in Table~\ref{tab:pred_cor}.
\begin{table}
\caption{Solving the differential equation~(\ref{eq:predcorrex}) using the predictor-corrector method.}
\begin{center}
\begin{tabular}{cccccccc}
$t_0$ & $x_0$ & $t_1$ & $m_0$ & $x^{*}_1$ & $m_1$ & $m=\frac{m_0+m_1}{2}$ & $\Delta x$\\ \hline
$0$ & $1$ & $0.5$ & $-1.0$ & $0.5$ & $-0.375$ & $-0.688$ & $-0.344$\\ \hline
$0.5$ & $0.656$ & $1.0$ & $-0.441$ & $0.435$ & $-0.246$ & $-0.344$ & $-0.172$\\ \hline
$1.0$ & $0.484$ &{}&{}&{}&{}&{}&{}\\\hline
\end{tabular}
\end{center}
\label{tab:pred_cor}
\end{table}
\end{example}
\markex
We could use $\mb{\bm$x$}(t_1)$ obtained from the corrector formula to construct a next approximation to $\mb{\bm$x$}(t_1)$. In general, we can use the formula
\[
\mb{\bm$x$}^{(k)}(t_1)=\mb{\bm$x$}_0 + \frac{h}{2} \left( \mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0)) + \mb{\bm$f$}(t_1,\mb{\bm$x$}^{(k-1)}(t_1),\mb{\bm$u$}(t_1) \right)
\]
to evaluate $\mb{\bm$x$}(t_1)$ until two successive iterations\index{successive iterations!predictor-corrector} agree to the desired accuracy. This general algorithm is referred to as a second-order predictor-corrector\index{predictor-corrector!second-order method} method.

An $n$-dimensional\index{predictor-corrector!n-dimensional} system of first-order differential equations can be worked analogously. We illustrate the procedure for $n=2$, where
\[
\left.\begin{array}{lll}
\frac{dx}{dt} &=& f(t,x,y), \quad x(t_0)=x_0\\
\frac{dy}{dt} &=& g(t,x,y), \quad y(t_0)=y_0
\end{array}\right\}
\]
Let
\begin{eqnarray*}
m_0 &=& f(t_0,x_0,y_0),\\
n_0 &=& g(t_0,x_0,y_0).
\end{eqnarray*}
Then, the predictor\index{predictor!n-dimensional} equations are
\begin{eqnarray*}
x^{*}_1 &=& x_0 +m_0\,h,\\
y^{*}_1 &=& y_0 + n_0\,h.
\end{eqnarray*}
Let
\begin{eqnarray*}
m_1 &=& f(t_1,x^{*}_1,y^{*}_1),\\
n_1 &=& g(t_1,x^{*}_1,y^{*}_1).
\end{eqnarray*}
Then, the corrector\index{corrector!n-dimensional} equations are
\[
\left.\begin{array}{lll}
x(t_1) &=& x_0 + \frac{m_0+m_1}{2}\,h\\
y(t_1) &=& y_0+\frac{n_0+n_1}{2}\,h
\end{array}\right\}
\]
Using the above, we can generate the corresponding formulas for $t_k$.

We now show that the predictor-corrector method yields a solution that agrees with the Taylor series solution through terms of degree two. We consider the scalar case, where $f$ is assumed to be differentiable with respect to its arguments and $\partial^2 f/\partial t \partial x = \partial^2 f/\partial x \partial t$. First observe that $m_1=f\left(t_1,x^{*}_1 \right)$ can be written as
\begin{eqnarray*}
f\left(t_1,x^{*}_1 \right) &=& f \left( t_0 + h,x_0+m_0h \right)\\
&\approx& f\left(t_0,x_0 \right) + \left[\begin{array}{cc}
f_{t} & f_{x}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right]+ \frac{1}{2} \left[\begin{array}{cc}
h & m_0h
\end{array}\right] \left[\begin{array}{cc}
f_{tt} & f_{tx}\\
f_{xt} & f_{xx}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right],
\end{eqnarray*}
where $f_t=\partial f/\partial t$, $f_x=\partial f/\partial x$, and so on. Substituting the above into the corrector equation yields
\begin{eqnarray*}
x(t_1) &=& x_0 + \frac{h}{2}\left( m_0 + m_1 \right)\\
&=& x_0 + \frac{h}{2}\left( m_0 + m_0\right)+ \frac{h}{2}\left(\left[\begin{array}{cc}
f_{t} & f_{x}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right]\right.\\
&{}&\quad +\left.\frac{1}{2} \left[\begin{array}{cc}
h & m_0h
\end{array}\right] \left[\begin{array}{cc}
f_{tt} & f_{tx}\\
f_{xt} & f_{xx}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right]\right)\\
&=& x_0+m_0h+\left(f_{t}+f_{x}m_0\right)\frac{h^2}{2}+\left(f_{tt}+2f_{tx}m_0 + f_{xx}m_0^2 \right)\frac{h^3}{4}.
\end{eqnarray*}
The Taylor series\index{Taylor series!predictor-corrector} expansion, on the other hand, gives
\begin{eqnarray*}
x\left(t_1 \right) &=& x_0 + \dot{x}\,h + \ddot{x}\, \frac{h^2}{2} + x^{(3)}\frac{h^3}{6} + \cdots\\
&=& x_0 + f(t_0,x_0) \,h + \left(f_{t} + f_{x}\dot{x} \right) \frac{h^2}{2}\\
{}&{}&+ \left(f_{tt}+2f_{tx} \dot{x} +f_{xx}\left(\dot{x}\right)^2+f_{x}\ddot{x} \right)\frac{h^3}{6} + \cdots,
\end{eqnarray*}
where $f_{t}$, $f_{x}$, $\dot{x}$, etc.\ are all evaluated at $t=t_0$. Comparing the right-hand sides of the last two equations, we see that the solution obtained using the predictor-corrector method agrees with the Taylor series expansion of the true solution through terms of degree two.

\subsection{Runge's Method}
Runge's method\index{Runge's method} is an improvement over the predictor-corrector method in the sense that it agrees with the Taylor series solution through terms of degree three. The method depends upon integration by means of Simpson's rule\index{Simpson's rule!Runge's method}, which says that the area under the curve $x=g(t)$ between $t_0$ and $t_0+h$ is given approximately by
\begin{equation}
\int_{t_0}^{t_0+h}g(t)dt \approx \frac{h}{6}\left(x_0+4x_1+x_2 \right),
\label{eq:simpsintru}
\end{equation}
where $x_0=g(t_0)$, $x_1=g\left( t_0 + h/2 \right)$, and $x_2=g(t_0+h)$. The above formula is derived by passing a parabola\index{parabola!through three points} through the three points $(t_0,x_0)$, $\left(t_0+h/2,x_1\right)$, and $(t_0+h,x_2)$, and assuming that the area under the given curve is approximately equal to the area under the parabola as it is illustrated in Figure~\ref{fig:simpson}.
\begin{figure}
\centering
\includegraphics[scale=0.55]{simpson.eps}
%\vskip3in
%\centerbmp{4in}{3in}{c:/mybook/simpson.bmp}
%\hskip1.25in\special{bmp:c:/mybook/simpson.bmp x=4in y=3in}
\caption{Integration using Simpson's rule.}
\label{fig:simpson}
\end{figure}
To verify Simpson's formula\index{Simpson's formula|see{Simpson's rule}}, we translate the axes so that in the new coordinates $\left( \tilde{t},\tilde{x} \right)$ we have $t=t_0+h/2+\tilde{t}$, $x=\tilde{x}$. In the new coordinates the three points become
\[
\left(-h/2,x_0 \right), \quad \left(0,x_1 \right),\quad \left(h/2,x_2 \right).
\]
The parabola through the above points is described by the equation
\[
\tilde{x}=a\tilde{t}^2 + b\tilde{t} + c.
\]
Hence
\begin{eqnarray*}
x_0 &=& a\frac{h^2}{4}-b\frac{h}{2}+c,\\
x_1 &=& c,\\
x_2 &=& a\frac{h^2}{4}+b\frac{h}{2}+c.
\end{eqnarray*}
Using the above three equations, we calculate
\begin{eqnarray*}
c &=& x_1,\\
b &=& \frac{x_2 -x_0}{h},\\
a &=& \frac{2}{h^2}\left( x_0 -2x_1 + x_2 \right).
\end{eqnarray*}
Then, the area under the parabola is
\begin{eqnarray*}
\int_{-h/2}^{h/2}\left( a\tilde{t}^2 + b\tilde{t} +c \right) d\tilde{t}
&=& \int_{-h/2}^{h/2}\left(\frac{2}{h^2}\left( x_0 -2x_1 + x_2 \right)\tilde{t}^2 + \frac{x_2 -x_0}{h}\tilde{t} +x_1 \right)d\tilde{t}\\
&=& \frac{h}{6}\left( x_0 +4x_1 + x_2 \right).
\end{eqnarray*}

If we are given the differential equation
\[
\frac{dx}{dt}=f(t,x),\quad x(t_0)=x_0,
\]
we integrate it to obtain
\[
x(t_1)=x(t_0)+\int_{t_0}^{t_1}f(s,x(s))ds.
\]
Then, we approximate the integral
\[
\int_{t_0}^{t_1}f(s,x(s))ds
\]
by
\[
\frac{h}{6}\left( m_0 + 4m_1 +m_3 \right),
\]
where
\begin{eqnarray*}
h &=& t_1-t_0,\\
m_0 &=& f(t_0,x_0),\\
m_1 &=& f \left( t_0+\frac{h}{2},x_0+m_0\frac{h}{2} \right),\\
m_2 &=& f\left( t_0+h,x_0+m_0h \right),\\
m_3 &=& f\left( t_0+h,x_0 + m_2h \right).
\end{eqnarray*}
Thus, we may write
\[
x(t_1)=x_0+ \Delta x,
\]
where $\Delta x = mh$, and $m=\frac{1}{6}\left( m_0 + 4m_1 + m_3\right)$.
\begin{example} \rm
\label{ex:runge}
We use Runge's method\index{Runge's method} to find $x(1)$, where the step length $h=0.5$, for the differential equation~(\ref{eq:predcorrex}). We found $x(1)$ using the predictor-corrector method in Example~\ref{ex:pred_cor}. The calculations for Runge's method are summarized in Table~\ref{tab:runge}.
\begin{table}
\caption{Solving the differential equation~(\ref{eq:predcorrex}) using Runge's method.}
\begin{center}
\begin{tabular}{cccccccc}
$t_0$ & $x_0$ & $m_0$ & $m_1$ & $m_2$ & $m_3$ & $m$ & $\Delta x$\\ \hline
$0$ & $1$ & $-1$ & $-0.610$ & $-0.375$ & $-0.482$ & $-0.654$ & $-0.327$\\ \hline
$0.5$ & $0.673$ & $-0.446$ & $-0.325$ & $-0.248$ & $-0.248$ & $-0.332$ & $-0.166$\\ \hline
$1.0$ & $0.507$ &{}&{}&{}&{}&{}&{}\\\hline
\end{tabular}
\end{center}
\label{tab:runge}
\end{table}
\end{example}
\markex
Runge's method can, of course, be used for systems of first-order differential equations\index{Runge's method!first-order differential equations}. If, for example,
\begin{eqnarray*}
\frac{dx}{dt} &=& f(t,x,y),\quad x(t_0)=x_0,\\
\frac{dy}{dt} &=& g(t,x,y),\quad y(t_0)=y_0,
\end{eqnarray*}
then
\begin{eqnarray*}
x(t_1) &=& x_0 + \Delta x,\\
y(t_1) &=& y_0 + \Delta y,
\end{eqnarray*}
where
\begin{eqnarray*}
\Delta x &=& \frac{h}{6} \left( m_0 + 4m_1 + m_3 \right),\\
\Delta y &=& \frac{h}{6} \left( n_0 + 4n_1 + n_3 \right).
\end{eqnarray*}
In the above expressions,
\begin{eqnarray*}
m_0 &=& f(t_0,x_0,y_0),\\
m_1 &=& f\left( t_0 + h/2,x_0+m_0\,h/2,y_0+n_0\,h/2 \right),\\
m_2 &=& f \left( t_0+h,x_0+m_0h,y_0+n_0h \right),\\
m_3 &=& f \left( t_0 +h,x_0+m_2h,y_0+n_2h \right),\\
n_0 &=& g(t_0,x_0,y_0),\\
n_1 &=& g\left( t_0 + h/2,x_0+m_0\,h/2,y_0+n_0\,h/2 \right),\\
n_2 &=& g\left( t_0+h,x_0+m_0h,y_0+n_0h \right),\\
n_3 &=& g \left( t_0 +h,x_0+m_2h,y_0+n_2h \right).
\end{eqnarray*}
\subsection{Runge-Kutta Method}
This method is an improvement over the Runge method, and agrees with the Taylor series solution through terms of degree four---for this reason it is sometimes called the fourth-order Runge-Kutta algorithm\index{Runge-Kutta algorithm|see{Runge-Kutta method}}\index{Runge-Kutta method}. It is based on the use of the formula
\begin{equation}
x(t_1)=x(t_0)+\frac{h}{6}\left( m_0+2m_1 + 2m_2 + m_3 \right),
\label{eq:runkuttintru}
\end{equation}
where
\begin{eqnarray*}
m_0 &=& f(t_0,x_0),\\
m_1 &=& f\left( t_0 + h/2, x_0 + m_0\,h/2 \right),\\
m_2 &=& f \left( t_0 + h/2, x_0 + m_1\,h/2 \right),\\
m_3 &=& f \left( t_0 + h, x_0 +h\,m_2 \right).
\end{eqnarray*}
Thus, we may write
\[
x(t_1)=x_0+ \Delta x,
\]
where $\Delta x = mh$, and $m=\frac{1}{6}\left( m_0 + 2m_1 + 2m_2 + m_3\right)$.
\begin{example} \rm
We apply the Runge-Kutta method\index{Runge-Kutta method} to find $x(1)$, where the step length $h=0.5$, for the differential equation~(\ref{eq:predcorrex}) that we used Examples~\ref{ex:pred_cor} and~\ref{ex:runge}. Our calculations for the Runge-Kutta method are summarized in Table~\ref{tab:runge_kutta}.
\begin{table}
\caption{Solving the differential equation~(\ref{eq:predcorrex}) using the Runge-Kutta method.}
\begin{center}
\begin{tabular}{cccccccc}
$t_0$ & $x_0$ & $m_0$ & $m_1$ & $m_2$ & $m_3$ & $m$ & $\Delta x$\\ \hline
$0$ & $1$ & $-1.0$ & $-0.609$ & $-0.668$ & $-0.444$ & $-0.666$ & $-0.333$\\ \hline
$0.5$ & $0.667$ & $-0.444$ & $-0.324$ & $-0.328$ & $-0.250$ & $-0.333$ & $-0.167$\\ \hline
$1.0$ & $0.500$ &{}&{}&{}&{}&{}&{}\\\hline
\end{tabular}
\end{center}
\label{tab:runge_kutta}
\end{table}
\end{example}
\markex
The above problem has been worked by the predictor-corrector\index{predictor-corrector}, Runge's, and the Runge-Kutta methods. This problem can be solved exactly by means of the Taylor series. The Taylor series solution is
\[
x(t) = 1-t + t^2 - t^3+\cdots
\]
The series converges to $1/(1+t)$ for $-1 < t < 1$. A check on the differential equation shows that
\[
x(t) = \frac{1}{1+t}
\]
is indeed a solution for all $t> -1$. Thus, $x(1)=0.500$.

The algorithms presented above require only one input point, $x(t_k)$, at each step to compute $x(t_{k+1})$. Such algorithms are referred to as single-step algorithms\index{single-step algorithms}. In the multi-step algorithms\index{multi-step algorithms}, a number of previous points are used, and the corresponding values of $f$, to evaluate $x(t_{k+1})$, that is, the multi-step algorithms reuse past information about the trajectory when evaluating its new point. In general, it may be difficult to decide which type of an algorithm is more efficient because the performance of a particular algorithm is problem dependent. Well known multi-step algorithms are: Adams-Bashforth\index{Adams-Bashforth!multi-step algorithms}, Adams-Moulton\index{Adams-Moulton!multi-step algorithms}, and Gear\index{Gear!multi-step algorithms}. More information about the above mentioned algorithms can be found in~\cite[Section~4.1]{parker89}.

%\section{Principles of Linearization}
%
%The replacement of a nonlinear system model by its linear approximation is called {\em linearization\/}. The motivation for linearization is that the dynamical behavior of many nonlinear system models can be well approximated within some range of variables by linear system models. Then, we can use well developed techniques for analysis and synthesis of linear systems to analyze a nonlinear system at hand. However, the results of analysis of nonlinear systems using their linearized models should be carefully interpreted in order to avoid unacceptably large errors due to the approximation in the process of linearization.
%
%We begin our discussion of linear approximations of nonlinear systems by considering a nonlinear element with a state variable $x$ and an output variable $y$ that are related by the equation
%\[
%y=h(x),
%\]
%where the function $h:\mb{\Bb R} \to \mb{\Bb R}$ is continuously differentiable, that is, $h \in {\cal C}^1$.
%Let $x_0$ be an operating point. We expand $h$ into a Taylor series about the operating point $x_0$ to obtain
%\begin{eqnarray*}
%y&=&h(x) \\
%&=&h(x_0)+\frac{dh(x_0)}{dx}(x-x_0)+\mbox{ higher order terms}.
%\end{eqnarray*}
%Linearization of $h(x)$ about the point $x_0$  consists of replacing $h$ by the linear approximation
%\begin{eqnarray}
%y&=&h(x_0)+\frac{dh(x_0)}{dx}(x-x_0)\nonumber\\
%&=&y_0+\frac{dh(x_0)}{dx}(x-x_0),
%\label{eq:aproxline}
%\end{eqnarray}
%where $y_0=h(x_0)$. Let $\Delta y=y-y_0$ and $\Delta x=x-x_0$. Then, we can represent~(\ref{eq:aproxline}) as
%\[
%\Delta y=\frac{dh(x_0)}{dx} \Delta x.
%\]
%Over a small range of $\Delta x$, the line~(\ref{eq:aproxline}) is a good approximation to the curve $y=h(x)$ in a neighborhood of the operating point $x_0$---see Figure~\ref{fig:lapprox} for an illustration of the above statement.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{lapprox.eps}
%%\vspace{3in}
%%\centerbmp{3.33in}{2.5in}{c:/mybook/lapprox.bmp}
%%\hskip1in\special{bmp:c:/mybook/lapprox.bmp x=4in y=3in}
%\caption{Linear approximation to a function $y=h(x)$.}
%\label{fig:lapprox}
%\end{figure}
%We now illustrate the process of linearization on the simple pendulum.
%\begin{example}\rm
%Consider the simple pendulum shown in Figure~\ref{fig:simpend}. Let $g$ be the gravity constant.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{simpend.eps}
%%\vspace{3.5in}
%%\centerbmp{3.33in}{2.5in}{c:/mybook/simpend.bmp}
%%\hskip1.2in\special{bmp:c:/mybook/simpend.bmp x=4.67in y=3.5in}
%\caption{Forces acting on the simple pendulum.}
%\label{fig:simpend}
%\end{figure}
%The tangential force component, $mg \sin \theta$, acting on the mass $m$ returns the pendulum to its equilibrium position. By Newton's second law
%\[
%\tau = -mgl \sin\theta,
%\]
%where $\tau = J\ddot{\theta}$ is the torque acting on the mass $m$ and $J=ml^2$ is the moment of inertia of the mass $m$. The equilibrium point for the pendulum is $\theta=0^{\circ}$. Note that $\tau\left(0^{\circ}\right)=0$. Therefore,
%\[
%\Delta \tau=\tau\;\mb{and}\; \Delta \theta=\theta.
%\]
%We also have
%\[
%\frac{d \tau}{d\theta}\left(0^{\circ}\right) = -mgl \cos 0^{\circ}=-mgl.
%\]
%Hence, the linearized model, about the equilibrium position $\theta=0^{\circ}$ has the form
%\[
%\tau=-mg l \theta.
%\]
%Thus, for small angular displacements $\theta$, the torque $\tau$ is proportional to the angular displacements. This approximation is quite accurate for $-\pi/4 \le \theta \le \pi/4$ as can be seen in Figure~\ref{fig:linapp}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.4]{linapp.eps}
%%\vspace{3.5in}
%%\hskip.6in\special{bmp:c:/mybook/linapp.bmp x=4.67in y=3.5in}
%\caption{Linearizing the equation modeling the simple pendulum.}
%\label{fig:linapp}
%\end{figure}
%\end{example}\markex
%If $h:\mb{\Bb R}^n \to \mb{\Bb R}$, that is, $y=h(x_1,x_2,\ldots,x_n)$ which means that the dependent variable depends upon several variables, the same principle applies. Let
%\[
%\b x_0=\left[\begin{array}{cccc}
%x_{10} & x_{20}&\cdots&x_{n0}
%\end{array}\right]^T
%\]
%be the operating point. The Taylor series expansion of $h$ about the operating point $x_0$ yields
%\[
%y-h(\b x_0)=\nabla h(\b x_0)^T(\b x-\b x_0)+\mb{higher order terms},
%\]
%where
%\[
%\nabla h(\b x_0)^T=\left[\begin{array}{cccc}
%\left.\frac{\partial h}{\partial x_1}\right\arrowvert_{\b x=\b x_0}&\left.\frac{\partial h}{\partial x_2}\right\arrowvert_{\b x=\b x_0}&\cdots&\left.\frac{\partial h}{\partial x_n}\right\arrowvert_{\b x=\b x_0}
%\end{array}\right].
%\]
%Geometrically, the linearization of $h$ about $\b x_0$ can be thought of as placing a tangent plane onto the nonlinear surface at the operating point $x_0$ as illustrated in Figure~\ref{fig:lintan}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{lintan.eps}
%%\vspace{3in}
%%\hskip1in\special{bmp:c:/mybook/lintan.bmp x=4in y=3in}
%\caption{Tangent plane as a linear approximation for a function of two variables.}
%\label{fig:lintan}
%\end{figure}
%
%\section{Linearizing Differential Equations}
%We now consider a dynamical system modeled by a set of nonlinear differential equations
%\[
%\left.\begin{array}{lll}
%\frac{dx_1}{dt} &=& f_1(x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
%\frac{dx_2}{dt} &=& f_2(x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
%{}&\vdots&{}\\
%\frac{dx_n}{dt} &=& f_n(x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)
%\end{array}\right\}
%\]
%We assume that the functions $f_i$, $i=1,2,\ldots,n$ are continuously differentiable. The above set of differential equations can be represented in vector form as
%\begin{equation}
%\dot{\b x}=\b f(\b x,\b u).
%\label{eq:xdoteqfxu}
%\end{equation}
%Let $\b u_e=\left[\begin{array}{cccc}
%u_{1e} & u_{2e}&\cdots&u_{me}
%\end{array}\right]^T$ be a constant input that forces the system~(\ref{eq:xdoteqfxu}) to settle into a constant, equilibrium, state $\b x_e=\left[\begin{array}{cccc}
%x_{1e} & x_{2e}&\cdots&x_{ne}
%\end{array}\right]^T$, that is, $\b u_e$ and $\b x_e$ satisfy
%\[
%\b f(\b x_e,\b u_e)=\b 0.
%\]
%We now perturb the equilibrium state by allowing
%\[
%\b x=\b x_e+\Delta \b x,\quad \b u=\b u_e+\Delta \b u.
%\]
%Taylor's expansion yields
%\begin{eqnarray*}
%\frac{d}{dt}\b x &=& \b f(\b x_e+\Delta \b x,\b u_e+\Delta \b u)\\
%&=& \b f(\b x_e,\b u_e)+\frac{\partial \b f}{\partial \b x}(\b x_e, \b u_e)\,\Delta \b x+\frac{\partial \b f}{\partial \b u}(\b x_e, \b u_e)\,\Delta \b u+\mb{higher order terms},
%\end{eqnarray*}
%where
%\[
%\frac{\partial \b f}{\partial \b x}(\b x_e,\b u_e)=\left.\left[\begin{array}{ccc}
%\frac{\partial f_1}{\partial x_1}&\cdots&\frac{\partial f_1}{\partial x_n}\\
%\vdots & {}&\vdots\\
%\frac{\partial f_n}{\partial x_1}&\cdots&\frac{\partial f_n}{\partial x_n}
%\end{array}\right]\right\arrowvert_{\b x=\b x_e \atop
%\b u=\b u_e}
%\!\!\mb{and}\quad\frac{\partial \b f}{\partial \b u}(\b x_e,\b u_e)=\left.\left[\begin{array}{ccc}
%\frac{\partial f_1}{\partial u_1}&\cdots&\frac{\partial f_1}{\partial u_m}\\
%\vdots & {}&\vdots\\
%\frac{\partial f_n}{\partial u_1}&\cdots&\frac{\partial f_n}{\partial u_m}
%\end{array}\right]\right\arrowvert_{\b x=\b x_e \atop
%\b u=\b u_e}
%\]
%are the Jacobian matrices of $\b f$ with respect $\b x$ and $\b u$ respectively, evaluated at the equilibrium point, $\left[\begin{array}{cc}\b x_e^T& \b u_e^T
%\end{array}\right]^T$. Note that
%\[
%\frac{d}{dt}\b x = \frac{d}{dt}\b x_e+\frac{d}{dt}\Delta \b x=\frac{d}{dt}\Delta \b x
%\]
%because $\b x_e$ is constant. Furthermore $\b f(\b x_e,\b u_e)=\b 0$. Let
%\[
%\b A=\frac{\partial \b f}{\partial \b x}(\b x_e,\b u_e)\quad\mb{and}\quad \b B= \frac{\partial \b f}{\partial \b u}(\b x_e,\b u_e).
%\]
%Neglecting higher order terms, we arrive at the linear approximation
%\[
%\frac{d}{dt}\Delta \b x=\b A\Delta \b x+\b B\Delta \b u.
%\]
%Similarly, if the outputs of the nonlinear system model are of the form
%\[
%\left.\begin{array}{lll}
%y_1 &=& h_1(x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
%y_2 &=& h_2(x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
%{}&\vdots&{}\\
%y_p &=& h_p(x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)
%\end{array}\right\}
%\]
%or in vector notation
%\[
%\b y=\b h(\b x,\b u),
%\]
%then Taylor's series expansion can again be used to yield the linear approximation of the above output equations. Indeed, let
%\[
%\b y=\b y_e+\Delta \b y,
%\]
%then we obtain
%\[
%\Delta \b y=\b C \Delta \b x+\b D \Delta \b u,
%\]
%where
%\[
%\b C=\frac{\partial \b h}{\partial \b x}(\b x_e,\b u_e)=\left.\left[\begin{array}{ccc}
%\frac{\partial h_1}{\partial x_1}&\cdots&\frac{\partial h_1}{\partial x_n}\\
%\vdots & {}&\vdots\\
%\frac{\partial h_p}{\partial x_1}&\cdots&\frac{\partial h_p}{\partial x_n}
%\end{array}\right]\right\arrowvert_{\b x=\b x_e \atop
%\b u=\b u_e}\!\in \mb{\Bb R}^{p\times n},
%\]
%and
%\[
%\b D=\frac{\partial \b h}{\partial \b u}(\b x_e,\b u_e)=\left.\left[\begin{array}{ccc}
%\frac{\partial h_1}{\partial u_1}&\cdots&\frac{\partial h_1}{\partial u_m}\\
%\vdots & {}&\vdots\\
%\frac{\partial h_p}{\partial u_1}&\cdots&\frac{\partial h_p}{\partial u_m}
%\end{array}\right]\right\arrowvert_{\b x=\b x_e \atop
%\b u=\b u_e}\!\in \mb{\Bb R}^{p\times m}
%\]
%are the Jacobian matrices of $\b h$ with respect $\b x$ and $\b u$ respectively, evaluated at the equilibrium point $\left[\begin{array}{cc}\b x_e^T& \b u_e^T
%\end{array}\right]^T$.
%\begin{example}\rm
%In this example, we linearize the model of Watt's governor that we derived in Section~\ref{subsec:watt}. Recall that the equations modeling the governor have the form:
%\[
%\left.\begin{array}{lll}
%\dot{x}_1 &=& x_2\\
%\dot{x}_2 &=& \frac{1}{2}N^2x_3^2\sin 2x_1-\frac{g}{L}\sin x_1 - \frac{b}{m L}x_2\\
%\dot{x}_3 &=& \frac{\kappa}{I}\cos x_1-\frac{\tau}{I}
%\end{array}\right\}
%\]
%We linearize the above model about an equilibrium state of the form
%\[
%\b x_e=\left[\begin{array}{ccc}
%x_{1e}&0&x_{3e}
%\end{array}\right]^T
%\]
%corresponding to the equilibrium control $u_e=\tau$.
%Equating the right-hand sides of the above engine-governor equations to zero yields
%\[
%\left.\begin{array}{lll}
%x_{2e} &=& 0\\
%N^2\,x_{3e}^2 &=& \frac{g}{L\cos x_{1e}}\\
%\cos x_{1e} &=& \frac{\tau}{\kappa}
%\end{array}\right\}
%\]
%The Jacobian matrix of the engine-governor model evaluated at the equilibrium point is
%\[
%\left.\frac{\partial \b f}{\partial \b x}\right\arrowvert_{\b x=\b x_e\atop
%u=u_e}=\left[\begin{array}{ccc}
%0 & 1 & 0\\
%N^2x_{3e}^2\cos 2x_{1e}-\frac{g}{L}\cos x_{1e} & -\frac{b}{m L} & N^2x_{3e}\sin 2x_{1e}\\
%-\frac{\kappa}{I}\sin x_{1e} & 0 & 0
%\end{array}\right].
%\]
%We use the following trigonometric identities,
%\[
%\cos 2\alpha=\cos ^2\alpha - \sin^2 \alpha,\quad \sin 2\alpha=2\sin \alpha \,\cos \alpha,
%\]
%and the previously obtained relation, $N^2\,x_{3e}^2 = \frac{g}{L\cos x_{1e}}$, to eliminate $N^2\,x_{3e}^2$ from the Jacobian matrix. Performing needed manipulations, we obtain
%\[
%\left.\frac{\partial \b f}{\partial \b x}\right\arrowvert_{\b x=\b x_e\atop
%u=u_e}=\left[\begin{array}{ccc}
%0 & 1 & 0\\
%-\frac{g \sin^2 x_{1e}}{L\cos x_{1e}} & -\frac{b}{m L} & \frac{2g\sin\,x_{1e}}{L x_{3e}}\\
%-\frac{\kappa}{I}\sin x_{1e} & 0 & 0
%\end{array}\right].
%\]
%Hence, the linearized model of the engine-governor system can be represented as
%\begin{equation}
%\frac{d}{dt}\left[\begin{array}{c}
%\Delta x_1\\
%\Delta x_2\\
%\Delta x_3
%\end{array}\right]=\left[\begin{array}{ccc}
%0 & 1 & 0\\
%-\frac{g \sin^2 x_{1e}}{L\cos x_{1e}} & -\frac{b}{m L} & \frac{2g\sin\,x_{1e}}{L x_{3e}}\\
%-\frac{\kappa}{I}\sin x_{1e} & 0 & 0
%\end{array}\right]\left[\begin{array}{c}
%\Delta {x_1}\\
%\Delta {x_2}\\
%\Delta {x_3}
%\end{array}\right] + \left[\begin{array}{c}
%0\\
%0\\
%-\frac{1}{I}
%\end{array}\right] \Delta u.
%\label{eq:linmoofwgo}
%\end{equation}
%\end{example}\markex
%
%\section{Describing Function Method}
%\label{sec:descfunme}
%The describing function method allows us to apply familiar frequency domain techniques, used in the linear system analysis, to the analysis of a class of nonlinear dynamical systems. The method can be used to predict limit cycles in nonlinear systems. The describing function method can be viewed as a ``harmonic linearization'' of a nonlinear element. The method provides a ``linear approximation'' to the nonlinear element based on the assumption that the input to the nonlinear element is a sinusoid of known, constant amplitude. The fundamental harmonic of the element's output is compared with the input sinusoid to determine the steady-state amplitude and phase relation. This relation is the describing function for the nonlinear element. The describing function method is based on the Fourier series. In our review of Fourier series, we use notion of a scalar product over a function space that we discuss next. Our presentation of the scalar product in the context of spaces of functions follows that of Lang~\cite{lang87}. The reader familiar with the Fourier series may go directly to Subsection~\ref{subsec:desfunana} where the describing function method analysis is presented.
%
%\subsection{Scalar Product of Functions}
%Let ${\cal C}([a,b])$ be the set of continuous functions on the interval $[a,b]\subset \mb{\Bb R}$.
%\begin{defn}
%The scalar product of two real-valued functions $u, v \in {\cal C}([a,b])$ is
%\[
%\langle u,v \rangle = \int_{a}^{b}u(t)v(t)dt.
%\]
%\end{defn}\markdefn
%
%Using simple properties of the integral, we can verify that the above scalar product satisfies the following conditions:
%\begin{enumerate}
%\item for all $u, v \in {\cal C}([a,b])$, we have $\langle u,v \rangle = \langle v,u \rangle$,
%\item if $u, v, w$ are elements of ${\cal C}([a,b])$, then $\langle u,v+w \rangle = \langle u,v \rangle + \langle u,w \rangle$,
%\item if $x$ is a number, then $\langle xu,v \rangle = x\langle u,v \rangle = \langle u,xv \rangle$,
%\item for all $v \in {\cal C}([a,b])$, we have $\langle v,v \rangle \ge 0$, and $\langle v,v \rangle > 0$ if $v \ne 0$.
%\end{enumerate}
%The functions $u, v \in {\cal C}([a,b])$ are said to be {\em orthogonal\/} if $\langle u,v \rangle = 0$. A set of functions from ${\cal C}([a,b])$ is said to be {\em mutually orthogonal\/} if each distinct pair of functions in the set is orthogonal.
%\begin{example}\rm
%The functions $\psi_n(t)=\sin\,n\omega t$, $\phi_n(t)=\cos\,n\omega t$, $n=1,2,\ldots$, and $\phi_0(t)=1$ form a mutually orthogonal set of functions on the interval $[t_0,t_0+{2\pi}/\omega]$. This can be verified by direct integration. For example, for positive $m$ and $n$ such that $m\ne n$, we have
%\begin{eqnarray}
%\langle \psi_m,\psi_n \rangle&=&\int_{t_0}^{t_0+{2\pi}/\omega}\sin\,m\omega t\,\sin\,n\omega t\, dt\nonumber\\
%&=&\frac{1}{2}\int_{t_0}^{t_0+{2\pi}/\omega}\left( \cos\,(m-n)\omega t-\cos\,(m+n)\omega
%t\right)dt\nonumber\\
%&=&\frac{1}{2\omega} \left.\left(\frac{ \sin\,(m-n)\omega t}{m-n}-\frac{\sin\,(m+n)\omega t}{m+n}\right)\right\arrowvert_{t_0}^{t_0+{2\pi}\omega}\nonumber\\
%&=&0
%\label{eq:psimpsin}
%\end{eqnarray}
%Note that $m+n\ne 0$ because $m$ and $n$ are positive. Using similar computations we can show that for $m\ne n$
%\begin{equation}
%\langle \phi_m,\phi_n \rangle = 0,
%\label{eq:himphin}
%\end{equation}
%and for all non-negative $m$ and positive $n$,
%\begin{equation}
%\langle \phi_m,\psi_n \rangle = 0.
%\label{eq:phimpsin}
%\end{equation}
%\end{example}\markex
%We define the {\em norm\/} of $v\in {\cal C}([a,b])$ to be
%\[
%\|v\|=\sqrt{\langle v,v \rangle}.
%\]
%\begin{example}\rm
%Consider $\psi_n(t)=\sin\,n\omega t$ on the interval $[t_0,t_0+{2\pi}/\omega]$. Then,
%\begin{eqnarray}
%\|\psi_n\|^2&=&\int_{t_0}^{t_0+{2\pi}/\omega}\sin\,n\omega t\,\sin\,n\omega t\, dt\nonumber\\
%&=& \int_{t_0}^{t_0+{2\pi}/\omega}\sin^2n\omega t\, dt\nonumber\\
%&=& \frac{1}{2}\int_{t_0}^{t_0+{2\pi}/\omega}\left(1-\cos\,2n\omega t\right)\, dt\nonumber\\
%&=& \frac{1}{2}\left.\left( t-\frac{\sin\,2n\omega t}{2n\omega}\right)\right\arrowvert_{t_0}^{t_0+{2\pi}/\omega}\nonumber\\
%&=&\frac{\pi}{\omega}.
%\label{eq:nsqofpsi}
%\end{eqnarray}
%Similarly,
%\begin{equation}
%\|\phi_n\|^2=\frac{\pi}{\omega},\quad n=1,2\ldots,
%\label{eq:nphisq}
%\end{equation}
%and
%\begin{equation}
%\|\phi_0\|^2=\frac{2\pi}{\omega}.
%\label{eq:nsforphio}
%\end{equation}
%\end{example}\markex
%
%So far we discussed only continuous functions. In many applications we have to work with more general functions, specifically with piecewise continuous functions.
%\begin{defn}
%A function $f$ is said to be {\em piecewise continuous\/} on an interval $[a, b]$ if the interval can be partitioned by a finite number of points, say $t_i$, so that
%\begin{enumerate}
%\item $a=t_0<t_1<\ldots<t_n=b$;
%\item the function $f$ is continuous on each open subinterval $(t_{i-1},t_i)$,
%\item the function $f$ approaches a finite limit as the end points of each subinterval are approached from within the subinterval, that is, the limits
%\[
%\lim_{{h \to 0}\atop {h>0}}f(t_i-h)\quad \mb{and}\quad \lim_{{h \to 0}\atop {h>0}}f(t_i+h)
%\]
%both exist.
%\end{enumerate}
%\end{defn}\markdefn
%
%The graph of a piecewise continuous function is shown in Figure~\ref{fig_piececon}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{piececon.eps}
%%\vspace{3in}
%%\centerbmp{4in}{3in}{c:/mybook/piececon.bmp}
%%\hskip1in\special{bmp:c:/mybook/piececon.bmp x=4in y=3in}
%\caption{A graph of a piecewise continuous function.}
%\label{fig_piececon}
%\end{figure}
%In further discussion, we assume that two piecewise functions $f$ and $g$ are equal if they have the same values at the points where they are continuous. Thus, the functions shown in Figure~\ref{fig_disfuns} are considered to be equal, though they differ at the points of discontinuity.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{disfuns.eps}
%%\vspace{3in}
%%\centerbmp{4in}{3in}{c:/mybook/disfuns.bmp}
%%\hskip1in\special{bmp:c:/mybook/disfuns.bmp x=4in y=3in}
%\caption{Two piecewise continuous functions that we consider to be equal.}
%\label{fig_disfuns}
%\end{figure}
%\begin{thm}
%Let $V$ be the space of functions that are piecewise continuous on the interval $[a,b]$. Let $f \in V$. Then, $\|f\|=0$ if and only if $f(x)=0$ for all but finite number of points $x$ in the interval $[a,b]$.
%\end{thm}
%{\bf Proof}\quad $(\Leftarrow)$ It is clear that if $f(x)=0$ except for a finite number of $x \in [a,b]$, then
%\[
%\|f\|^2=\int_a^b f^2(x)dx =0.
%\]
%
%
%$(\Rightarrow)$ Suppose $f$ is piecewise continuous on $[a,b]$, and let $a=t_0<t_1<\ldots<t_n=b$ be a partition of the interval $[a,b]$ such that $f$ is continuous on each subinterval $[t_{i-1},t_i]$ except possibly at the end points. Suppose that $\| f \|=0$. Then, also $\| f \|^2= \langle f,f \rangle =0$. The above means that
%\[
%\int_a^b f^2(x)dx=0,
%\]
%where the integral is the sum of the integrals over the intervals $[t_{i-1},t_i]$, that is,
%\[
%\int_a^b f^2(x)dx=\sum_{i=1}^n \int_{t_{i-1}}^{t_i}f^2(x)dx=0.
%\]
%Each integral satisfies
%\[
%\int_{t_{i-1}}^{t_i}f^2(x)dx \ge 0.
%\]
%Hence each such integral must be equal to $0$. The function $f$ is continuous on $(t_{i-1},t_i)$. Therefore,
%\[
%f^2(x)=0 \quad \mb{for} \quad t_{i-1} < x < t_i,
%\]
%which means that
%\[
%f(x)=0 \quad \mb{for} \quad t_{i-1} < x < t_i.
%\]
%Hence, $f(x)=0$ except at a finite number of points. The proof is complete.\\
%\hspace*{5.5in} $\Box$
%
%The above theorem implies that the scalar product on the space $V$ of piecewise continuous functions is positive definite. In other words, for a piecewise continuous function $f \in V$, we have $\langle f,f \rangle \ge 0$, and $\langle f,f \rangle =0$ if and only if $f = 0$ for all but a finite number of points in the interval $[a,b]$.
%
%
%
%
%\subsection{Fourier Series}
%We begin with the definition of a periodic function.
%\begin{defn}
%A function $f:\mb{\Bb R} \to \mb{\Bb R}$ is said to be periodic with the period $T$ if the domain of $f$ contains $t+T$ whenever it contains $t$, and if
%\[
%f(t+T)=f(t)
%\]
%for every value of $t$. The smallest positive value of $T$ for which $f(t+T)=f(t)$ is called the fundamental period of $f$.
%\end{defn}\markdefn
%
%It follows from the above definition that if $T$ is a period of $f$, then $2T$ is also a period, and so is any integral multiple of $T$.
%
%A periodic function, satisfying certain assumptions spelled out later, may be represented by the series
%\begin{equation}
%f(t)=A_0+\sum_{k=1}^{\infty}\left(A_k \cos k\omega t +B_k \sin k \omega t \right),
%\label{eq:fouriersero}
%\end{equation}
%where $\omega=2\pi/T$, with $T$ being the fundamental period of $f$. On the set of points where the series~(\ref{eq:fouriersero}) converges, it defines $f$, whose value at each point is the sum of the series for that value of $t$, and we say that the series~(\ref{eq:fouriersero}) is the Fourier series for $f$.
%
%Suppose that the series~(\ref{eq:fouriersero}) converges. We use the results of the previous subsection to find expressions for the coefficients $A_k$ and $B_k$. We first compute $A_k$ coefficients. Multiplying~(\ref{eq:fouriersero}) by $\phi_m=\cos m\omega t$, where $m>0$ is a fixed positive integer, and then integrating both sides with respect to $t$ from $t_0$ to $t_0+2\pi/\omega$ yields
%\begin{equation}
%\langle \phi_m, f \rangle = A_0\langle \phi_m,\phi_0 \rangle + \sum_{k=1}^{\infty} A_k\langle \phi_m,\phi_k\rangle + \sum_{k=1}^{\infty}B_k\langle \phi_m,\psi_k\rangle ,
%\label{eq:amceff}
%\end{equation}
%where $\phi_0=1$. It follows from~(\ref{eq:psimpsin}), (\ref{eq:himphin}), and~(\ref{eq:phimpsin}) that the only nonvanishing term on the right of~(\ref{eq:amceff}) is the one for which $k=m$ in the first summation. Using~(\ref{eq:nphisq}), we obtain
%\begin{eqnarray}
%A_m &=& \frac{\langle \phi_m, f \rangle}{\langle \phi_m,\phi_m\rangle}\nonumber\\
%&=& \frac{\omega}{\pi}\int_{t_0}^{t_0+2\pi/\omega}f(t)\cos m\omega t\, dt,\quad m=1,2,\ldots
%\label{eq:exforam}
%\end{eqnarray}
%To obtain $A_0$, we integrate~(\ref{eq:fouriersero}) and using~(\ref{eq:nsforphio}), we get
%\begin{equation}
%A_0=\frac{\langle \phi_0, f \rangle}{\langle \phi_0,\phi_0\rangle}
%= \frac{\omega}{2\pi}\int_{t_0}^{t_0+2\pi/\omega}f(t) dt.
%\label{eq:exforaz}
%\end{equation}
%An expression for $B_m$ may be obtained by multiplying~(\ref{eq:fouriersero}) by $\sin m\omega t$ and then integrating term by term from $t_0$ to $t_0+2\pi/\omega$. Using the orthogonality relations~(\ref{eq:psimpsin})--(\ref{eq:phimpsin}) yields
%\begin{equation}
%B_m=\frac{\langle \psi_m,f\rangle}{\langle \psi_m,\psi_m \rangle}=\frac{\omega}{\pi}\int_{t_0}^{t_0+2\pi/\omega}f(t)\sin m\omega t \,dt,\quad m=1,2,\ldots
%\label{eq:exforbm}
%\end{equation}
%Suppose now that a function $f$ is given that is periodic with period $2\pi/\omega$ and integrable on the interval $[t_0,t_0+2\pi/\omega]$. Then, coefficients $A_k$ and $B_k$ can be computed using~(\ref{eq:exforam})--(\ref{eq:exforbm}) and series of the form~(\ref{eq:fouriersero}) can be formally constructed. However, we do not know whether this series converges for each value of $t$, and if so, whether its sum is $f(t)$. The answer is given by the theorem that we state after explaining the notation used. We denote the limit of $f(t)$ as $t$ approaches $c$ from the right by $f(c+)$, that is,
%\[
%f(c+)=\lim_{t \to c+}f(t).
%\]
%Similarly, $f(c-)=\lim_{t \to c-}f(t)$ denotes the limit of $f(t)$ as $t \to c$ from the left. The mean value of the right and left hand limits at the point $c$ is
%\[
%\left(f(c+)+f(c-) \right)/2.
%\]
%Note that at any point $c$ where $f$ is continuous, we have
%\[
%f(c+)=f(c-)=f(c).
%\]
%\begin{thm}
%Assume that $f$ and its derivative are piecewise continuous on the interval $\left[t_0, t_0+2\pi/\omega\right]$. Furthermore, suppose that $f$ is defined outside the interval $[t_0, t_0+2\pi/\omega]$ so that it is periodic with period $2\pi/\omega$. Then, $f$ has a Fourier series of the form~(\ref{eq:fouriersero}) whose coefficients are given by~(\ref{eq:exforam})--(\ref{eq:exforbm}). The Fourier series converges to $\left(f(t+)+f(t-) \right)/2$ at all points of the interval $\left[t_0, t_0+2\pi/\omega\right]$.
%\label{thm:fourierb}
%\end{thm}
%
%With this knowledge of the Fourier series it is now possible to discuss the describing function method for analyzing a class of nonlinear control systems.
%
%\subsection{Describing Function in the Analysis of Nonlinear Systems}
%\label{subsec:desfunana}
%The describing function method can be considered as an extension of the Nyquist stability criterion to nonlinear systems. In our further discussion, we sometimes use a double box to represent a nonlinear element as illustrated in Figure~\ref{fig_non_el}.
%\begin{figure}
%\vspace{1in}
%\setlength{\unitlength}{.05in}
%\begin{center}
%\begin{picture}(40,20)
%\thicklines
%\put(10,5){\framebox(20,10){}}
%\put(11,6){\framebox(18,8){}}
%\put(4,10){\vector(1,0){6}}
%\put(30,10){\vector(1,0){6}}
%\end{picture}
%\end{center}
%\caption{A symbol for a nonlinear element.}
%\label{fig_non_el}
%\end{figure}
%
%Consider a feedback system that contains a nonlinear element as shown in Figure~\ref{fig:fbknonl}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{fbknonl.eps}
%%\vspace{3in}
%%\centerbmp{4in}{3in}{c:/mybook/fbknonl.bmp}
%%\hskip1.2in\special{bmp:c:/mybook/fbknonl.bmp x=4in y=3in}
%\caption{A nonlinear feedback system with a nonlinear element in the forward path.}
%\label{fig:fbknonl}
%\end{figure}
%Assume that the input to the nonlinear element is a sinusoidal signal $x(t)=X\sin \omega t$, where $X$ is the amplitude of the input sinusoid. The output of the nonlinear element is, in general, not sinusoidal in response to a sinusoidal input. Suppose that the nonlinear element output is periodic with the same period as its input and that it may be expanded in the Fourier series,
%\begin{eqnarray}
%y(t)&=&A_0+\sum_{k=1}^{\infty}\left(A_k \cos k\omega t +B_k \sin k \omega t \right)\nonumber\\
%&=& A_0+\sum_{k=1}^{\infty} Y_k \sin \left(k\omega t +\phi_k \right).
%\label{eq:nonoutfour}
%\end{eqnarray}
%Thus $y(t)$ consists of a fundamental component
%\[
%Y_1\sin \left(\omega t + \phi_1 \right),
%\]
%together with a mean level $A_0$ and harmonic components at frequencies $2\omega, 3\omega,\ldots$. The describing function method relies on the assumption that only the fundamental harmonic is significant. This assumption is called the {\em filtering hypothesis\/}. The filtering hypothesis can be expressed as
%\[
%\left|G(j\omega )\right| \gg \left|G(k j\omega )\right|\quad\mb{for}\quad k=2,3,\ldots,
%\]
%which requires that the linear element, modeled by the transfer function $G(s)$, be a low-pass filter. The filtering hypothesis is often valid because many, if not most, control systems are low-pass filters. This results in higher harmonics being more attenuated compared with the fundamental harmonic component. The nonlinear element is then approximately modeled by its describing function
%\begin{equation}
%N(X)=\frac{Y_1}{X}\angle \phi_1=\frac{\mb{fundamental component of output}}{\mb{input}}.
%\label{eq:descfunofnnel}
%\end{equation}
%The describing function is, in general, a complex-valued function of the amplitude $X$ of the input sinusoid.
%
%The stability of the nonlinear closed-loop system depicted in Figure~\ref{fig:fbknonl} can be analyzed using Nyquist's criterion applied to the associated linear system shown in Figure~\ref{fig:fbkdefu}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{fbkdefu.eps}
%%\vspace{3in}
%%\centerbmp{3.33in}{2.5in}{c:/mybook/fbkdefu.bmp}
%%\hskip1.2in\special{bmp:c:/mybook/fbkdefu.bmp x=4in y=3in}
%\caption{Associated linear system used in the describing function analysis of the nonlinear system shown in Figure~\ref{fig:fbknonl}.}
%\label{fig:fbkdefu}
%\end{figure}
%The associated linear closed-loop system characteristic equation is
%\begin{equation}
%1+N(X)G(s)=0.
%\label{eq:aslclce}
%\end{equation}
%In the $G(j\omega)$ complex plane, the locus traced out by $-1/N$ as $X$ is varied, is being considered as a generalization of the $(-1,j0)$ point in the linear case. The stability of the nonlinear system is determined by computing the encirclements of the $-1/N$ locus by the $G(j\omega)$ plot. We might have situations shown in Figure~\ref{fig:limcyde}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{limcyde.eps}
%%\vspace{4in}
%%\centerbmp{5.33in}{4in}{c:/mybook/limcyde.bmp}
%%\hskip.4in\special{bmp:c:/mybook/limcyde.bmp x=5.33in y=4in}
%\caption{Stability analysis using describing function.}
%\label{fig:limcyde}
%\end{figure}
%In case~(i), the $-1/N$ locus is not encircled by the $G(j\omega)$ plot. Assuming that the poles of $G(s)$ are in the left-half complex plane, the closed-loop system is asymptotically stable. For a limit cycle to exist, we must have an intersection between the polar plot of $G(j\omega)$ and the $-1/N$ locus. We have two intersections between $-1/N$ and $G(j\omega)$ in case~(ii). At intersection point~1, any increase in $X$ will cause the operating point to be encircled and consequently $X$ to grow. If there is a decrease in $X$, then the oscillations will decay. Thus, point~1 represents an unstable limit cycle. At intersection point~2, an increase in $X$ would shift the operating point outside the $G(j\omega)$ plot resulting in the decaying oscillations. Point~2 represents a stable limit cycle. In case~(iii), any oscillation will grow without limit.
%\begin{example}\label{ex:dfexamp}\rm
%Consider the nonlinear control system shown in Figure~\ref{fig:a_96_f1} that was analyzed in~\cite[pages~339--341]{power78}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{a_96_f1.eps}
%%\vspace{3.5in}
%%\centerbmp{4.67in}{3.5in}{c:/mybook/a_96_f1.bmp}
%%\hskip.5in\special{bmp:c:/mybook/a_96_f1.bmp x=4.67in y=3.5in}
%\caption{Nonlinear feedback control system analyzed in Example~\ref{ex:dfexamp}.}
%\label{fig:a_96_f1}
%\end{figure}
%We first derive the describing function of the ideal relay. When the input of the ideal relay is greater than zero, its output is $M$. When the relay's input is less than zero, its output is $-M$. Thus, the ideal relay can be described as
%\begin{equation}
%y=\left\{\begin{array}{rcc}
%M & \mb{for}& x>0\\
%-M & \mb{for}& x<0.
%\end{array}\right.
%\label{eq:desidrel}
%\end{equation}
%We assume that the input to the relay is a sinusoid, $x(t)=X\sin \omega t$, drawn vertically down from the relay characteristic in Figure~\ref{fig:relayde}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{relayde.eps}
%%\vspace{4.5in}
%%\centerbmp{4in}{3in}{c:/mybook/relayde.bmp}
%%\hskip.2in\special{bmp:c:/mybook/relayde.bmp x=6in y=4.5in}
%\caption{Input and output waveforms of an ideal relay.}
%\label{fig:relayde}
%\end{figure}
%The output $y(t)$ of the relay is a rectangular wave. It is drawn horizontally to the right of the relay characteristic in Figure~\ref{fig:relayde}. This rectangular wave is an odd function. For an odd function, we have $A_k=0$, $k=0,1,\ldots$, in its Fourier series. Thus, the Fourier series of the relay's output has the form
%\[
%y(t)=\sum_{k=1}^{\infty}B_k \sin k\omega t.
%\]
%The fundamental harmonic of $y(t)$ is
%\[
%B_1\sin\omega t=Y_1\sin \omega t,
%\]
%where
%\begin{eqnarray}
%B_1 &=& \frac{\omega}{\pi}\int_{0}^{2\pi/\omega}y(t)\sin \omega t\, dt\nonumber\\
%&=& \frac{2\omega}{\pi}\int_0^{\pi/\omega}M\sin \omega t\, dt\nonumber\\
%&=& \frac{2M\omega}{\pi}\int_0^{\pi/\omega}\sin \omega t \,dt\nonumber\\
%&=& \left.-\frac{2M}{\pi}\cos\left( \omega t\right)\right\arrowvert_0^{\pi/\omega}\nonumber\\
%&=&\frac{4M}{\pi}.
%\label{eq:funharex}
%\end{eqnarray}
%Hence, the describing function of the ideal relay is
%\[
%N(X)=\frac{Y_1}{X}\angle 0^{\circ}=\frac{4M}{\pi X}.
%\]
%A normalized plot of the describing function of the ideal relay nonlinearity is depicted in Figure~\ref{fig:reldepl}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{reldepl.eps}
%%\vspace{3in}
%%\centerbmp{4in}{3in}{c:/mybook/reldepl.bmp}
%%\hskip1.2in\special{bmp:c:/mybook/reldepl.bmp x=4in y=3in}
%\caption{A plot of the describing function for the ideal relay.}
%\label{fig:reldepl}
%\end{figure}
%In our example, $M=1$ and thus the describing function of the relay is
%\[
%N(X)=\frac{4}{\pi X}.
%\]
%The sinusoidal transfer function of the linear plant model is
%\[
%G(j\omega)=\frac{16}{j\omega (j \omega +2)^2}=\frac{-64\omega+j16(\omega ^2-4)}{\omega ( 4+\omega ^2)^2}.
%\]
%At $\omega = 2$ rad/sec the polar plot of $G(j \omega)$ intersects the plot of $-1/N$, see Figure~\ref{fig:nyqdesf}, to form a stable limit cycle.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{nyqdesf.eps}
%%\vspace{3in}
%%\hskip1.2in\special{bmp:c:/mybook/nyqdesf.bmp x=4in y=3in}
%\caption{Plots of $G(j\omega)$ and $-1/N$ of Example~\ref{ex:dfexamp}.}
%\label{fig:nyqdesf}
%\end{figure}
%The amplitude of the limit cycle is obtained from the relation, $G(j2)=-1/N$, that is,
%\[
%\frac{-64}{(4+4)^2}=-\frac{\pi X}{4}.
%\]
%Hence, the amplitude of the limit cycle is
%\[
%X=\frac{4}{\pi},
%\]
%while the amplitude of the fundamental harmonic of the plant output, $c$, is
%\[
%\frac{4}{\pi}\left|G(j2) \right|=\frac{4}{\pi}.
%\]
%Note that the filtering hypothesis is satisfied in this example. Indeed, the Fourier series of the output of the relay is
%\[
%\frac{4}{\pi}\left(\sin \omega t + \frac{1}{3}\sin 3\omega t +\frac{1}{5}\sin 5 \omega t\right).
%\]
%Thus, for example, the amplitude of the third harmonic of the plant output is
%\[
%\frac{4}{3\pi}\left|G(j6) \right|=\frac{1}{45}\frac{4}{\pi}.
%\]
%The higher harmonics are even more attenuated compared with the fundamental harmonic.
%\end{example}\markex


\section{Notes}
The system of differential equations~(\ref{eq:nosywtlicyex}) in Exercise~\ref{ex:limicycfbdp} comes from Boyce and DiPrima~\cite[page~523]{boyce97}. Al-Khafaji and Tooley~\cite{alkhafaji86} provide an introductory treatment of numerical solutions of differential equations. For further discussion of solving differential equations using numerical methods\index{numerical methods} and analyses of their computational efficiency\index{computational efficiency}, the reader is referred to Salvadori and Baron~\cite{salvadori61}, Conte and de~Boor~\cite[Chapter 8]{conte80}, or Parker and Chua~\cite{parker89}.
%For an in-depth treatment of the describing function method, we refer to Graham and McRuer~\cite{graham61} and Chapters~6 and~7 of Hsu and Meyer~\cite{hsu68}. Graham and McRuer~\cite{graham61} make an interesting remark about the origins of the describing function method. They note on page~92 of their 1961 book: ``While the notion of representing a nonlinearity by an ``equivalent'' linear element is quite old and has been used by many writers, the first instance of a major systematic exploitation of the technique was probably by N. M. Kryloff and N. N. Bogoliuboff, {\em Introduction to Non-linear Mechanics\/} (a free translation of the Russian edition of 1937, by S. Lefschetz), Princeton University Press, Princeton, N. J., 1943. See also N. Minorsky, {\em Introduction to Non-linear Mechanics\/}, J. W. Edwards, Ann Arbor, Mich., 1947.'' The so-called sinusoidal describing function that we analyzed in the previous section, was developed almost simultaneously in several different countries during and just after World War II. The first to introduce the method in the open literature appears to be A.~Tustin of England who published his paper in the Journal of the IEE in 1947.
%
%
%The behavior of numerous dynamical systems is often modeled by a set of linear, first-order, differential equations. Frequently, however, a linear model is a result of linearization of a nonlinear model. Linear models seem to dominate the controls' literature. Yet nature is nonlinear. ``All physical systems are nonlinear and have time-varying parameters in some degree.''~\cite[page~1]{graham61}. How can one reconcile this paradox?  Graham and McRuer~\cite[page~1]{graham61} have the following answer: ``Where the effect of the nonlinearity is very small, or if the parameters vary only slowly with time, linear constant-parameter methods of analysis can be applied to give an approximate answer which is adequate for engineering purposes. The analysis and synthesis of physical systems, predicated on the study of linear constant-parameter mathematical models, has been, in fact, an outstandingly successful enterprise. A system represented as linear, however, is to some extent a mathematical abstraction that can never be encountered in a real world. Either by design or because of nature's ways it is often true that experimental facts do not, or would not, correspond with any prediction of linear constant-parameter theory. In this case nonlinear or time-varying-parameter theory is essential to the description and understanding of physical phenomena.''

\section{Exercises}
\begin{exercise}\label{ex:lin_algebra_test1_09}\rm
Suppose that $\b x=\left[\begin{array}{ccc}
1 & 2 & 3
\end{array}\right]^{\top}$ is the representation of the given vector in the natural basis $\{\b e_1, \b e_2, \b e_3\}$. Find the representation of this vector in the basis $\{\b e_1', \b e_2', \b e_3'\}$, where
\begin{eqnarray*}
\b e_1' &=& \b e_1 + \b e_3\\
\b e_2' &=& \b e_2-\b e_3\\
\b e_3' &=& \b e_3.
\end{eqnarray*}
\end{exercise}

\begin{exercise}\label{ex:dist_eigenval_pf}\rm Prove Theorem~\ref{thm:dist_eigenval}.
\end{exercise}

\begin{exercise}\label{ex:trans_mat to_get_J}\rm Find a similarity transformation\index{similarity transformation} that transforms
\[
\b A=\left[\begin{array}{cccc}
\lambda & 0 & 0 & 0\\
0 & \lambda & 0 & 1\\
1 & 0 & \lambda & 0\\
0 & 0 & 1 & \lambda
\end{array}\right]
\]
into its Jordan form. Do not fix $\lambda$; construct a transformation that works for any $\lambda$. Verify your calculations using MATLAB's Symbolic Toolbox\index{MATLAB!Symbolic Toolbox}.
\end{exercise}

\begin{exercise}\label{ex:complex_eig_in_A}\rm Calculate the solution for the system,
\[
\dot{\b x}=\left[\begin{array}{cc}
0 & -1\\
1 & -1
\end{array}\right]\b x + \left[\begin{array}{c}
1\\
0
\end{array}\right]u,\quad \b x(0)=\left[\begin{array}{c}
x_1(0)\\
x_2(0)
\end{array}\right],
\]
where $u(t)$ is the unit step, that is, $u(t)=1(t)$. Plot the system trajectories on the time interval $[0,\; 2]$ for $\b x(0)=\b 0$.
\end{exercise}


\begin{exercise}\label{ex:discretize_with_ZOH}\rm
Construct an exact linear discrete state-space model\index{exact linear discrete model} for the system,
\begin{eqnarray*}
\dot{\b x} &=& \left[\begin{array}{cc}
0 & -1\\
1 & -1
\end{array}\right]\b x + \left[\begin{array}{c}
1\\
0
\end{array}\right]u\\
y &=& \left[\begin{array}{cc}
0 & 1
\end{array}\right]\b x
\end{eqnarray*}
The sampling period\index{sampling period} is $h=0.2$~sec. Verify you calculations using MATLAB's function \verb!c2d!\index{c2d!MATLAB's function}. Find the transfer function of the discrete systems and its step response. Plot the step response of the discrete system on the time interval $[0,\; 2~\mb{sec}]$.
\end{exercise}
\begin{exercise}\rm
Sketch a phase-plane portrait\index{phase-plane portrait} of the differential equation
\[
\ddot{x}-\dot{x}+2x+2\mb{sign}\left(2x+\dot{x} \right)=0.
\]
\end{exercise}
\begin{exercise}\rm
For the system of differential equations
\[
\left.\begin{array}{lll}
\dot{x}_1 &=& -\frac{7}{4}x_1+\frac{1}{4}x_2\\
\dot{x}_2 &=& \frac{3}{4}x_1-\frac{5}{4}x_2
\end{array}\right\}
\]
find the phase-plane asymptotes\index{phase-plane asymptotes}.
\end{exercise}
\begin{exercise}\rm
For a dynamical system modeled by
\[
\left.\begin{array}{lll}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& -x_1
\end{array}\right\}
\]
find the time $T$ it takes the system trajectory to move from the state $\b x(0)=\left[\begin{array}{cc}
0 & 1
\end{array}\right]^T$ to the state $\b x(T)=\left[\begin{array}{cc}
1 & 0
\end{array}\right]^T$.
\end{exercise}
\begin{exercise}\label{exer:straitraj}\rm
Consider a trajectory in the state plane shown in Figure~\ref{fig:traject}. Find the time it takes the representative point to move from A to B.
\begin{figure}
\centering
\includegraphics[scale=0.5]{traject.eps}
%\vspace{2.5in}
%\hskip1.5in\special{bmp:c:/mybook/traject.bmp x=3.33in y=2.5in}
\label{fig:traject}
\caption{A trajectory for Exercise~\ref{exer:straitraj}.}
\end{figure}
\end{exercise}
\begin{exercise}\rm
Given the following initial value problem
\[
\dot{x}=f(t,x),\quad x(t_0)=x_0.
\]
\begin{itemize}
\item[(i)] Derive the Euler formula\index{Euler formula} for solving the above initial value problem using a Taylor series.
\item[(ii)] Make use of a Taylor series\index{Taylor series} with a remainder to compute the local error formula assuming that the data at the $n$-th step are correct.
\end{itemize}
\end{exercise}
\begin{exercise}\rm
Use Bendixson's theorem\index{Bendixson's theorem} to investigate the existence of limit cycles\index{limit cycle!existence} for the system
\begin{equation}
\left.\begin{array}{ccl}
\dot{x}_1&=& x_1+x_2-x_1\left(x_1^2+x_2^2 \right)\\
\dot{x}_2&=& -x_1+x_2-x_2\left(x_1^2+x_2^2 \right)
\end{array}\right\}
\label{eq:nosywtlicyex}
\end{equation}
Sketch a phase-plane portrait for this system.
\label{ex:limicycfbdp}
\end{exercise}
\begin{exercise}\rm
Use Bendixson's theorem to investigate the existence of limit cycles for the system
\begin{equation}
\left.\begin{array}{ccl}
\dot{x}_1&=& x_2-x_1\left(x_1^2+x_2^2 -c\right)\\
\dot{x}_2&=& -x_1-x_2\left(x_1^2+x_2^2 -c\right)
\end{array}\right\}
\label{eq:licyex_parker}
\end{equation}
where $c$ is a constant Sketch phase-plane portraits for this system for $c=-0.2$ and $c=0.2$. The above nonlinear system was analyzed in~\cite[page~204]{parker89}.
\label{ex:parker_limicyc}
\end{exercise}

\begin{exercise}\rm
Verify that the Runge method\index{Runge method} agrees with the Taylor series solution through terms of degree three\index{Taylor series!of degree three}.
\end{exercise}
\begin{exercise}\rm
Verify that the Runge-Kutta method\index{Runge-Kutta method} agrees with the Taylor series solution through terms of degree four\index{Taylor series!of degree four}.
\end{exercise}
%\begin{exercise}\label{ex:nonlindamp}\rm
%Consider the nonlinear mass-spring-damper system depicted in Figure~\ref{fig:j99qe}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{qe99jf1.eps}
%%\vspace{3in}
%%\hskip1in\special{bmp:c:/mybook/qe99jf1.bmp x=4in y=3in}
%\caption{The mass-spring-damper system of Exercise~\ref{ex:nonlindamp}.}
%\label{fig:j99qe}
%\end{figure}
%The spring coefficient $K=0.5 x^2$~N/m, the damping coefficient $B=0.1 \dot{x}^2$~N$\cdot$sec/m, and $M=1$~kg.
%\begin{itemize}
%\item[(i)] Write the second-order differential equation describing the motion of the system.
%\item[(ii)] Choose $x_1=x$ and $x_2=\dot{x}$. Represent the equation of motion derived in~(i) in state-space format.
%\item[(iii)] For the model from~(ii), find the equilibrium state corresponding to $u=0$. Then, linearize the model about the obtained equilibrium point.
%\end{itemize}
%
%\end{exercise}
%\begin{exercise}\label{ex:cyclotro}\rm {(based on~\cite[pages~29--31]{toptscheev77})}  A cyclotron, whose operation is illustrated in Figure~\ref{fig:cyclotron}, accelerates charged particles into high energies so that they can be used in atom-smashing experiments.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{cyclotrn.eps}
%%\vspace{3.5in}
%%\hskip.7in \special{bmp:c:/mybook/cyclotrn.bmp x=4.67in y=3.5in}
%\caption{Schematic of a cyclotron of Exercise~\ref{ex:cyclotro}.}
%\label{fig:cyclotron}
%\end{figure}
%The operation of the cyclotron can be described as follows. Positive ions enter the cyclotron from the ion source, S, and are available to be accelerated. An electric oscillator establishes an accelerating potential difference across the gap of the two D-shaped objects, called dees. An emerging ion from the ion source, S, finds the dee that is negative. It will accelerate toward this dee and will enter it. Once inside the dee, the ion is screened from electric fields by the metal walls of the dees. The dees are immersed in a magnetic field so that the ion's trajectory bends in a circular path. The radius of the path depends on the ion's velocity. After some time, the ion emerges from the dee on the other side of the ion source. In the meantime, the accelerating potential changes its sign. Thus the ion again faces a negative dee. The ion further accelerates and describes a semicircle of larger radius. This process goes on until the ion reaches the outer edge of one dee, where it is pulled out of the cyclotron by a negatively charged deflector plate. The moving ion is acted upon by a deflecting force of magnitude
%\[
%q(v+\Delta v)(B+\Delta B),
%\]
%where $q$ is the charge of the ion, $v$ is the ion's speed, $\Delta v$ is the increase in the ion's speed, $B$ is the magnetic field acting on the ion. The nominal radius of the ion's orbit is $r$. The change in the magnetic field as the result of a change in the orbit of the ion is $\Delta B$. As the ion speed increases, the relativistic mass $m$ increases by $\Delta m$. The ion moving in a circular path undergoes centrifugal force of magnitude
%\[
%\frac{(m+\Delta m)(v+\Delta v)^2}{r+\Delta r}.
%\]
%By Newton's second law, the equation describing the motion of the ion is
%\[
%(m+\Delta m)\frac{d^2}{dt^2}(r+\Delta r)=\frac{(m+\Delta m)(v+\Delta v)^2}{r+\Delta r}-q(v+\Delta v)(B+\Delta B).
%\]
%Linearize the above equation. The motion of the ion along a nominal path is described by
%\[
%m\frac{d^2 r}{dt^2}=\frac{mv^2}{r}-qvB.
%\]
%Find the transfer function relating the change in radius of motion, $\Delta r$, of the ion, and the change in the momentum of the ion
%\[
%\Delta p=m\Delta v+\Delta m v.
%\]
%In other words, find the transfer function ${\cal L}( \Delta r)/{\cal L}(\Delta p)$, where ${\cal L}$ is the Laplace transform operator.
%\end{exercise}
%\begin{exercise}\rm
%Linearize the permanent-magnet stepper motor model, analyzed in Subsection~\ref{subsec:pmstepmo}, about $x_{1eq}=\theta_{eq}=\mb{constant}$.
%\end{exercise}
%\begin{exercise}\rm
%Let $v_1,v_2,\ldots,v_n$ be non-zero, mutually orthogonal elements of the space $V$ and let $v$ be an arbitrary element of $V$.
%\begin{itemize}
%\item[(i)] Find the component of $v$ along $v_j$.
%\item[(ii)] Let $c_j$ denote the component of $v$ along $v_j$. Show that the element
%\[
%v-c_1v_1-\cdots c_nv_n
%\]
%is orthogonal to $v_1,v_2,\ldots, v_n$.
%\item[(iii)] Let $c_j$ be the component of $v\in V$ along $v_j$. Show that
%\[
%\sum_{i=1}^n c_i^2\le \|v\|^2.
%\]
%The above relation is called the Bessel inequality.
%\end{itemize}
%\end{exercise}
%\begin{exercise}\label{ex:limplusdeadzo}\rm
%Consider characteristics of two nonlinear elements that are shown in Figure~\ref{fig:non_els}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{non_els.eps}
%%\vspace{3in}
%%\hskip1.25in\special{bmp:c:/mybook/non_els.bmp x=4in y=3in}
%\caption{Characteristics of a limiter and dead zone nonlinear elements of Exercise~\ref{ex:limplusdeadzo}.}
%\label{fig:non_els}
%\end{figure}
%Let
%\begin{equation}
%N_S(x)=\frac{2}{\pi}\left[\arcsin \frac{1}{x}+\frac{1}{x}\cos\left(\arcsin \frac{1}{x}\right)\right].
%\label{eq:nsexpress}
%\end{equation}
%Then, the describing function of the limiter is
%\[
%N_{l}=\left\{ \begin{array}{lc}
%K, & M \le S\\
%KN_S \left(\frac{M}{S}\right), & M>S.
%\end{array} \right.
%\]
%The describing function of the dead zone is
%\[
%N_{d-z}=\left\{ \begin{array}{lc}
%0, & M \le A\\
%K\left[1-N_S\left(\frac{M}{A}\right)\right], & M>A.
%\end{array}\right.
%\]
%Derive the describing function of the nonlinear element whose characteristic is shown in Figure~\ref{fig:non_sat}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.55]{non_sat.eps}
%%\vspace{3in}
%%\hskip1.25in\special{bmp:c:/mybook/non_sat.bmp x=4in y=3in}
%\caption{Characteristic of a nonlinear element of Exercise~\ref{ex:limplusdeadzo}.}
%\label{fig:non_sat}
%\end{figure}
%\end{exercise}
%\begin{exercise}\label{ex:nonelenone}\rm
%The describing function of the nonlinear element shown in Figure~\ref{fig:fin_non}, is
%\[
%N_1=\left\{ \begin{array}{lc}
%0, & M \le A\\
%K\left[1-N_S\left(\frac{M}{A}\right)\right]+\frac{4B}{\pi M}\sqrt{1-\left( \frac{A}{M}\right)^2}, & M>A ,
%\end{array}\right.
%\]
%where $N_S(x)$ is given by~(\ref{eq:nsexpress}). The describing function of the dead-zone nonlinearity, is given in Exercise~\ref{ex:limplusdeadzo} and its characteristic is shown in Figure~\ref{fig:non_els}. Derive the describing function of the ideal relay with dead zone whose characteristic is shown in Figure~\ref{fig:fin_non3}.
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{fin_non.eps}
%%\vspace{3in}
%%\centerbmp{3.33in}{2.5in}{c:/mybook/fin_non.bmp}
%%\hskip1.25in\special{bmp:c:/ee675/fin_non.bmp x=4in y=3in}
%\caption{Nonlinear element with the describing function $N_1$ of Exercise~\ref{ex:nonelenone}.}
%\label{fig:fin_non}
%\end{figure}
%%\begin{figure}
%%\vspace{3in}
%%\hskip1.25in\special{bmp:c:/mybook/fin_non2.bmp x=4in y=3in}
%%\caption{Dead-zone nonlinearity.}
%%\label{fig:fin_non2}
%%\end{figure}
%\begin{figure}
%\centering
%\includegraphics[scale=0.5]{fin_non3.eps}
%%\vspace{3in}
%%\centerbmp{3.33in}{2.5in}{c:/mybook/fin_non3.bmp}
%%\hskip1.25in\special{bmp:c:/ee675/fin_non3.bmp x=4in y=3in}
%\caption{Characteristic of the ideal relay with dead zone.}
%\label{fig:fin_non3}
%\end{figure}
%\end{exercise} 