demo.m

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

    flipud(z2(1:l,1))],[z2(:,2);flipud(z1(1:l,2))],'y','linestyle',...
    'none','facealpha',.5);
fill3([10 10 10 10 10],[5 -5 -5 5 5],[-15 -15 15 15 -15],'y','linestyle',...
    'none','facealpha',.5);
text(33,20,'Press ENTER for \beta < 0 Case','color','y',...
    'horizontalAlignment','center','fontsize',14)
pause

%% 2.4. beta = -1, alpha = -10
cla reset, clf reset, text(1,1,'.'), pause(0.1)
beta = -1; alpha = -10;
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14), axis off
text(0.5,0.8,'Now consider negative \beta case (w.l.o.g., \beta = -1).',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.7,'Thus, for equilibria we again have two distinct cases.',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.5,...
    'Case 1: \alpha < 0. Only one equilibrum (x,dx/dt) = (0,0)',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.3,'Press ENTER to display motion on the Energy Surface',...
    'color','y','horizontalAlignment','center','fontsize',14)
pause

%% 2.4.1. motion on energy surface
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14)
[X,Y] = meshgrid(-3:0.05:3,-9:.1:9);
surface(X,Y,E(alpha,beta,X,Y)+90,'facealpha',.9), view([-55,25])
axis([-3 3 -9 9 0 150]), 
shading interp
xlabel('x'),ylabel('v'),zlabel('First Integral, E(x,v)')
title(['This is stiil the total energy: \beta = -1, \alpha = -10, '... 
    '\gamma = 0.25'],'Color','y')
% -------------------------------------------------------------------------
pause
title('Press ENTER key to see motion on this surface')
pause
yini = [-1 -1 -1 -1 1 1 1 1 -0.9 -0.75 -0.64 -0.5 -0.25 0.9 0.75 0.64 0.5 0.25;...
         3 2.5 2 1.5 -3 -2.5 -2 -1.5 3 3 3 3 3 -3 -3 -3 -3 -3]*3;
