Week 3 - Solving second order ODEs

% Program for solving second order Ordinary Differential Equations ODEs

close all;
clear all;
clc;

%{
ODEs are extremely useful for solving many Engineering Applications

An example of a second ODE is a simple pendulum oscillating at various
angles theta. Solving the equation for theta for different time interval
determines the velocity at which the pendulum oscillates

The pendulum concept can be extended to the mass, spring and damper system
which find applications in mechanical control system design, vehicle
suspension system design, earthquake modelling and in design of tall
buildings
%}

figure(1)
Pendulum = imread("penphoto.JPG");
image(Pendulum);
%{
The pendulum equation is represented as a second order equation as
follows:
dtheta2_dt2 + (b/m)*dtheta_dt + (g/l)*sin(theta) = 0 where

b = damping coefficient
m - mass of pendulum bulb
l - length of string
g - acceleration due to gravity
theta - angular displacement of bulb

Using the concept of suspended variables, the pendulum equation is reduced
to two first order equations as follows

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
The pendulum equation then becomes:

d_dt(dtheta1_dt) + (b/m)*dtheta1_dt + (g/l)*sin(theta1) = 0
dtheta2_dt + (b/m)*theta2 + (g/l)*sin(theta1) = 0

%}
b = 0.5;
g = 9.81;
l = 1;
m = 0.1;

%{
To use ode45, one need to define a function and then use the concept of function
handles in MATLAB to implement ode45. Please, refer to the MATLAB file Assignweek3ode45

The ode45 takes an initia condition of a two row vector of the variables of interest
and returns a two column vector of the variables of the initial conditions
%}

% Initial conditions
theta_0 = [0;5]; % Initial angular displacement is zero and velocity is 5 rad/sec

% time stamp
t1 = linspace(0, 10, 500);

%{
Solve ODE using ode45.
To solve a second order ODE, after writing the ODE as a matrix system of equations by reducing the system from
second order system to two first order couppled equations by using a method of suspended variable
ode45 which is a numerical integrator using Runke Kutta 4th order scheme to integrate the vector field

The vector is dtheta_dt = [dtheta1_dt; dtheta2_dt] defined in a function.
Refer to MATLAB file Assignweek3ode45.

The ode solver returns a column vector of y where y(:1) is displacement and y(:2) is velocity

%}
[t, y] = ode45(@(t, theta) Assignweek3ode45(t, theta, b, g, l, m), t1, theta_0);

figure(2)
plot(t, y(:,1)) % Plot of angular displacement with time
hold on
plot(t, y(:,2)) % Plot of angular velocity with time
ylabel('Angular Velocity & Angular Displacement')
xlabel('time')
title('Displacement-Velocity Time Series')
legend('Displacement', 'Velocity', 'Location','NorthEast')
grid on
set(gca,'FontSize', 12, 'LineWidth', 1)

%{
Pendulum Animation
To animate the pendulum, we have to plot the angular displacement vector returned
by ode45 which is y(:,1).
The initial position is (0,0) at the fixed end of the pendulum.
From trigonometry, the other position is l*sin(theta) and l*cos(theta) in
the negative y-direction where l is the length of pendulum string. Please,
refer to the file penphoto.jpg
ode45 returned a 500x2 vector where y(500,1) is the angular displacement vector.
Hence, y(1,1) is the first angle of theta and the corresponding velocity is y(1,2).
%}

Assignment = pwd ; % current directory
Pendulum = 'Pendulum' ; % Folder name to save all plots
folder = mkdir([Assignment,filesep,Pendulum]);
Assignment = [Assignment,filesep,Pendulum];
count=3; % Initialize count from 3 for the purpose of figures 1 & 2
%{
Likewise, the forloop need to start from 3, otherwise, the video will not
created due to uninitialized frames.
%}
for i = 3:length(y([1:100,1]))
x0 = 0; % initial position in x direction at the fixed end.
y0 = 0; % initial position in y direction at the fixed end
x1 = l*sin(y(i,1)); % first angular position in the x-direction returned by ode45
y1 = -l*cos(y(i,1));% first angular position in the negative y-direction returned by ode45

figure(count)
plot([-1,1],[0,0],'linewidth',3,'color',[rand rand rand]);
line([x0 x1],[y0 y1],'linewidth',2,'color',[rand rand rand]);
hold on;

plot(x1,y1,'o','markers',20,'markerfacecolor',[rand rand rand]);
grid on ;
axis([-1.5 1.5 -1.5 1.5]);
pause(0.1);
hold off;
count = count + 1;
M(count-3) = getframe(gca);
temp=[Assignment,filesep,'fig',num2str(count),'.jpeg'];
saveas(gca,temp);
end

% Video Creation
% The Videowriter command create a video file from an array or MATLAB Movie
movie(M)
v = VideoWriter('pendulum.avi', 'uncompressed AVI');
open(v)
writeVideo(v,M)
close(v)

%{
Google Link to video
https://drive.google.com/file/d/1GluzOd9dJJBDGUT5ZhKaKOFolwiGnM3a/view?usp=sharing

%}