%\setcounter{chapter}{2}
\chapter{Linearization}
\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 dynamic behavior of many nonlinear system models can be accurately approximated by linear system models within some range of variables. Then, we can use well developed techniques for analysis and synthesis of linear systems. 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) \\
&=&h(x_0)+\frac{dh(x_0)}{dx}(x-x_0) \\
&=&y_0+\frac{dh(x_0)}{dx}(x-x_0),
\end{eqnarray*}
where $y_0=h(x_0)$.

Over a small range of $\Delta x=x-x_0$, the line $y=y_0+\frac{dh(x_0)}{dx}(x-x_0)$ is a good approximation to the curve $y=h(x)$ in a neighborhood of the operating point $x_0$---see Figure~\ref{fig:lapprox}.
\begin{figure}
\vspace{3.5in}
\hskip.7in\special{bmp:c:/mybook/lapprox.bmp x=4.67in y=3.5in}
\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}.
\begin{figure}
\vspace{3.5in}
\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$, where $g$ is the gravity constant, acting on the mass $m$ returns the pendulum to its equilibrium position. By Newton's second law \[
F=-mg \sin\theta.
\]
The equilibrium point for the pendulum is $\theta=0^{\circ}$. Note that $F(0^{\circ})=0$. Therefore,
\[
\Delta F=F,
\]
and
\[
\Delta \theta=\theta.
\]
We have
\begin{eqnarray*}
\frac{dF}{d\theta}(0) &=& -mg \cos 0^{\circ}\\
&=&-mg.
\end{eqnarray*}
Hence, the linearized model, about the equilibrium position $\theta=0^{\circ}$ has the form
\[
F=-mg \theta.
\]
Thus, for small angular displacements $\theta$, the force $F$ is proportional to the 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}
\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 $y:\mb{\Bb R}^n \to \mb{\Bb R}$, that is, $y=g(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 $g$ about the operating point $x_0$ yields
\[
y-g(\b x_0)=\nabla^Tg(\b x_0)(\b x-\b x_0)+\mb{higher order terms},
\]
where $\nabla^Tg(x_0)=\left[\begin{array}{cccc}
\left.\frac{\partial g}{\partial x_1}\right\arrowvert_{\b x=\b x_0}&\left.\frac{\partial g}{\partial x_2}\right\arrowvert_{\b x=\b x_0}&\cdots&\left.\frac{\partial g}{\partial x_n}\right\arrowvert_{\b x=\b x_0}
\end{array}\right]$. Geometrically, the linearization of $g$ 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}
\vspace{3.5in}
\hskip.6in\special{bmp:c:/mybook/lintan.bmp x=4.67in y=3.5in}
\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 dynamic system modeled by a set of nonlinear differential equations
\begin{eqnarray*}
\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{eqnarray*}
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
\[
\dot{\b x}=\b f(\b x,\b u).
\]
Let $\b u_0=\left[\begin{array}{cccc}
u_{10} & u_{20}&\cdots&u_{m0}
\end{array}\right]^T$ be a constant input that forces the system to settle into a constant, equilibrium, state $\b x_0=\left[\begin{array}{cccc}
x_{10} & x_{20}&\cdots&x_{n0}
\end{array}\right]^T$, that is,
\[
\b f(\b x_0,\b u_0)=\b 0.
\]
We now perturb the equilibrium state by allowing
\[
\b x=\b x_0+\Delta \b x,\quad \b u=\b u_0+\Delta \b u.
\]
Taylor's expansion yields
\begin{eqnarray*}
\frac{d}{dt}\b x &=& \frac{d}{dt}\b x_0+\frac{d}{dt}\Delta \b x\\
&=& \b f(\b x_0+\Delta x,\b u_0+\Delta u)\\
&=& \b f(\b x_0,\b u_0)+\frac{\partial \b f}{\partial \b x}(\b x_0, \b u_0)\,\Delta \b x+\frac{\partial \b f}{\partial \b u}(\b x_0, \b u_0)\,\Delta \b u+\mb{higher order terms},
\end{eqnarray*}
where
\[
\frac{\partial \b f}{\partial \b x}(\b x_0,\b u_0)=\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_0 \atop
\b u=\b u_0}
\!\!\mb{and}\quad\frac{\partial \b f}{\partial \b u}(\b x_0,\b u_0)=\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_0 \atop
\b u=\b u_0}
\]
are the Jacobian matrices of $\b f$ with respect $\b x$ and $\b u$ respectively, evaluated at the equilibrium point. Note that
\[
\frac{d}{dt}\b x = \frac{d}{dt}\b x_0+\frac{d}{dt}\Delta \b x=\frac{d}{dt}\Delta \b x
\]
because $\b x_0$ is constant. Furthermore $\b f(\b x_0,\b u_0)=\b 0$. Let
\[
\b A=\frac{\partial \b f}{\partial \b x}(\b x_0,\b u_0)\quad\mb{and}\quad \b B= \frac{\partial \b f}{\partial \b u}(\b x_0,\b u_0).
\]
By 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
\begin{eqnarray*}
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{eqnarray*}
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_0+\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_0,\b u_0)=\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_0 \atop
\b u=\b u_0}\!\in \mb{\Bb R}^{p\times n},
\]
and
\[
\b D=\frac{\partial \b h}{\partial \b u}(\b x_0,\b u_0)=\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_0 \atop
\b u=\b u_0}\!\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 $[\b x_0^T\,\b u_0^T]^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:
\begin{eqnarray*}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& \frac{1}{2}N^2x_3^2\sin 2x_1-g\sin x_1 - \frac{b}{m}x_2\\
\dot{x}_3 &=& \frac{\kappa}{I}\cos x_1-\frac{\tau}{I}.
\end{eqnarray*}
We assume that the equilibrium state has the form
\[
\b x_e=\left[\begin{array}{ccc}
x_{1e}&0&x_{3e}
\end{array}\right]^T.
\]
Equating the right-hand sides of the engine-governor equations to zero yields
\begin{eqnarray*}
x_{2e} &=& 0\\
N^2\,x_{3e}^2 &=& \frac{g}{\cos x_{1e}}\\
\cos x_{1e} &=& \frac{\tau}{\kappa}.
\end{eqnarray*}
The Jacobian matrix of the engine-governor model evaluated at the equilibrium state is
\[
\left.\frac{\partial \b f}{\partial \b x}\right\arrowvert_{\b x=\b x_e}=\left[\begin{array}{ccc}
0 & 1 & 0\\
N^2x_{3e}^2\cos 2x_{1e}-g\cos x_{1e} & -\frac{b}{m} & N^2x_{3e}\sin 2x_{1e}\\
-\frac{\kappa}{I}\sin x_{1e} & 0 & 0
\end{array}\right].
\]
We use the 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}{\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}=\left[\begin{array}{ccc}
0 & 1 & 0\\
-\frac{g \sin^2 x_{1e}}{\cos x_{1e}} & -\frac{b}{m} & \frac{2g\sin\,x_{1e}}{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}
\left[\begin{array}{c}
\Delta \dot{x}_1\\
\Delta \dot{x}_2\\
\Delta \dot{x}_3
\end{array}\right]=\left[\begin{array}{ccc}
0 & 1 & 0\\
-\frac{g \sin^2 x_{1e}}{\cos x_{1e}} & -\frac{b}{m} & \frac{2g\sin\,x_{1e}}{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].
\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 was developed during World War II. The first to introduce the method in the open literature seems to be A.~Tustin of England.

The describing function method is based on the Fourier series. In our review of Fourier series we use the notion and properties of general scalar product that we discuss next. The reader familiar with the Fourier series may go directly to Subsection~\ref{subsec:desfunana} where the describing function method analysis is presented.

\subsection{General Scalar Product}
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}
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 $V$ 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$ and $\phi_n(t)=\cos\,n\omega t$, $n=1,2,\ldots$, 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 positive $m$, $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
Let $\psi_n(t)=\sin\,n\omega t$ be defined 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}
The graph of a piecewise continuous function is shown in Figure~\ref{fig_piececon}.
\begin{figure}
\vspace{3in}
\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}
\vspace{3in}
\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=0= \langle f,f \rangle$. 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 $x$ 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}
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 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 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 periodia with period $2\pi/\omega$ and integrable on the interval $[t_0,t_0+2\pi/\omega]$. Then, coefficients $A_n$ and $B_n$ 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, if 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 $t$ is
\[
\left(f(t+)+f(t-) \right)/2.
\]
Note that at any point $t$ where $f$ is continuous, we have
\[
f(t+)=f(t-)=f(t).
\]
\begin{thm}
Assume that $f$ and its derivative are piecewise continuous on the interval $[t_0, t_0+2\pi/\omega]$. 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.
\label{thm:fourierb}
\end{thm}

