Week 3 - Solving second order ODEs

%2nd Order ODE

clear all
close all
clc

%b = damping coefficient;
%g = accn due to gravity;
%l = length of the pendulum string;
%m = mass of the pendulum sphere;

b = 0.05;
g = 9.81;
l = 1;
m = 1;

% initial condition
% initial displacement is 0
% initial velocity is 3

theta_0 = [0;3];



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

% solve ODE
[t, results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span, theta_0);
%t, theta are main inputs
%define function with input parameters = ode_func(t,theta,b,g,l,m)
%define time limits = t_span
%defeine initial conditions = t_span


%Plot results ;
%plot partition into quadrant

%plot displacement & velocity curve in lower left quadrant
subplot(5,4,[13,14,17,18]);
plot(t,results(:,1));
hold on;
plot(t,results(:,2));
hold off;
%plot angular displacement & angular velocity curve in lower right quadrant
subplot(5,4,[15,16,19,20]);
plot(results(:,1),results(:,2));
ylabel('Angular velocity');
xlabel('position radians');

ct=1;

%for loop for the pendulum animation
%define length of the loop time
for i=1:length(results(:,1))
    %define initial conditions
    x0=0;
    y0=0;
    
    %define the sweep coordinates
    x1=1*sin(results(i,1));
    y1=-1*cos(results(i,1));
    
    %define plot quadrant
    subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
    plot([1 1],[0 0],'linewidth',2);
    hold on
    line([x0 x1],[y0 y1],'linewidth',2);
    hold on;
    
    %define the sphere for the pendulum
    plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
    
    %display grid
    grid on;
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.0003);
    hold off
    G(ct)=getframe(gcf);
    ct=ct+1;
end
    
%generate a movie
movie(G)

%Create an avi file using all the frames, give a file name and the
%compression details
videofile = VideoWriter('second_order_ode.avi','uncompressed AVI');

%Open the video file
open(videofile)

%write the video file with captured frames
writeVideo(videofile,G)

%close the video file window
close(videofile)  

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

%define theta1 & theta2
theta1 = theta(1)
theta2 = theta(2)

%define first ODE
dtheta1_dt = theta2;

%define second ODE
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);

%Combine the equation into an array of two rows
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end  

https://youtu.be/i0LdzI1kXFk

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