Python Program to solve the 2nd order ODE and to display the Damped oscillations of a simple pendulum

Aim:- A python program to solve this 2nd order ODE as well as to crate the damped oscillations of a simple pendulum.

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.

(d²θ/dt²) +(b/m)(dθ/dt) + (g/l)sinθ =0

It is the equation for the damped simple harmonic motion.

Where

b=damping coefficient of the string material

m=mass of the bob in kg

g=gravity in m/s^2

l=length of the string in mtr

the python program to solve this 2nd order ODE:

\"\"\"
A program to simulate the Damped Simple pedulum
By Rupesh
\"\"\"
import math
import numpy as np 
from scipy.integrate import odeint
import matplotlib.pyplot as plt


#Function that returns dz/dt 
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

# pendulum Inputs
b = 0.05 	# damping coefficient 
g = 9.81	# force due to gravity in m/s^2
l = 1 		# length of the pedullum in mtr
m = 1   	# mass of the bob kg

# Initial conditions
theta_0 =[0,3]  # initial displacement=0 and the angular velocity=3rad/sec

# time points
t = np.linspace(0,20,150)   # time span ranging between 0-20sec 

# calculating the 2nd order ODE with the help of funtion
theta = odeint(model, theta_0, t, args=(b,g,l,m))

# 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(\'Angular_Velocity & Anular_displacement\')
plt.xlabel(\'time\')
plt.legend(loc=\'best\')
plt.show()

# creating the loop for the oscillations of pendulum
ct=1
for i in range(len(theta)): #as the results contains both time and angular velocity
        x0 = 0
        y0 = 0
        x1 = l*np.sin(theta[i,0]) #assigning the direction of oscillation of the bob
        y1 = -l*np.cos(theta[i,0])
    

        fn=str(ct)+\'.png\' # fn is the filename

        plt.figure()
        #creating the regid support
        plt.plot([-1,1],[0,0], lw=4, color=\'k\')
        
        #creating the rope
        plt.plot([x0,x1],[y0,y1], lw=2, color=\'b\')
        
        #creating the shape and size of the bob using marker  command
        plt.plot(x1,y1,\'-o\',markersize=20, color=\'r\')
        plt.grid ()
        plt.axis([-1.5,1.5,-1.5,0.5])  #assaigning the plot area
        plt.savefig(fn) #capturing the each frame for loop
        plt.close()
        ct = ct+1

The out put graph as shown below:

The Damped simple pendulum animation obtained as below:

the link to the youtube uploaded video is https://www.youtube.com/watch?v=MkGuHZk5FwE

the errors obtained during the programming are :

un bounded pedulum oscillaions as shown below

thus after enalizing the mistake as shown below it was rectified

here the \'1\' inside the theta is replace with \'0\' thus the problem is solved.