Week 3 - Solving second order ODEs

%Program to solve 2nd order ODEs

clear all
close all
clc

%Inputs
L=1;
m=1;
b=0.05;
g=9.81;

%Initial Condition
theta_0 = [0;3];

%Time points
t_span = linspace(0,20,500);

%Solving ODE
[t , results] = ode45(@(t,theta) ode_fun(t,theta,b,g,L,m), t_span,theta_0);
figure(1)
hold on
plot(t,results(:,1));
plot(t,results(:,2));
xlabel('Time')
ylabel('Displacement/Velocity')
hold off
legend('Angular Displacement/Time','Angular Velocity/Time')
displacement = results(:,1);

%Counter
ct = 1;

%For Loop
for i=1:length(displacement)

x0 = 0;
y0 = 0;
pos = displacement(i);

x1 = L*sin(pos);
y1 = -L*cos(pos);

%Plotting
plot([x0 x1],[y0 y1],x1,y1,'g.','markersize',20,'linewidth',3)
axis([-1.5 1.5 -1.5 1.5])
xlabel('Horizontal Displacement')
ylabel('Vertical Displacement')
pause(0.03)
A(ct) = getframe(gcf);
ct = ct+1;
end
%Movie Generation
movie(A)
videofile=VideoWriter('Second_order_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,A)
close(videofile)

%ODE function
function [detheta_dt] = ode_fun(t,theta,b,g,L,m)
theta1 = theta(1);
theta2 = theta(2);
detheta1_dt = theta2;
detheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1);
[detheta_dt] = [detheta1_dt ; detheta2_dt];
end

%Explanation
%The program starts with clear all (clears the workspace), close all
%(closes the figure windows) and clc (clears the command window). Then we
%define the input parameters L,m,b,g from line 8-11. Then we form an array
%of initial conditions (displacement and velocity) in store in theta_0. We
%have to plot the values between 0 and 20 sec so we use linspace command
%to generate 500 values between 0 and 20. Now we broke the 2nd order
%equation into two first order equations and write a funtion for it. The
%function name is ode_fun and use t,theta,l,m,b,g as inputs. Then we use
%the function in the main program in the ode45 command. We declare that
%the main inputs will be t,theta and the output will be stored in variables t
% and results. The results have 2 columns one is displacement and other is
%the velocity. We then declare a counter ct and start a for loop. In for
%loop we store value of each displacement in variable called pos and use
%this variable in our calculation. Then as similar to the 2R-Robotic arm
%kinematics program we generate a movie called Second_order_pendulum.avi

Link for the file:- https://drive.google.com/file/d/10FOyX1Ofcb1MkPbQatCEe7wIOxIfTtJS/view?usp=sharing