\documentclass{article}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphics}
\usepackage{listings}
\usepackage{tcolorbox}
\usepackage{verbatim}

\DeclareMathOperator*{\argmin}{arg\,min}

\marginparwidth 0pt
\oddsidemargin  0pt
\evensidemargin  0pt
\marginparsep 0pt
\topmargin   0pt
\textwidth   6.5in
\textheight  8.5in
\setlength{\parindent}{0em}

\newenvironment{packed_itemize}{
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{itemize}}

\newenvironment{packed_enumerate}{
\begin{enumerate}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{enumerate}}

\definecolor{codeback}{RGB}{235,255,255}

\lstset{ %
language=Matlab,                % choose the language of the code
basicstyle=\footnotesize,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\footnotesize,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=10pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{codeback},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,           % adds a frame around the code
tabsize=2,          % sets default tabsize to 2 spaces
captionpos=b,           % sets the caption-position to bottom
breaklines=true,        % sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={(*@}{@*)}          % if you want to add a comment within your code
}

\begin{document}

\begin{center}
  \Huge{Ordinary Differential Equations:}\\
  \Huge{Pendulum Problems}\\
  \Large{MATH2070: Numerical Methods in Scientific Computing I}
\end{center}
Location: http://people.sc.fsu.edu/$\sim$jburkardt/classes/math2070\_2019/odes\_pendulum/odes\_pendulum.pdf
\vskip 0.1in
\hrule
\vskip 0.1in
\begin{center}
  \includegraphics[scale=0.40]{wrecking_ball.png}
  \vskip 0.1in
  {\it{Pendulum problems can involve many modeling choices.}}
\end{center}

\begin{tcolorbox}[colback=red!5,title=Pendulum problems]
{\it{Consider various models of the pendulum problem, the effect of 
stepsize and initial conditions, the modeling on nonlinear effects and
damping, and the use of phase plots to analyze the solution.}}
\end{tcolorbox}

\section {Flash back: Phase plot analysis for the predator problem}

The predator problem was introduced in the previous lab as a test
problem and you were asked to plot the rabbit and fox populations
over time.
\begin{align*}
  r'(t) &= 2 r - 0.001 * r * f\\
  f'(t) &=- 10 f + 0.002 * r * f\\
  r(0.0) &= 5000\\
  f(0.0) &= 100
\end{align*}  
We can investigate this problem with the code {\it{predator\_euler(n)}},
where we are allowed to specify the number of time steps as an input argument.  
\vskip 0.1in
We are aware that the results produced by the Euler method are approximations.
For our experiments, we take $n$ equal steps between a starting and final
time.  We know that accuracy depends in part on the stepsize, and so we
want to feel confident that we have picked a number of steps large enough
that our approximation is reasonable.  
\vskip 0.1in
How can we do this when we are
solving a problem whose exact answer is unknown?  Sometimes, even if we
don't know the true solution, we have certain expectations of the result,
such as that it must stay positive, or be bounded, or maintain a fixed
energy.  In such cases, comparing the results of two values of $n$,
or looking at plots, can be useful in judging our results.
\vskip 0.1in
For our predator problem, using a small value of $n=100$, the 
computation blows up.  Increasing $n$ to 500 or 1000 can give a 
bounded solution that seems to oscillate and grow.  Is this growth a 
true property of the system, or is it really just a symptom of a weak
approximation?  As we did for the pendulum problem, we can try an even 
larger value of $n$, at which point it seems that we are getting a periodic
solution.
\vskip 0.1in
Because the right hand sides do not depend on time, we can get a lot
of information from a phase plot, which in this case maps $r(t)$ against
$f(t)$, or in general, $y(t)$ against $y'(t)$.
The phase plot of the data is more convincing, and should suggest that
an exact solution of this problem would be periodic.  
\vskip 0.1in
In fact, this kind of predator-prey problem (also called a Lotka-Volterra
equation) has a conserved quantity, somewhat like the energy in a physical
system.  For our problem, this quantity has the form
\begin{displaymath}
  E(r,f) = 0.002 f + 0.001 r - 10 \ln(r) - 2 \ln(f)
\end{displaymath}
and you can convince yourself, using a little algebra, that
\begin{displaymath}
  \frac{dE(r,f)}{dt} = \frac{\partial E}{\partial r} \frac{dr}{dt} + \frac{\partial E}{\partial f} \frac{df}{dt} = 0
\end{displaymath}
Hence, the phase plot is tracing out contour lines $E=\text{constant}$.


