Week 3 - Solving second order ODEs

https://youtu.be/UcAJvRuRb_Y




clear all
close all
clc

% Solving a second order ODE which defines the motion of a pendulum

% Inputs

b = 0.005;
g = 9.81;
l = 1;
m = 1;
ct = 1;

% Initial condition

theta_0 = [0;3];

% Time points

t_span = linspace(0,20,500);

% Solving ODE

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

% Plotting

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Variation of displacement and velocity (angular) w.r.t')
ylabel('time')

% Plotting the pendulum

position_theta = [results(:,1)*(180/pi)]


for i =1:length(position_theta)
    POS_THETA = position_theta(i)


    x0 = 0;
    y0 = 0;

    x1 = l*sind(POS_THETA);
    y1 = -l*cosd(POS_THETA);


figure(2)
plot([-1 1],[0 0],'color','y','LineWidth',5)
hold on
plot([x0 x1],[y0 y1],'-o','color','r','linewidth',1)
hold on
plot(x1,y1,'o','Color','b','LineWidth',10)
hold off
axis([-2 2 -2 2])
pause(0.02)
M(ct) = getframe(gcf);
ct = ct+1;
end

movie(M);
videofile = VideoWriter('simple_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

% Program for defining function 

function[dtheta_dt] = ode_func(t,theta,b,g,l,m,t_span,theta_0);

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