Week 3 - Solving second order ODEs
clear all
close all
clc
b = 0.05;
g = 9.81;
l = 1;
m = 1;
%initial condition
theta_0 = [0;5];
%time points
t_span = linspace(0,10,500);
%solve ODE
[t,r] = ode45(@(t,theta)odefunc(t,theta,b,g,l,m),t_span,theta_0);
plot(t,r(:,1))
hold on
plot(t,r(:,2))
xlabel('time')
ylabel('plot')
%animation pendulum
x0 = 0;
y0 = 0;
displacement_values = r(:,1);
ct = 1;
for i = 1:length(displacement_values)
x1 = 1*sin(displacement_values(i));
y1 = -1*cos(displacement_values(i));
figure(2)
plot([-3 3],[0 0])
hold on
plot([x0 x1],[y0 y1])
plot(x1,y1,'.','MarkerSize',40)
axis([-3 3 -3 3])
M(ct) = getframe(gcf);
ct = ct+1;
hold off
end
movie(M)
videofile=VideoWriter('ode_main.AVI','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
function [dtheta_dt] = odefunc(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