Week 3 - Solving second order ODEs
AIM:-
1.To obtain a numerical solution for the ordinary differential equation which represents the damped motion of a simple pendulum.
2. Simulating the transient behavior of the pendulum.
THEORY:-
A simple pendulum consists of a massless rod/string of length 'L', which is fixed at one end and a bob of mass 'm'is attached to the other end of the rod.
Fig(1): Simple Pendulum
Its position concerning time t can be described merely by the angle `theta` measured against a reference line, usually taken as the vertical line straight down.
Mathematical Formulation:-
Differential equations are widely observed in physics for describing physical phenomena as they describe how anything changes as it involves a derivative.
In the case of a simple pendulum; the equation describes how its position (angle '`theta`') is varying with the time while also taking into account the other factors like mass, damping co-efficient etc.
The differential equation which represents the motion of the simple pendulum is given by
`(d^2 theta)/dt^2 +(b/m) * ((d theta)/dt) + (g/l) * sin theta = 0`
where m = mass of the bob,
b = damping co-efficient,
g = gravitational acceleration = 9.8 `m/(s^2)`,
l = Length of the rod/string.
Now for solving the second order linear differential equation in MATLAB we break it down into first order ODEs.
Now lets say initially `theta = theta_1`
Therefore;
`(d^2 theta_1)/dt^2 +(b/m) * ((d theta_1)/dt) + (g/l) * sin theta_1 = 0 -> (I)`
Now if there is a change from position `theta_1` in time t. We can say that:
`((d theta_1)/dt) = theta_2 ->(II)`
Now;
`(d^2 theta_1)/dt^2 = (d/dt)*((d theta_1)/dt) = (d theta_2)/dt -> (III)`
Substituting (II) & (III) IN EQUATION (I)
`(d theta_2)/dt +(b/m) * theta_2 + (g/l) * sin theta_1 = 0`
`(d theta_2)/dt = -b/m *theta_2 - g/l * sin theta_1 -> (IV)`
Equation (2) and (4) are the Two first order ODE's which will help to find the solution of the ODE.
ASSUMPTIONS:-
1. The string/rod is considered to be massless.
2. The bob is assumed to be initially in the horizontal straight position.
3. Air resistanc eis negligible
SOLUTION:-
Python was used to obtain the numerical solution to the diffferential equation by using the odeint function from the scipy module.
CODE:-
# PROGRAM: Numerical Solution to the ODE for damped motion of the pendulum.
# Author: HARSH SHAH
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
# Defining Function
def pendulum(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
#Input Parameters
b = 0.6 # DAMPING CO-EFFICIENT
g = 9.8 # GRAVITATIONAL ACCELERATION
l = 2 # lENGTH OF THE PENDULUM
m = 2 # MASS OF THE BOB
# Initial Conditions
theta_0 = [math.pi/2,0] # Angular displacement & Angular velocity
# Time points
t = np.linspace(0,10,101)
# Solving using the odeint command
theta = odeint(model, theta_0, t, args = (b,g,l,m))
# Plotting
plt.plot(t, theta[:,0], 'b', label = r'$\frac{d\theta_1}{dt} = \theta_2$')
plt.plot(t, theta[:,1],'g--', label = r'$\frac{d\theta2}{dt} = -\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.xlabel('Time')
plt.ylabel('Plot')
plt.legend(loc = "best")
plt.grid()
plt.show()
# Generating images to create animation
# Counter for plots
ct = 1
# Array for accessing angle values
disp = theta[:,0]
for angle in disp:
# Initial position of the bob
x0 = 0
y0 = 0
# Computing positions
x1 = l * math.sin(angle)
y1 = - l * math.cos(angle)
filename = ' %04d.png' %ct # Generates filename as 0001.png
ct = ct + 1
# Generating plot
plt.figure()
plt.plot(0,0, marker = 'o', markersize = '4')
plt.plot([x0, x1], [y0,y1], 'b-')
plt.plot (x1, y1, marker = 'o', markersize = '9')
plt.xlim((-5,5))
plt.ylim((-5,5))
plt.grid()
plt.savefig(filename)CODE ANALYSIS:-
1. Importing required libraries like math, matplotlib, numpy. Here we directly import the command odeint from scipy.integrate module.
2. Let us look at the code by breaking it down into different sections:
- def pendulum(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 - b = 0.6
g = 9.8
l = 2
m = 2 - theta_0 = [math.pi/2,0]
t = np.linspace(0,10,101) - theta = odeint(pendulum, theta_0, t, args = (b,g,l,m))
Section D:- Solving using odeint.
ODEINT integrates a system of differential equations. It solves first order equations, which is the primary reason for breaking down the original equation into a system of two first order ODE's.
In the program, we have passed certain arguments to the odeint function like:-
- pendulum which is a function specified by us containing the mathematical formulations of the two first order ODEs that will be used for calculating theta at time t.
- theta_0 represents the initial conditions like displacement and velocity specified by us.
- t is the time array. The values of time for which the value for `theta` is to be computed.
- args are some extra arguments that need to be provided. Like in our case equation has constants like b,g,l,m
It outputs the values of angular dispacement and angular velocity.
Section A:- Defining the function.
Initially; we define a function pendulum which shall return us an array dtheta/dt which consists of values dtheta1/dt and dtheta2/dt which are represented by equations (II) and (IV) from theory section.
Now from equations (II) and (IV) we know parameters like 'theta', 't' (time), 'b' (damping co-efficient), 'g' (gravitational acceleration), 'l' (length of the rod), 'm' (mass of the bob) will be neccesary to be specified for the calculation.
So we pass on this parameters as arguements for the function pendulum. theta[0] extract the first value from the initial conditions array i.e. theta_0 (Look C.) and similarly theta[1] indexes the second value from theta_0 array.
Then we represent the two ODEs in a code format and store them in an array dtheta/dt.
Section B, C:- Specifying constant parameters, initial conditions and time array.
The pendulum function requires certain parameters to be specified. So we assign the constant parameters. We generate a time array as we want to compute for `theta` varying with time t.
Section E:- We generate plots for variation of computed values of `(d theta_1)/dt` and `(d theta_2)/dt` varying with respect to time 't'.
Labels and legends are used to specify the plot denoting the parameter.
Section F:- Creating images which shall be used for simulating the motion of the pendulum.
ct = 1
disp = theta[:,0]
for angle in disp: x0 = 0
y0 = 0 x1 = l * math.sin(angle)
y1 = - l * math.cos(angle) filename = ' %04d.png' %ct
ct = ct + 1 plt.figure()
plt.plot(0,0, marker = 'o', markersize = '4')
plt.plot([x0, x1], [y0,y1], 'b-')
plt.plot (x1, y1, marker = 'o', markersize = '9')
plt.xlim((-5,5))
plt.ylim((-5,5))
plt.grid()
plt.savefig(filename)
A counter was used for storing the count of the generated images. The counter increments with the number of iterations processed by the for loop.
An array disp was used to store the values of theta i.e the angular displacement.
Then using a for loop; the cartesian positions (X and Y) were calculated for different angle values. Likewise, the plots for varying positons calculated in the for loop were stored as different images was stored.
Later the images were stitched together to create the animation file for the pendulum.
PLOTS:-
The plot for `(d theta_1)/dt` and `(d theta_2)/dt` varying with respect to time 't' is as follows:-
The video of the simulated motion of the pendulum can be found at the following link:-
From the plot and video; it is clearly observed that the amplitude of the pendulum decreases with time.
This decrease is due to the damping factor which is clearly evident from the mathematical formulation.
By varying the parmeters we can observe the changes in plots and various iterations can be achieved like that of a overdamped system or critically damped system.
It has applications in designing suspensions, or areas where mechanical vibrations play a major role.