\vskip 0.1in
Note that every point in this plot represents a pair of values that can
be used to specify the initial value for an ODE problem.  In this case,
these would be the initial number of rabbits and foxes.  
\vskip 0.1in
Also, note that if we compute several solutions with differing initial
conditions, we will get a nested sequence of curves, one inside the other.
These curves can't cross, because that would mean that two different
solutions share a common point.  If we used that common point as the
initial condition, we can't get two different solution curves.
\vskip 0.1in
An interesting feature of this predator problem is that it has a fixed
point, that is, there is a pair of positive values of $r$ and $f$ for
which the solution over time does not change at all.  We can determine
these values by assuming that neither $r$ nor $f$ is zero, and then
solving
\begin{align*}
\frac{dr}{dt} &= 2 r - 0.001 * r * f = 0 &\rightarrow f = 2000\\
\frac{df}{dt} &= -10 f + 0.002 * r * f &\rightarrow r = 500
\end{align*}

\section{The pendulum problem}

We have previously mentioned the linear pendulum problem, which
models the variation in $\theta$, the angle made by a pendulum, as measured 
counterclockwise from the downward vertical position.  Here, $\theta'$
is the angular velocity, and $\theta"$ the angular acceleration.
The pendulum has a length $l$, and a mass $m$ and is subject to a 
gravitational force vector which points downward with magnitude $g$.
\vskip 0.1in
Newton's law relates force and acceleration by $F=ma$.  Because the
pendulum is constrained on a string, it can only move
in the angular direction.  Therefore, its motion is only in response
to the projection of gravity in that direction, which is $-g \sin(\theta)$.
The position $s$ of the pendulum along the circle of radius $l$ 
is related to $\theta$ by $s=l \, \theta$.  Thus, $F=ma$ for this 
problem becomes:
\begin{align*}
  F &= m \, a\\
  -m \, g \, \sin(\theta) &= m \, \frac{\partial^2 s}{\partial t^2}\\
  -m \, g \, \sin(\theta) &= m \, l \, \frac{\partial^2 \theta}{\partial t^2}\\
  \theta"(t) &= - \frac{g}{l}\sin(\theta(t))
\end{align*}
For the linear pendulum problem, we assume that we have chosen a scale
for our units so that $\frac{g}{l}=1$, and we assume that $\theta$ is
so small that we can take $\sin{\theta} \approx \theta$.  Then our
model system is:
\begin{align*}
  \theta"(t) &= - \theta(t)\\
  \theta(t_0) &= u\\
  \theta'(t_0) &= v\\
  \theta_{exact}(t-t_0) &= u \cos(t)+v\sin(t)
\end{align*}
\vskip 0.1in
For convenience, we might take $t_0=0$, $u=1$, $v=0$.  Note that
the linear pendulum model is a good approximation for small
angles $\theta$, so using $u=1$ is stretching the model beyond
the limits where it is a good approximation!
\vskip 0.1in
Here $\theta$ is the angle made with the vertical axis, measured
in the counterclockwise direction, so that $\theta=0$ corresponds
to the pendulum hanging straight down, and a value 
$\theta=\frac{\pi}{2}\approx 1.5$
would have the pendulum pointing in the direction we might think of
as 3 o'clock.
\vskip 0.1in
The quantity $\theta'$ measures the angular velocity, that is, the rate
at which the angle theta is increasing.  A positive value means the
pendulum is traveling counterclockwise.
\vskip 0.1in
The initial conditions specify the angle and angular velocity at the 
initial time.

\section{Euler method solution}

As outlined in lst week's lab, we can use the Euler method to
estimate the solution of a pendulum problem.  First, of course,
we transform the second order system into a pair of coupled
first order equations.  Then, using the default data, and specifing
a fairly small number of steps $n$ as an input quantity, we can compute with:
\begin{lstlisting}
function pendulum_euler ( n )

  a = 0.0;
  b = 6.0 * pi;
  u = 1.0;
  v = 0.0;
  y0 = [ u, v ];
  
  [ t, y ] = euler ( @(t,y)pendulum_prime(t,y), n, a, b, y0 );

  return
end
\end{lstlisting}
with the derivative specified by:
\begin{lstlisting}
function value = pendulum_prime ( t, y )

  u = y(1);
  v = y(2);

  dudt = v;
  dvdt = u;

  value = [ dudt; dvdt ];

  return
end
\end{lstlisting}

\section { What happens to the solution as we change $n$?}  

Run the code {\it{pendulum\_euler(n)}} using $n=100$ steps.  We get a 
solution, but are we comfortable with it?  If we didn't know the exact 
solution, how could we judge our result?
\vskip 0.1in
Certainly, if our approximation is good, then using $n=200$ steps
over the same time interval should not change it appreciably.  
But when we plot the new solution, it seems noticeably different.
Jumping ahead, we try $n=1000$ steps and see that the approximate
solution seems to reach roughly the same peak value at each period,
rather than growing.  This gives us some reason to trust the accuracy
of the result.
\vskip 0.1in
Since we have the exact solution for this problem, we can run the
code {\it{pendulum\_compare\_euler(n)}}, starting with $n=100$, to
see how the stepsize is affecting our accuracy.  But in general,
we won't know the solution of our problem, and so we need to rely
on comparing the solution at two stepsizes to judge whether we are
converging to a reasonable solution.
\vskip 0.1in
Another way to judge the solution arises because we happen to know that
the solution ought to be periodic.  This means that, instead of plotting
the two solution components against $t$,  we can consider plotting 
$\theta(t)$ versus $\theta'(t)$.  This is known as a {\bf{phase plot}}.
A periodic solution will repeatedly trace the same closed curve.
Let's go back to $n=100$ and look at the corresponding phase plot.
Now we can see that the trace of our approximate solution is not close to 
a closed curve.  The gradual growth in size is obvious.
{\bf{phase plot}}

