%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
We start by reviewing how to solve the simple pendulum. In equations, the problem is
\begin{align} \frac{d\omega}{dt} &= - \frac{g}{L} \theta\;, \\ \frac{d\theta}{dt} &= \omega\;, \end{align}Let's solve this first with Euler's method.
# Euler calculation of motion of simple pendulum
# by Kevin Berwick,
# based on 'Computational Physics' book by N Giordano and H Nakanishi,
# section 3.1
#pendulum length in metres
length = 1
#acceleration due to gravity
g = 9.8
#Discretize time into 250 intervals
N = 250
#Time step in seconds
dt = 0.04
#initializes omega, a vector of dimension N,to being all zeros
omega = np.zeros(N)
#initializes theta, a vector of dimensionN,to being all zeros
theta = np.zeros(N)
#this initializes the vector time to being all zeros
t = np.zeros(N);
#you need to have some initial displacement, otherwise the pendulum will not swing
theta[0]=0.2
#loop over the timesteps
for i in range(N-1):
omega[i+1] = omega[i] - (g/length)*theta[i]*dt
theta[i+1] = theta[i] + omega[i]*dt
t[i+1] = t[i] + dt
fig, axes = plt.subplots()
#plots the numerical solution in red
plt.plot(t, theta, 'r-');
# label axes
plt.xlabel('time (seconds) ');
plt.ylabel('theta (radians)');
plt.title('Simple Pendulum - Euler Method')
plt.show()
This problem with growing oscillations is addressed by performing the solution using the Euler - Cromer method. The code is below
# Euler_cromer calculation of motion of simple pendulum
# by Kevin Berwick,
# based on 'Computational Physics' book by N Giordano and H Nakanishi,
# section 3.1
#pendulum length in metres
length = 1
#acceleration due to gravity
g = 9.8
#Discretize time into 250 intervals
N = 250
#Time step in seconds
dt = 0.04
#initializes omega, a vector of dimension N,to being all zeros
omega = np.zeros(N)
#initializes theta, a vector of dimensionN,to being all zeros
theta = np.zeros(N)
#this initializes the vector time to being all zeros
t = np.zeros(N);
#you need to have some initial displacement, otherwise the pendulum will not swing
theta[0]=0.2
#loop over the timesteps
for i in range(N-1):
omega[i+1] = omega[i] - (g/length)*theta[i]*dt
theta[i+1] = theta[i] + omega[i+1]*dt #This line is the only change between this program and standard Euler
t[i+1] = t[i] + dt
fig, axes = plt.subplots()
#plots the numerical solution in red
plt.plot(t, theta, 'r-');
# label axes
plt.xlabel('time (seconds) ');
plt.ylabel('theta (radians)');
# axes scales
plt.xlim(0, 10);
plt.ylim(-0.2, 0.2);
# Figure 3. Simple Pendulum: Euler - Cromer method
plt.show()
In this example I use a variety of approaches in order to solve the following, very simple, equation of motion. This is the general equation for all simple harmonic oscillators which all have the same mathematical structure.
$$ \frac{\mathrm{d}^2y}{\mathrm{d}t^2} = -y $$
Here, four approaches are used to solving the equation, illustrating the use of the Euler, Euler Cromer, Second order Runge-Kutta and finally the python ODE solver integrate.odeint from SciPy. We will need to import this function.
The first step is to take the second order ODE equation and split it into 2 first order ODE equations that can be solved using the mentioned methods. These are
$$ \frac{\mathrm{d}y}{\mathrm{d}t} = v $$and
$$ \frac{\mathrm{d}v}{\mathrm{d}t} = -y $$These can then be solved using the preferred method, by dividing them into individual python functions that return the values of interest, in this case the time and y. Here are the functions to do the individual calculations and the needed code to plot them.
# Simple harmonic motion - Euler, Euler Cromer and Runge Kutta methods
# by Kevin Berwick,
# based on 'Computational Physics' book by N Giordano and H Nakanishi,
# section 3.1
# Equation is d²y/dt² = -y
def Euler(initial_displacement):
N = 2500 #Discretize time into 2500 intervals
dt = 0.04 #time step in seconds
#Initialize v,y and time (t) to all zeros
v = np.zeros(N)
y = np.zeros(N)
t = np.zeros(N)
y[0] = initial_displacement #Add initial displacement
#Euler Solution
for i in range(N-1):
v[i+1] = v[i] - y[i]*dt
y[i+1] = y[i] + v[i]*dt
t[i+1] = t[i] + dt
return t,y
def Euler_Cromer(initial_displacement):
N = 2500 #Discretize time into 2500 intervals
dt = 0.04 #time step in seconds
#Initialize v,y and time (t) to all zeros
v = np.zeros(N)
y = np.zeros(N)
t = np.zeros(N)
y[0] = initial_displacement #Add initial displacement
#Euler Cromer Solution
for i in range(N-1):
v[i+1] = v[i] - y[i]*dt
y[i+1] = y[i] + v[i+1]*dt
t[i+1] = t[i] + dt
return t,y
def Runge_Kutta(initial_displacement):
N = 2500 #Discretize time into 2500 intervals
dt = 0.04 #time step in seconds
#Initialize v,y along with v',y' and time (t) to all zeros
v = np.zeros(N)
vp = 0.0
y = np.zeros(N)
yp = 0.0
t = np.zeros(N)
y[0] = initial_displacement #Add initial displacement
#Runger Kutta Solution
for i in range(N-1):
vp = v[i]-0.5*y[i]*dt
yp = y[i]+0.5*v[i]*dt
v[i+1] = v[i] - yp*dt
y[i+1] = y[i] + vp*dt
t[i+1] = t[i] + dt
return t,y
#With the functions written, they need to be called to obtain the results
#Set the initial displacement
initial_displacement = 10
#Obtain the results for the three first methods
[t_euler,y_euler] = Euler(initial_displacement)
[t_euler_cromer,y_euler_cromer] = Euler_Cromer(initial_displacement)
[t_runge_kutta,y_runge_kutta] = Runge_Kutta(initial_displacement)
#Python Integration method
#To use the the python method, we need to import the SciPy integrate
# module
from scipy import integrate
#A function is created that has information on the ODEs
def Python_ODEint(variables, t):
y, v = variables #Seperate the two variables
dydt = [v, -y] #Create an array that holds the factors of the ODEs
return dydt
#For the python method, first create an array for the time variable
t_integrate_odeint = np.arange(0, 100, 0.04)
#Call the integrate.odeint method with the function that was written,
#giving initial conditions and the time array
sol = integrate.odeint(Python_ODEint, [initial_displacement, 0],
t_integrate_odeint)
#Save the arrays that can be used to plot
[t_integrate_odeint, y_integrate_odeint] = [t_integrate_odeint, sol[:,0]]
#Plot the four figures next to each other in a 2x2 grid.
fig, plt_array = plt.subplots(2, 2, figsize=(9, 9))
#Plot the Euler values
plt_array[0, 0].plot(t_euler, y_euler,'r-')
plt_array[0, 0].set_xlabel('Time')
plt_array[0, 0].set_ylabel('Displacement')
plt_array[0, 0].set_ylim(-100, 100) #Set the y range to be higher for Euler
plt_array[0, 0].set_title('Euler')
#Plot the Euler Cromer values
plt_array[0, 1].plot(t_euler_cromer, y_euler_cromer,'b-')
plt_array[0, 1].set_xlabel('Time')
plt_array[0, 1].set_ylabel('Displacement')
plt_array[0, 1].set_ylim(-15, 15)
plt_array[0, 1].set_title('Euler Cromer')
#Plot the Runge Kutta values
plt_array[1, 0].plot(t_runge_kutta, y_runge_kutta,'g-')
plt_array[1, 0].set_xlabel('Time')
plt_array[1, 0].set_ylabel('Displacement')
plt_array[1, 0].set_ylim(-15, 15)
plt_array[1, 0].set_title('Runge Kutta')
#Plot the integrate.odeint values
plt_array[1, 1].plot(t_integrate_odeint,y_integrate_odeint,'k-')
plt_array[1, 1].set_xlabel('Time')
plt_array[1, 1].set_ylabel('Displacement')
plt_array[1, 1].set_ylim(-15, 15)
plt_array[1, 1].set_title('integrate.odeint')
fig.suptitle('Simple pendulum solution using Euler, Euler Cromer, Runge Kutta and Python ODE solver.')
plt.show()
fig
If we add some damping to our model, we can begin to include the effect of friction into the pendulum.
The friction is assumed to be linearly proportional to the velocity:
$$ F_{fric} \propto - \frac{d \theta}{dt} $$We can model this with a damping coefficient, $q$, that represents the amount of friction in the system. The equation of motion then becomes:
\begin{align} \frac{d^2\theta}{dt^2} &= - \frac{g}{L} \theta -q \frac{d \theta}{dt} \;, \\ \end{align}or written as a pair of first-order equations
\begin{align} \frac{d\omega}{dt} &= - \frac{g}{L} \theta -q \frac{d \theta}{dt} \;, \\ \frac{d\theta}{dt} &= \omega\;, \end{align}The analytical/theoretical solution to such differential equations are typically studied in a course on Ordinary Differential Equations.
\begin{align} \frac{d^2\theta}{dt^2} &= - \frac{g}{L} \theta -q \frac{d \theta}{dt} \;, \\ \end{align}Formally, we call this a linear, second-order differential equation with constant coefficients.
Without assuming knowledge about how to solve such an equation with analytical methods, let's try a numerical approach.
This solution uses q=1
# Euler_cromer calculation of motion of simple pendulum with damping
# by Kevin Berwick,
# based on 'Computational Physics' book by N Giordano and H Nakanishi,
# section 3.2
#pendulum length in metres
length = 1
#acceleration due to gravity
g = 9.8
#damping strength
q = 1
#Discretize time into 250 intervals
N = 250
#Time step in seconds
dt = 0.04
#initializes omega, a vector of dimension N,to being all zeros
omega = np.zeros(N)
#initializes theta, a vector of dimensionN,to being all zeros
theta = np.zeros(N)
#this initializes the vector time to being all zeros
t = np.zeros(N);
#you need to have some initial displacement, otherwise the pendulum will not swing
theta[0]=0.2
#loop over the timesteps
for i in range(N-1):
omega[i+1] = omega[i] - (g/length)*theta[i]*dt - q*omega[i]*dt
theta[i+1] = theta[i] + omega[i+1]*dt #This line is the only change between this program and standard Euler
t[i+1] = t[i] + dt
fig, axes = plt.subplots()
#plots the numerical solution in red
plt.plot(t, theta, 'r-')
# label axes
plt.xlabel('time (seconds) ')
plt.ylabel('theta (radians)')
# axes scales
plt.xlim(0, 10)
plt.ylim(-0.2, 0.2)
plt.title('The damped pendulum using the Euler-Cromer method')
plt.show()
- Write code using a loop to calculate and plot the motion of damped pendulum for several different values of
q- Damped pendulums exhibit three regimes that depend on
q. What do you think this means?
- underdamped
- critically damped
- overdamped
- Try to estimate the value
qwhen the system transitions from underdamped to overdamped; this is what critically damped means.
# in class exercise
In the damped pendulum, the motion decay from an initial perturbation until it comes to rest.
To keep the pendulum swinging, we need to provide a source of energy to drive the motion.
One way of modeling such a system is to drive it with a driving amplitude $F_D$ and driving angular frequency $\Omega_D$
\begin{align} \frac{d^2\theta}{dt^2} &= - \frac{g}{L} \theta -q \frac{d \theta}{dt} + F_D \sin(\Omega_D t)\;, \\ \end{align}This driving force will pump energy into (or out of) the system. The external forcing frequency, $\Omega_D$, "competes" with the natural frequency of the pendulum, $\Omega = \sqrt{g/L}$.
- Modify the code below so that it includes a driving force term in the equation of motion.
# Euler Cromer Solution for damped, driven pendulum
#pendulum length in metres
length = 9.8
#acceleration due to gravity
g = 9.8
#damping strength
q = 0.5
#The following is created as a function to be able to easily change values.
def Euler_Cromer_Pendulum(N, dt, F_D, Omega_D):
#initialize array storage
omega = np.zeros(N)
theta = np.zeros(N)
t = np.zeros(N)
#need to have some initial displacement, otherwise the pendulum will not swing
theta[0]=0.2
omega[0]=0
#loop over the timesteps.
for i in range(N-1):
omega[i+1] = omega[i] + (-(g/length)*theta[i] - q*omega[i])*dt ## CHANGE THIS LINE
theta[i+1] = theta[i] + omega[i+1]*dt
t[i+1] = t[i] + dt
return t, omega, theta
F_D = 1.2 # Driving force ampltiude
Omega_D = 0.6 # Driving force frequency
N = 15000 # Discretize time
dt = 0.04 # Time step in seconds
t, omega, theta = Euler_Cromer_Pendulum(N, dt, F_D, Omega_D)
fig, axes = plt.subplots()
plt.plot(t, theta, 'b-');
# label axes
plt.xlabel('time (seconds) ');
plt.ylabel('theta (radians)');
# axes scales
plt.xlim(0, 60);
plt.ylim(-2.0, 2.0);
plt.title('Results from Physical pendulum,\nusing the Euler-Cromer method\n' +
'$F_D = {:.1f}, \Omega_D = {:.1f}$'.format(F_D, Omega_D))
plt.show()
- Explore the damped, driven pendulum for different forcing amplitudes and frequencies.