Week 3 - Solving second order ODEs

close all
clc

%inputs
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

%initial condition
theta_0 = [0,3];

%time points
t_span = linspace(0,10,250);

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

%plotting results
figure(1)
plot(t, results(:, 1),'linewidth',1)
hold on
plot(t, results(:, 2),'linewidth',1)
ylabel('displacement')
xlabel('time')
grid on
legend('displacement', 'angular velocity')

%simulation

%inputs
x0 = 0;
y0 = 0;
ct = 1;

for i = 1:length(t)
    theta = results(i,1);
    x1 = l*sin(theta);
    y1 = l*(-cos(theta));
    
    figure(2)
    plot([-1 1],[0 0],'linewidth',3,'color','b');
    hold on
    p = plot([x0 x1], [y0 y1], 'b');
    hold on
    plot(x1,y1,'marker','o','markersize',20,'markerfacecolor','b')
    title('motion of the pendulum')
    hold off
    axis([-2 2 -2 0])
    grid on
    pause(0.03)
    M(ct) = getframe(gcf);
    ct = ct+1;
end

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


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

https://youtu.be/6zvVlCHt7r8