经过近期的研究发现，目前对于系统单参数分岔图的计算共有以下的几种方法：

1)最大值法

即对系统微分方程(组)进行求解，对求解的结果用getmax函数进行取点，并绘图。

2)Poincare截面法

对系统参数的每一次取值，绘制其Poincare截面，进而得到其分岔图。

这种方法需要注意的是，自治系统的Poincare截面是选取一超平面，平面上点的分布即构成一Poincare截面，非自治系统的Poincare截面则是根据系统激励的频率进行取点并绘图。

本帖将以Lorenz系统为例，对这两种方法进行比较

首先对第二种方法进行阐述。

编程如下(matlab)

Lorenz系统：

function dy = Lorenz(t,y)

% Lorenz系统

% 系统微分方程：

%        dx/dt = -a(x-y)

%        dy/dt = x(r-z)-y

%        dz/dt = xy-bz

%    a=y(4)

%    r=y(5)

%    b=y(6)

dy=zeros(6,1);

dy(1)=-y(4)*(y(1)-y(2));

dy(2)=y(1)*(y(5)-y(3))-y(2);

dy(3)=y(1)*y(2)-y(6)*y(3);

dy(4)=0;

dy(5)=0;

dy(6)=0;

随r的分岔图求解程序：——按照x=y平面取截面

function Lorenz_bifur_r

Z=[];

for r=linspace(1,500,1000);

% 舍弃前面迭带的结果，用后面的结果画图

[T,Y]=ode45('Lorenz',[0,1],[1;1;1;16;r;4]);

[T,Y]=ode45('Lorenz',[0,50],Y(length(Y),:));

Y(:,1)=Y(:,2)-Y(:,1);

% 对计算结果进行判断，如果点满足x=y，则取点

for k=2:length(Y)

f=k-1;

if Y(k,1)<0

if Y(f,1)>0

y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));

Z=[Z r+abs(y)*i];

end

else

if Y(f,1)<0

y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));

Z=[Z r+abs(y)*i];

end

end

end

end

plot(Z,'.','markersize',1)

title('Lorenz映射分岔图')

xlabel('r'),ylabel('|y| where x=y')

结果如图1所示。

getmax法——取最大值法

function [Xmax] = getmax(y)

a=length(y);

j=1;

for i=(a-1)/2:a

b=(y(i,1)-y(i-2,1))/2;

c=(y(i,1)+y(i-2,1))/2-y(i-1,1);

if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0

Xmax(j)=y(i-1,1);

j=j+1;

elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)

Xmax(j)=y(i-1,1)-b^2/(4*c);

j=j+1;

end

end

function Lorenz_bifur_r_getmax

% 最大值法求解分岔图

clear all

t0=[0 100];%积分时间

%bifurcation

for r=linspace(1,500,1000);  %r的变化精度

[t,y]=ode45('Lorenz',t0,[1;1;1;16;r;4]);

[Xmax]=getmax(y(:,1));

plot(r,Xmax,'b','markersize',1)

hold on

clear Xmax

end

计算结果如图2所示！

最后上传一下参考计算机仿真第22卷第12期上一篇文章“李雅普诺夫指数的研究与仿真”中Lorenz系统的分岔图计算结果，大家比较一下即可看出孰优孰劣了！