\section {The pendulum with a friction term}

The pendulum model can include a term for
wind resistance or friction as a force that opposes the velocity,
with a proportionality constant of $\alpha$:
\begin{align*}
  \theta"(t) &= - \alpha \theta'(t) - \theta(t)\\
  \theta(t_0) &= u\\
  \theta'(t_0) &= v
\end{align*}
Friction means that energy is lost, and that therefore, the true
physical motion cannot be  periodic.  Instead, it must slowly decay 
towards zero.  Does our computation also reveal this behavior?
\vskip 0.1in
The function {\tt{pendulum\_friction\_euler(n,alpha)}} allows us to vary
the number of steps $n$ and the friction coefficient $\alpha$ in our code.
To begin our exploration, we can use $n=100$ and $\alpha = 0.05$.
\vskip 0.1in
While the time plots are useful, it is the phase plot that will
be easiest to interpret.
\begin{verbatim}
  n          alpha         Remarks
  100          0.05        ____________________
  500          0.05        ____________________
 1000          0.05        ____________________
 2000          0.05        ____________________
\end{verbatim}

\section {Controlling the pendulum with a friction term}

For the linear pendulum with no friction, the exact solution is
\begin{align*}
  \theta(t) &= u \cos(t) + v \sin(t) \\
  \theta'(t) &= -u \sin(t) + v \cos(t)
\end{align*}
from which we can conclude that something we can think of as
the energy $E$ satisfies:
\begin{displaymath}
  E(t) = \theta(t)^2 + \theta'(t)^2 = u^2 + v^2 = \text{constant}
\end{displaymath}
By contrast, we expect that energy of the friction pendulum will
steadily decline towards zero from its initial value of $u^2+v^2$.
Could we monitor this decline, and give the pendulum a ``kick'' from
time to time to boost its energy and keep it from running down?  
This is a simple version of a control problem.  
\vskip 0.1in
To code it is simple.  After each Euler step, we can compute the
current energy $E(t)$.  If it falls below, say, 90\% of its initial
value, we can simulate a jolt from an energy source that multiplies
$\theta$ and $\theta'$ by a factor that restores the initial energy level.
\vskip 0.1in
We can use the function {\tt{pendulum\_control\_euler(n,alpha)}}
to do these experiments.
This should make for an interesting phase plot!

\section {The nonlinear pendulum model}

A simple version of the nonlinear pendulum model has the form
\begin{align*}
  \theta"(t) &= - \sin(\theta(t))\\
  \theta(t_0) &= u\\
  \theta'(t_0) &= v
\end{align*}
Again, we might take by default the values $t_0=0$, $u=1$ and $v=0$.
We can run these experiments with {\tt{pendulum\_nonlinear\_euler(n)}}.
\vskip 0.1in
For large values of $\theta$, we might expect the linear and nonlinear
to differ in some way, but we may have no idea of how this will happen.
\vskip 0.1in
But we can easily experiment, running both models and then comparing
their plots, setting $v=0$ and choosing a sequence of $u$ values of
increasing size.  We can do this using the code
{\tt{pendulum\_angle\_euler(n, theta0)}}, where $n$ is the number of
time steps, and {\tt{theta0}} is the initial value for the angular
displacement .  What can we observe?

\section {Checking results against MATLAB's ode23() }

The solution of differential equations is so important that, of course,
MATLAB includes a number of ODE solvers.  A simple one, which we will
discuss in the next lab, is called {\tt{ode23()}}.  It has the advantage
that it automatically estimates the appropriate stepsize to get good
accuracy, so we don't need to worry about picking that ourselves.
\vskip 0.1in
The input is arranged a little differently than what we have used for
our {\tt{euler()}} code.  In particular, the starting and stopping times
are input as a vector.  Here is an example of how we might solve
the pendulum problem, as coded in {\tt{pendulum\_ode23()}}:
\begin{lstlisting}
  [ t, y ] = ode23 ( @(t,y)pendulum_prime(t,y), [tstart, tstop], y0 ]; 
\end{lstlisting}
We may find it convenient to refer to {\tt{ode23()}} to check the answers
we get with our own ODE solvers.

\section { No computing assignment for this lab!}  

\end{document}
