Week 3 - Solving second order ODEs

AIM : To solve second order ordinary differential equation using MATLAB.

Theory- Ordinary Differential Equations (ODEs) - An Ordinary differential equation (ODE) is a differential equation contating one or more functions of one independent variable and the derivatives of those funtions.The term oridnary is used in contrast with the term parital differential equation which may be with respect to more than one independent variable.

Simple Pendulum - The simple pendulum is an idealized mathmatical model of a pendulum.This is a weight (or bob)on the end of massless cord suspended from a pivot,without friction.When given an initial push,it will swing back and forth at constant amplitude. Real pendulums are subject to friction and air drag,so the amplitude of their swings declines.

The point from which the bob is suspended is known as point of suspension.Wegiht of the bob acts vertically downwards and the tension acts equillibrium position and the bob starts oscillacting if even smaller external force is applied.

Governing Equations -

`theta` is the angle made by the movable link with vertical axis`

'L' is the length of the movable rod

'm' is the mass of bob,`beta` is damping coefficient of the system.

1. `d^2theta/dt^2+b/m*d(theta)/(d t)+g/L sin theta =0` second order ODE -------(1)

To solve this ODE we need to reduce the order of equation into lower order

let `theta`= `theta 1`

so this implies `d(theta)/(dt) = d(theta1)/(dt)=theta2`---------(2)

therefore `d^2(theta)/(dt^2) = d^2(theta1)/(dt^2)=d(theta2)/dt`-----------(3)

Putting the value of (2) and (3) in (1),we get-

`d^2 (theta1)/dt^2+b/m*d(theta1)/(d t)+g/L sin theta1 =0`

therefore,

`d (theta1)/dt^2+b/m*(theta2)+g/L sin theta1 =0`

we got,

`d(theta1)/dt^2=-b/m*(theta2)-g/L sin theta1`

Objectives-

1. To solve the second order ODE with graph of Velocity v/s Displacement.
2. To create the animation of the simple pendulum.

Procedure-

Given details-

simulate the motion between 0-20 sec for angular displacement =0

Angular velocity =3 rad/sec

at time t = 0

L = 1 meter

M = 1 kg

b = 0.05

g = 9.81 m/s^2

MATLAB code file -

1. Code functioof ode45 for main program -

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
theta1 = theta(1);
theta2 = theta(2);

dtheta1_dt = theta2;
dtheta2_dt = (-(b/m)*theta2 - (g/l)*sin(theta1));
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

2.Main program -

clear all
close all
clc

%sloving the second order ODE and showing the animation

b = 0.05;           %damping coefficient
l = 1;              %length of pendulum in meter
m = 1;              %mass of bob in kg
g = 9.81;           %gravity in m/s^2

%initial condition 
theta_0 = [0;3];

%Time points 
t_span = linspace(0,20,500);

%solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Velocity')
legend('Displacement','Angular velocity')
grid on
hold off

displacement = results(:,1);
velocity = results(:,2);

%fixed coordinates

x0=1;
y0=1;
ct=1;

for i=1:length(results(:,1))
    x_0=0;
    y_0=0;
    x_1=l*sin(results(i,1));
    y_1=-l*cos(results(i,1));
   figure(2) 
    plot([-1.5 1.5 ],[0 0])
    hold on
    plot([x_0 x_1],[y_0 y_1],'linewidth',1,'color','b')
    axis([-1.5 1.5 -1.5 1.5])
    plot(x_1,y_1,'o','MarkerSize',20,'color','r')
    pause(0.0001)
    M(ct)=getframe(gcf);
    hold off
    ct=ct+1;
end
movie(M)
videofile = VideoWriter('ode_main.m','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Code Explaination-

1. Firstly,the given details are assigned as variables in the code.Then time span of the simple pendulum is set by using linspce command then,theta value is given as a column matrix. After this we use ode45 function and corresponding 'ode_func' function call in main program.
2. Then used plot command to plot graph corresponding to time and velocity,legend command is used for identiy each values in the graph.
3. For getting the accurate simulation of the simple pendulum correspoding to the second order ODE with damping we need to use the looping statements,here for loop is used to create this.
4. The pendulum has fixed coordinates mentioned as x0,y0 and x1,y1 are the other coordinates.
5. Marker command is used for marking bob of the pendulum,The axis is fixed using axis command. The counter command is used to make the value changes in each for loop. Then videoWriter command is used to create the video file.

Results -

Outputs -

1. Graph for the velocity v/s time which shows Displacement and Angular velocity-

The Displacement (blue color) and Angular velocity (red color) is decreasing by the increased time,this is due to the introduction of damper i.e. Damping coefficient which absorbs the vibrations to make the body in equillibrium position.

2. Position of the setup -

3.Animated file -

https://youtu.be/0oT3benxigc

Conclusion -

1.The graph of Displacement and Angular velocity is plotted in velocity v/s Time diagram.

2.Animation of simple pendulum is created.