Week 3 - Solving second order ODEs

Explanation

ODE Function Program

function [dtheta_dt] = ode_function(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

Program to simulate pendulum motion and creating animation

% Program to solve the ODE and simulate pendulum motion
clear all
close all
clc

%Inputs
b = 0.05; %damping coefficient value
g = 9.81; %gravitational value in m/s^2
l = 1; %length of the pendulum in m
m = 1; %mass of the body in kg

theta_0 = [0;3]; %Initial condition for displacement and velocity
t_span = linspace(0,20,400); %Time span points

%Solving ODE
[t,results] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0);

%Plotting
figure(1)
plot(t,results(:,1))
hold on
grid on
plot(t,results(:,2))
ylabel('Plot')
xlabel('Time')
legend('Displacement','Velocity')

figure(2)
ct = 1;
for i = 1:length(results(:,1))
THETA_1 = results(i,1);
%Initial coordinates
x0 = 0;
y0 = 0;

x1 = l*sin(THETA_1);
y1 = -l*cos(THETA_1);

%Plot
plot([-0.5 0.5],[0,0],'color','g')
line([x0 x1],[y0 y1],'color','r')
hold on
plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','b')
grid on
axis([-1.5 1.5 -1.5 1])
hold off
pause(0.0001)
M(ct) = getframe(gcf);
ct = ct + 1;
end
movie(M)
videofile = VideoWriter('Second_ODE_Pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Results

The displacement v/s time and velocity v/s time graphs were plotted.

Attach below

The animation video for the pendulum motion was created and uploaded in youtube.

Link for video:https://www.youtube.com/watch?v=gbOlm5hSx94