\documentclass{simiode}

%*********** Toggle the following line to turn the 
%            teacher version on (1) or off (0). ***********

\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}

\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{enumitem}
\usepackage{hyperref}
%\hypersetup{
%  colorlinks   = true, % Color links instead of boxes
%  urlcolor     = blue, % Color for external hyperlinks
%  linkcolor    = blue, % Color of internal links
%  citecolor    = blue  % Color of citations
%}
\usepackage{url}  
\usepackage{caption}
\usepackage{subcaption}
\captionsetup[figure]{labelfont=bf,labelsep=period}
\captionsetup[subfigure]{labelformat=parens}

\renewcommand{\baselinestretch}{1.3}


\begin{document}

\newpage

\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\ 
AN INTERACTIVE ILLUSTRATION OF THE \\ VAN DER POL OSCILLATOR}



\markboth{Van der Pol Oscillator}{Van der Pol Oscillator}

\author{Mark A. Lau Kwan\\
University of Turabo\\
Gurabo  00777 PUERTO RICO\\
mlau@suagm.edu
}

\makeStitlePDFLaTex


\begin{abstract}
This paper presents an electronic spreadsheet 
model of the Van der Pol oscillator, a well-known 
nonlinear second-order ordinary differential 
equation. 
The spreadsheet features a number of dynamic 
controls that permit the user to alter the 
parameters of the Van der Pol equation and 
explore the effects of such changes on the 
solution. 
Furthermore, the spreadsheet exploits its 
graphics capabilities to produce a rendition 
of the limit cycle associated with the 
(periodic) steady-state solution to this 
nonlinear equation. 
The Van der Pol equation models many 
phenomena in which periodic behavior occurs 
such as the oscillations encountered in electronic 
circuits and biological systems. 
After investigating the behavior of the Van der Pol 
oscillator with the accompanying spreadsheet, the 
reader will have gained a better appreciation of 
limit cycles and nonlinear dynamics. 
\end{abstract}

\section*{SCENARIO DESCRIPTION}

The purpose of this paper is to present an 
interactive spreadsheet application that 
allows the reader to investigate the Van der 
Pol oscillator. 
Through a number of suggested activities, 
the reader will gain a better understanding 
of the qualitative behavior of solutions as 
the parameters of the governing differential 
equation are varied. 

\subsection*{A Brief Description of the 
Van der Pol Oscillator} 

The Van der Pol oscillator is described by 
the nonlinear, second-order ordinary 
differential equation 
\begin{equation} 
  \ddot{x} + \mu(x^2 - 1)\dot{x} 
    + \omega^2 x = 0, 
\label{vdp}
\end{equation} 
where $\mu$ and $\omega$ are positive parameters, 
and $x$ (a function of time~$t$) is the quantity 
of interest (e.g., position or voltage). 
The quantities $\dot{x}$ and $\ddot{x}$ 
denote the first- and second-order derivative 
with respect to time~$t$, respectively. 

This celebrated equation~\eqref{vdp} is 
attributed to the Dutch electrical engineer 
Balthasar Van der Pol, who pioneered 
experimental investigations in nonlinear 
dynamics in the early days of radio and 
telecommunications. 
He proposed his equation in order to 
describe the nonlinear oscillations 
produced by self-sustained oscillating 
triode circuits~\cite{vanderpol}. 
His equation has been extensively studied 
and finds numerous applications in 
engineering, physics, biology, sociology, 
and economics. 
 
The Van der Pol equation~\eqref{vdp} 
introduces a nonlinear damping term, 
namely $\mu(x^2 - 1)\dot{x}$. 
Unlike the familiar linearly damped 
harmonic oscillator, the presence of the 
nonlinear damping term can introduce 
negative damping depending on the magnitude 
of~$x$. 
More specifically, if $|x| > 1$, there is 
\textit{positive} damping and energy is 
dissipated from the system resulting 
in decaying motion (as in the linearly 
damped harmonic oscillator). 
On the other hand, if $|x| < 1$, there is 
\textit{negative} damping and energy is 
pumped into the system resulting in 
amplification of motion. 
Therefore, it is reasonable to expect 
that the steady-state motion consists of 
sustained periodic oscillations (e.g., the 
system may start at a small value of~$x$, 
then be driven to a large value of~$x$ by 
the amplification, then damped back to a 
small~$x$ by the decay, and repeat this 
process {\it ad infinitum\/}). 
These periodic oscillations are referred 
to as \textit{limit cycles}. 
The interested reader can learn more about 
nonlinear dynamics from these excellent 
sources~\cite{khalil}, \cite{strogatz}. 


