Homework 1: Swinging Pendulum

[Notebook]

Due: 1/19/2017

Jupyter Notebooks

As discussed, Jupyter Notebooks provide a powerful tool for us to incorporate code, data, analysis, notes, and documentation in one single place that can be reviewed by anyone (or ourselves) to understand our experiment. In terms of performing the coding and tests together, this is known as Literate Programming as dubbed by Dr. Donald Knuth. It is the idea that your code and results should be readable by humans rather than computers. Let the computer do the job of interpretation. We will extend this idea to include our data, results, and analysis for our homeworks.

For each homework, the assignment, code, results, and your writeup will all be together. When you are finished, you can do File -> Download as -> HTML (.html) and turn this in via [Laulima].

Background

A swinging pendulum is merely a system of forces that forms a second-order, ordinary differential equation when using the small-angle approximation. We neglected the effects of friction and forcing such that we have a basic pendulum purely under the force of gravity. Splitting the second-order equations into two, first-order equations, we derived the equations of motion for a swinging pendulum as:

\[\begin{split}\frac{d\theta}{dt} = \omega,\\ \frac{d\omega}{dt} = -\frac{g}{l}\sin\left(\theta\right).\end{split}\]

where \(\theta\) is the angle the pendulum has rotated in radians (its position), \(\omega\) is the rotation rate in rad sec (its angular velocity), \(g\) is the gravitational acceleration in m sec, and \(l\) is the length of the pendulum arm in m. The closed-form solution to this set of equations is given by,

\[\theta(t)=\theta(0)\cos\left(\sqrt{\frac{g}{l}}t\right).\]

The natural period of the system is given by:

\[2\pi\sqrt{l/g}\]

Employing the forward Euler, we were able to derive numerical approximations for these equations as given by:

\[\begin{split}\theta_{n+1} = \theta_n + \Delta t\left(\omega_n\right), \\ \omega_{n+1} = \omega_n + \Delta t\left(-\frac{g}{l}\sin\theta_n\right).\end{split}\]

These equations are trivial to integrate forward with a computer and this notebook does just this.

Recall that the explicit, forward Euler is very unstable for periodic functions (further recall that periodic functions are defined by \(e^{-i\omega t}\)) as numerical noise will grow exponentially. This requires an extremely small timestep. The way to avoid this problem is to move to a semi-implicit algorithm, the Euler-Cromer, given by:

\[\begin{split}\omega_{n+1} = \omega_n + \Delta t\left(-\frac{g}{l}\sin\theta_n\right),\\ \theta_{n+1} = \theta_n + \Delta t\left(\omega_{n+1}\right).\end{split}\]

Exercise

In this exercise, you will use the this Jupyter Notebook to experiment with various model parameters and settings. First, download the notebook to your computer and read over both the code and assignment (this document). You should understand the code and how it works. Notice that it is organized into distinct sections: preamble, parameter declaration, and model code. I strongly recommend that you follow this structure and format for your work this semester.

To analyze, interpret, and experiment with the model, we will use this Jupyter Notebook.

NOTE: If you do not wish to use the notebook, you could save the code into a file hw1.py. This would allow you to run the model within IPython or SPYder and plot the results manually.

Model Code

Here we define the code for our model in the homework exercise. This is primarily composed of a single function, which we will call pendulum. We must:

  1. Load any modules required
  2. Define any fixed Parameters that would not change.
  3. Define a list of arguments that can be specified by the user. We must set default values for these parameters.
  4. Define a list of initial conditions for our model. We must also specify these in the calling arguments so that the user can change them. We must set default values for these parameters.
  5. Write documentation for our function that describes what it does and what the arguments are.
  6. Create the structures that are required for our model to run
  7. Set the initial conditions
  8. Set the boundary conditions (if any)
  9. Integrate our model

Modules

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np

Parameters

  1. The acceleration of gravity, g [m s\(^{-2}\)]
  2. The start time of the simulation, tstart [s]
  3. The end time of the simulation, tend [s]
g = 9.81
tstart = 0
tend = 10

Model

Implement both explicit and semi-implicit integration of a simple pendulum.

Input

l      : length of moment arm in m [0.5]
dt     : time-step of integration in s [0.05]
stable : [True/False] use explicit [False] or semi-implicit (True)
omega0 : Initial angle of pendulum in rad [0]
theta0 : Initial angular velocity of pendulum in rad/s [pi/6]

Output

