Solving second-order ODEs

Problem statement

Create animation of a simple pendulum by solving ode. The ode is given along with the necessary constants:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s
and initial conditions as angular displacement=0,angular velocity=3 rad/sec at time t=0.

Objective

To write a python code to solve the given ode and with the result plot the simple pendulum in question.

Solution procedure

code:

#importing the necessary modules

from scipy.integrate import odeint

import matplotlib.pyplot as plt

import numpy as np

import math

#defining the constants

L=1 #m

m=1 #Kg

b=.05 #damping coefficient

g=9.81 #m/s^2

#initail condition

theta_0=[0,3] #angular displacement = 0 and angular veloxity = 3 rad/s

#time points

t=np.linspace(0,20,200) #time from 0 to 20 sec

#defining the function

def func(theta,t,L,m,b,g):

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

#calling odeint function

theta=odeint(func,theta_0,t,args=(L,m,b,g))

#for pendulum animation

th_=theta[:,0]

i=0

R=np.array([[1,1],[1,1]]) #orthogonal transformation

for t in th_:

x=L*math.cos(t)

y=L*math.sin(t)

p=np.array([x,y])

p_rot=np.matmul(R,p)

plt.figure()

plt.plot(np.linspace(-1.5,1.5,100),np.zeros(100),'r')

plt.plot(np.zeros(100),np.linspace(-1.5,0,100),'--')

plt.xlim([-1.5,1.5])

plt.ylim([-1.5,1])

plt.plot(p_rot[0],p_rot[1],marker='o',markersize=12,color='b') #using marker to create the sphere

plt.plot([0,p_rot[0]],[0,p_rot[1]],'b')

plt.savefig('test%.05d.png'%i)

i+=1

plt.close()

1. Import required modules, odeint to solve differential equations, matplotlib.pyplot for plotting, numpy for numerical operations, and math for basic trigonometric functions.

1. Define constants, L = 1 m (length of the pendulum), m = 1 kg (mass), b = 0.05 (damping coefficient), and g = 9.81 m/s² (acceleration due to gravity).

2. Set initial conditions: theta_0 = [0, 3], where 0 is the initial angular displacement (in radians) and 3 is the initial angular velocity (rad/s).

3. Create a time array t from 0 to 20 seconds with 200 equally spaced points using np.linspace.

4. Define a function func(theta, t, L, m, b, g) that represents the system of ODEs, where theta[0] is the angular position and theta[1] is the angular velocity.

5. Inside the function, calculate the derivatives: dtheta1_dt = theta2 and dtheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1), then return them as a list.

6. Use odeint to numerically solve the system defined by func.

7. Extract angular displacement from the result using theta[:,0] and store it in th_.

8. Define a rotation matrix R = [[0,1],[-1,0]] to rotate coordinates by 90 degrees counterclockwise.

9. Loop over each value in th_, which represents the angle of the pendulum at each time step.

10. For each angle, compute the Cartesian position of the pendulum bob using x = L*cos(t) and y = L*sin(t).

11. Form the position vector p = [x, y] and rotate it using matrix multiplication p_rot = R * p to adjust the orientation.

12. Create a new plot: draw a fixed horizontal red line (ground), a dashed vertical reference line, and plot the pendulum bob as a blue circle.

13. Draw a blue line from the origin to the pendulum bob representing the rod.

14. Save each frame as a .png image with filename format test00000.png, test00001.png, etc.

15. Close the plot using plt.close() to prevent memory overload during the loop.
    

16. About 200 image files are saved, I have used MATLAB live script to create the animation, link for the video will be provided in the result. A screenshot of the live script along with the plot is provided as a file.

Errors faced

I haven't face any major error, but since x and y is plotted w.r.t the x axis the orientation of the plot was not correct. So I used that concept that I have learned form linear algebra called orthogonal transformation to rotate the pendulum by 90⁰

Result

Animation is created successfully and the link is provided below

https://youtu.be/oPI2Krvf8-c

(code is also provided as .py file)