\subsection*{The Interactive Spreadsheet Model 
for Investigating the Van der Pol Oscillator} 

The accompanying electronic spreadsheet 
model has been created using Microsoft Excel. 
The user interface is shown in 
Figure~\ref{limitcycle}. 
The top left section of the spreadsheet shows 
a number of active controls for changing the 
parameters and initial conditions for the 
Van der Pol equation~(\ref{vdp}). 
To change the damping $\mu$, the frequency 
$\omega$, the initial position $x(0)$, or the 
initial velocity $\dot{x}(0)$, the user simply 
slides the thumb on the corresponding scroll 
bar; the values selected by the user will be 
automatically adjusted and displayed in the 
cell range \textsf{K3:K6}. 
As the user varies the parameters and/or the 
initial conditions, the phase plane portrait 
of the limit cycle is also automatically 
adjusted to reflect the changes introduced by 
the user. 

\begin{figure} 
\centering{\includegraphics[width=15 cm]{3-4-VanderPol-JPG-Figure1.jpg}} 
\caption{Main user interface of the Microsoft 
Excel spreadsheet for investigating the behavior of the 
Van der Pol oscillator. 
The screen capture displays the limit cycle for $\mu = 2$ 
(only the first $28$~rows of the spreadsheet are shown here).} 
\label{limitcycle}
\end{figure} 

To obtain the limit cycle, the Van der Pol 
equation~(\ref{vdp}) is solved numerically 
within the spreadsheet. 
The numerical solution requires 
rewriting~\eqref{vdp} as a set of first-order 
differential equations.  
Defining $x_1 = x$ and $x_2 = \dot{x}$, we 
obtain the following system: 
\begin{equation} 
  \begin{array}{rcl} 
    \dot{x}_1 &\!\!\!=\!\!\!& x_2, \\ 
    \dot{x}_2 &\!\!\!=\!\!\!& - \mu(x_1^2 - 1)x_2 
      - \omega^2 x_1. 
  \end{array} 
\label{vdpset} 
\end{equation} 
The set of differential equations~\eqref{vdpset} 
is solved in the spreadsheet using the Runge-Kutta 
method. 
The tabulated solutions for the state variables 
$x_1(t)$ and $x_2(t)$, as well as the Runge-Kutta 
increments, are shown in the spreadsheet. 
If necessary, the user can manually alter the 
time step size, $\Delta t$, stored in cell 
\textsf{D8} to obtain more refined solutions; 
making $\Delta t$ too small may require the 
user to copy the pre-programmed formulas 
occupying the cell range \textsf{O10:AA10} 
(or any complete set of formulas below row 
\textsf{10} spanning columns \textsf{O} 
through \textsf{AA}) in order to generate 
sufficient points to plot the limit cycle. 


\subsubsection*{Activity 1. Harmonic Oscillator}

Assume the initial conditions 
$x_1(0) = x(0) = x_0$ and 
$x_2(0) = \dot{x}(0) = y_0$ (not both zero), 
and set $\mu = 0$ in \eqref{vdp}, 
which will then become the \textit{harmonic 
oscillator} 
\begin{equation}
  \ddot{x} + \omega^2 x = 0.  
\label{harmonic}
\end{equation} 
Show that the trajectory (solution to \eqref{harmonic}) is an 
ellipse described by 
\begin{equation}
  \frac{x_1^2}{\displaystyle{x_0^2 
  + \frac{y_0^2}{\omega^2}}}
  + \frac{x_2^2}{\omega^2 x_0^2 + y_0^2} 
  = 1. 
\end{equation}
Confirm this with the spreadsheet 
application by setting the parameters 
$\mu = 0$, $\omega \ne 0$, and the 
initial conditions to any allowable values 
(not both zero). 

For what value of $\omega$ does the 
trajectory become a circle? 
For this value of $\omega$, what is 
the radius of the circle? 



\subsubsection*{Activity 2. Rigid Body}

Assume the initial conditions 
$x_1(0) = x(0) = x_0$ and 
$x_2(0) = \dot{x}(0) = y_0$ (not both zero), 
and set both $\mu = 0$ and $\omega = 0$ 
in \eqref{vdp}, which will result 
in the \textit{rigid body} described by 
\begin{equation}
  \ddot{x} = 0.  
