Week 3 - Solving second order ODEs

SIMPLE PENDULUM

%inputs
b = 0.05;
g=9.81;
l=1;
m=0.1;

% initial condition
in_ang_disp = 0;
in_ang_v = 3;

theta_0 = [in_ang_disp;in_ang_v];


% time points
t_span = 20; 
fps = 30; 
t_pts = linspace(0,t_span,fps*t_span);


% solve ODE
[t, results] = ode45(@(t,theta)simp_pend_ode_func(l,theta,b,g,m),t_pts,theta_0);
ang_disps = results(:,1); 
ang_vs = results(:,2); 

% Fixed point of the pendulum
x0 = 0; 
y0 = 0; 
ct = 1; 

% Pendulum ball dimensions
r = 0.1; 

% Angular displacement & Angular velocity v/s time plot

figure()
hold on
title("Angular displacement & Angular Velocity v/s time of a simple damped pendulum")
plot(t,ang_disps)
text(3.052,0.3727,'leftarrow Angular Displacement')
plot(t,ang_vs)
text(2.088,2.83,'leftarrow Angular Velocity ')
ylabel('Plots')
xlabel('time (s)')

% The animation loop

figure()
for i = 1:length(t_pts)
    pos = ang_disps(i); 
    
    y = -l*cos(pos); 
    x = -l*sin(pos); 
        
    % Animation plots
    
    plot([x0 x],[y0 y],'linewidth',3,'Marker','O','Markersize',15,'MarkerFaceColor','r','MarkerEdgeColor','r')
       axis([-1.5 1.5 -1.5 1.5])
    
    rectangle('Position',[-1.5 0 3 0.5],'Facecolor',[0.2 0 0],'Edgecolor',[0 0 0])

    pause(t(2)) 
    M(ct) = getframe(gcf);
    ct = ct+1; 
end

% The movie

videofile = VideoWriter('simple pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
 

% Input values

L = 1; % Radius of Pendulum arc in m

m = 1; % Mass of pendulum in kg

b = 0.05; % Damping coefficient

g = 9.81; % Acceleration due to gravity m/s^2

% Conditions

in_ang_disp = 0; % Initial angular displacement

in_ang_v = 3; % Initial angular velocity in radians/s

theta_0 = [in_ang_disp;in_ang_v];

t_span = 20; % Total time of video in seconds

fps = 30; % Frames per second

t_pts = linspace(0,t_span,fps*t_span); % Time points for t_span seconds

% ODE Solver

[t, results] = ode45(@(t,theta) simp_pend_ode_func(L,theta,b,g,m),t_pts,theta_0);

ang_disps = results(:,1); % Angular displacement array

ang_vs = results(:,2); % Angular velocity array

% Fixed point of the pendulum

x0 = 0; % x origin coordinate

y0 = 0; % y origin coordinate

ct = 1; % Counter for counting the total number of frames

% Pendulum ball dimensions

r = 0.1; % Radius in m

% Angular displacement & Angular velocity v/s time plot

figure()

hold on

title("Angular displacement & Angular Velocity v/s time of a simple damped pendulum")

plot(t,ang_disps)

text(3.052,0.3727,'leftarrow Angular Displacement')

plot(t,ang_vs)

text(2.088,2.83,'leftarrow Angular Velocity ')

ylabel('Plots')

xlabel('time (s)')

% The animation loop

figure()

for i = 1:length(t_pts)

pos = ang_disps(i); % Angular displacement of pendulum at a point of time

y = -L*cos(pos); % x coordinate of pendulum at a point of time

x = -L*sin(pos); % y coordinate of pendulum at a point of time

% Animation plots

plot([x0 x],[y0 y],'linewidth',3,'Marker','O','Markersize',15,'MarkerFaceColor','r','MarkerEdgeColor','r') % The pendulum plot

axis([-1.5 1.5 -1.5 1.5]) % Axis of the plot

rectangle('Position',[-1.5 0 3 0.5],'Facecolor',[0.2 0 0],'Edgecolor',[0 0 0])

pause(t(2)) % Determines the time b/w the frames Also, t(2) = 1/fps

M(ct) = getframe(gcf);

ct = ct+1; % Counter gets incremented after the frame is collected

end

% The movie

movie(M);

videofile = VideoWriter('simple_pendulum.avi','Uncompressed AVI'); % Naming and formating the video

open(videofile) % Open the video editor

writeVideo(videofile,M) % Convert the image to frame

close(videofile) % Close the video writer

And it has the calling of ODE solve function

 function [dtheta_dt] = simp_pend_ode_func(l,theta,b,g,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