Week 4.1 - Solving second order ODEs

OUTPUT ANIMATION LINK:

https://youtu.be/VeUvBUIMoXM

MAIN PROGRAM:

close all %delete all plots  
clear all %delete all variables 
clc %clear command window 
 
%inputs
b = 0.05;
g = 9.81;
l = 1;
m = 1;

%intial condition 
theta_0= [0;5];

%time period
t_span = linspace(0,20,500);

%solving partal differential equation 
[t, results] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Plot')
ylabel('Time')
legend('Dispacement','Angular velocity')

% Pendulum Animation
ct = 1
for i = 1:length(results(:,1))
    %coordinates for the fixed end
    x0 = 0;
    y0 = 0;
    %coordinates
    x1 = 1*sind(results(i,1));     
    y1 = -1*cosd(results(i,1));
    
    figure(2)
    plot([-0.9 0.9],[0 0],'linewidth',2);
    hold on
    plot([x0 x1],[y0 y1],'linewidth',2,'color','g');
    plot(x1,y1,'o','Markersize',25,'markerfacecolor','r');
    hold off
    axis([-0.2 0.2 -1.5 0]);
    pause(0.01)
    M(ct)=getframe(gcf);
    grid on
    ct = ct+1;
    
end
movie(M)
videofile = VideoWriter('Pendulum.avi','uncompressed avi');
open(videofile)
writeVideo(videofile, M)
close(videofile)

FUNCTION CODE:

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

OUTPUT: