Week 3 - Solving second order ODEs

Project title: Python programme for pendulum motion

Objective: 1. generate MatLab code to solve the second-order partial differential equation.

2. generate a pendulum simulation video.

theory:

In many real-life modelling situations, a differential equation for a variable of interest won’t just depend on the first derivative, but on higher ones as well. Naturally then, higher-order differential equations arise in STEP and other advanced mathematics examinations. For anything more than a second derivative, the question will almost certainly be guiding you through some particular trick very specific to the problem at hand. For second-order differential equations though, you need to know how to tackle them in general. Fortunately, the technique involved is straightforward, and this article guides you through all you need to know, with a helpful example as well.

Homogeneous Second Order Differential Equations

The first major type of second-order differential equations you’ll have to learn to solve are ones that can be written for our dependent variable y and independent variable t as:

Here a, b and c are just constants. In general, the coefficients next to our derivatives may not be constant, but fortunately, you don’t need to worry about how to approach such problems like that in general for STEP.

Inhomogeneous Second Order Differential Equations

One extension to the above that we must tackle is the case when the RHS in our DE is non-zero, i.e. when we have:

Simple Pendulum

Assumption

1. mass of the rod is not considered.

2. mass of the ball is concentrated on a single point.

Pendulum follows below second order differential equation.

So, in the end, we solve the matrix.

Problem Statement:

In Engineering, ODE is used to describe the transient behaviour 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 these 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.

Solution:

Programme Logic:

First second order differential equation is converted into the first-order differential equation.

Second Create a function in python which solve the ODE

Third create time array initial condition array

Forth call the ODE function with input arguments

Fifth Use plotting command to plot the graphs and animation.

Code

"""
Programme to solve the Ordanary Differential Equation
By Rahul Date
"""
from math import *
from numpy import * 
from matplotlib.pyplot import * 
from scipy.integrate import *
def model(theta,t,b,g,l,m):
	theta1 = theta[0]
	theta2 = theta[1]
	dtheta1_dt = theta2
	dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1)
	dtheta_dt = [dtheta1_dt,dtheta2_dt]
	return dtheta_dt
b = 0.05
g = 9.81
l = 1
m = 0.1
# Initial condition 
theta_0 = [0,3]
# Time points
t = linspace(0,20,150)
# solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m))
# Plot results
figure(1)
plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta2-\frac{g}{l}sin\theta_1$')
ylabel('Displacement & Velocity')
xlabel('time')
legend(loc='best')
show()
# Animation 
ct = 1
for i in range(len(theta)):
	figure()
	x0 = 0
	y0 = 0
	x1 = sin(theta[i,0])
	y1 = -cos(theta[i,0])
	filename = str(ct) +'.png' 
	cord = plot([x0,x1], [y0,y1],color = 'k', lw = 2)
	bob = plot(x1,y1,'o',markersize=15, color = 'r')
	xlim(-1, 1)
	ylim(-1.5, 0)
	ct = ct+1
	savefig(filename)

Code description:

"""
Programme to solve the Ordinary Differential Equation
By Rahul Date
"""
from math import * ....................................................math library import for math operation
from numpy import * ..................................................numpy library import for trinometric operations
from matplotlib.pyplot import * ...................................matplotlib.pyplot library import for ploting command
from scipy.integrate import * ...................................... Scipy.intergrate library import for solving ordinary differential equation
def model(theta,t,b,g,l,m): .......................................... Calling ordinary differential equation function
theta1 = theta[0] ...................................................... theta1
theta2 = theta[1] .......................................................theta2
dtheta1_dt = theta2 ...................................................derivative of theta1
dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1) ............... derivative of theta2
dtheta_dt = [dtheta1_dt,dtheta2_dt] ............................ matrix formation od theta1 and theta2
return dtheta_dt ......................................................... return down the theta value
b = 0.05 ..................................................................... Dammping coefficient
g = 9.81 ..................................................................... Gravitational force
l = 1 .......................................................................... Length of pendulum
m = 0.1 ...................................................................... Mass of bob
# Initial condition
theta_0 = [0,3] ............................................................Initial condition for displacement and velocity
# Time points
t = linspace(0,20,150) ...................................................Time array from 0 to 20 with 150 time steps
# solve ODE
theta = odeint(model,theta_0,t,args = (b,g,l,m)) ...............Solve ordinary differential equation
# Plot results
figure(1)
plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plot(t,theta[:,1],'r--',label=r'$\frac{d\theta_2}{dt}=-\frac{b}{m}\theta2-\frac{g}{l}sin\theta_1$')
ylabel('Plot')
xlabel('time')
legend(loc='best')
show()
# Animation
ct = 1
for i in range(len(theta)):
figure()
x0 = 0
y0 = 0
x1 = sin(theta[i,0])
y1 = -cos(theta[i,0])
filename = str(ct) +'.png'
cord = plot([x0,x1], [y0,y1],color = 'k', lw = 2)
bob = plot(x1,y1,'o',markersize=15, color = 'r')
xlim(-1, 1)
ylim(-1.5, 0)
ct = ct+1
savefig(filename)

Output

Fig shows Displacement and velocity Vs Time graph

Fig determine the behaviour of second-order differential equation with input argument for 20 sec

Animation

Animation of the simple pendulum