twowell.m

双势井duffing方程分叉演示程序。演示Duffing方程在不同系数选择下单周期

plot3(X,Y,Z,'r','linewidth',1) % trajectories in phase space
title('An Actual P-1 Trajectory in the Phase Space','color','y')
pause
title('Press ENTER to see the projection of trajectory onto the PS',...
    'color','b')
pause
Xp=(5+y(:,1));
Yp=zeros(size(X))-10*eps;
% projection of trajectories onto poincare section
plot3(Xp,Yp,Z,'m','linewidth',1) 
title('Projection of the P-1 onto the Poincare Plane','color','y')
pause
title('Press ENTER for the Intersection Point','color','b')
pause
% poincare section
plot3(X(1:36:end),Y(1:36:end),Z(1:36:end),'k.','markersize',12)
title('Press ENTER to see the motion on the trajectory','color','b')
pause
title('Evolution of the Trajectory in the Phase Space','color','y')
plotting([1.19718071928 0.08451532560 0],1000)
title('For larger amplitudes P-n orbits are possible: press ENTER',...
    'color','b')
pause

%% 3.2 Now look at P-10 motion

clf reset, cla reset, set(gcf,'color','k')
f = .259;
[t,y]=ode45(@duffing,[0:2*pi/omega/36:400*pi/omega],...
   [-1.1712132945 -0.3137585346 0]);
subplot(2,1,1)
plot(y(1:700,1),'color','r','linewidth',1), grid
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on')
xlabel('Time'),ylabel('Amplitude')
title('P-10 motion observed for f = 0.259','color','y')
subplot(2,1,2)
psd(y(:,1)-mean(y(:,1)),1024)
set(get(gca,'children'),'color','r','linewidth',1)
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on')
axis tight
pause(2)
text(0.5,40,'Press ENTER for Corresponding Poincere Section',...
    'HorizontalAlignment','center','color','b','fontsize',14)
pause

%% 3.2.1 Now look at the corresponding poincare sections
clf reset,set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on')
plot(y(1:36:end,1),y(1:36:end,2),'ro')
axis([-2 2 -2 2]), xlabel('u'), ylabel('v')
title(['Poincare Section of P-10 Orbits, '...
    '\omega = 1, f = 0.259'],'color','y'), hold on
pause
title('Press ENTER to see P-10 Projection onto the PS','color','b')
pause
plot(y(:,1),y(:,2),'y')
title('Press ENTER to visualize the P-10 in the Phase Space','color','b')
pause

%% 3.2.2 Now visualize all this in the Phase Space
clf reset,axis off,hold on
set(gca,'color','k'),set(gcf,'color','k')
% plot time line C^1
plot3(5*cos(y(:,3)),5*sin(y(:,3)),zeros(size(y(:,3))),'y','linewidth',1)
view([23,73])
axis([-8.0000   8.0000   -6.0000    6.0000   -1    1])
title('Time Lives on a Compact Manifold (Circle)','color','y')
pause
title('Press ENTER to Fix a Poincare Section in Time','color','b')
pause
fill3([3 7 7 3 3],[0 0 0 0 0],[-1 -1 1 1 -1],'c',...
    'linestyle','none','facealpha',.5) % poincare plane
title('Poincare Plane for Stroboscopic Sampling','color','y')
pause
title('Press ENTER to view one of the trajectories','color','b')
pause
X=(5+y(:,1)).*cos(y(:,3));
Y=(5+y(:,1)).*sin(y(:,3));
Z=y(:,2);
plot3(X,Y,Z,'r','linewidth',1) % trajectories in phase space
title('An Actual P-10 Trajectory in the Phase Space','color','y')
pause
title('Press ENTER to see the projection of trajectory onto the PS',...
    'color','b')
pause
Xp=(5+y(:,1));
Yp=zeros(size(X))-10*eps;
% projection of trajectories onto poincare section
plot3(Xp,Yp,Z,'m','linewidth',1) 
title('Projection of the P-10 onto the Poincare Plane','color','y')
pause
title('Press ENTER for the Intersection Points','color','b')
pause
% poincare section
plot3(X(1:36:end),Y(1:36:end),Z(1:36:end),'k.','markersize',12)
title('Press ENTER to see the motion on the trajectory','color','b')
pause
title('Evolution of the Trajectory in the Phase Space','color','y')
plotting([-1.1712132945 -0.3137585346 0],2000)
title('For larger amplitudes CHAOTIC orbits are possible: press ENTER',...
    'color','b')
pause

%% 3.3 Now look at chaotic motion

clf reset, cla reset, set(gcf,'color','k')
f = .3;
[t,y]=ode45(@duffing,[0:2*pi/omega/36:2000*pi/omega],...
   [-0.67079530230 -0.36933661878 0]);
subplot(2,1,1)
plot(y(1:2000,1),'color','r','linewidth',1), grid
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on')
xlabel('Time'),ylabel('Amplitude')
title('Chaotic motion observed for f = 0.3','color','y')
subplot(2,1,2)
psd(y(:,1)-mean(y(:,1)),1024)
set(get(gca,'children'),'color','r','linewidth',1)
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on')
axis tight
pause(2)
text(0.5,40,'Press ENTER for Corresponding Poincere Section',...
    'HorizontalAlignment','center','color','b','fontsize',14)
pause

%% 3.3.1 Now look at the corresponding poincare sections
clf reset,set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on')
plot(y(1:36:end,1),y(1:36:end,2),'r.','markersize',4)
axis([-2 2 -2 2]), xlabel('u'), ylabel('v')
title(['Poincare Section of Chaotic Orbit (Strange Attractor), '...
    '\omega = 1, f = 0.3'],'color','y'), hold on
