matlab编程分叉与混沌,[转载]Matlab——混沌分岔图

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

1)最大值法

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

2)Poincare截面法

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

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

本帖将以Lorenz系统为例&＃xff0c;对这两种方法进行比较

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

编程如下(matlab)

Lorenz系统&＃xff1a;

function dy &＃61; Lorenz(t,y)

% Lorenz系统

% 系统微分方程&＃xff1a;

%

dx/dt &＃61; -a(x-y)

%

dy/dt &＃61; x(r-z)-y

%

dz/dt &＃61; xy-bz

%

a&＃61;y(4)

%    r&＃61;y(5)

%

b&＃61;y(6)

dy&＃61;zeros(6,1);

dy(1)&＃61;-y(4)*(y(1)-y(2));

dy(2)&＃61;y(1)*(y(5)-y(3))-y(2);

dy(3)&＃61;y(1)*y(2)-y(6)*y(3);

dy(4)&＃61;0;

dy(5)&＃61;0;

dy(6)&＃61;0;​

随r的分岔图求解程序&＃xff1a;——按照x&＃61;y平面取截面

function Lorenz_bifur_r

Z&＃61;[];

for r&＃61;linspace(1,500,1000); %

舍弃前面迭带的结果&＃xff0c;用后面的结果画图

[T,Y]&＃61;ode45(&＃39;Lorenz&＃39;,[0,1],[1;1;1;16;r;4]);

[T,Y]&＃61;ode45(&＃39;Lorenz&＃39;,[0,50],Y(length(Y),:));

Y(:,1)&＃61;Y(:,2)-Y(:,1); %

对计算结果进行判断&＃xff0c;如果点满足x&＃61;y&＃xff0c;则取点

for k&＃61;2:length(Y)

f&＃61;k-1;

if

Y(k,1)<0

if Y(f,1)>0

y&＃61;Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));

Z&＃61;[Z r&＃43;abs(y)*i];

end

else     if

Y(f,1)<0  y&＃61;Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));

Z&＃61;[Z r&＃43;abs(y)*i];

end

end

end

end

plot(Z,&＃39;.&＃39;,&＃39;markersize&＃39;,1)title(&＃39;Lorenz映射分岔图&＃39;)xlabel(&＃39;r&＃39;),ylabel(&＃39;|y|

where x&＃61;y&＃39;)

getmax法——取最大值法

function [Xmax] &＃61; getmax(y)

a&＃61;length(y);

j&＃61;1;

for i&＃61;(a-1)/2:a

b&＃61;(y(i,1)-y(i-2,1))/2;

c&＃61;(y(i,1)&＃43;y(i-2,1))/2-y(i-1,1);

if

y(i-2,1)<&＃61;y(i-1,1)&y(i-1,1)>&＃61;y(i,1)&c&＃61;&＃61;0

Xmax(j)&＃61;y(i-1,1);

j&＃61;j&＃43;1;

elseif

y(i-2,1)<&＃61;y(i-1,1)&y(i-1,1)>&＃61;y(i,1)Xmax(j)&＃61;y(i-1,1)-b^2/(4*c);

j&＃61;j&＃43;1;

end

end

function Lorenz_bifur_r_getmax% 最大值法求解分岔图

clear all

t0&＃61;[0 100];%积分时间

%bifurcation

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

[t,y]&＃61;ode45(&＃39;Lorenz&＃39;,t0,[1;1;1;16;r;4]);

[Xmax]&＃61;getmax(y(:,1));

plot(r,Xmax,&＃39;b&＃39;,&＃39;markersize&＃39;,1)

hold on

clear Xmax

end