Solving second order ODE

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

Theory -

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

Simple Pendulum - The simple pendulum is an idealized mathematical model of a pendulum.This is a weight (or bob) on the end of a massless cord suspended from a pivot, without friction. When given an initial push, it will swing back and forth at a 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. Weight of the bob acts vertically downwards and the tension acts vertically upwards, so, the bob comes under equillibrium position and the bob starts oscillating 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 , `b` is damping coefficient of the system.

1. `(d^2theta)/dt^2 +b/m*(d(theta))/dt+g/Lsintheta = 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(theta_1))/dt` = `theta_2` ---------(2)

therefore `(d^2theta)/dt^2 = (d^2theta_1)/dt^2 = (d(theta_2))/dt`-----(3)

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

`(d^2theta_1)/dt +b/m*(d(theta_1))/dt +g/L *sintheta_1 = 0`

therefore , `(d(theta_2))/dt +b/m *(theta_2) + g/L*sintheta_1 = 0`

we got, `(d(theta_2))/dt = -b/m*(theta_2)-g/Lsintheta_1`

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 metre

M = 1 kg

b = 0.05

g = 9.81 `m/s^2`

MATLAB code file -

1.Code function of 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

% Solving the Second order ODE and show the animation

b = 0.05;                     % damping coefficient
l = 1;                        % length of pendulum in metre
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(displacement)
    
    %coordinates
    x1=1-l*sin(displacement(i));
    y1=1-l*cos(displacement(i));
    
    figure(2)
     %plot fixed link
     plot([0 2],[1 1],'linewidth',15,'color','b')
     hold on                              
     plot([1 x1],[1 y1],[x1 x1],[y1 y1],'-O','markersize',25,'markerfacecolor','r','linewidth',4)
     plot(x0,y0,'-O','markersize',3,'markerfacecolor','b')
     %axis setting
     axis([0 2 -0.5 2])
     
     M(ct)= getframe(gcf);                 % to change the frame of the movie
     ct=ct+1;
     hold off
end
movie(M)                                                                %for creating animation video
videofile=VideoWriter('simple pendulum.avi','uncompressed AVI');
open(videofile)                                                         %open the videofile while running the code
writeVideo(videofile,M)                                                 %for writing the video file
close(videofile)                                                        %for closing the videofile after code run

Code Explanation -

1.Firstly, the given details are assigned as variables in the code. Then time span of the simple pendulum is set by using linspace 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 identify each values in the graph.

3.For getting the accurate simulation of the simple pendulum corresponding 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 making 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/CeoGsPBqUuY

Conclusion -

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

2. Animation of Simple pendulum is created.