Animating Behaviour of Pendulum

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
import math

\"\"\"
Challenge - WEEK 3 : 2nd Order ODE
\"\"\"

# 2nd oder ODE to system of first order ODE
def model(theta_init, t, b, g, l, m ):
    theta1     = theta_init[0]
    theta2     = theta_init[1]
    dtheta1_dt = theta2
    dtheta2_dt = -(b/m)*theta2 -(g/l)*math.sin(theta1)
    dtheta_dt  = [dtheta1_dt, dtheta2_dt]

    return dtheta_dt


# given data
length     = 1  # m
mass       = 1  # kg
gravity    = 9.81  # m/s2
damp_coeff = 0.005

# time to simulate
n = 300 # number of time steps to take
t = np.linspace(0,20,n)

# initial conditions
initial = [0, 3] #angular displacement, angular velocity(rad/sec)

# solve ode
theta = odeint(model, initial, t, args=(damp_coeff,gravity,length,mass))

ct = 0

#here theta is the angle -  from this we can get the x and y coordinates
#plot that point and animate
#x = lsin(theta) y = -lcos(theta)

for i in range(n):
    filename = \'loc%04d.png\' % ct
    x = length*math.sin(theta[i][0])
    y = -1*length*math.cos(theta[i][0])
    #print(\"%0.2f and %0.2f\" %(x,y))

    plt.figure()
    plt.plot([0, x], [0, y], linewidth = 2)
    plt.plot(x, y, marker = \"o\", ms = 15)
    plt.xlabel(\'x axis\')
    plt.ylabel(\'y axis\')
    plt.title(\'Simple Pendulum Behaviour\')
    plt.xlim([-1.5, 1.5])
    plt.ylim([-1.5, 1.5])
    plt.savefig(filename)
    ct += 1
    #plt.show()


\"\"\"
#plot results
plt.plot(t,theta[:,0],\'b-\',label=r\'$\\frac{d\\theta_1}{dt}=\\theta2$\')
plt.plot(t,theta[:,1],\'r--\', label = r\'$\\frac{d\\theta2}{dt}=-\\frac{b}{m}\\theta_2-\\frac{g}{l}sin\\theta_1$\')
plt.ylabel(\'plot-position and velocity\')
plt.xlabel(\'time\')
plt.legend(loc=\'best\')
plt.show()
\"\"\"

In this code, we initially convert the 2nd order ODE to a system of 1st order ODE and put that into a function (model). This is used later by odeint to integrate and get the system of solution to the system of 1st order ODEs given.

The given data, time divided into number of points at which we wish to solve the ODE, initial conditions of position and velocity are defined.

Using all of the above data and odeint from scipy library we can solve the ODEs to integrate and get the output of [position, velocity] at predefined time steps between 0-20s. We get output of position and velocity as our model function has ODEs defined for velocity and acceleration.

Further to plot the motion of the pendulum, we have the data for position. Since this is the angular displacement, we can convert it to cartesian coordinates using x=Lsin(theta) and y=-Lcos(theta). Each of these points are plotted for each time step. And each figure is saved using plt.savefig() as previously done in Week 1 challenge. Using these .png files, it can be stitched together to animate the behaviour of the pendulum.

Errors faced - solved

Warning regarding more than 20 images being saved

Initially started looping with theta[:][0] data and plotted with entire array of time which gave an error. Corrected it by looping with number of timesteps(n) and plotting correct position and time.

Only one file loc000.png saved. filename was defined outside for loop, ct was being incremented inside the for loop. Probably, it was not updating the filename as it was not in the for loop. Did not give an error in the command window though. Corrected it by defining filename within the loop

However, plot of position v/s time was not required. So corrected it by finding x and y with the theta data and plotting x,y instead

Plotted point was not easily understood, introduced xlim and ylim which gave a clear picture of the simple pendulum motion

Animation Link: https://youtu.be/91LTWFLZ2VU