Week 3 - Solving second order ODEs

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time. 

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

When you do this you will get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum, as shown in the below figure.

Your animation does not have to look exactly like this and here is the logic.

• For each value in time
• Create an image that represents the position of the pendulum
• You can assume the length of the string and the size of the ball
• Use the marker command to create a sphere
• One end of the pendulum is fixed
• The other end moves according to the solution of the ODE you had solved.
• If you are using using Matlab, use the movie command create an animation
• If you are using Octave you need to generate images first and then stitch them to form the animation (https://skilllync1.freshdesk.com/support/solutions/articles/43000012529-creating-animations-in-octave)

• Upload an animation on You tube and provide the You tube link for an animation in your report.

2. Submit a report that has the following content:

•  Step wise explanation of your code.
•  Errors faced while programming.
•  How have you overcome those errors.
•  Screenshots showing errors thrown in command window.

The presentation of your code will be taken into consideration while grading. Labelling sections of your code using comments is appreciated.

Upload the code using "INSERT SAMPLE CODE" option and then post screenshot of the plots generated.

In case if you are not able to finish the challenge, upload until what you have done and explain the procedure and plot. Deserved marks will be allotted. 

GENERAL INSTRUCTIONS:

• The codes, plots, images, GIFs, and videos should be submitted only through the answer box. 
• Use the Google Drive to share only the cad models or other answers which couldn't be done through the answer box.

PROGRAM:

#Second order ODE
#by Sallavoddin

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

# calling function
def model(theta,t,b,g,l,m):
          theta1=theta[0]
          theta2=theta[1]
          dtheta1_dt=theta2
          dtheta2_dt=-(b/m)*theta2-(g/l)*math.sin(theta1)
          dtheta_dt=[dtheta1_dt,dtheta2_dt]
          return dtheta_dt
#input parameters
b=0.05
g=9.81
l=1
m=1
#intial condition
theta_0=[0,3] #angular displacement abd angular velcity

#time point
t=np.linspace(0,20,200) 

#Solve ODE
theta=odeint(model,theta_0,t,args=(b,g,l,m))

#point results
plt.plot(t,theta[:,0],'g-', label=r'$frac{dtheta_1}{dt}=theta2$')
plt.plot(t,theta[:,1],'b--', label=r'$frac{dtheta2}{dt}=-frac{b}{m}theta_2-frac{g}{l}sintheta_1$')
plt.xlabel('Displacement')
plt.ylabel('Time')
plt.legend(loc='best')
plt.show()
#Animation of Simple pendulum
x0=0 #this is fixed point at the wall
y0=0 #this is fixed point at the wall
ct=1
for ang in theta[:,0]: #angle
    x1=l*math.sin(ang)
    y1=-l*math.cos(ang)
    filename=str(ct)+'.png'
    ct=ct+1
    plt.figure()
    plt.plot([-1,1],[0,0],'red',linewidth=10)
    plt.plot([x0,x1],[y0,y1],'black')
    plt.plot(x1,y1,'o',markersize=20)
    plt.xlim(-1.2,1.2)
    plt.ylim(-1.2,1.2)
    plt.savefig(filename)
    plt.show()