Week 3 - Solving second order ODEs

%To solve the two order ODE

%This is very important inorder to understand Vehicle dynamics like
%suspensions

clear all
close all
clc

%Inputs
b=0.05; %Damping co-efficient
g=9.81; %accleration due to gravity m/s^2
m=1; %Mass kg
l=1; %Length of the pendulum

%Since we are solving ODE, we know that there are two 1st order ODE which
%are velocity and displacement

%Initial Condition
theta_0=[0;3]; %1st value is the initial angular displacement (rad) and the 2nd value is the initial angular velocity (rad/s)

%Now we need to integrate the ODE system to obtain the solution. So we need
%to define the Time.

%Time points
t_span = linspace(0,50,500); %Though you've asked for 20 seconds, I've considered it for 50 Seconds, to clearly witness that the angular displacement decrease along with the time and also the velocity reduces. This can be witnessed in the Animation

%Solving ODE
%We make use of the inbuilt solver which is 'ODE45', its syntax is
%ode45(@(odefun, timespan).
%t_span is basically the time array which tells us how far the integration
%needs to be done.
%theta_0 is the initial condition corresponding to the starting time.
%ode45(@(t, theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);
%There are many other constants but t and theta are the dependent and
%independent variables of our ODE system
[t, theta_d_v] = ode45(@(t, theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);

AD = theta_d_v(:,1); %Seperating and assigning the angular displacement values to a variable
AV = theta_d_v(:,2); %Seperating and assigning the angular velocity values to a variable

%Giving counter
ct=1;

for i = 1:length(AD) %We consider only this because we need only the angular displacement to define the angle made by the pendulum
W=AD(i);

x0 = 0; %Initial position of the Pendulum string, which is attached to the support
y0 = 0;

% [G] = l.*a; %This gives the arc angles covered by the ball.
x1 = l*cos(W); %Not taking cosd or sind because the angular displacement is already given in radians.
y1 = l*sin(W); %This gives the position of the ball travelled

%Plotting
plot([y0 y1], -[x0 x1], 'linewidth', 2, 'color', 'b') %Have purposefully taken y 1st and -x 2nd in order to reverse and orient the axis to our view
hold on
plot(y1,-x1,'-o','MarkerSize',10,'MarkerFaceColor','r') %There are multiple markers but I've chosen '-o' type and have increased it's size and colored the marker face.
hold on
plot([-1 1], [0 0], 'linewidth', 5, 'color', 'k') %This can be considered as a wall or the support holding the Pendulum
hold off %This is given in order to remove any traces left by 'hold on' in the previous lines
axis([-2 2 -2 2])
pause(0.03)
M(ct)=getframe(gcf);
ct=ct+1;

end

movie(M)
Videofile = VideoWriter('Pendulum_Animation.avi', 'uncompressed AVI');
open(Videofile) %This is the 1st operation. to open the video file
writeVideo(Videofile,M) %We are going to write the frames(Which is stored in M) into this video file
close(Videofile)

%Plotting
figure
plot(t,theta_d_v(:,1))
hold on
plot(t,theta_d_v(:,2))
xlabel('Time')
ylabel('Angular Displacement and Angular velocity')

%Results, we can clearly see the motion of the pendulum decrease along the
%time. AS, we know that the Pendulum swing time depends on the length
%and accl. due to gravity, which inturn depends on Mass. It also depends
%inversely on the damping coefficient.

%Errors faced during the program
%1. I tried to theta_d_v in the for loop, for which I got some uneven
%motion. It is due to two different arrays of disp and vel.
%2. I used cosd and sind, later relaising that the angles are given in
%radians
%3. The trace made by pendulum were seen clearly and had to use hold off to
%remove any traces made by the swing
%4. This is not an error but an orientation problem. In plot I used x 1st
%and later y. But the pendulum was drawn on x axis line. So had to shift
%the axis to make it look properly like a pendulum.
%5. Initially I used x marker, then later in help section found different
%kinds of markers and how to colour the outerline and the innerspace.

%6. Was claculating the arc lenth to know the ball motion. But I understood it

%is just the length covered by the ball and not the position of the ball.

%Function file

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

%we are defining and assigning the thetha's of the array
theta1=theta(1)
theta2=theta(2)

%The ode function which we have created
dtheta1_dt=theta2;
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);

%Finally assigning the above value in theta array, which is having 2
%rows(; inbetween)
[dtheta_dt] = [dtheta1_dt; dtheta2_dt];

end

https://youtu.be/tPclsb_yc3o - Pendulum Animation

https://youtu.be/9qWyyk1JQ-I - Pendulum Animation with increased Damping