twowell.m

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

function [] = twowell()
% this is just a demo script fot two-well duffing equation. see
% corresponding static bifurcation script.
% Call: > twowell
% copirught (c) david chelidze 2008
% last modified 02/26/2008

close all, clear all

%% 0. initialize duffing global parameters
global alpha beta gamma f omega % define global variables

warning off all

V = inline('alpha.*x.^2/2+beta.*x.^4/4'); % define a potential
E = inline('alpha.*x.^2/2+beta.*x.^4/4+y.^2/2'); % define a first integral
dv = inline('-gamma.*y-alpha.*x-beta.*x.^3'); % define vector field

clf reset, cla reset, set(gcf,'color','k'), hold on
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'FontSize',14), axis off
txt = [{'This demos two-well Duffing Equation:'};
       {[]};
       {'(d^2x/dt^2) + \gamma (dx/dt) + \alpha x + \beta x^3'};
       {[]};
       {'We start by considering unforced conservative case,'};
       {[]};
       {'that is, Hamiltonian system (\gamma = 0; f = 0; & \omega = 0):'};
       {[]};
       {'(d^2x/dt^2) + \alpha x + \beta x^3 = 0 .'};
       {[]};
       {'PRESS ENTER TO CONTINUE...'}]; 
text(0.5,0.5,txt,'HorizontalAlignment','center','color','y','fontsize',14)
pause

%% 1. free hamiltonian system
gamma = 0;

%% 1.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
txt = [{'For two-well system we need: \alpha < 0.'};
       {[]};
       {'There are three equilibria: (0,0) and (+/-(-\alpha)^{1/2},0).'};
       {[]};
       {'PRESS ENTER TO DISPLAY THE ASSOCIATED POTENTIAL'}]; 
text(0.5,0.5,txt,'HorizontalAlignment','center','color','y','fontsize',14)
pause

%% 1.1.1. potential function
X = -5:0.02: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) = \alpha x^2/2 + \beta x^4/4')
title('Hamiltonian Duffing: \beta = 1, \alpha = -10','color','y')
pause
title('Press ENTER key to display the Conserved Energy Surface',...
    'color','b')
pause

%% 1.1.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([-33,34]), 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,[-24 -20 -15 -10 -5 0 5 12 18 25]+50)
h = findobj('Type','patch');
set(h,'LineWidth',1,'EdgeColor','w')
title('Trajectories are the Energy Level Sets')
contour(X,Y,E(alpha,beta,X,Y),[-24 -20 -15 -10 -5 0 5 12 18 25],'linecolor','b')
pause
title('Press ENTER to view the corresponding Phase Portrait',...
    'color','b')
pause

%% 1.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(-5:0.5:5,-10:1:10);
quiver(x,y,y,dv(alpha,beta,gamma,x,y),1.8,'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.05:5,-10:.1:10);
contour(x,y,E(alpha,beta,x,y),[-24 -20 -15 -10 -5 0 5 12 18 24],'y','linewidth',1);
pause
title('Press ENTER to consider dissipative case',...
    'color','b')
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
txt = [{'Now consider unforced dissipative case,'};
       {[]};
       {'that is, (\gamma = 0.25; f = 0; & \omega = 0):'};
       {[]};
       {'(d^2x/dt^2) + \gamma (dx/dt) + \alpha x + \beta x^3 = 0 .'};
       {[]};
       {'PRESS ENTER TO CONTINUE'}]; 
text(0.5,0.5,txt,'HorizontalAlignment','center','color','y','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
txt = [{'Again, for two-well we need \alpha < 0.'};
       {[]};
       {'There are the same three equilibria: (0,0) and (+/-(-\alpha)^{1/2},0).'};
       {[]};
       {'PRESS ENTER TO DISPLAY THE ASSOCIATED POTENTIAL'}]; 
text(0.5,0.5,txt,'HorizontalAlignment','center','color','y','fontsize',14)
pause

%% 2.1.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','color','b')
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','color','b')
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(-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','color','b')
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','y')
pause
title('Press ENTER to consider forced two-well Duffing','color','b')
pause
% End UNFORCED

%% 3. Begin FORCED

clf reset, cla reset, set(gcf,'color','k'), set(gca,'color','k'), hold on
txt = [{'Now consider FORCED Two-Well Duffing.'};
       {[]};
       {'The time is now another state variable.'};
       {[]};
       {'Thus, STROBOSCIPIC POINCARE SECTION is appropriate.'};
       {[]};
       {'For small amplitudes, sinks become limit cycles.'};
       {[]};
       {'We first look at the time and frequency domain data.'};
       {[]};
       {'Press ENTER to see time history and PSD of P-1 motion.'}]; 
text(0.5,0.5,txt,'HorizontalAlignment','center','color','y','fontsize',14)
pause,

%% 3.1 Look at the small P-1 motion in one energy well
clf reset, cla reset, set(gcf,'color','k')
gamma = .25; omega = 1; f = .2; beta = 1; alpha = -1;
[t,y]=ode45(@duffing,[0:2*pi/omega/36:400*pi/omega],[1.19718071928 0.08451532560 0]);
subplot(2,1,1)
plot(y(1:300,1),'color','r','linewidth',1), grid
set(gca,'XColor','y','YColor','y','ZColor','y','Color','k',...
    'Box','on')
xlabel('Time'),ylabel('Amplitude')
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,240,'Press ENTER for Poincere Sections',...
    'HorizontalAlignment','center','color','b','fontsize',14)
pause

%% 3.1.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',15)
axis([-2 2 -2 2]), xlabel('u'), ylabel('v')
title(['Poincare Sections of P-1 orbits in Each Energy Well, '...
    '\omega = 1, f = 0.2'],'color','y'), hold on
[t,yn]=ode45(@duffing,[0:2*pi/omega/36:400*pi/omega],...
    [-0.79007534131 0.04362085397 0]);
plot(yn(1:36:end,1),yn(1:36:end,2),'r.','markersize',15)
pause(2)
title('Press ENTER to see P-1 Projections onto Poincare Sections',...
    'color','b')
pause
plot(y(:,1),y(:,2),'y',yn(:,1),yn(:,2),'y','linewidth',1)
title('Press ENTER to visualize the trajectory in the Phase Space',...
    'color','b')
pause

%% 3.1.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([27,67])
axis([-6.0000   6.0000   -5.0000    5.0000   -0.3000    0.3000])
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],[-.3 -.3 .3 .3 -.3],'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);