matlab绘制庞加莱截面_matlab庞加莱截面法画Lorenz系统分岔图(附图).doc

利用庞加莱截面法 画的Lorenz 系统的分岔图，复制改成其他系统即可运行，Matlab12a可以运行，附分岔图，见下页！

function Lorenz_bifur_r

Z=[];

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

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

[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')

function dy = Lorenz(t,y)

% Lorenz系统

% 系统微分方程：

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

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

% dz/dt = xy-bz

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;