Week 3 - Solving second order ODEs

https://drive.google.com/file/d/1HVmJVJqO45UQNiMaivDt0V8QZ8SiRdzZ/view?usp=sharing

%Function to be called in the main program
function [dtheta_dt]=ode_func(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
%main program to solve second order ODE
clear all
close all
clc
%Inputs
%b being damping co-efficient
%g being accelaration due to gravity
%l being length of the thread
%m being the mass of bob
b=0.05;
g=9.81;
l=1;
m=1;
%Initial condition
theta_0=[0;3];
%Time points
%Pendulum oscillation time is 20 seconds
t_span=linspace(0,20,400);
%Solving ODE
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
%Plot
figure(1)
plot(t, results(:,1))
hold on
plot(t,results(:,2))
title('Damped harmonic motion of pendulum')
ylabel('Plot')
xlabel('Time in seconds')
figure(2)
ct=1;
for i=1:length(results(:,1))
x_start=0;
y_start=0;
x_end=l*sin(results(i,l));
y_end=-l*cos(results(i,l));
%Plotting the animation
plot([-1 1],[0 0])
hold on
plot([x_start x_end],[y_start y_end],'linewidth',3)
hold on
plot(x_end,y_end,'o','markersize',20,'color','r')
title('animation of transient behaviour of Simple pendulum')
axis([-1.5 1.5 -1.5 0.5])
pause(0.03)
hold off
M(ct)=getframe(gcf);
ct=ct+1;
end
%Animation of the pendulum
movie(M)
videofile=VideoWriter('Pendulummotion.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)