Week 3 - Solving second order ODEs

Aim:

To Write a program in Python that will simulate the pendulum motion. ODE is used to describe the transient behaviour of a system.

Objective:

To simulate the transient behaviour of a simple pendulum and to create an animation of its motion.

Problem statement:

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 the Newtons 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 theses systems.

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

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.

Write a program in Python that will simulate the pendulum motion, just like the one shown in the start of this challenge.

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

When you do this you will get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum, as shown in the below figure.

Program code:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
# function that return d2/dt
def model1(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
#inputs
b=0.05
m=1
g=9.81
l=1
#initial condition
theta_0=[0,3]
#time point
t=np.linspace(0,20,200)
#solve ODE
theta=odeint(model1,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('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.savefig('ODE Solution')
plt.show()

# for simple pendulum animation
x0 = 0
y0 = 0
ct = 0
theta_1 = theta[:,0]
for time in theta_1:
     x1 = l*math.sin(time)
     y1 = -l*math.cos(time)
     filename = 'test%05d.png' % ct
     ct = ct + 1
     plt.figure()
     plt.plot([-0.5,0.5],[0,0,],linewidth = 5)
     plt.plot([x0,x1],[y0,y1],linewidth = 2)
     plt.plot(x1,y1,'o',markersize = 20)
     plt.xlim(-1.5,1.5)
     plt.ylim(-1.5,0.5)
     plt.savefig(filename)

Solution procedure:

• Importing the libraries such as math, matplotlib.pyplot, numpy & scipy.
• dtheta_dt is the variable which contains dtheta1_dt & dtheta2_dt

where

dtheta1_dt = dθ1/dt=θ2 & dtheta2_dt = −b/mθ2−g/Lsinθ1
dtheta1_dt represent the 1st ode while dtheta2_dt represent the 2nd ode.
After that both the ode's have been put under dtheta_dt.

• The constants has been defined b=damping coefficient, g = gravity in m/s^2, l = length of the pendulum in m, m = mass of the ball in kg.
• As the 2nd ode is broken down into two ode, the initial condition has been defined for both of the ode's, theta_0 = [0,3] the 1st condition is for the displacement while the 2nd condition is for the velocity.
• t = np.linspace(0,10,200), creating the time array between 0-10, with the help of numpy 200 values are created between 0-10 seconds.
• Solving the ode with the help of the odeint command which is the command in the scipy module. It will integrate the ode.
• theta = odeint(model,theta_0,t, args=(b,g,l,m)) defining the system, initial condition, time array, and the arguments.
• Plotting the results, theta is the solution vector with 2 columns & multiple rows.
• The theta_0 = [0,3], displacement is zero & velocity is 3 m/s.

Output:

Plot,

Result:

• The ode solved and plot for the simple pendulum generated.
• The animation of simple pendulum was created. Link given below.

https://drive.google.com/file/d/1stijIvhXDd0bnCdbOpbQz1XJnlg7bmyT/view?usp=sharing

Errors faced while programming:

• Unexpected indent, out of function.
• giving correct spacing, get into function.