for i = 1:size(yini,2),     
    [t,z]=ode45(@duffing,0:0.01:3,[yini(:,i);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot3(z(:,1),z(:,2),E(-10,-1,z(:,1),z(:,2))+90,'w')
    plot(z(:,1),z(:,2),'b')
end
[evecs2,evals] = eig([0,1;-alpha,-gamma]);
for i = 1:2,     
    [t,z]=ode45(@duffing,0:0.005:7,[(-1)^i*0.01*evecs2(:,1);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot3(z(:,1),z(:,2),E(-10,-1,z(:,1),z(:,2))+90,'w')
    plot(z(:,1),z(:,2),'b')
    [t,z]=ode45(@duffing,0:-0.005:-7,[(-1)^i*0.01*evecs2(:,2);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot3(z(:,1),z(:,2),E(-10,-1,z(:,1),z(:,2))+90,'w')
    plot(z(:,1),z(:,2),'b')
end
title('Trajectories are still confined to the energy surface','Color','y')
pause
title('Press ENTER to plot the Phase Portrait')
pause

%% 2.4.2. phase portrait
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14,'Box','on')
[x,y] = meshgrid(-3:0.5:3,-9:1:9);
quiver(x,y,y,dv(alpha,beta,gamma,x,y),1.3,'w','linewidth',1), %vector field
xlabel('position, x'),ylabel('velocity, v'),
title(['Dissipative Phase Portrait: \beta = -1, \alpha = -10,'...
    '\gamma = 0.25'],'color','y')
axis([-3 3 -9 9])
for i = 1:size(yini,2),     
    [t,z]=ode45(@duffing,0:0.01:3,[yini(:,i);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot(z(:,1),z(:,2),'y')
end
[evecs2,evals] = eig([0,1;-alpha,-gamma]);
for i = 1:2,     
    [t,z]=ode45(@duffing,0:0.005:7,[(-1)^i*0.01*evecs2(:,1);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot(z(:,1),z(:,2),'r','linewidth',2)
    [t,z]=ode45(@duffing,0:-0.005:-7,[(-1)^i*0.01*evecs2(:,2);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot(z(:,1),z(:,2),'b','linewidth',2)
end
pause
title('Press ENTER to considere Case 2: \alpha > 0.')
pause

%% 2.5. beta = -1, alpha = 10
beta = -1; alpha = 10;
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',12), axis off
text(0.5,0.7,'Now consider Case 2: \alpha > 0.',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.5,'There are three equilibria: (0,0) and (+/-(\alpha)^{1/2},0).',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.3,'PRESS ENTER TO DISPLAY THE ENERGY SURFACE',...
    'color','y','horizontalAlignment','center','fontsize',14)
pause

%% 2.5.1. first integral
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14)
warning off
[X,Y] = meshgrid(-5:0.05:5,-10:.1:10);
surface(X,Y,E(alpha,beta,X,Y)+50,'facealpha',.9), 
view([-32,32]), shading interp
xlabel('x'),ylabel('v'),zlabel('First Integral, E(x,v)')
title('This is still the total energy: \beta = -1, \alpha = 10','color','y')
% -------------------------------------------------------------------------
pause
title('Press ENTER key to see motion on this surface')
pause
yini = [4.1:0.5:5.1 5.2,-4.1:-0.5:-5.1 -5.2;...
        -10*ones(1,4),10*ones(1,4)];
for i = 1:size(yini,2),     
    [t,z]=ode45(@duffing,0:0.005:2,[yini(:,i);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot3(z(:,1),z(:,2),E(alpha,beta,z(:,1),z(:,2))+50,'w')
    plot(z(:,1),z(:,2),'b')
end
xe = sqrt(-alpha/beta);
[evecs3,evals] = eig([0,1;-alpha-3*beta*xe^2,-gamma]);
for i = 1:2, 
    for j=1:2,     
    [t,z]=ode45(@duffing,0:0.01:25,(-1)^j*[xe;0;0]+[(-1)^i*0.01*evecs3(:,1);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot3(z(:,1),z(:,2),E(alpha,beta,z(:,1),z(:,2))+50,'w')
    plot(z(:,1),z(:,2),'b')
    [t,z]=ode45(@duffing,0:-0.01:-24,(-1)^j*[xe;0;0]+[(-1)^i*0.01*evecs3(:,2);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot3(z(:,1),z(:,2),E(alpha,beta,z(:,1),z(:,2))+50,'w')
    plot(z(:,1),z(:,2),'b')
    end, 
end
title('Trajectories are still confined to the energy surface','Color','y')
pause
title('Press ENTER to plot the Phase Portrait')
pause

%% 2.5.2. phase portrait
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14,'Box','on')
[x,y] = meshgrid(-5:0.5:5,-10:1:10);
quiver(x,y,y,dv(alpha,beta,gamma,x,y),2,'w'), % show vector field
xlabel('position, x'),ylabel('velocity, v'),
title(['Dissipative Phase Portrait: \beta = -1, \alpha = 10, '...
    '\gamma = 0.25'],'color','y')
axis([-5 5 -10 10])
for i = 1:size(yini,2),     
    [t,z]=ode45(@duffing,0:0.005:2,[yini(:,i);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot(z(:,1),z(:,2),'y','linewidth',1)
end
for i = 1:2, 
    for j=1:2,     
    [t,z]=ode45(@duffing,0:0.01:50,(-1)^j*[xe;0;0]+[(-1)^i*0.01*evecs3(:,1);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot(z(:,1),z(:,2),'r','linewidth',2)
    [t,z]=ode45(@duffing,0:-0.01:-24,(-1)^j*[xe;0;0]+[(-1)^i*0.01*evecs3(:,2);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot(z(:,1),z(:,2),'b','linewidth',2)
    end, 
end
pause
title('Press ENTER to plot the Bifurcation Diagram')
pause

%% 2.5.3. bifurcation diagram
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',12,'Projection','perspective')
alpha = 10;
yini = [4.1:0.5:5.1 5.2,-4.1:-0.5:-5.1 -5.2;...
        -10*ones(1,4),10*ones(1,4)];
for i = 1:size(yini,2),     
    [t,z]=ode45(@duffing,0:0.005:2,[yini(:,i);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10.5|z(:,2)<-10.5); z(indx,:) = NaN;
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'y','linewidth',1)
end
for i = 1:2, 
    for j=1:2,     
    [t,z]=ode45(@duffing,0:0.01:50,(-1)^j*[xe;0;0]+[(-1)^i*0.01*evecs3(:,1);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'r','linewidth',2)
    [t,z]=ode45(@duffing,0:-0.01:-24,(-1)^j*[xe;0;0]+[(-1)^i*0.01*evecs3(:,2);0]);
    indx = find(z(:,1)>5|z(:,1)<-5); z(indx,:) = NaN;
    indx = find(z(:,2)>10|z(:,2)<-10); z(indx,:) = NaN;
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'b','linewidth',2)
    end, 
end
alpha = -10; 
yini = [-1 -1 -1 -1 1 1 1 1 -0.9 -0.75 -0.64 -0.5 -0.25 0.9 0.75 0.64 0.5 0.25;...
         3 2.5 2 1.5 -3 -2.5 -2 -1.5 3 3 3 3 3 -3 -3 -3 -3 -3]*3;
for i = 1:size(yini,2),     
    [t,z]=ode45(@duffing,[0:0.01:3],[yini(:,i);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'y')
end
for i = 1:2,     
    [t,z]=ode45(@duffing,[0:0.005:7],[(-1)^i*0.01*evecs2(:,1);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'r','linewidth',2)
    [t,z]=ode45(@duffing,[0:-0.005:-7],[(-1)^i*0.01*evecs2(:,2);0]);
    indx = find(z(:,1)>3|z(:,1)<-3); z(indx,:) = NaN;
    indx = find(z(:,2)>9|z(:,2)<-9); z(indx,:) = NaN;
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'b','linewidth',2)
end
a = 0:.1:15; aa = -15:.1:0;
x11 = zeros(size(a)); x12 = sqrt(a); x2 = x11;
plot(a,x11,'b','linewidth',2)
plot(aa,x2,'r','linewidth',2)
plot(a,x12,'r','linewidth',2)
plot(a,-x12,'r','linewidth',2)
view(-40,40), axis([-10 10 -5 5 -10 10])
xlabel('bifurcation parameter, \alpha'),ylabel('x'),zlabel('v')
title('Dissipative Bifurcation Diagram (\beta < 0)','color','y')


% -------------------------------------------------------------------------
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];