\label{rigid}
\end{equation} 
Observe the phase plane portrait displayed 
by the spreadsheet application. 
Provide an analytical justification to your 
observations. 



\subsubsection*{Activity 3. Empty Phase Portrait}

Set both initial conditions to zero 
$\big(x_1(0) = x(0) = 0$ and 
$x_2(0) = \dot{x}(0) = 0\big)$ on the 
spreadsheet application; any values 
of $\mu$ and $\omega$ allowed by the 
spreadsheet can be used. 
In this case, the phase portrait will 
appear empty. 
Give an explanation. 



\subsubsection*{Activity 4. Behavior 
of Van der Pol Oscillator for Small 
$\boldsymbol{\mu}$}

Describe the trajectory and the shape of the 
limit cycle for the Van der Pol oscillator 
when $\mu$ is small. 
To do this using the spreadsheet model, set 
$\mu$ to a small value, say $\mu = 0.2$ or 
less; set the other parameters $\omega$, 
$x(0)$, and $\dot{x}(0)$ to any nonzero 
values permitted by the spreadsheet 
(remember not to set both 
initial conditions to zero). 
If desired, use the \textsf{Scattered Chart} 
feature of the spreadsheet to obtain the time 
history plots for $x_1(t)$ and $x_2(t)$ from 
their tabulated values. 



\subsubsection*{Activity 5. Behavior 
of Van der Pol Oscillator for Large 
$\boldsymbol{\mu}$}

Describe the shape of the limit cycle for the 
Van der Pol oscillator when $\mu$ is large. 
To see this on the spreadsheet model, set 
$\mu$ to a large value, say $\mu = 5$ or 
greater; set the other parameters $\omega$, 
$x(0)$, and $\dot{x}(0)$ to any nonzero 
values permitted by the spreadsheet 
(remember not to set both 
initial conditions to zero). 
If desired, use the \textsf{Scattered Chart} 
feature of the spreadsheet to obtain the time 
history plots for $x_1(t)$ and $x_2(t)$ from 
their tabulated values. 



\subsubsection*{Activity 6. Existence of Limit Cycle}

In this section, the reader will be able to 
show the existence of a stable limit cycle 
for the Van der Pol equation~\eqref{vdp}. 
To do this, consider the \textit{Li\'{e}nard 
equation} 
\begin{equation} 
  \ddot{x} + u(x)\dot{x} + v(x) = 0, 
\label{lienard} 
\end{equation} 
which can be written as the system 
\begin{equation} 
  \begin{array}{rcl} 
    \dot{x}_1 &\!\!\!=\!\!\!& x_2, \\ 
    \dot{x}_2 &\!\!\!=\!\!\!& -u(x)x_2 - v(x). 
  \end{array} 
\label{lienardset}
\end{equation}
Then apply the \textit{Levinson-Smith 
theorem}~\cite{mattuck} given below to 
show the existence of a limit cycle for 
\eqref{vdp}.  

\textbf{Levinson-Smith Theorem.}\footnote{The 
Levinson-Smith theorem is reproduced 
from~\cite{mattuck} and is provided here 
for the convenience of the reader in showing 
the existence of a limit cycle for the Van 
der Pol equation~\eqref{vdp}.} 
Given equation~\eqref{lienard}, suppose the 
following conditions are satisfied: 
\begin{enumerate}[label=(\alph*)]
  \item $u(x)$ is even and continuous. 
  \item $v(x)$ is odd, $v(x) > 0$ if $x > 0$, 
        and $v(x)$ is continuous for all $x$. 
  \item $V(x) \rightarrow \infty$ as 
        $x \rightarrow \infty$, where 
        $$ V(x) = \int_{0}^{x} v(t)\,dt. $$ 
  \item For some $k > 0$, the following holds: 
        \begin{itemize} 
          \item $U(x) < 0$, for $0 < x < k$, 
          \item $U(x) > 0$ and increasing, 
                for $x > k$, 
          \item $U(x) \rightarrow \infty$ as 
                $x \rightarrow \infty$, 
        \end{itemize} 
        where 
        $$ U(x) = \int_{0}^{x} u(t)\,dt. $$ 
