Week 3 - Solving second order ODEs

clear all

close all

clc

%All units in the SI system

b = 0.05;

g= 9.81;

l=1;

m=1;

%initial conditions

theta_r= [0;3];

%time points

t_span = linspace(0,20,300);

%calling the ODE function

[t,results]= ode45(@(t,theta) ode_funct(t,theta,b,g,l,m), t_span, theta_r);

%hold on

%plot(t,results(:,1))

%plot(t,results(:,2))

x label('angular displacement')

y label('time')

ct = 1;

for i=1:length(t_span)

x0 = 0;%fixed end of the pendlum

y0=0;

x1 = l*sin(results(i,1));%position of the bob for the time value

y1=-1*cos(results(i,1));

grid on

plot([-0.5 0.5],[0 0], 'linewidth',3,'color','black');%plotting the fixed end

hold on

plot([x0 x1],[y0 y1],'linewidth',3);%plotting the pendulum rod

plot(x1,y1,'marker','0',markersize',15,'markersizecolor','b')%plotting the bob

hold off

axis([-2 2 -2 2])

pause(0.001);

title('Simple pendulum');

m(ct)= getframe(gcf);

ct= ct+1;

end

%to create a movie

movie(M)

videofile= videowriter('pendulum.avi','Uncompressed AVI');

open(videofile)

writevideo(videofile,M)

close(videofile)

ode_funct code:

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