====== Math 445 lab 11: tick tock, the pendulum clock ======

==== Mathematical models of the plane pendulum ====

In class we derived four different mathematical models of the plane pendulum, listed here in order of decreasing mathematical complexity. 

** Nonlinear damped pendulum.** This is the most physically realistic model. It includes a linear model of air resistance the $\alpha/m \; d\theta/dt$ term and is accurate for large displacement angles $\theta$.
\begin{eqnarray*}
\frac{d^2\theta}{dt^2} + \frac{\alpha}{m} \frac{d\theta}{dt} + \frac{g}{\ell} \sin \theta = 0
\end{eqnarray*}

** Nonlinear undamped pendulum.** This model neglects air resistance. It is derived from the previous model by setting the air resistance constant $\alpha$ to zero.
\begin{eqnarray*}
\frac{d^2\theta}{dt^2} + \frac{g}{\ell} \sin \theta = 0
\end{eqnarray*}

** Linear damped pendulum.** This is valid only for small oscillations. You can derive it from the nonlinear damped pendulum model by substituting the //small angle approximation// $\sin \theta \approx \theta$ for small $\theta$. 
\begin{eqnarray*}
\frac{d^2\theta}{dt^2} + \frac{\alpha}{m} \frac{d\theta}{dt} + \frac{g}{\ell} \theta = 0
\end{eqnarray*}

** Linear undamped pendulum.** The simplest pendulum model. It neglects air resistance and employs the small-angle approximation. 
\begin{eqnarray*}
\frac{d^2\theta}{dt^2} + \frac{g}{\ell} \theta = 0
\end{eqnarray*}

In each of these, the variables are

  * $\theta(t)$: the unknown function that describes the angle of the pendulum away from vertical as a function of time.
  * $m$: the mass of the pendulum bob.
  * $\ell$: the length of the pendulum string
  * $g$: the acceleration of gravity
  * $\alpha$: an air resistance constant


In class we showed that the linear undamped pendulum has a solution of the form $\theta(t) = \theta_0 \cos \omega t$ where $\omega = \sqrt{g/\ell}$ and $\theta_0$ is an arbitrary initial angular displacement. (Note that $\theta_0$ must be small for the small-angle approximation to be valid!)

We also showed that the 2nd-order equation in the scalar variable $\theta(t)$ can be transformed into a 1st order equation in the vector variable $\vec{x}(t)$, using the substitution

\begin{eqnarray*}
\vec{x} = \left(\begin{array}{l} x_1\\ x_2 \end{array} \right) = \left(\begin{array}{l} \theta \\ d\theta/dt \end{array} \right)
\end{eqnarray*}

If you differentiate the above equation in time and perform a few substitutions between $\theta, d\theta/dt$ and $x_1, x_2$, you get 

\begin{eqnarray*}
\frac{d\vec{x}}{dt} = \left(\begin{array}{l} dx_1/dt\\ dx_2/dt \end{array} \right) = \left(\begin{array}{cc} 0 & 1 \\ -g/\ell & 0 \end{array} \right) \left(\begin{array}{c} x_1 \\ x_2 \end{array} \right)
\end{eqnarray*}

This equation is now of the form $d\vec{x}/dt = \vec{f}(\vec{x}) = A\vec{x}$, so it can be solved numerically with Matlab's ''ode45'' function.

----

====Problem 1: the linear undamped pendulum====

Compare the analytic solution $\theta(t) = \theta_0 \cos \omega t$ of the linear undamped pendulum to a numerical solution computed with Matlab's ''ode45'' function. Use constants

  * $g = 9.8$ (meters per second^2)
  * $\ell = 1.0$ (meters)
  * $\theta_0 = 0.1$ (radians)
  * $m=1$ (kilogram)