\end{enumerate} 
Then, the system given in 
\eqref{lienardset}, and 
hence~\eqref{lienard}, has 
\begin{enumerate}[label=\roman*.] 
  \item a unique critical point at the 
        origin; 
  \item a unique non-zero closed 
        trajectory~$C$, which is a stable 
        limit cycle around the origin; 
  \item all other non-zero trajectories 
        spiraling toward~$C$ as 
        $t \rightarrow \infty$. 
\end{enumerate} 


\subsubsection*{Activity 7. Model of Aortic Blood Flow}

One application of the Van der Pol oscillator is the 
modeling of biological systems in which periodic (or 
nearly periodic) phenomena occur. 
The work presented in~\cite{atlas} uses precisely 
\eqref{vdp} to model blood flow in the human 
aorta. 
To run a simulation on the accompanying spreadsheet, 
set $\mu = 1$, $\omega = 2.37\ (\omega^2 = 5.6)$, 
$x(0) = 1$, and $\dot{x}(0) = -3$. 
The reader can use the tabulated values to obtain 
a pattern of blood flow $(x_1\ \mbox{vs.}\ t)$ using 
the \textsf{Scattered Chart} feature of Excel. 
The pattern will not match that of~\cite{atlas} 
because the solution to the Van der Pol equation 
requires some post-processing in order to produce 
a more meaningful physiological signal. 
Nevertheless, this example serves to illustrate 
the usefulness of the Van der Pol oscillator in 
real-life applications. 



\section*{CONCLUDING REMARKS} 

This paper presented the Van der Pol oscillator 
and how its qualitative behavior can be 
investigated using an electronic spreadsheet 
model. 
The spreadsheet features a number of dynamic 
controls that permit the user to alter the 
parameters of the Van der Pol equation and 
interactively explore the effect of such 
changes on the solution and the shape of the 
limit cycle as the damping and other parameters 
change. 
The accompanying spreadsheet can be used to 
investigate other problems arising from physics, 
engineering, and biology, among others, in which periodic 
(or nearly periodic) phenomena occur and the Van der 
Pol oscillator can be profitably used to gain 
further insight or to obtain approximate solutions. 


\begin{thebibliography}{98}

\bibitem{atlas}
Atlas, G. M. and Desiderio, M. C. 2006. 
Solutions to the Van der Pol Equation: 
A Model of Aortic Blood Flow. 
textit{Proceedings of the IEEE 32nd Annual Northeast 
Bioengineering Conference}. 143--144.

\bibitem{khalil}
Khalil, H. K. 2001. 
\textit{Nonlinear Systems, Third Edition}. 
Upper Saddle River  NJ: Pearson Education. 