With this knowledge of the Fourier series it is now possible to discuss the describing function method of analyzing a class of nonlinear control systems.

\subsection{Describing Function 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, sometimes we 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}
\vspace{3in}
\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 since 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}
\vspace{3in}
\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}
\vspace{4in}
\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)$ is 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}\rm
Consider the nonlinear control system shown in Figure~\ref{fig:a_96_f1}.
\begin{figure}
\vspace{3.5in}
\hskip.5in\special{bmp:c:/mybook/a_96_f1.bmp x=4.67in y=3.5in}
\caption{Nonlinear feedback control system.}
\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,
\]
that is drawn vertically down from the relay characteristic in Figure~\ref{fig:relayde}.
\begin{figure}
\vspace{4.5in}
\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}
\vspace{3in}
\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}
\vspace{3in}
\hskip1.2in\special{bmp:c:/mybook/nyqdesf.bmp x=4in y=3in}
\caption{Plots of $G(j\omega)$ and $-1/N$.}
\label{fig:nyqdesf}
\end{figure}
The amplitude of the limit cycle is obtained from the relation, $G(j\omega)=-1/N$. We have
\[
\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(j3) \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 behavior of numerous dynamical systems can be modeled by a set of first-order differential equations. Almost always, however, a linear model is an approximation to a nonlinear model. This point was illustrated with an example of flyball governor. Linear models dominate the literature. Yet nature seems to be all but linear. Modeling of natural objects and systems is the subject of Chapter~\ref{chp:chaos}.
\section{Exercises}
\begin{exercise}\rm
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}
\vspace{3.5in}
\hskip.7in \special{bmp:c:/mybook/cyclotrn.bmp x=4.67in y=3.5in}
\caption{Schematic of a cyclotron.}
\label{fig:cyclotron}
\end{figure}
The operation of a cyclotron can be described briefly 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 denoted by $\Delta B$. As the ion speed increases, the relativistic mass $m$ also 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 t}{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)$.
\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}