Week 4.1 - Solving second order ODEs

https://drive.google.com/file/d/1J3BnGoUm0o-joDnuvQDoaQltZbuvokzZ/view?usp=drive_linkAim :

To simulate the pendulum motion

Code :

clear all
close all
clc


%% Main Code
L = 1;      % Length of the pendulum (m)
m = 1;      % Mass of the ball (kg)
b = 0.05;   % Damping coefficient
g = 9.81;   % Gravity (m/s^2)


% Initial conditions or inputs
 % Time span
t = linspace(0,50,500); % Simulation time span (s)
x0=[0 3];

% Define the ODE function
[t,result]=ode45(@(t,theta) pendulum_2nd_order(t,theta),t,x0);



% Plot for angular velocity and angular displacement

plot(t,result(:,1),'color','g','linewidth',1);
xlabel('Time (s)');
ylabel('Angular displacement(rad) & Angular Velocity(rad/s)');
hold on;
plot(t,result(:,2),'color','r','linewidth',1);
legend('time vs displacement','time vs velocity')
grid on;



%Animation

X0=0;
Y0=0;
ct=1;%counter represents each and every frame
for i=1:length(t)
    X1=L*sin(result(i,1));
    Y1=-L*cos(result(i,1));
   
figure(2)
plot([-1 1],[0 0],'linewidth',6,'color','b');
axis([-2 2 -2 2]);
line([X0 X1], [Y0 Y1],'linewidth',2,'color','r');
hold on;
plot(X0,Y0,'o','markersize',2,'markerfacecolor','g');
plot(X1,Y1,'o','markersize',15,'markerfacecolor','b');
title('Damping of simple pendulum');
legend('Support rod','length of pendulum','fixed point','Bob of pendulum');
grid on;
hold off;
M(ct)=getframe(gcf);
    ct=ct+1;
end

movie(M)
videofile=VideoWriter('simplependulum.avi','uncompressed AVI');
open(videofile)
writerVideo(videofile,M)
close(videofile)


plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')

function code:

%% Function code
function [dtheta_dt]=pendulum_2nd_order(t,theta)
% Parameters for damping co-efficient
L = 1;      % Length of the pendulum (m)
m = 1;      % Mass of the ball (kg)
b = 0.05;   % Damping coefficient
g = 9.81;   % Gravity (m/s^2)

% Calculation
theta1 = theta(1);
theta2 = theta(2);
%Angular dispalcement
dtheta1_dt=theta2;

dtheta2_dt=(-b/m)*theta2-(g/L)*sin(theta1);
%Angular Velocity
dtheta_dt = [dtheta1_dt;dtheta2_dt]