Make three plots: 
 
   - $\theta(t)$ versus $t$ (a //time series// plot) 
   - $d\theta/dt$ versus $t$ (a //time series// plot)
   - $d\theta/dt$ versus $\theta$ (a //phase portrait//)

For the first two plots, show the analytic solution in solid lines and the numerical solution with solid lines with superimposed dots. 

For the $d\theta/dt$ versus $\theta$ plot, plot just the numerical solution.  Show a number of numerical solutions, using initial conditions $\vec{x}_0 = [\theta_0, 0]$ with $\theta_0$ varying from $0.01$ to $0.10$ in steps of $0.01$. Superimpose a quiver plot that shows the vector field $\vec{f}(\vec{x})$ which governs the time-evolution of the pendulum. 

The plots should look something like this:

{{:gibson:teaching:spring-2016:math445:lab11:problem1a.png?400|}} {{:gibson:teaching:spring-2016:math445:lab11:problem1b.png?400|}}

----

====Problem 2: the linear damped pendulum==== 

If we include the damping of air resistance, the equations of motion of the linear damped pendulum are 

\begin{eqnarray*}
\frac{d}{dt} \left( \begin{array}{c} x_1 \\ x_2 \end{array} \right) = \left[ \begin{array}{cc} 0 & 1 \\ -g/l & -\alpha/m \end{array} \right] \left( \begin{array}{c} x_1 \\ x_2 \end{array} \right)
\end{eqnarray*}

and the time series and phase portrait for (for $g=9.8$, $\ell=1$, $m=1$, and $\alpha=1$) look like

{{:gibson:teaching:spring-2016:math445:lab11:problem2a.png?direct&400|}} {{:gibson:teaching:spring-2016:math445:problem2b.png?direct&400|}}

Now the temporal oscillations get smaller and smaller as time goes on. The phase portrait shows that all initial conditions eventually spiral into the origin, i.e. the pendulum hangs straight down ($\theta = 0$) and doesn't move $d\theta/dt = 0$).

Write Matlab code to simulate the linear damped pendulum and reproduce the above plots.

----
====Problem 3: nonlinear undamped pendulum====

If $\theta$ gets large the approximation $\sin \theta \approx \theta$ is no longer valid. This lab problem is an investigation of how the frequency of oscillation of the pendulum changes when the amplitude of the oscillation is large. Specifically, you'll determine if a more realistic nonlinear pendulum clock runs faster or slower than the idealized linear clock, and by how much, depending on amplitude.

**(a)** Starting from the nonlinear undamped pendulum equation, revise the derivation from class to develop a nonlinear system of differential equations $d\vec{x}/dt= \vec{f}(\vec{x})$ that is valid for large $\theta$.

**(b)** Set up Matlab code to integrate this system of equations numerically, using Matlab's ''ode45'' function.

**%%(c)%%** Determine the frequency of oscillation of the nonlinear pendulum for 
for the constants $g=9.8$ and $\ell=1.0$ and a variety of oscillation amplitudes.
The oscillation amplitude is given by the angle $\theta_0$ at which the pendulum is released at time $t=0$, with no initial velocity. Determine the frequency $\omega$ for a variety of values in the range
$0 < \theta_0 < 2 \approx 115^{\circ}$ and plot $\omega$ versus $\theta_0$. 

**(d)** For this range of amplitudes, the frequency $\omega$ of the nonlinear pendulum should vary with $\theta_0$ as

<latex>
\omega = \sqrt{g/l} + c \; \theta_0^2
</latex>

for some value of $c$.  Determine what the value of $c$ is from your graph in %%(c)%%. Do this by adding a curve of the above form to your plot from %%(c)%%, and adjusting the constant $c$ until this curve passes through the data points computed in %%(c)%%. What is the value of $c$? What does this mean for the accuracy of a pendulum clock? Would a real pendulum clock run fast or slow compared to the idealized linear clock if the amplitude of oscillation is too big? 

----

====Problem 4: nonlinear damped pendulum==== 

For this lab problem, you are to recreate the time series and phase portrait for the nonlinear damped pendulum, whose equations of motion are

\begin{eqnarray*}
\frac{d}{dt} \left( \begin{array}{c} x_1 \\ x_2 \end{array} \right) = \left( \begin{array}{c} x_2 \\ -g/l \, \sin x_1 - \alpha/m \, x_2 \end{array} \right)
\end{eqnarray*}

Use parameter values $g=9.8, \ell=1$, $m=1$, and $\alpha = 1$. For the time series plots, initiate the pendulum with $x_1 = \theta=0$ and $x_2 = d\theta/dt = 2$. For the phase portraits, show the range $-8 \leq \theta \leq 8$ on the horizontal axis and $-10 \leq d\theta/dt \leq 10$ on the vertical. On top of the quiver plots, show trajectories with initial conditions $\theta=0$ and a variety of $d\theta/dt$ ranging from -10 to 10 in steps of 1. 

----
Turn in your code, your plots, and answer the following questions
====Questions (to be answered at the end of your lab)==== 

**(a)** Describe the differences that you see in the phase portraits of the nonlinear pendulum compared to the linear pendulum.

**(b)** For the undamped nonlinear pendulum, you should see trajectories with three distinctly different types of behavior. Give a physical interpretation of each of those three types of trajectories in terms of the motion of the pendulum.

**%%(c)%%** You should also see three distinct type of trajectories for damped nonlinear pendulum. Give a physical interpretation of each of those types of trajectories.

**(d)** The damped nonlinear pendulum has quite a few more distinct groups of trajectories than the ones shown in shown in your phase portrait. Based on your physical interpretation, what would those be? Revise your nonlinear damped phase portrait to show some of these other trajectories.