%
%  duff_per.m -- djm -- 17 jan 2007
%
%  - period
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - -

%  parameters
global epsi

%  ODE, timing & accuracy (events on)
Xode = 'Xduff';  T = 8*pi;  tol = 1e-6;
options = odeset('reltol',tol*1e-3,'abstol',tol*1e-6, ...
		'maxstep',T/5,'stats','off','events','on');

%  loop of ODE solves (A=1)  
period = [];  

%for epsi = 0:0.005:0.2
for epsi = 0:0.02:0.5
    %  correct IV by O(eps) term
    IV = 1 + epsi/32;
    [t y] = ode45(Xode,[0 T],[IV 0],options);
    
    period = [period ; epsi t(end)];
    disp(['epsilon = ' num2str(epsi,3) ' ; period = ' num2str(t(end))])
end

%  last solution plotted here
figure(3);  clf;  hold on
plot(t,y(:,1))
omega = sqrt(1 + epsi*(3/4));
plot(t,cos(omega*t) + (epsi/32)*cos(3*omega*t),'r')
axis equal;  axis([0 2*pi -1.2 1.2])

title('\bf duffing solutions','fontsize',12)
xlabel('\bf t','fontsize',12)
ylabel('\bf y_{comp}(t) & y_{thy}(t)','fontsize',12)

%  final plot
figure(4);  clf;  hold on
plot(period(:,1),period(:,2),'.','markersize',10)
plot(0,2*pi,'r.')
theory = 2*pi./sqrt(1 + (3/4)*period(:,1));
plot(period(:,1),theory,'kx','markersize',10)

title('\bf period vs \epsilon','fontsize',12)
xlabel('\bf \epsilon','fontsize',12)
ylabel('\bf period','fontsize',12)

%  log-log error plot
figure(5);  clf;  hold on
plot(log(period(:,1)),log(abs(period(:,2)-theory)),'b.','markersize',10)
plot([-4.5 -3],[-11 -8],'k--','linewidth',2)
axis([-5.5 -1.5 -14 -5])

title('\bf log-log error plot','fontsize',12)
ylabel('\bf log |period_{comp}-period_{thy}|','fontsize',12)
xlabel('\bf log |\epsilon|','fontsize',12)
text(-4.7,-9,'\bf slope = 2','fontsize',14)