Week 3 - Solving second order ODEs

clear all
close all
clc
%%
%input parameters (b-damping coefficient, g-gravity,l-length,m-mass)
b = 0.05;
g = 9.81;
l = 1;
m = 1;

% initial condition :- at intial postion speed is equal to 3
theta_0 = [0;3];

% time points
t_span = linspace(0,20,600);

% solving ODE using ode45

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

%%
%plotting angular velocity
%creating separate graphs for angular displacement and angular velocity
tiledlayout(2,1) 
nexttile
plot(t, results(:,1))
ylabel('angular displacement')
xlabel('time')
grid on;
%plotting angular acceleration
nexttile
plot (t, results(:,2))
ylabel('angular velocity')
xlabel('time')
grid on;

%%
%animation
x0 = 0;
y0 = 0;
%counter represents each frame
ct=1; 
for i=1:length(t)
    %x coordinates as function of t and theta
    x1 = l*sin(results(i,1));
    %y coordinates as function of t and theta
    y1 = -l*cos(results(i,1));
    
    figure(2)
    plot([-.5 .5],[0,0], 'linewidth',3,'color','b');
    axis([-1.5 1.5 -1.5 .5]);
    line([x0 x1],[y0 y1], 'linewidth',2,'color','g');
    hold on;
    plot(x0,y0,'o','markersize',6,'markerfacecolor','w');
    plot(x1,y1,'o','markersize',20,'markerfacecolor','r');
    title('Damped simple pendulum');
    legend('Support rod','length of pendulum','fixed point', 'ball of pendulum');
    grid on;
    hold off;
    M(ct) = getframe(gcf);
    ct = ct+1;
end
%creating movie file
movie(M);
videofile = VideoWriter('pendulum');
open(videofile);
writeVideo(videofile,M);
close(videofile);

%%
%code for the 2nd order ODE function
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
%theta1 is first derivative or angular velocity
%theta2 is second derivative or angular accelaration
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