Week 3 - Solving second order ODEs

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

% Program to solve second order ODE
b = 0.05;% in Ns/m
g = 9.81;%in m/s^2
l = 1; % in m
m = 1; % in kg 

% initial condition
theta_0 = [0;3];

% time span of simulation
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);
hold on
plot(t,results(:,1), '-b','linewidth',3,'color', 'r')
plot(t,results (:,2),'-k','linewidth',3,'color','b')
xlabel('Time t in s')
ylabel('Angular displacement')
legend('Angular displacement','Angular velocity in rad')
title('simple pendulum with damping')
grid on

%creating animation of simple pendulum
X0 = 0;
Y0 = 0;
ct=1;
for i = 1:length(t)
    X1 = l*sin(results(i,1));
    Y1 = -l*cos(results(i,1));
    figure(2)
    plot([-1 1],[0  0 ],'linewidth',6,'Color','black');
    axis([-2 2 -2 2]);
    line([X0 X1],[Y0 Y1],'linewidth',2,'color','b');
    hold on
    plot(X0, X0,'O','markersize',20,'MarkerFaceColor','y')
    plot(X1,Y1,'O','Markersize',20,'MarkerFaceColor','r')
    title('Damping of single pendulum');
    legend('Support rod','length of pendulum','Fixed point','bob of pendulum');
    grid on;
    pause(0.03)
    hold off;
    M(ct) = getframe(gcf);
    ct = ct+1;
   
end
% creating a movie for the animation
movie(M);
videofile = VideoWriter('Simple_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile);




% To write a function program which will be used in 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