\bibitem{mattuck}
Mattuck, A. and M.I.T. 2011. 
\textit{Limit Cycles}. 
\href{http://math.mit.edu/\~jorloff/suppnotes/suppnotes03/lc.pdf}{http://math.mit.edu/\~jorloff/suppnotes/suppnotes03/lc.pdf}. 
Accessed March 1, 2017. 

\bibitem{strogatz}
Strogatz, S. H. 2014. 
\textit{Nonlinear Dynamics and Chaos, 
Second Edition}. 
Boulder  CO: Westview Press.  

\bibitem{vanderpol} 
Van der Pol, B. 1927. 
On Relaxation-Oscillations. 
\textit{The London, Edinburgh and Dublin Philosophical \& 
Journal of Science}. 2(7): 978--992.

\end{thebibliography}


\teacher{
\section*{NOTES FOR TEACHERS}

\subsection*{\textit{Solution to Activity 1. 
Harmonic Oscillator}} 

When $\mu = 0$, \eqref{vdp} reduces 
to the harmonic oscillator 
$$ \ddot{x} + \omega^2 x = 0. $$ 
With the initial conditions $x(0) = x_0$ and 
$\dot{x}(0) = y_0$, the solution becomes  
$$ x(t) = x_0\cos(\omega t) 
     + \frac{y_0}{\omega}\sin(\omega t). $$ 
Since the state variables are defined as 
$x_1 = x$ and $x_2 = \dot{x}$, then 
\begin{eqnarray*} 
  x_1(t) &\!\!\!=\!\!\!& x_0\cos(\omega t) 
     + \frac{y_0}{\omega}\sin(\omega t), \\[2mm] 
  x_2(t) &\!\!\!=\!\!\!& -\omega x_0\sin(\omega t) 
     + y_0\cos(\omega t). 
\end{eqnarray*} 
Computing $\omega^2 x_1^2(t)$ and $x_2^2(t)$, 
and adding them, after some simplification,   leads to
$$ \omega^2 x_1^2 + x_2^2 
     = \omega^2 x_0^2 + y_0^2\, , $$ 
or, after dividing through by $\omega^2 x_0^2 
+ y_0^2$, 
$$ \frac{x_1^2}{\displaystyle{x_0^2 
  + \frac{y_0^2}{\omega^2}}}
  + \frac{x_2^2}{\omega^2 x_0^2 + y_0^2} 
  = 1, $$  
which corresponds to the equation of an 
ellipse with center at the origin and 
semiaxes given by 
$$ a = \sqrt{x_0^2 
  + \frac{y_0^2}{\omega^2}} 
  \qquad \mbox{and} \qquad 
  b = \sqrt{\omega^2 x_0^2 + y_0^2}. $$ 
Notice that the equation describing the 
ellipse does not depend on time~$t$. 

The ellipse becomes a circle when 
$\omega = 1$. 
The equation of the circle will be 
$$  x_1^2 + x_2^2 = x_0^2 + y_0^2\, , $$ 
with center at $(0,0)$ and radius 
$r = \sqrt{x_0^2 + y_0^2}$. 



\subsection*{\textit{Solution to Activity 2. 
Rigid Body}} 

The solution to the differential 
equation~\eqref{rigid}, subject 
to the initial conditions $x_1(0) = x(0) = x_0$ 
and $x_2(0) = \dot{x}(0) = y_0$ 
(not both zero), is 
$$ x(t) = y_0 t + x_0. $$ 
Now observe that 
$$ x_2(t) = \dot{x}(t) = y_0, $$ 
which is a constant, independent of~$t$ 
and independent of~$x_1 = x$. 
Hence, the phase portrait is a horizontal 
line whose equation is given by $x_2 = y_0$. 



\subsection*{\textit{Solution to Activity 3. 
Empty Phase Portrait}} 

With zero initial conditions, the point 
$(x_1,x_2) = (0,0)$ is an equilibrium point 
of the system given in \eqref{vdpset}. 
Therefore, if unperturbed, the trajectory will 
consist of the origin (a single point). 
The origin is an unstable equilibrium point, 
however. 



\subsection*{\textit{Solution to Activity 4. 
Behavior of Van der Pol Oscillator for 
Small $\boldsymbol{\mu}$}}

Figure~\ref{fig2} shows the behavior of the 
Van der Pol oscillator for $\mu = 0.2$, 
$\omega = 1$, $x_1(0) = 1$, and $x_2(0) = 1$. 

\begin{figure}[!ht]	
  \centering
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure2a.jpg}
    \caption{Limit cycle.}
    \label{fig2a}
  \end{subfigure}
  \newline
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure2b.jpg}
    \caption{Position $x_1(t)$.}
    \label{fig2b}
  \end{subfigure}
  \quad
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure2c.jpg}
    \caption{Velocity $x_2(t)$.}
    \label{fig2c}
  \end{subfigure} 
\caption{Behavior of Van der Pol oscillator for small $\mu$ 
$\big(\mu = 0.2,\ \omega = 1,\ x_1(0) = 1, x_2(0) = 1\big)$.}
\label{fig2}
\end{figure}

As can be seen from Figure~\ref{fig2}(a), 
all trajectories emanating from 
$\big(x_1(0),x_2(0)\big) \neq \big(0,0\big)$ 
converge to a limit cycle that resembles the 
periodic orbit exhibited by a harmonic oscillator. 
Figures~\ref{fig2}(b) and~\ref{fig2}(c) confirm 
this observation since the time history plots 
for $x_1$ and $x_2$ are oscillatory, with 
slightly time-varying amplitudes and phases. 



\subsection*{\textit{Solution to Activity 5. 
Behavior of Van der Pol Oscillator for 
Large $\boldsymbol{\mu}$}}  

Figure~\ref{fig3} shows the behavior of the 
Van der Pol oscillator for the choice $\mu = 5$, 
$\omega = 1$, $x_1(0) = 1$, and $x_2(0) = 1$. 

