Week 3 - Solving second order ODEs

1) Main program for plotting the graph and creating the animations for the pendulum motion.

%To solve the ODE, that represents the equation of motion of a simple pendulum with damping.

clear all
close all
clc

%% Input required for solving ODE equation of motion of a simple pendulum with damping

%Damping coefficient
b=0.05;

%Acceleration due to gravity (Unit: m/s^2)
g=9.81; 

%Length of the pendulum (Unit: m)
l=1;

%mass of the ball (Unit: kg)
m=1;

%The initial values of angular Displacement= 0 radian and angular Velocity=3 rad/s is given in 2X1 matrix 
theta_0 = [0;3];

%The motion with in the given time i.e from 0 to 20 sec is calculated.
t_span = linspace(0,20,500);

%% The ODE function is solved in the matlab using an inbuilt function called as ODE45
%Syntax of ODE45: [t,y] = ode45(odefun,tspan,y0), where tspan = [t0 tf], integrates the system of differential equations 
%y′=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

%% Plotting the results from solving the ODE function with the given conditions

%The figure is assigned and named using the below command
f = figure('name','Motion of a simple pendulum','numbertitle','off');
%The below-used command is for positioning the figure on the screen and to enlarge the size of the figure for better visibility
f.Position = [100 100 800 700];

%% Subplot for angular Velocity & Displacement (VS) time in the same figure

% The subplot would divide the figure into a specified MXN matrix, where the
% individual elements of matrix denote each division of the figure
% In the below case the figure is divided into 6 equal parts, forming a 3X2 matrix
% This graph is plotted at the 5 division of figure

subplot(3,2,5)
plot(t,results(:,1),'color','b')
hold on
plot(t,results(:,2),'color','r')
hold on

%Provides the title and labels for the subplot
title('Ang-Velocity & Displacement (VS) time')
ylabel('Ang-Velocity & Displacement')
xlabel('Time (sec)')

%% Subplot for angular Velocity vs Position of the pendulum in the same figure

% The subplot would divide the figure into a specified MXN matrix, where the
% individual elements of matrix denote each division of the figure
% In the below case the figure is divided into 6 equal parts, forming a 3X2 matrix
% This graph is plotted at the 6 division of figure

subplot(3,2,6)
plot(results(:,1),results(:,2))
hold on

%Provides the title and labels for the subplot
title('Angular velocity (VS) Position')
ylabel('Angular velocity (rad/s)')
xlabel('Angle (radian)')


%% Subplot for creating the animation for the motion of the pendulum in the same figure

%One end of the pendulum must be fixed in order for correct animation of the motion of the pendulum
x0=0;
y0=0;

%ct variable contains the frame number for creating a movie from the obtained frames
ct=1;

for i = 1:length(results(:,1))
    
    %The values obtained from X1 and Y1 are the positions of the endpoint of the pendulum (or) the center point of the ball with mass m.
    x1 = l*sin(results(i,1));
    %Without the negative sign the pendulum will be oscillating at the top of fixed point which is incorrect
    y1 = -l*cos(results(i,1));
    
    % The subplot would divide the figure into specified MXN matrix, where the
    % individual elements of matrix denote each division of the figure
    % In the below case the figure is divided into 6 equal parts, forming a 3X2 matrix
    % This graph is plotted at the 1,2,3 and 4 divisions of figure

    subplot(3,2,[1,2,3,4])
    plot([-1 1],[0 0],'linewidth',3,'color','black')
    hold on
    
    %The line command is used to represent the connecting rod of the pendulum with the suspended mass
    line([x0 x1],[y0 y1],'linewidth',2,'color','b')
    
    %The marker feature is used for representing the ball at the end of the pendulum
    plot(x1,y1,'-o','markersize',25,'markerfacecolor','r')
    
    %This command is used to display the graph in the figure
    grid on
    
    %Provides the title for the subplot
    title('Pendulum Animation')
    
    %Axis command is used for setting the limits of X and Y-axis
    axis([-1.5 1.5 -1.5 0.5])
    hold off
    
    %Getframe command Captures the axes or figure as movie frame
    M(ct)=getframe(gcf);
    ct=ct+1;
end

%% These commands are used for creating the movie file from each frame 
%This command is used to create a movie file using variable M
movie(M)
%A video named Pendulum_ODE.avi is created and assigned to variable video_2
video_2 = VideoWriter('Pendulum_ODE.avi','Uncompressed AVI');
%Video_2 is opened using this command
open(video_2)
%Every frame in the variable M is written to this video_2
%This command is case sensitive and found that V in writevideo should be in uppercase
writeVideo(video_2,M);
%After writing every frame to this variable, the file is closed
close(video_2)

%% Result

% The animation representing the motion of the pendulum is successfully generated using Matlab

% The Second-order differential equation is reduced into sets of two first-order differential equations.

% After solving the ODE equation we would get different values of angular velocity and angular displacement
% for the given time limit of 0 to 20 secs for 500 iterations

% With the values obtained from solving this ODE function, we can plot different graphs as required 

% As per the obtained values from solving the ODE, the oscillation of the pendulum is getting damped, 
% but not fully within the given time limit within 0 to 20 secs for the given inputs.

% So this oscillation of the pendulum can be completely damped by varying the value of the Damping coefficient

2) Subfunction for defining the ODE:

%% Function for defining the ODE 

% The variables required for solving the ODE function is obtained from the main function

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

% The Second-order diffrential equation is reduced into sets of two first-order differential equations.
theta1 = theta(1)
theta2 = theta(2)

% The dtheta1_dt represents angular displacement
dtheta1_dt = theta2;

% The dtheta2_dt represents angular velocity
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);

% These values of angular displacement and angular velocity
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

The output video link:

https://youtu.be/wC9fasCN5Wg