% Compare Euler and Euler-Cromer

clear
g=9.8;
L=0.5;
m=0.005;
per=2*pi*sqrt(L/g);

a=0;
%b=5*per;
b=8;
N=500;
h=(b-a)/N;

% Initial Condition
t(1)=0;
theta(1)=pi/6;
omega(1)=0;

% Euler's Method
for i=2:N+1
    omega(i)=omega(i-1)-g/L*h*sin(theta(i-1));
    theta(i)=theta(i-1)+h*omega(i-1);
    t(i)=t(i-1)+h;
end

figure(1)
plot(t,omega,'k-')
xlabel('t (s)')
ylabel('\theta (\pi rad)')
set(gca,'YTick',-3*pi/2:pi/2:3*pi/2)
set(gca,'YTickLabel',{'-3/2','-1','-1/2','0','1/2','1','3/2'})
hold on

figure(2)
E1 = 1/2*m*L^2*omega.^2+m*g*L*(1-cos(theta)); 
plot(t,E1,'k-')
xlabel('t (s)')
ylabel('E (J)')
hold on

% Initial Condition
t2(1)=0;
theta2(1)=pi/6;
omega2(1)=0;

% Euler-Cromer Method
for i=2:N+1
    omega2(i)=omega2(i-1)-g/L*h*sin(theta2(i-1));
    theta2(i)=theta2(i-1)+h*omega2(i);
    t2(i)=t2(i-1)+h;
end

figure(1)
plot(t2,omega2,'k--')
legend('Euler Method','Euler-Cromer','location','southwest')
axis([0,8,-2*pi,2*pi])
hold off

figure(2)
E2 = 1/2*m*L^2*omega2.^2+m*g*L*(1-cos(theta2)); 
plot(t2,E2,'k--')
xlabel('t (s)')
ylabel('E (J)')
axis([0,8,0,0.03])
legend('Euler Method','Euler-Cromer','location','northwest')
hold off