clc;
close all;
clear all;
theta0 = pi/3; g = 9.81; L = 1;
ti = 0; tf = 20; dt = 0.005;
nt = (tf-ti)/dt;
t = linspace(ti,tf,nt);
theta = zeros(1,nt);
theta(1) = theta0;
omega = zeros(1,nt);
alpha = zeros(1,nt);
Etotal = zeros(1,nt);
H = L - L*cos(theta);
f = 'Time'; z = 'Omega'; h = 'Theta'; p = 'Alpha'; q = 'Energy';
fprintf('Time \t Omega \t Theta \t Alpha \t Energy\n',f, z, h, p, q);
for k = 1:nt-1
omega(k+1) = (-g/L)*sin(theta(k))*dt + omega(k);
theta(k+1) = omega(k)*dt + theta(k);
alpha(k+1) = (omega(k+1) - omega(k))/dt;
Etotal = g*H + (1/2)*(L*omega).^2;
end
fprintf('%5.2f \t %5.2f \t %5.2f \t %5.2f \t %5.2f\n',[t,omega(k),theta(k),alpha(k),Etotal]');