Week 4.1 - Solving second order ODEs
CODE:
Main
close all
clear variables
clc
b = 0.05;
g = 9.81;
l = 1;
m = 1;
%intial conditions
theta_0 = [0;5];
%time points
t_span = linspace(0,10,150);
%solve ODE
[t, results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span, theta_0);
%ploting
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('plot')
ylabel('time')
legend('plot vs time')
figure;
ct=1;
for i = 1:length(t)
% Pendulum coordinates
% Using the results(:,1) for angular displacement of pendulum
x = l * sin(results(i, 1));
y = -l * cos(results(i, 1));
% Plot pendulum
plot([0, x], [0, y],'color', 'r', 'linewidth', 2); % Rod
hold on;
% pendulum ball
plot(x, y, 'bo', 'markersize', 8, 'markerfacecolor', 'k');
% Pendulum pivot
plot(0, 0, 'ro', 'markersize', 4, 'markerfacecolor', 'b');
axis equal;
axis([-1.5 1.5 -1.5 0.3]);
grid on;
title('Pendulum Motion');
pause(0.03);
hold off;
M(ct) = getframe(gcf);
ct=ct+1;
end
movie(M)
videofile = VideoWriter('Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Function
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
CODE:
Function
Animation movie link:
https://drive.google.com/file/d/1cauBYBjDwp3V2nqmwZ-7drzsRMjQmIvU/view?usp=sharing
also attached drive link in attach link