Week 3 - Solving second order ODEs

Q = theta;

Q = Q₁ ; dQ/dt = dQ₁/dt= Q₂ ;

d²Q/dt² = d²Q₁/dt² = d/dt (dQ₁/dt) = d/dt (Q₂) => d²Q/dt² = d/dt (Q₂);

d²Q/dt² + b/m dQ/dt + g/l sinQ = 0;

substituting Q = Q₁ ; dQ₁/dt= Q₂ ;

d/dt (Q₂) + b/m Q₂ + g/l sinQ₁ = 0;

d/dt |Q₁ | |Q₂ |

|Q₂ | = |- b/m Q₂ - g/l sinQ₁|

I dont have a player in my computer so i cannot see the video file is written or not , therfore I did it with for loop and pause command for the simulation

clear all
close all
clc
%% input parameters
b = 0.05;   %damping coefficient
g = 9.81;   %acceleration due to gravity(m/s^2)
l = 1;          %length of pendulam
m = 0.1;    %mass of the bob
%the simple pendulam motion has two 1st order ODEs one for defining the
%displacement and 2nd for defining the velocity of the pendulam 
%so there two initial conditions need to be defined to solve
%initial conditions at time t =0 ==> displacement  =0 & velocity = 5 m/sec^2 (instantaneous) 
theta_0 = [0;5];    %1st row  in the theta_0 [2x1]matrix represents the initial angular displacement 0, and 2nd row initial angular velocity 5 
% time span
t_span = linspace(0,10,500);%taking 500 sampling points to ingerate for 0 to 10secs
%% solving the 2nd order ODE using inbuilt ODE45
% ode45(@(t,theta) this term calls the inbuilt function by passing the
% variables contained in the ODES to be solved; 2nd part is to pass the
% varibles contained in the two 1st order ODE as indepent var(t)and and
% depandant var(theta)
% ODES  we get them  as a output of the function ode_funand its parameters
% as defined in the term  ode_fun(t, theta, b, g, l, m)
%the inputs to solve the ODEs such as t and the initial conditionst_span, theta_0
[t, results] = ode45(@(t,theta) ode_fun(t, theta, b, g, l, m),t_span, theta_0);% results is a 2 column matrix it contains 500 rows 1st column contains the results of the displacement and second columnof 500 rows contains the velocity of the pendulam
%plot the s and velocity
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Displacement in meters & Velocity in m/sec')
xlabel('time in seconds')
legend('displacement','velocity')
hold off

%%
%creating the simulation of the pendulum motion
X0 = 0; Y0 = 0; % pendulum fixed top position
S = results(: ,1); % displacement values from the first column in all rows of the results matrix.
S_point = 1;

figure(2)
hold on
for i = 1: length(S)
    X1 = 1*sin(S(i));
    Y1 = -1*cos(S(i));
    plot(X1,Y1, 'Marker','.','markersize', 50, 'Color','r')
    plot([X0 X1],[Y0 Y1],'Color','b');
    plot([-1 1], [0 0], 'LineWidth', 2,'Color','k')
    axis([-3 3 -3 3])
    pause(0.05);
end
hold off
 %%




function [dtheta_dt] = ode_fun(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