Week 3 - Solving second order ODEs

AIM: To write program to plot and animate the motion of pendulum by solving the 2nd order ODE equation for pendulum . Please find the equation below .

Explanation : The Pendulum we consider is assumed to be suspended from a frictionless static surface by a massless rod of length of 'l' and the susupended pendulum is having a mass of ~0.1 . The Theta (angle with the y axis ) is varied for varioius values to plot the graph. Please find below the pendulum specification below :

FIG1 pendulum picture

Assumptions : [Note : these are the values considered in the program also .Please find the assumptions below]

damping co-efficient - b=0.05

gravity - g=9.81`m/s^2`

Length of rod - l=1m

mass of pendulum - m= 0.1

Formulae simplification to 1st order

-- Main formulae

When ,

`theta=theta1` ---> Assumpion 1

`d(theta1)/dt = theta2` --> Assumption 2

`d(theta2)/dt = -(b/m)(theta2) -(g/l)*sin(theta1)`---> Simplified formulae by applying assumption 1&2 in main formulae

Program :

close all
clear all
clc

%Solving 2nd Order ODE45 Equation 
%Declaration of variables 
b=0.05; %damping co-efficient 
g=9.81; %gravitational pull 
l=1; %length of pendulum 
m=0.1; %mass of pendulum 

%initial values declaration  
theta_0=[0;5];

%Time Values
time=linspace(0,10,200);

%ODE45 solution 
[t,results]=ode45(@(t,theta) odefunc(t,theta,b,g,l,m),time,theta_0);

%plotting 
figure(1)
plot(t,results(:,1));
hold on
plot(t,results(:,2));
grid on 
ylabel 'plot'
xlabel 'time'

figure(2)
ct=1;
for i=1:length(time) 
    x=0;
    y=0;
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
    plot([x x1], [y y1],'LineWidth',2)
    xlabel 'time'
    ylabel 'plot'
    grid on
    hold on
    plot(x1,y1,'o','MarkerSize',25,'MarkerFaceColor','r')
    hold off 
    axis([-1.5 1.5 -2 1])
    M(ct)=getframe(gcf);
    ct=ct+1;
end

%animation creation 
Videofile=VideoWriter('pendulum.avi');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)

FIG 2 Main Program for ploting

Program for Function :

function [dtheta_dt]=odefunc(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

FIG 3 Function fragment

Explanation of code :

close all
clear all
clc

//Solving 2nd Order ODE45 Equation
//Declaration of variables
b=0.05; //damping co-efficient
g=9.81; //gravitational pull
l=1; //length of pendulum
m=0.1; //mass of pendulum

%initial values declaration
theta_0=[0;5];

%Time Values
time=linspace(0,10,200); //Time is takes for 0-10s

%ODE45 solution
[t,results]=ode45(@(t,theta) odefunc(t,theta,b,g,l,m),time,theta_0);

// The syntax for ODE45 is appllied for the current scenario

%plotting
figure(1)
plot(t,results(:,1)); // The plotting for the necessary values is done using plot command
hold on
plot(t,results(:,2));
grid on
ylabel 'plot'
xlabel 'time'

figure(2)
ct=1; // to animate the required motion of pendulum the procedure is initiated by ct=1
for i=1:length(time) // for loop implimenred to iterate each time to get the motion of pendulum in animation
x=0;
y=0;
x1=l*sin(results(i,1));
y1=-l*cos(results(i,1));

// animation for pendulum was in the upward direction ; to get the pendulum facing downward the y value is condidered in -ve
plot([x x1], [y y1],'LineWidth',2)
xlabel 'time'
ylabel 'plot'
grid on
hold on
plot(x1,y1,'o','MarkerSize',25,'MarkerFaceColor','r')
hold off
axis([-1.5 1.5 -2 1])
M(ct)=getframe(gcf);
ct=ct+1;
end

%animation creation // Trying to save the file in a avi format
Videofile=VideoWriter('pendulum.avi');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)

Errors Encountered:

• The markersize command syntax was implemented worngly , was rectified in a 1-on-1 session . Thanks FIG4 Error snapshot

*****************************************END*********************************************

OUTPUT :

1.Figure 1

ANIMATION :
Link : https://drive.google.com/file/d/1dySObILRWfj9vYhzIkc-PlFSU-Qj9e1V/view?usp=sharing

Thanks ,

Vigneshwar Viswanathan