demo.m

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

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 ASSOCIATED POTENTIAL',...
    'color','y','horizontalAlignment','center','fontsize',14)
pause

%% 1.5.1. potential function
X = -5:0.05:5;
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14,'Box','on')
plot(X,V(alpha,beta,X),'r','linewidth',3)
xlabel('x'),ylabel('Potential, V(x)')
title('Hamiltonian Duffing: \beta = -1, \alpha = 10','color','y')
pause
title('Press ENTER to display the Conserved Energy Surface')
pause

%% 1.5.2. first integral
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(-5:0.05:5,-10:.1:10);
surface(X,Y,E(alpha,beta,X,Y)+50,'facealpha',.9), view([-30,30]), shading interp
xlabel('x'),ylabel('v'),zlabel('First Integral, E(x,v)')
title('Hamiltonian Duffing: \beta = -1, \alpha = 10','color','y')
pause
% -------------------------------------------------------------------------
contour3(X,Y,E(alpha,beta,X,Y)+50,[-20 -10 0 2 6 12 19 25 30 36 42 50]+50)
h = findobj('Type','patch');
set(h,'LineWidth',1,'EdgeColor','w')
title('Trajectories are the Energy Level Sets')
contour(X,Y,E(10,-1,X,Y),[-20 -10 0 2 6 12 19 25 30 36 42 50],'linecolor','b')
pause
title('Press ENTER to view the corresponding Phase Plane')
pause

%% 1.5.3. 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),3,'r','linewidth',1), % show vector field
xlabel('position, x'),ylabel('velocity, v'),
title('Hamiltonian Phase Portrait: \beta = -1, \alpha = 10','color','y')
axis([-5 5 -10 10])
[x,y] = meshgrid(-5:0.02:5,-10:.01:10);
contour(x,y,E(alpha,beta,x,y),[-20 -10 0 2 6 12 19 25 30 36 42 50],'y','linewidth',1);
pause
title('Press ENTER to view the corresponding Bifurcation Diagram')
pause

%% 1.5.4. bifurcation diagram
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14,'Projection','perspective')
[A,X,Y] = meshgrid(-10:20:10,-5:0.05:5,-10:.1:10);
cont1 = contourslice(A,X,Y,E(A,-1,X,Y),-10,[],[],...
    [-35 -30 -25 -20 -15 -10 -5 0 5 15 25 35]);
set(cont1,'EdgeColor','y','LineWidth',1)
cont2 = contourslice(A,X,Y,E(A,-1,X,Y),10,[],[],...
    [-20 -10 0 2 6 12 19 25 30 36 42 50]);
set(cont2,'EdgeColor','y','LineWidth',1)
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('Bifurcation Diagram (\beta < 0)','color','y')
pause
title('Bifurcation Diagram (\beta < 0). Press ENTER for dissipative system.')
pause

%% 2. free dissipative system
gamma = 0.25; f = 0; omega = 0;
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.5,'Now consider unforced dissipative case,',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.45,...
    'that is, (\gamma = 0.25; f = 0; & \omega = 0):',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.35,...
    '(d^2x/dt^2) + \gamma (dx/dt) + \alpha x + \beta x^3 = 0 .',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.1,'PRESS ENTER TO CONTINUE',...
    'color','y','horizontalAlignment','center','fontsize',14)
pause

%% 2.1. 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.8,'First consider positive \beta case (w.l.o.g., \beta = 1).',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.7,'Thus, for equilibria wehave the same two 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 see the Energy Surface',...
    'color','y','horizontalAlignment','center','fontsize',14)
pause

%% 2.1.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(-1:0.05:1,-3:.1:3);
surface(X,Y,E(alpha,beta,X,Y)+3,'facealpha',.9), view([-35,22]), shading interp
xlabel('x'),ylabel('v'),zlabel('First Integral, E(x,v)')
title(['This is still the Total Energy: \beta = 1, \alpha = 10,'...
       '\gamma = 0.25'],'Color','y')
% -------------------------------------------------------------------------
pause
title('Press ENTER to see the motion on this surface')
pause
yini = [0 0; -3 3];
for i = 1:2,     
    [t,z]=ode45(@duffing,[0 40],[yini(:,i);0]);
    plot3(z(:,1),z(:,2),E(10,1,z(:,1),z(:,2))+3,'w')
    plot(z(:,1),z(:,2),'b'), axis([-1 1 -3 3 0 13])
end
title('Trajectories are still confined to the energy surface','Color','y')
pause
title('Press ENTER to plot the Phase Portrait')
pause

