\chapter{Linearization}

\begin{flushright}
\parbox{5.5in}{Early amplifiers were made of crude materials that tended to disintegrate over use, causing the amp to ``run away.'' Not only would an aging relay amplify the phone signal, it would mistakenly compound any tiny deviation from the range it expected until the mushrooming error filled and killed the system. What was needed was Heron's {\it regula}, a counter signal to rein in the chief signal, to dampen the effect of the perpetual recycling. Black came up with a {\it negative} feedback loop, which was designated negative in contrast to the snowballing positive loop of the amplifier. Conceptually, the electrical negative feedback loop is a toilet flusher or thermostat. This braking circuit keeps the amplifier honed in on a steady amplification in the same way a thermostat hones in on a steady temperature. But instead of metallic levers, a weak train of electrons talks to itself.~\cite[page~116]{kelly94}}
\end{flushright}




The replacement of a nonlinear system\index{nonlinear system!linearization} model by its linear approximation\index{linear approximation} is called {\em linearization\/}\index{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\index{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 present two different approaches  to linear approximation of nonlinear systems.

\section{Taylor 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\index{differentiable!continuously}, that is, $h \in {\cal C}^1$.
Let $x_0$ be a point of interest that we refer to as an \emph{operating point}\index{operating point}\index{operating point|see{equilibrium}}. 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\index{simple pendulum|see{pendulum}} nonlinear model\index{nonlinear model!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\index{force!tangential component} component, $mg \sin \theta$, acting on the mass $m$ returns the pendulum to its rest, or equilibrium\index{equilibrium} position. By Newton's second law\index{Newton's second law}
\begin{equation}
\tau = -mgl \sin\theta,
\label{eq:mod_simp_pendul}
\end{equation}
where $\tau = I\ddot{\theta}$ is the torque acting on the mass $m$ and $I=ml^2$ is the moment of inertia of the mass $m$.

We are interested in linearizing the above pendulum model about its rest position,\index{rest position}\index{rest position|see{equilibrium}} or equilibrium\index{equilibrium}; call it $\theta_e$. At such a point, $\tau(\theta_e)=0$.


There are many equilibrium points of the pendulum model~\eqref{eq:mod_simp_pendul}. We select $\theta_e=0^{\circ}$. Note that $\tau_e=\tau\left(0^{\circ}\right)=0$. Therefore,
\[
\Delta \tau=\tau-\tau_e=\tau\;\mb{and}\; \Delta \theta=\theta-\theta_e=\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\index{several variables}, the same principle applies. Let
\[
\b x_0=\left[\begin{array}{cccc}
x_{10} & x_{20}&\cdots&x_{n0}
\end{array}\right]^{\top}
\]
be the operating point. The Taylor series\index{Taylor series!several variables} expansion of $h$ about the operating point\index{operating point!several variables} $\b x_0$ yields
\[
y-h(\b x_0)=\nabla h(\b x_0)^{\top}(\b x-\b x_0)+\mb{higher order terms},
\]
where
\[
\nabla h(\b x_0)^{\top}=\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\index{tangent plane} onto the nonlinear surface at the operating point $\b 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}
\subsection{Time-Invariant Unforced System}
We consider a dynamical system modeled by a set of first-order differential equations,
\[
\left.\begin{array}{lll}
\frac{dx_1}{dt} &=& f_1(x_1,x_2,\ldots,x_n)\\
\frac{dx_2}{dt} &=& f_2(x_1,x_2,\ldots,x_n)\\
{}&\vdots&{}\\
\frac{dx_n}{dt} &=& f_n(x_1,x_2,\ldots,x_n)
\end{array}\right\}
\]
We assume that the functions $f_i$, $i=1,2,\ldots,n$ are continuously differentiable. We represent the above system of differential equations in vector form as
\begin{equation}
\dot{\b x}=\b f(\b x).
\label{eq:time_inv_no_u}
\end{equation}
Note that the right-hand side of the above modeling equations does not depend explicitly on time. Thus~\eqref{eq:time_inv_no_u} is a time-invariant system model\index{system model!time-invariant}. We are interested in constructing a linear model approximation of the above nonlinear model\index{nonlinear model!linear approximation} in a vicinity\index{vicinity} of a point of interest, which is an equilibrium point\index{equilibrium point} of the model.
\begin{defn}
A point $\b x_e$ is an equilibrium\index{equilibrium state} state, or an equilibrium point\index{equilibrium point|see{equilibrium state}}, of the system modeled by~\eqref{eq:time_inv_no_u} if it has the property that when the system starts at $\b x_e$, it will stay at $\b x_e$ for all future time. The equilibrium points are obtained by solving the algebraic equation\index{algebraic equation!solving for equilibrium points}
\[
\b f(\b x)=\b 0.
\]
\end{defn}\markdefn

A nonlinear system may have a number of equilibrium states.
\begin{example}\rm
The nonlinear system model
\[
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2
\end{array}\right]=\left[\begin{array}{c}
{x}_2\\
x_1 - {x}_2 - x_1^2
\end{array}\right],
\]
has two equilibrium states
\[
\b x_e^{(1)}=\left[\begin{array}{c}
0\\
0
\end{array}\right]\quad\mb{and}\quad \b x_e^{(2)}=\left[\begin{array}{c}
1\\
0
\end{array}\right].
\]
\end{example}\markex

We are interested here in isolated equilibrium states\index{isolated equilibrium}\index{equilibrium!isolated}; these are the equilibrium states that have no other equilibrium states in some neighborhood. In the above example, the equilibrium states are isolated. We can formally define an isolated equilibrium state using notion of the $r$-neighborhood\index{r-neighborhood!isolated equilibrium}.
\begin{defn}
The $r$-neighborhood is a set of points
\[
\{\b x: \|\b x-\b x_e\| < r \},
\]
where $\|\cdot\|$ can be any $p$-norm on $\mb{\Bb R}^n$---see Section~\ref{sec_app:vectors} in the Appendix.
\end{defn}\markdefn

Having the above definition, we formally define an isolated equilibrium state.

\begin{defn}
An equilibrium state $\b x_e$ is isolated if there is an $r>0$ such that the $r$-neighborhood of $\b x_e$ contains no equilibrium states other than $\b x_e$.
\end{defn}\markdefn

We note that a linear system, $\dot{\b x}=\b A\b x$, can have only one isolated equilibrium state. Note also that if there two different points, $\b x_e^{(1)}$ and $\b x_e^{(2)}$ in $\mb{\Bb R}^n$, that are equilibrium states of the linear system, then so is any point on the line connecting these two points.

To proceed, we expand the right-hand side of~\eqref{eq:time_inv_no_u} into a Taylor series\index{Taylor series!linearization} about an isolated equilibrium state $\b x_e$ to obtain
\begin{equation}
\frac{d}{dt}\b x=\b f(\b x_e)+\frac{\partial \b f}{\partial \b x}(\b x_e)\,\Delta \b x + \mb{higher-order terms},
\label{eq:t_exapansion_lti}
\end{equation}
where $\frac{\partial \b f}{\partial \b x}(\b x_e)$ is the Jacobian matrix\index{Jacobian matrix!linearization} of $\b f$ with respect $\b x$, evaluated at the equilibrium state $\b x_e$ and $\Delta \b x = \b x-\b x_e$. 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. Taking the above into account and the fact that $\b f(\b x_e)=\b 0$, we obtain a linear approximation of~\eqref{eq:t_exapansion_lti},
\[
\frac{d}{dt}\Delta \b x = \frac{\partial \b f}{\partial \b x}(\b x_e)\,\Delta \b x
\]
Let
\[
\b A= \frac{\partial \b f}{\partial \b x}(\b x_e).
\]
Then, a linear model of~\eqref{eq:time_inv_no_u} is written as
\begin{equation}
\frac{d}{dt}\Delta \b x= \b A\, \Delta \b x .
\label{eq:t_lti_lin_mo}
\end{equation}
\begin{example}\rm For the nonlinear dynamical system model,
\[
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2
\end{array}\right]=\left[\begin{array}{l}
2-x_1x_2\\
-4 + x_1^2
\end{array}\right],
\]
\begin{itemize}
\item find the equilibrium states;
\item construct the corresponding linear models about the equilibrium\index{equilibrium!linearization about} states found above.
\end{itemize}
\markex

\begin{itemize}
\item To find equilibrium states, we solve the algebraic equations, $\dot{x}_1=\dot{x}_2=0$, that is,
\[
2-x_1x_2=0\quad\mb{and}\quad -4 + x_1^2=0.
\]
We have two equilibriums states,
\[
\b x_e^{(1)}=\left[\begin{array}{c}
2\\
1
\end{array}\right]\quad\mb{and}\quad \b x_e^{(2)}=\left[\begin{array}{c}
-2\\
-1
\end{array}\right].
\]
\item A linear model has the form,
\[
\frac{d}{dt}\Delta \b x = \b A\: \Delta \b x,
\]
where
\[
\b A=\left[\begin{array}{cc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2}\\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2}
\end{array}\right]=\left[\begin{array}{cc}
-x_2 & -x_1\\
2x_1 & 0
\end{array}\right]
\]
evaluated at the equilibrium state about which we linearize the nonlinear system. We have two linear models:
\begin{enumerate}
\item about $\b x_e^{(1)}$,
\[
\frac{d}{dt}\Delta\b x = \left[\begin{array}{cc}
-1 & -2\\
4 & 0
\end{array}\right]\Delta \b x,
\]
\item about $\b x_e^{(2)}$,
\[
\frac{d}{dt}\Delta\b x = \left[\begin{array}{cc}
1 & 2\\
-4 & 0
\end{array}\right]\Delta \b x.
\]
\end{enumerate}
\end{itemize}


\end{example}\markex

\subsection{Time-Varying Forced System}

We now consider a dynamical system modeled by a set of nonlinear differential equations\index{differential equations!nonlinear}
\[
\left.\begin{array}{lll}
\frac{dx_1}{dt} &=& f_1(t, x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
\frac{dx_2}{dt} &=& f_2(t, x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
{}&\vdots&{}\\
\frac{dx_n}{dt} &=& f_n(t, x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)
\end{array}\right\}
\]
The above set of differential equations can be represented in vector form as
\begin{equation}
\dot{\b x}=\b f(t, \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]^{\top}$ 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]^{\top}$ at time $t_0$, that is, $\b u_e$ and $\b x_e$ satisfy
\[
\b f(t, \b x_e,\b u_e)=\b 0\quad\mb{for all}\quad t\ge t_0 .
\]
We say that the pair $(\b x_e, \b u_e)$ is the equilibrium pair\index{equilibrium pair} at $t_0$. We note that if the pair $(\b x_e, \b u_e)$ is an equilibrium pair at $t_0$, then it is also an equilibrium pair for all $\tilde{t}_0 \ge t_0$.

We perturb\index{perturb!equilibrium} 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(t, \b x_e+\Delta \b x,\b u_e+\Delta \b u)\\
&=& \b f(t, \b x_e,\b u_e)+\frac{\partial \b f}{\partial \b x}(t, \b x_e, \b u_e)\,\Delta \b x+\frac{\partial \b f}{\partial \b u}(t, \b x_e, \b u_e)\,\Delta \b u+\mb{higher order terms},
\end{eqnarray*}
where
\[
\frac{\partial \b f}{\partial \b x}(t, \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}(t, \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^{\top}& \b u_e^{\top}
\end{array}\right]^{\top}$. 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(t, \b x_e,\b u_e)=\b 0$ for all $t\ge t_0$. Let
\[
\b A(t)=\frac{\partial \b f}{\partial \b x}(t, \b x_e,\b u_e)\quad\mb{and}\quad \b B(t) = \frac{\partial \b f}{\partial \b u}(t, \b x_e,\b u_e).
\]
Neglecting higher order terms, we arrive at the linear approximation\index{approximation!linear}
\[
\frac{d}{dt}\Delta \b x=\b A(t) \Delta \b x+\b B(t) \Delta \b u.
\]
Similarly, if the outputs of the nonlinear system model are of the form
\[
\left.\begin{array}{lll}
y_1 &=& h_1(t, x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
y_2 &=& h_2(t, x_1,x_2,\ldots,x_n,u_1,u_2,\ldots,u_m)\\
{}&\vdots&{}\\
y_p &=& h_p(t, 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(t, \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,
\]
where $\b y_e= \b h(t, \b x_e,\b u_e)$. Then we obtain
\[
\Delta \b y=\b C(t) \Delta \b x+\b D(t) \Delta \b u,
\]
where
\[
\b C(t)=\frac{\partial \b h}{\partial \b x}(t, \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(t) =\frac{\partial \b h}{\partial \b u}(t, \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 pair $\left[\begin{array}{cc}\b x_e^{\top}& \b u_e^{\top}
\end{array}\right]^{\top}$.
\begin{example}\rm
In this example, we linearize the model of Watt's governor\index{Watt's governor!linearize the model} 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]^{\top}
\]
corresponding to the equilibrium control $u_e=\tau$.

Equating the right-hand sides of the above engine-governor\index{engine-governor!Watt's governor}\index{Watt's governor!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\index{Jacobian matrix} of the engine-governor\index{engine-governor!Jacobian matrix} 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\index{engine-governor!linearized model}\index{linearized model!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\index{describing function} method allows us to apply familiar frequency domain\index{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\index{limit cycles!predict} in nonlinear systems. The describing function method can be viewed as a ``harmonic linearization''\index{harmonic linearization!describing function} 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\index{sinusoid input!describing function} of known, constant amplitude. The fundamental harmonic\index{fundamental harmonic} of the element's output is compared with the input sinusoid to determine the steady-state\index{steady-state} amplitude\index{steady-state!amplitude} and phase\index{steady-state!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\index{Fourier series}, we use notion of a scalar product\index{scalar product!over a function space} over a function space that we discuss next. Our presentation of the scalar product in the context of spaces of functions\index{space!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\index{continuous!set of functions} functions on the interval $[a,b]\subset \mb{\Bb R}$.
\begin{defn}
The scalar product\index{scalar product!of two functions} of two real-valued functions\index{functions!real-valued} $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\/}\index{orthogonal!functions} if $\langle u,v \rangle = 0$. A set of functions from ${\cal C}([a,b])$ is said to be {\em mutually orthogonal\/}\index{mutually orthogonal!functions} 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\/}\index{norm!of a function} 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\/}\index{piecewise continuous!function}\index{function!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\index{subinterval!open} $(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\index{discontinuity!points of}.
\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\index{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\index{piecewise continuous} on $[a,b]$, and let $a=t_0<t_1<\ldots<t_n=b$ be a partition\index{partition!of the interval} 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\index{integrals!the sum of} 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\index{scalar product} on the space $V$ of piecewise continuous functions is positive definite\index{scalar product!positive definite}\index{positive definite!scalar product}. 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\index{function!periodic}.
\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\index{fundamental period}\index{period!fundamental} 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\index{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\index{integrable!on the interval} 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\index{left hand!limits} 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\index{derivative!piecewise continuous} are piecewise continuous\index{piecewise continuous!on the interval} 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 the above background of the Fourier series we are ready 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\index{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\index{sinusoidal signal} $x(t)=X\sin \omega t$, where $X$ is the amplitude of the input sinusoid\index{sinusoid!amplitude}. 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}
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}\rm
%Verify that the Runge method agrees with the Taylor series solution through terms of degree three.
%\end{exercise}
%\begin{exercise}\rm
%Verify that the Runge-Kutta method agrees with the Taylor series solution through terms 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} 