Aim:- Solving second order ODEs

• In mathematics, an ordinary differential equation(ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. Mathematical descriptions of change use differentials and derivatives. Various differentials, derivatives, and functions become related via equations, such that a differential equation is a result that describes dynamically changing phenomena.

In engineering, ODE is used to describe the transient behaviour of a system. A example is a simple pendulum. The position(theta) of the bob with respect to time of a simple pendulum can be determined by using the second order differential equations.

• Governing equations:-

The second order ode for simple pendulum will be as follows,

Programme written to solve 2nd order ordinary differential equation (ODE)

To solve 2nd order ODE converted to 1st order such as expaline below:

`1st order ODE is [(dθ1/dt)=0] ,

2nd order ODE is [ (dθ1^2/dt + (b/m.(dθ1/dt))+(g/l.(sinθ1))=0 ]`

Now,

(dθ^2/dt) = (dθ2/dt),

hence dθ/dt=θ2

hence equation becomes [(dθ1^2/dt + (b/m.(dθ1/dt))+(g/l.(sinθ1))= (dθ2/dt)+ b/m.(θ2)+ g/l.(sinθ2)]

Where, b = damping coefficient

m = mass of the bob ( kg),

g = acceleration due to gravity( m/s^2),

l = length of the string.(m)

The simple pendulum problem can be solved by using the second order ode functions and by solving these we will get the theta1 an theta2 values of the simple pendulum

3.Objectives of the project:

1. to write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping.

4.Solution procedure:

Introduction:

Steps in code

1.Programme wtitten in python using spyder notebook

2.Reuired libraries are imported such as: math, matplotlib,numpy,scipy.integrate

2.All the Input parameters values are given to programme such as: b,l,g,m

3.function model is defined that represent parameters of ODE.

3.Initial values of theta, time,angular velocity are entered

4. Use of odeint function by importing it from scipy.integrate library is done to solve problem.

3. Explanation of ‘odeint’ function syntax
4. model-defines function of ode
5. model(theta,t,b,g,l,m)-defining function,
6. args are arguments i.e( b,g,l,m)
7. t-time array to decide integration limits,
8. theta_0-initial condition of pendulum.

5.Commands to generate plot and position of pendulum and video file are written.

# -*- coding: utf-8 -*-
"""
 Programme for ordinary differential equation(ODE)
 1st order ODE is [(d(theta1)/dt)=0]  
 2nd order ODE is [ (d(theta1)^2/dt + (b/m.(d(theta1)/dt))+(g/l.(sin(theta1)))=0 ]
 to solve 2nd order ODE converted to 1st order such as considered (d(theta)^2/dt)=(d(theta2)/dt), hence d(theta)/dt=(theta2)
hence equation becomes [(d(theta1)^2/dt + (b/m.(d(theta1)/dt))+(g/l.(sin(theta1))= (d(theta2)/dt)+ b/m.(theta2)+ g/l.(sin(theta2))]
"""

import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint  ## importing ode function
import numpy as np



## argument Parameters
###  b-damping coefficient, m- mass, g-gravity, l-length of pendulum

## function for ode
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
    

b=0.05     ## dampig coefficient
g=9.81     ## gravitational force [m/s^2)
l=1       ##lengh of pendulum [m]
m=1       ##  mass of pendulum bob [kg]

# Intial condition for pendulum
theta_0 = [0,5]  ## intial displacemnt of pendulum is 0 with velocity 5 m/s

# time points
t= np.linspace(0,10,150)   ## time array i.e time increseaed from 0 to 10 sec in 150 divisions

theta= odeint(model,theta_0,t,args=(b,g,l,m))
'''
odeint is function to integrate ordinary diffferential equation reperesented by model.
'''
## plotting graphs
plt.plot(t,theta[:,0],label=r'Angular displacement')
plt.plot(t,theta[:,1],label=r'Angular Velocity')
plt.xlabel('Time')
plt.ylabel('Plot')

plt.show()
plt.plot(theta[:,0],theta[:,1])
plt.xlabel('time')
plt.ylabel('velocity')
plt.show()

## animation creation
ct=1   
x0=0;
y0=0;
for i in range (0,len(t)):
    x1=x0+1*math.sin(theta[i,0])
    y1=y0-1*math.cos(theta[i,0])
    x2=-1
    y2=0
    x3=1
    y3=0
    filename= 'pendulum'+str(ct)+'.png'   ## defining how to assign name to generating image files for each theta value
    plt.plot([x0,x1],[y0,y1])
    plt.plot([x2,x3],[y2,y3])
    plt.xlim([-3,3])     ## defining limits for axis of graph
    plt.ylim([-3,3])
    plt.savefig(filename)
    ct=ct+1

Gif/Animation of 2nd order ODE(Pendulum):

Error faced:

Error solution: syntax is wrong corrected as

x1=x0+1*math.sin(theta[i,0])

Plots:

Variables:

final pendulum all position tracing graph:

screenshot of graphs generated:

Result and Discussion:

Total 150 plots are generted for pendulum position.Starting position of pendulum is 0deg. with velocity of 5m/s.

and position changes in between 0deg. to 90 deg.Timespan is of 0 to 10 divided in 150 unit interval.so plots are 150 in no. Animation shows how the posoition of pendulum changes.