%% 2.1.3. 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(-1:0.1:1,-3:.3:3);
quiver(x,y,y,dv(alpha,beta,gamma,x,y),1,'r'), % show vector field
xlabel('position, x'),ylabel('velocity, v'),
title('Dissipative Phase Portrait (\beta = 1, \alpha = 10, \gamma = 0.25)','color','y')
axis([-1 1 -3 3])
for i = 1:2,     
    [t,z]=ode45(@duffing,[0 40],[yini(:,i);0]);
    plot(z(:,1),z(:,2),'y','linewidth',1), axis([-1 1 -3 3 0 13])
end
pause
title('Press ENTER to considere Case 2: \alpha < 0.')
pause

%% 2.2. 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 the same 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 ASSOCIATED POTENTIAL',...
    'color','y','horizontalAlignment','center','fontsize',14)
pause

%% 2.2.2. 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(-5:0.05:5,-10:.1:10);
surface(X,Y,E(alpha,beta,X,Y)+50,'facealpha',.9), view([-30,30]), shading interp
xlabel('x'),ylabel('v'),zlabel('First Integral, E(x,v)')
title('The Same Energy Surface: \beta = 1, \alpha = -10','color','y')
% -------------------------------------------------------------------------
pause
title('Press ENTER key to see the motion on this surface'), pause
[evecs1,evals] = eig([0,1;-alpha,-gamma]); % linearized eigen value problem
for i = 1:2,     
    [t,z]=ode45(@duffing,[0 40],[(-1)^i*0.01*evecs1(:,1);0]);
    plot3(z(:,1),z(:,2),E(alpha,beta,z(:,1),z(:,2))+50,'w')
    plot(z(:,1),z(:,2),'r'), axis([-5 5 -10 10 0 130])
    [t,z]=ode45(@duffing,[0 -5.1],[(-1)^i*0.01*evecs1(:,2);0]);
    plot3(z(:,1),z(:,2),E(alpha,beta,z(:,1),z(:,2))+50,'w')
    plot(z(:,1),z(:,2),'b'), axis([-5.2 5.2 -12 12 0 130])
end
title('Trajectories traverse the energy surface','Color','y')
pause
title('Press ENTER to plot the corresponding Phase Portrait')
pause

%% 2.2.3. 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.5:0.5:5.5,-12:1:12);
quiver(x,y,y,dv(-10,1,0,x,y),1.8,'r','linewidth',1), % show vector field
xlabel('position, x'),ylabel('velocity, v'),
title('Dissipative Phase Portrait (\beta = 1, \alpha = -10, \gamma = 0.25)','color','y')
axis([-6 6 -12 12])
for i = 1:2,     
    [t,z]=ode45(@duffing,[0 40],[(-1)^i*0.01*evecs1(:,1);0]);
    plot(z(:,1),z(:,2),'r','linewidth',2)
    [t,z]=ode45(@duffing,[0 -12],[(-1)^i*0.01*evecs1(:,2);0]);
    plot(z(:,1),z(:,2),'b','linewidth',2)
    eval(['z' num2str(i) ' = z;']);
end
pause
title('Press ENTER key and observe BASINS OF ATTRACTION')
pause
l = 320;
fill(-[z1(:,1);flipud(z2(1:l,1))],[z2(:,2);flipud(z1(1:l,2))],'y',...
    'linestyle','none','facealpha',.7)
title('Unforced Two-Well Duffing''s Basins of Attraction','color','w')
pause
title('Press ENTER to view the Bifurcation Diagram')
pause

%% 2.3. bifurcation diagram
clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14,'Projection','perspective')
alpha = -10;
for i = 1:2,     
    [t,z]=ode45(@duffing,[0 30],[(-1)^i*0.01*evecs1(:,1);0]);
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'r','linewidth',1),
    [t,z]=ode45(@duffing,[0 -7.5],[(-1)^i*0.01*evecs1(:,2);0]);
    plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'b','linewidth',1),
    eval(['z' num2str(i) ' = z;']);
end
alpha = 10; 
[t,z]=ode45(@duffing,[0 40],[0;11;0]);
plot3(alpha*ones(size(z(:,1))),z(:,1),z(:,2),'r','linewidth',1)
a = 0:.1:15; aa = -15:.1:0;
x11 = zeros(size(a)); x12 = sqrt(-aa); x2 = x11;
plot(a,x11,'b','linewidth',2)
plot(aa,x2,'r','linewidth',2)
plot(aa,x12,'b','linewidth',2)
plot(aa,-x12,'b','linewidth',2)
view(-45,25), axis([-10 10 -6 6 -18 18])
xlabel('bifurcation parameter, \alpha'),ylabel('x'),zlabel('v')
title('Bifurcation Diagram (\beta > 0)','color','y')
pause
title('Press ENTER to see Basins of Attraction')
pause
l = 157;
fill3(-10*(ones(size([z2(:,2);flipud(z1(1:l,2))]))),-[z1(:,1);...