Solving second order ODEs using PYTHON

Aim:

Write a program in PYTHON that will simulate the pendulum motion.

Theory:

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 Newton's 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 these systems.

This ODE represents the equation of motion of a simple pendulum with damping.

`(d ^2 θ )/(dt ^2 ) + b/m (dθ )/dt + g/l sinθ = 0`

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.

From the above equation, we create two ODEs’ i.e.,

`(dθ )/dt = θ2,`

`(dθ2 )/dt = -(b/m (θ2) + g/l sin⁡(θ1))`

Now we write a function code in which we are going to use some parameters that are going to use in the main program. We basically call this function in our program to solve our problem.

And, after the function code, we run the function code and call its value in code. And in code we will first write all the inputs,

• L=1 meter,
• m=1 kg,
• b=0.05.
• g=9.81 m/s2.

Explanation of the main for solving ODE and simulation of the pendulum: -

1. Now we are solving the ODE as mentioned above and we split that ode into two single order ODE’s one for velocity and one for displacement.
2. So, when we need to solve the first ODE, we need to provide initial conditions for the displacement and velocities.
3. Our ODE depends upon the time, theta and other parameters i.e.,

• Damping coefficient (b)
• Acceleration due to gravity (g)
• Length (l)
• Mass of bob (m)

5. For simulation, we created an image for each value in the time that represents the position of the pendulum.
6. For making bob I use the marker command. And keeping the one end of the pendulum fixed so that we can oscillate our pendulum.

PYTHON Code for Solving Second-order ODE

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


def model(theta,t,b,g,m,l):

	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


b= 0.05
g = 9.81
l = 1
m = 0.1


theta_0 = [0,3]
t= np.linspace(0,20,100)


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


#Plot Generation
plt.plot(t,theta[:,0],'g' ,label=r'Displacement')               # g -> green 
plt.plot(t,theta[:,1],'b' ,label=r'Velocity')                   # b -> blue
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.legend()
plt.savefig('output.png')
plt.show()


#Creating Animation
ct = 1
for theta in theta[:,0]:
	x0 = 0
	y0 = 0                                                        # x0 y0 fixed end of pendulum
	x1 =  l*math.sin(theta)
	y1 = -l*math.cos(theta)                                       # using sin and cos function to locate the cordinates of mass
	plt.figure()
	plt.plot([-1,1],[0,0],"black",linewidth=5)                    # horizontal support
	plt.plot([x0,x1],[y0,y1],"red")                               # pendulum string by connecting 2 points
	plt.plot(x1,y1,'o',markersize=8)                              # x1 , y1 = center of mass / pendulum representation  !
	plt.xlim([-1,1])
	plt.ylim([-1.5,1])
	filename = str(ct) + 'Pendulum.png'
	plt.savefig(filename)
	ct = ct +1

OUTPUTS:

PLOT:

VIDEO:

Click here to Watch Simulation Video