demo.m

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

function [] = demo()
% This is a small graphical demonstration of the Duffing equation and its
% dynamical characteristics.
% CALL: > demo
%
% script by David Chelidze
% (c) 2008 David Chelidze
% last modified 02/25/2008

%% 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
text(0.5,0.7,'This demonstration introduces the Duffing Equation:',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.6,...
    ['(d^2x/dt^2) + \gamma (dx/dt) + \alpha x + \beta x^3'...
    ' = f cos(\omega t) .'],...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.5,'We start by considering unforced conservative case,',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.45,...
    'that is, Hamiltonian system (\gamma = 0; f = 0; & \omega = 0):',...
    'color','y','horizontalAlignment','center','fontsize',14)
text(0.5,0.35,...
    '(d^2x/dt^2) + \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

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

%% 1.1.1. potential function
X = -1:0.05:1;
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 TO SEE THE CORRESPONDING ENERGY SURFACE','color','y')
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(-1:0.05:1,-3:.1:3);
surface(X,Y,E(alpha,beta,X,Y)+3,'facealpha',.9), view([-30,36]), 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)+3,10)
h = findobj('Type','patch');
set(h,'LineWidth',1,'EdgeColor','w')
title('Trajectories Are the Energy Level Sets','Color','y')
contour(X,Y,E(alpha,beta,X,Y),[0.2 0.5 1 1.5 2 2.5 3 3.5 4],'linecolor','b')
pause
title('Press ENTER to plot the corresponding Phase Portrait')
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(-1:0.1:1,-3:.3:3);
quiver(x,y,y,dv(alpha,beta,gamma,x,y),1,'r','linewidth',1), % show vector field
xlabel('position, x'),ylabel('velocity, v'),
title('Hamiltonian Phase Portrait: \beta = 1, \alpha = 10','color','y')
axis([-1 1 -3 3])
[x,y] = meshgrid(-1:0.01:1,-3:.03:3);
contour(x,y,E(alpha,beta,x,y),[0.15 0.5 1.1 1.9 3 4.3],'y','linewidth',1);
pause
title('Press ENTER to considere Case 2: \alpha < 0.')
pause

%% 1.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 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.2.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')
pause

%% 1.2.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')
pause

%% 1.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: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 view the corresponding Bifurcation Diagram')
pause

%% 1.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')
[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,[],[],...
    [-24 -20 -15 -10 -5 0 5 12 18 24]);
set(cont1,'EdgeColor','y','LineWidth',1)
cont2 = contourslice(A,X,Y,E(A,1,X,Y),10,[],[],[0.15 0.5 1.1 1.9 3 4.3]*2);
set(cont2,'EdgeColor','y','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(-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 to consider \beta < 0 case.')
pause

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

%% 1.4.1. potential function
X = -1:0.05:1;
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 to display the conserved Energy Surface')
pause

%% 1.4.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(X,-3:.1:3);
surface(X,Y,E(alpha,beta,X,Y)+10,'facealpha',.9), view([-55,25])
axis([-1 1 -3 3 0 15]), 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)+10,[-4 -3 -2 -1 0 1 2 3 4]+10)
h = findobj('Type','patch');
set(h,'LineWidth',1,'EdgeColor','w')
title('Trajectories are the Energy Level Sets','Color','y')
contour(X,Y,E(alpha,beta,X,Y),[-4 -3 -2 -1 0 1 2 3 4],'linecolor','b')
pause
title('Press ENTER to view the corresponding Phase Space')
pause

%% 1.4.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.8,'r','linewidth',1), % show vector field
xlabel('position, x'),ylabel('velocity, v'),
title('Hamiltonian Phase Portrait: \beta = -1, \alpha = -10','color','y')
axis([-1 1 -3 3])
contour(x,y,E(alpha,beta,x,y),[-4 -3 -2 -1 0 1 2 3 4],'y','linewidth',1);
pause
title('Press ENTER to considere Case 2: \alpha > 0.')
pause


%% 1.5. beta = -1, alpha = 10