Week 3 - Solving second order ODEs
SOLVING 2ND ORDER ODE
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.
SOLUTION:
1) CREATING ODE FUNCTION
Thus ODE function is created which can be used to code and it will look like below
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
2) CODING INSTRUCTIONS:
a) As per example problem solved in lecture, follow the code format : library import, define given parameters, initial condition, time points and the function to solve ODE.
b) Plot the figures with proper title, label,legend.
c) To animate the result, Use previous idea explained in previous lectures. Create images for each stage and stitch them together as animation.
3) PYTHON CODE
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def ode_func(theta, t, b, g, l, m):
theta0 = theta[0]
theta_dot = theta[1]
theta_ddot = -b/m * theta[1] - g/l * math.sin(theta[0])
return [theta_dot, theta_ddot]
# Inputs
b = 0.05 # Damping Coefficient
g = 9.81 # Gravity m/s^2
m = 1 # Mass of the Pendulum in kg
l = 1 # Length of the string attached to pendulum in meters
# Initial conditions
theta_0 = [0, 3] # [angular displacement=0; angular velocity=3]
# Time points
t_span = np.linspace(0, 20, 500) # Pendulum oscillates between 0-20 seconds
# Solve ODE
result = odeint(ode_func, theta_0, t_span, args=(b, g, l, m))
# Plotting the graph
plt.figure()
plt.plot(t_span, result[:, 0], 'b')
plt.xlabel('Time')
plt.ylabel('Theta')
plt.title('Angular Displacement (theta) vs Time')
plt.figure()
plt.plot(t_span, result[:, 1], 'r')
plt.xlabel('Time')
plt.ylabel('Velocity')
plt.title('Angular Velocity vs Time')
plt.show()
# Pendulum simulation using ODE
theta = result[:, 0]
# Creating images for each time step
for i in range(len(result)):
THETA = result[i, 0]
# Initial boundary condition
x0 = 0
y0 = 0
# Defining x and y coordinates
x1 = l * np.sin(THETA)
y1 = -l * np.cos(THETA)
plt.figure(figsize=(6, 6))
plt.plot([x0, x1], [y0, y1], 'k-', linewidth=4)
plt.plot(x1, y1, 'go')
plt.plot(x0, y0, 'ro')
plt.axis([-2, 2, -2, 0])
plt.title('Pendulum at t = {:.2f}'.format(t_span[i]))
plt.savefig('pendulum_{:03d}.png'.format(i))
plt.close()
4) RESULTS
Angular Velocity vs Time
as per given values angular velocity starts from 3m/s and displacement starts from 0.
Angular Displacement (theta) vs Time
ANIMATION:
Animation file has been uploaded here and in Youtube as well.