Week 3 - Solving second order ODEs

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

Theory:

Pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. When a pendulum is displaced sideways
from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the
restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The mathematics of pendulums is in
general quite complicated. Simplifying assumptions can be made, which in the case of a simple pendulum allow the equations of motion to be solved analytically
for small-angle oscillations.
Damped Simple Harmonic Motion
When the motion of an oscillator reduces due to an external force, the oscillator and its motion are damped. These periodic motions of gradually decreasing
amplitude are damped simple harmonic motion. An example of a damped simple harmonic motion is a simple pendulum.
In the damped simple harmonic motion, the energy of the oscillator dissipates continuously. But for a small damping, the oscillations remain approximately
periodic. The forces which dissipate the energy are generally frictional forces.

code:

import math
import numpy as n1
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def model(initial,t,b,l,g,m):
    t1=initial[0]
    t2=initial[1]
    dt1=t2
    dt2=(-b*t2/m)-(g*math.sin(t1)/l)
    dt=[dt1,dt2]
    return dt
b=0.05
l=1
g=9.81
m=1

initial=[0,3]
t=n1.linspace(0,20,200)
disvelo=odeint(model,initial,t,args=(b,l,g,m))

plt.plot(t,disvelo[:,0],'-',label='displacement')
plt.plot(t,disvelo[:,1],'--',label='velocity')
plt.legend()
plt.xlabel('time')
plt.ylabel('displacement')
plt.show()

y=disvelo[:,0]
ct=1
for angle in y:
    x1=0;y1=0
    x2=l*math.sin(angle)
    y2=-l*math.cos(angle)
    file='test%05d.png'%ct
    plt.figure()
    plt.plot([-1.5,1.5],[0,0],linewidth=8)
    plt.plot([x1,x2],[y1,y2],linewidth=3)
    plt.plot(x2,y2,'o',markersize=15)
    plt.xlim([-1.5,1.5])
    plt.ylim([-1.5,1.5])
    plt.savefig(file)
    ct=ct+1

codeexplanation:

The equation to be solved is a 2nd order ODE. For solving second order ODEs, we first have to express it as a set of 1st order ODEs and then solve for theta.
ii. This is done by assuming θ = θ1 and . Hence, the second order ODE becomes :
iii. Thus we can express the ODE in a matrix form as : .
iv. A function called dθ_dt is defined which would take in the values of θ,t,b,g,l and m as arguments and return a matrix, dtheta_dt 
v. After getting this matrix, we just need to integrate this matrix between the time limits to get the solution. To get the time limits, we define time points such
that we have, for example, 200-time points between 0 sec and 10 sec. This is done using the linspace command from the numpy module.
vi. To integrate, we can use the odeint fuction after importing it from the scipy.itegrate module. This odeint function is a solver which will integrate the matrix
returned by the function dθ_dt and give us the solution matrix: , where θ1 is the angular displacement and θ2 is the angular velocity.
vii. To obtain the plots, the displacement and velocity arrays are plotted against the time points on x-axis.
viii. To animate the motion of the pendulum, 1st we have to plot a pivot point at origin. Then to obtain the position of the bob at each time point, a for loop is
used with t increasing from lower limit to upper limit. To form the thread, a line is plotted to join the bob position with the pivot point at each time point.
The plot obtained for every time point is saved in a folder and stitched together in the correct sequence to generate the animations a gif file

errorsfaced:

i forget to use dot operator in linspacespace linecommand after n so i faced error

2.i forget to use mathmodle befor sine and cos functions

fianal output:i. The output shows angular displacement and angular velocity plotted against time.
ii. The animation of the motion of the pendulum.

youtubelink:

https://youtu.be/bsoVpIBsS3A