times : array of times in s
theta : array of angles in rad
omega : array of angular velocities in rad/s
truth : array of true angles from closed-form solution in raad
def pendulum(l=0.5, dt=0.05, stable=False, omega0=0, theta0=np.pi/6.0):

    # Create the variables to store our results
    times=np.arange(tstart, tend+dt, dt)  # Create the array of model times
    Nt=len(times)                         # Number of time-steps
    omega=np.zeros(Nt)                    # Create the array for our velocity with time
    theta=np.zeros(Nt)                    # Create the array for our angle with time
    truth=np.zeros(Nt)                    # Create the array for the true solution with time

    # Set our initial conditions
    omega[0]=omega0
    theta[0]=theta0
    truth[0]=theta0

    # Loop over all of our times and integrate the model
    for n in range(Nt-1):
        # Time-step the velocity
        omega[n+1]=omega[n] + dt*(-g/l*np.sin(theta[n]))

        # Time-step the angle. If we are using 'stable', then use the semi-implicit time-step;
        # otherwise, use the standard Euler time-step.
        if stable:
            theta[n+1]=theta[n] + dt*omega[n+1]
        else:
            theta[n+1]=theta[n] + dt*omega[n]

        # Since we know the closed-form solution, we can compute the true pendulum solution
        # for comparison against our numerical solution.
        truth[n+1]=truth[0]*np.cos(np.sqrt(g/l)*times[n+1])

    return times, theta, omega, truth

Exercise 1

Using the default parameters, run the model and make three plots:

  1. Time-series of both \(\theta\) and the true \(\theta\) versus time
  2. Time-series of \(\omega\) versus time
  3. Phase-space plot of \(\theta\) versus \(\omega\)
times, theta, omega, truth = pendulum()

Exercise 1a)

Write about the results. Do they make sense? What could be wrong?

Results

Write up your analysis of the results above here.

Exercise 1b)

Repeat using stable=True and compare with the default results.

Results

Write up your analysis of the results above here.

Exercise 2

Explore the stability of the two integrators. First, use stable=True in the parameters and execute. Calculate the root-mean-square (RMS) difference between the modeled and the true solution. The RMS is a measure of how the model deviates from the truth. It is defined by:

\[RMS=\sqrt{\frac{1}{N}\sum_i^N\left(x_i-y_i\right)^2},\]

where \(N\) are the number of values, \(x\) is dataset one, and \(y\) is dataset two. It is what the name says, it is the square-root of the mean of the squared differences.

You can compute it with our arrays using:

np.sqrt( np.mean((theta-tru)**2) )

And display the value by:

print( "RMS=", np.sqrt( np.mean(( theta-tru )**2 )))

This will display the result of the calculation.

Exercise 2a)

Run the model with the parameters of: \(l=0.5\), \(dt\) as 20 steps of the natural frequency, Stable=True, initial angle \(\theta_0=\pi/6\), and initial velocity \(\omega=0\). Print the RMS value. You should get a value of \(\approx 0.0835\).

Exercise 2b

Next, use the same parameters as above, but revert to the unstable integrator by setting stable=False. Integrate and calculate the RMS. This RMS will be three orders of magnitude larger than the stable version.

Experiment with the timestep. Can you get the RMS of the unstable solution to within a factor of 10 of the RMS from the stable solution? What time-step was required? How long did the program take to run relative to the stable?

Results

Write up your analysis of the results above here.

Problem 3

Using the stable solution with the default parameters, experiment with the initial angle from which the pendulum is released. Recall that in our derivation, we made the small-angle approximation. What happens if you release it from large angles? Try a number of large angles and compare. Compare plots of each one along with the RMS difference.

Results

Write up your analysis of the results above here.

Problem 4

Using the stable solution with the default parameters, experiment with the length of the pendulum arm. Inspect the resulting dt value for each case. A system that requires very small \(\Delta t\) values is considered to be numerically sensitive to \(\Delta t\), and it can go unstable easily.

Plot your choices of \(l\) and the resulting \(\Delta t\) values. Recall that the stability is related to \(\Delta t\), so what do the results say about the stability?

Create a figure that plots the natural frequency \((2\pi\sqrt{l/g})\) for a range of \(l\) values.

Results

Write up your analysis of the results above here.

Problem 5 [OCN 681 ONLY]

Try adding an initial velocity (by default it is set to 0). What happens to the period and amplitude of the solution?

Results

Write up your analysis of the results above here.