Solving second order ODEs Using Python

Solving Secind Order ODE's in Python: Simple Pendulum

1. AIM: To Solve a second order ordinary differnential eqaution in Python.

2. Goverining Eqautions:

`(d^2theta)/(dt^2)+b/m (d theta)/dt +g/Lsintheta = 0`

Varables:

`theta` -Angle made by rod with vertical

B - damping coefficient

m and L are mass of the bob and length of the rod/string

3. Objective of the Project:

- To solve the seocnd order ODE in python with graphs for velocity and displacement

- To simulate the motion of simple pendulum (animation)

4. Procedure

- Simplify the ODE using substituition method. Prepare a function for the same code

- Specify the geometrical and physical parametrs like pendulum length, mass, damping coefficient and gravity.

- Set initial conditions and timespan for which you want to run the simulation.

- Use the inbuilt ODE solver "odeint" to solve the ODE which is reduced by calling our function, to obtain teh values of displacement and velocity and store in array.

- Plot the graph of displcement v/s time and velocity v/s time.

- Create Simulation of the positions of the pendulum bob and arm using sticing of images by saving the images.

5. Python Code

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

#Function that returns dz/dt
def model (theta,t,b,g,l,m):

    theta1 = theta[0] # Refer governing equation section
    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 #Damping Coefficient
g = 9.81
l = 1 #length of figure
m = 1 #mass of the bob

#initial conditions
theta_0 = [0,3] #initial condition where angle displacement = 0 and velocity =3m/s

#line plots
t = np.linspace(0,20,150) #Running Time for 2 seconds and t array to have 150 equally spaced between 0 - 20

#Solving ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m)) # equation solved by odeint function.

#plot results
plt.figure() #Plot as a figure

#PLotting time and angle displcements
plt.plot(t,theta[:,0],'-b',label = r'$frac{dtheta_1}{dt}=theta2$')

#plotting time and velocity
plt.plot(t,theta[:1],'r--',label = r'$frac{dtheta_2}{dt}=frac{b}{m}theta_2-frac{g}{l}sintheta_1$')
plt.ylabel('Plot')
plt.xlabel('Time')
plt.legend(loc='best')
plt.show()

val = theta[:,0] #Assigning all rows and colums to val
ct = 1 #Counter for file names
for i in range(0,len(val)):
    plt.figure() #open new figure
    x1 = 0.5 - l*np.sin(val[i])
    y1 = 0.5 - l*np.cos(val[i])
    plt.plot([0.5,x1],[0.5,y1])
    plt.plot(x1,y1,'o',markersize = 10)
    plt.plot([-1,2],[0.5,0.5],linewidth = 4)
    plt.xlim([-1,2])
    plt.ylim([-1,2])
    filename = str(ct) + '.png'
    plt.savefig(filename)
    ct = ct + 1

Explanation

1. Set the parameters of teh pendulum like mass,length, gravity and damping coefficient.

2. Set the initial conditions and timespan for which you want ti run the simulation

3. Use ode solve for solving the equation

4. Results are stored in an array of theta

5. Since results contain 2 Components, angle displcement and velocity, to plot both the graphs, use the (:,0) adn (:,1) property of python and plot these graphs.

6. To Create animation, use a for loop so that we can get X any Y co-ordinaltes using trignomotery.

7. the ct counter to make up of for the file name using str(ct) and a'.png' file extension.

8. use image sticher to stich to stich all the images to create the animation.

6. Errors

- Changing the val=[:,2] to val=theta[:,0]. Indexing starts from 0.

- plt.figure() to be added inside the for loop than outside the for loop.

7. Graphical Results

- Displacement (blue Color) and Velocity (Red) is plotted against time. We can see the velocity and angle of displacement is reducing as the time increases.

8. Conclusion

- The reduction in velocity and angle diaplcement is due to damping coefficient from the second order ODE of a simple pendulum.

- Animation of the pendulum.