混沌时间序列的几个例子

最近在研究混沌时间序列，想了解基本知识，看的一本书叫做《混沌时间序列分析及其应用》（作者：吕金虎）（是中文书籍，不是翻译的书籍）

书中第二章 混沌动力学 提到了几个典型的混沌的例子

1.虫口模型 Logistics

1.1 原理

虫口模型---昆虫变化的简单数学模型。（昆虫一次产卵的数量极大，子代出生后，其亲代的数量可以忽略不计）认为昆虫的种群是由单一的世代构成，在历代之间没有重叠。因此，昆虫种群的增长（虫口数）是分步进行的，例如，每个3年、5年或7年出现一次成虫期的蝉就是这样。

1.2 代码

说明：首先进行一些迭代，使之大概达到一种有迹可循的状态；然后就记录结果。

横轴可以看做是迭代次数，横轴上的一点，对应的纵轴上很多个点的叠加。

%% -------------------虫口模型-分岔图-----------

% 虫口模型 x(n+1) = lambda * x(n) * (1 - x(n)) 

% 当 lambda 从 3 到 4 的过渡图像

% 参见<<混沌动力学初步>>陈士华，陆君安编著 P46

clc

clear

close all

lambda = 0:5e-4:4;%原来是3:4

x = 0.4*ones(1,length(lambda));

N1 = 800;                   % 前面的迭代点数(原来是400)

N2 = 200;                   % 后面的迭代点数（原来是100）

f = zeros(N1+N2,length(lambda));

for i = 1:N1+N2

    x = lambda .* x .* (1 - x);

    f(i,:) = x;

end

f = f(N1+1:end,:);

plot(lambda,f,'k.','MarkerSize',1)

xlabel('\lambda')

ylabel('x');

title('Logistic虫口模型-虫口随着λ的变化')

1.3 结果

2.洛伦兹方程-天气对流模型

2.1 原理

1963年，美国著名气象学家洛伦兹在研究天气预报问题时，发现确定性方程中出现混沌现象。

振荡曲线与完全随机的噪声不一样，它似乎有一定的规律：无论时间多长，它都不会完全重复以前某时刻的振荡，即此时振荡具有非周期性和随机性。这些轨迹被限制在相平面（空间）的有限区域内，这样的有限区域称为动力系统的吸引子（attractor）。对于周期运动，吸引子是一条简单的闭曲线。但形如洛伦兹系统这样的具有随机性的费周期运动，其吸引子远非什么复杂的闭曲线，其轨线往往是反复折叠和相互交叉而形成的密集的带。这种具有复杂结构的吸引子称为奇怪吸引子（strange attractor）。具有奇怪吸引子的运动就是混沌的。

2.2 代码

2.2.1 主函数

%本程序的作用：求解&绘图
[t,y]=ode45('Lorenz',[0 30],[12,2,9]);

（可以写作： [t,y]=ode45(@Lorenz,[0 30],[12,2,9]) ）
figure(1)

plot(t,y(:,1));
title('x关于t的变化图');
xlabel('t');
ylabel('x(t)');

figure(2)
plot(t,y(:,2));
title('y关于t的变化图');
xlabel('t');
ylabel('y(t)');

figure(3)
plot(t,y(:,3));
title('z关于t的变化图');
xlabel('t');
ylabel('z(t)');

figure(4)
plot3(y(:,1),y(:,2),y(:,3))
title({'x,y,z的空间关系图','Lorenz方程在三维相空间中的轨迹图形'}');
xlabel('x');ylabel('y'),zlabel('z');
view([20,32]);


%% ----验证“蝴蝶效应”----------
%解释：洛伦兹方程的解对初始值十分敏感，现在对y的初始值稍加修改，将2改为2.01和1.99，
%           然后求解z的数值解。
%结果说明：随着时间的推移，三条曲线的吻合程度越来越差，差距越来越大，变化也越来越不明显，成为混沌状态、。
%                   前面一部分（t<5），三条曲线几乎完全重合。

% clf %清楚当前figure内的图形
figure(5)
hold on
[t,u]=ode45('Lorenz',[0 15],[12,2,9]);
plot(t,u(:,3),'Color','r');
[t,v]=ode45('Lorenz',[0 15],[12,2.01,9]);
plot(t,v(:,3),'Color','b');
[t,w]=ode45('Lorenz',[0 15],[12,1.99,9]);
plot(t,w(:,3),'Color','k');
hold off
legend('y(t=0)=2','y(t=0)=2.01','y(t=0)=1.99');
xlabel('t');
ylabel('z(t)');
title('蝴蝶效应-y的不同的初始值下，z的取值变化');

2.2.2 函数 Lorenz.m

% （1）洛伦兹方程是一组常微分方程，Runge-Kutta方法是求解常微分方程的经典方法，
% Matlab提供了多个采用该方法的函数命令：ode23 ，ode45 等。
% （2）下面采用ode45命令求解Lorenz在研究大气对流模型时得到的洛伦兹方程：
% dx/dt=A*(-x+y)
% dy/dt=B*x-y-x*z
% dz/dt=x*y-C*z
% 经常引用的参数是： A=10，B=2，C=8/3   相空间维数为3
% （3）为了用Matlab求解，将x,y,z表示为y(1),y(2),y(3)，即为列向量中的三个分量
function dy=Lorenz(t,y)
A=10;
B=28;
C=8/3;
dy=zeros(3,1);
dy(1)=A*(-y(1)+y(2));
dy(2)=B*y(1)-y(2)-y(1)*y(3);
dy(3)=y(1)*y(2)-C*y(3);
end

3.Chen's吸引子

3.1 微分方程求解

3.1.1定义函数Chens.m

function dy=Chens(t,y)
dy=zeros(3,1);%一个三行的矢量
A=35;
B=3;
C=28;
dy(1)=A*(y(2)-y(1));
dy(2)=(C-A)*y(1)-y(1)*y(3)+C*y(2);
dy(3)=y(1)*y(2)-B*y(3);
end

3.1.2求解方程

%二阶、三阶Runge-Kutta法积分，采用函数：ode23
%四阶、五阶Runge-Kutta法积分，采用函数：ode45


[T,Y]= ode45(@Chens,[0 200],[0 1 0]);
plot3(Y(:,1),Y(:,2),Y(:,3));
xlabel('x axis');ylabel('y axis');zlabel('z axis');


view(0,90);%在xoy面的投影
view(90,0);%在yoz面的投影
view(0,0);%在xoz面的投影

3.2 线性方程组求解

%% ---------求解平衡点---------------
syms x y z
[solx,soly,solz]= solve(35*(y-x) = =0,-7*x-x*z+28*y==0,x*y-3*z==0,x,y,z)
solutions=[solx,soly,solz]

求得的解如下：

solutions =
 
[          0,          0,  0]
[ -3*7^(1/2), -3*7^(1/2), 21]
[  3*7^(1/2),  3*7^(1/2), 21]