\begin{figure}[!ht]	
  \centering
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure3a.jpg}
    \caption{Limit cycle.}
    \label{fig3a}
  \end{subfigure}
  \newline
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure3b.jpg}
    \caption{Position $x_1(t)$.}
    \label{fig3b}
  \end{subfigure}
  \quad
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure3c.jpg}
    \caption{Velocity $x_2(t)$.}
    \label{fig3c}
  \end{subfigure} 
\caption{Behavior of Van der Pol oscillator for large $\mu$ 
$\big(\mu = 5,\ \omega = 1,\ x_1(0) = 1, x_2(0) = 1\big)$.}
\label{fig3}
\end{figure}

As can be seen from Figure~\ref{fig3}(a), 
the limit cycle has markedly different 
characteristics from those corresponding 
to small~$\mu$. 
For large~$\mu$, the position $x_1$ varies 
very slowly for an extended period of time, 
where the velocity $x_2$ is very small; 
the position $x_1$ then changes very rapidly 
over a short period of time, where the 
velocity $x_2$ reaches a large value until 
the position starts to vary slowly and the 
process repeats ad infinitum, as shown in 
Figures~\ref{fig3}(b) and~\ref{fig3}(c). 
This type of behavior is referred to as  
\textit{relaxation oscillation}~\cite{vanderpol}. 



\subsection*{\textit{Solution to Activity 6. 
Existence of Limit Cycle}} 

In applying the Levinson-Smith theorem 
to the Van der Pol equation~\eqref{vdp}, 
the following functions can be identified 
$$ u(x) = \mu(x^2 - 1), \qquad 
   v(x) = \omega^2 x, \qquad 
   U(x) = \mu(\frac{x^3}{3} - x), 
   \qquad \mbox{and} \qquad 
   V(x) = \frac{\omega^2}{2}x^2. $$ 
Clearly, $u(x)$ is even 
$\big(u(x) = u(-x)\big)$ and 
continuous $\big(u(x)$ is a polynomial 
function$\big)$. 
The function $v(x)$ is odd 
$\big(v(x) = -v(-x)\big)$, $v(x) > 0$ 
for $x > 0$, and continuous $\big($the 
graph of $v(x)$ is a straight line 
passing through the origin with a positive 
slope of~$\omega^2\big)$. 
The function $V(x)$ becomes unbounded when 
$x$ grows very large. 
As for the function~$U(x)$, it can be 
said that 
\begin{eqnarray*} 
  && U(x) < 0, \ \mbox{for}\ 
     0 < x < \sqrt{3}, \\ 
  && U(x) > 0\ \mbox{and increasing for}\ 
       x > \sqrt{3}, \ \mbox{and}\\ 
  && U(x) \ \mbox{becomes unbounded as}\ 
       x\ \mbox{grows very large}. 
\end{eqnarray*} 
Since all of the conditions of the 
Levinson-Smith theorem are fulfilled, 
it can be concluded that the Van der Pol 
equation~\eqref{vdp} -- or the 
system~\eqref{vdpset} -- has a stable 
limit cycle. 



\subsection*{\textit{Solution to Activity 7. 
Model of Aortic Blood Flow}} 

This is not a solution per se. 
Rather, the comment emphasizes the limitations 
of the solution obtained with the spreadsheet. 
The solution shown in Figure~\ref{fig4} requires 
post-processing (e.g., restore the time and 
amplitude scales, restore a baseline, and so on) 
so as to produce a more realistic physiological 
signal~\cite{atlas}. 
The solution in Figure~\ref{fig4} should be taken 
as a first approximation; further refinements are 
required to obtain a reasonable model of aortic 
blood flow. 


\begin{figure}[!t]	
  \centering
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure4a.jpg}
    \caption{Limit cycle.}
    \label{fig4a}
  \end{subfigure}
  \newline
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure4b.jpg}
    \caption{Flow $x_1(t)$.}
    \label{fig4b}
  \end{subfigure}
  \quad
  \begin{subfigure}[b]{2.9in}
    \centering
    \includegraphics[width=2.9in]{3-4-VanderPol-JPG-Figure4c.jpg}
    \caption{Flow velocity $x_2(t)$.}
    \label{fig4c}
  \end{subfigure} 
\caption{Aortic blood flow 
$\big(\mu = 1,\ \omega = 2.37,\ x_1(0) = 1, x_2(0) = -3\big)$. 
The solution provided by the spreadsheet is only an 
approximation and requires post-processing so as to 
produce a more realistic physiological signal.}
\label{fig4}
\end{figure}


} 

\end{document}