pause
title('Press ENTER to See Trajectory Projection onto the PS','color','b')
pause
plot(y(1:2000,1),y(1:2000,2),'y','linewidth',0.25)
title('Press ENTER to visualize this Chaos in the Phase Space','color','b')
pause

%% 3.3.2 Now visualize all this in the Phase Space
clf reset,axis off,hold on
set(gca,'color','k'),set(gcf,'color','k')
% plot time line C^1
plot3(5*cos(y(:,3)),5*sin(y(:,3)),zeros(size(y(:,3))),'y','linewidth',1)
view([23,73])
axis([-8.0000   8.0000   -6.0000    6.0000   -1    1])
title('Time Lives on a Compact Manifold (Circle)','color','y')
pause
title('Press ENTER to Fix a Poincare Section in Time','color','b')
pause
fill3([3 7 7 3 3],[0 0 0 0 0],[-1 -1 1 1 -1],'c',...
    'linestyle','none','facealpha',.5) % poincare plane
title('Poincare Plane for Stroboscopic Sampling','color','y')
pause
title('Press ENTER to view one of the trajectories','color','b')
pause
X=(5+y(:,1)).*cos(y(:,3));
Y=(5+y(:,1)).*sin(y(:,3));
Z=y(:,2);
plot3(X(1:2000),Y(1:2000),Z(1:2000),'r','linewidth',1) % trajectories in phase space
title('An Actual Chaotic Trajectory in the Phase Space','color','y')
pause
title('Press ENTER to see the projection of trajectory onto the PS',...
    'color','b')
pause
Xp=(5+y(:,1));
Yp=zeros(size(X))-10*eps;
% projection of trajectories onto poincare section
plot3(Xp(1:2000),Yp(1:2000),Z(1:2000),'m','linewidth',1) 
title('Projection of the Trajectory onto the Poincare Plane','color','y')
pause
title('Press ENTER for the Intersection Points','color','b')
pause
% poincare section
plot3(X(1:36:end),Y(1:36:end),Z(1:36:end),'k.','markersize',12)
title('Press ENTER to see the motion on the trajectory','color','b')
pause
title('Evolution of the Trajectory in the Phase Space','color','y')
plotting([-1.1712132945 -0.3137585346 0],2000)
title('To Look at the Evolution of Strange Attractor press ENTER',...
    'color','b')
[t,y]=ode45(@duffing,[0:2*pi/omega/360:2000*pi/omega],...
   [-0.67079530230 -0.36933661878 0]);
pause
clf reset,
set(gcf,'color','k')
for i=[1:360 1:360],
    plot(y(i:360:end,1),y(i:360:end,2),'y.','markersize',4)
    set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on'), 
    axis([-1.6 1.6 -1 1])
    xlabel('u'), ylabel('v'), 
    title('Evolution of STRANGE ATTRACTOR','color','y')
    pause(0.02)
end
title('This is it: Go Home!','color','b')

%% ------------------------------------------------------------------------
function dy=duffing(t,y)
% duffing equation for integration
% revised 10/04/2006 by DC
global alpha beta gamma f omega
dy=[y(2);-alpha*y(1)-beta*y(1)^3-gamma*y(2)+f*cos(omega*t);1];

%% ------------------------------------------------------------------------
function plotting(y,time)
global alpha beta gamma f omega axis
% y -- initial condition
% t -- lenght of iteration (300~3000)

t=0;
head = line( ...
        'color','y', ...
        'Marker','.', ...
        'markersize',18, ...
        'erase','xor', ...
        'xdata',y(1),'ydata',y(2),'zdata',y(3));
body = line( ...
        'color','y', ...
        'LineStyle','-', ...
        'erase','none', ...
        'xdata',[],'ydata',[],'zdata',[]);
tail=line( ...
        'color','r', ...
        'LineStyle','-', ...
        'erase','none', ...
        'xdata',[],'ydata',[],'zdata',[]);
proj = line( ...
        'color','y', ...
        'Marker','.', ...
        'markersize',15, ...
        'erase','xor', ...
        'xdata',y(1),'ydata',y(2),'zdata',y(3));

% Save L steps and plot like a comet tail.
L = 15;
Y = ones(L,1)*[(5+y(1)).*cos(omega*y(3)),(5+y(1)).*sin(omega*y(3)),y(2)];

for i=1:time,
    [t,y] = ode45(@duffing,[t t+0.05],y);
    y = y(end,:);
    t = t(end);
    yn = [(5+y(1)).*cos(omega*y(3)),(5+y(1)).*sin(omega*y(3)),y(2)];
    yp = [5+y(1),0,y(2)];
    Y = [yn; Y(1:L-1,:)];
    % Update the plot
    set(head,'xdata',Y(1,1),'ydata',Y(1,2),'zdata',Y(1,3))
    set(body,'xdata',Y(1:2,1),'ydata',Y(1:2,2),'zdata',Y(1:2,3))
    set(tail,'xdata',Y(L-1:L,1),'ydata',Y(L-1:L,2),'zdata',Y(L-1:L,3))
    set(proj,'xdata',yp(1,1),'ydata',yp(1,2),'zdata',yp(1,3))
    drawnow;
end

X=(5+y(:,1)).*cos(omega*y(:,3));
Y=(5+y(:,1)).*sin(omega*y(:,3));
Z=y(:,2);