Week 3 - Solving second order ODEs

%animation link

https://drive.google.com/file/d/1Ym_nKDo4ONAhgdGepzztyjHaR1UxLmND/view?usp=sharing

Step 1: commands to clear previous data

clear all %to clear workspace
close all %to close plots
clc %clear command window

Step 2: providing the inputs to solve the ODE

b= 0.05; %damping coefficient
g= 9.81; %gravitational constant(m/s^2)
l= 1; %length of the pendulum(m)
m= 0.1; %mass of the body(kg)

% initial condition
theta_0= [0;3];

% time points
t_span= linspace(0,20,1000);

Step 3: to solve the ODE by calling the ode_func

%solve ODE by the function call ode_func
[t, results]= ode45(@(t,theta) ode_func(t, theta,b,g,l,m),t_span, theta_0);

Step 4: to plot angilar displacement and angular velocity

%plot showing the variation of angular displacement and angular velocity
%with respect to time.
figure(1)
hold on
plot(t,results(:,1),'o','MarkerSize',1) %to plot angular displacement
plot(t,results(:,2),'o','MarkerSize',1) %to plot angular velocity
ylabel('parameter')
xlabel('time(s)')
legend('angular displacement((rad)/s)','angular velocity(rad/s^2)')
grid on

Step 5: to create the animation of the motion of the simple pendulum

%to create the movie
ct=1;
%for loop to inculcate t_span
for i= 1:length(t_span)
x0=0; %initial x-coordinate of the pendulum
y0=0; %initial y-coordinate of the pendulum

x1= sin(results(i,1))*l; %end x-coordinate of the pendulum
y1= -cos(results(i,1))*l; %end y-coordinate of the pendulum
hold off
%plot to create animation
figure(2)
plot([-1.5 1.5],[0 0],'linewidth',3,'color','b') %support for the pendulum
hold on
axis([-1.5 1.5 -3 0]) %axis of the pendulum
plot([x0 x1],[y0 y1],'linewidth',3,'color','b') %initial and final coordinates of the string
plot(x1, y1, 'o', 'MarkerSize',20,'MarkerFaceColor','g') %to define the bob position

plot(x0, y0, 'o', 'MarkerSize',10,'MarkerFaceColor','c') %to define the support of the string
s=['time=',num2str(t(i)),'s'] %to plot time in the animation
text(0,-2,s) %shoe position of the time on axis
M(ct)=getframe(gcf); %array to create movie
ct=ct+1;
end

%movie(M)
%to create pendulum animation
videofile = VideoWriter('ode_main.avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

program to define function ode_func

%to define the function for motion of simple pendulum
function [dtheta_dt] = ode_func(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