混沌matlab模拟

1.Lorenz模型

洛仑兹在研究天气的不可预测性时，从流体的运动方程出发，通过简化方程获得了具有三个自由度的系统

dx dt =−10x+10y

dy dt =28x−y−xz

dz dt =−

8

3

z+xy

其中x、y、z为无量纲量，分别表征对流强度，对流中升流与降流间的温差和竖直方向温度分布的非线性度。任意给定初值，系统最终都会回到状态空间的特定区域内

若状态轨迹经过一段时间之后停在一个不动点上，那么意味着系统进入了一个稳定的状态，这相轨迹将是一个平庸吸引子。然而，事实上，相轨迹在两片上“随机”地跳来跳去，说明系统的状态演变着有某种规律性，这种相图不对应任何一种定常状态，因此，被称为奇异吸引子，又称洛仑兹吸引子。

m.function dx=Lorenz(t,x)

dx=[10*(x(2)-x(1));28*x(1)-x(1)*x(3)-x(2);x(1)*x(2)-(8/3)*x(3)];

end

[T,X]=ode45('Lorenz',[0,30],[12;4;0]);

vx=10*(X(:,2)-X(:,1));

vy=28*(X(:,1)-X(:,1).*X(:,3)-X(:,2));vz=X(:,1).*X(:,2)-(8 /3)*X(:,3);

>> subplot(2,2,1); plot3(X(:,1),X(:,2),X(:,3))

>>subplot(2,2,2); plot(X(:,1),vx)

>>subplot(2,2,3);plot(X(:,2),vy)

>>subplot(2,2,4);plot(X(:,3),vz)

2.虫口

x n = g n x0

x0是开始计算的那一代人口数。只要g>1，x n很快就趋向无穷大，发生“人口爆炸”。这样的线性模型，不能完全反应人口的变化规律，但是稍加修正，就可以称为描述某些没有世代交叠的昆虫数目的虫口方程。

虫口数目太多时，由于争夺有限的食物和生存空间发生咬斗，由于接触传染而导致疾病蔓延，争斗使虫口数目减少的事件，这些事件的数目比例于x n2，于是方程可以修正为：

x n+1 = gx n (1-x n)

取g =2, x0=0.9, x1=0.18,…, x n=0.5，它停在那儿不动了。即在x n=0.5处有一个点吸引子，一个稳定定态。若追踪这个种群，则会发现种群数目随着时间的演化而保持稳定的数值。

振荡称为周期2循环，即若跟踪种群，会发现种群数目每隔一年，数目重复循环一次，就象有些果树有大年小年一样，x1n和x2n也是定点吸引子。

继续增加g值，还可得周期4循环，周期8循环，周期16循环等等。每一次解的周期都增加一倍。当g 达到某一临界值时，继续增加，迭代结果再也不循环了，而是疯狂地振荡，永远也不会稳定下来，我们称为混沌态。

若以g 为横坐标，迭代结果为纵坐标，可得分岔图。从临界值开始，逻辑斯蒂映射进入了混沌区，在这种情况下，种群的数目就完全不能预测了。

（1）g=1,2,4时对应迭代次数和最后X n的关系

syms x0

for k=1:3

t=0;x0=0.8;g=2^(k-1);

while(t<20)

t=t+1;

x0=g*x0*(1-x0);

u(t,:)=x0;

v(t,:)=t;

end

subplot(1,3,k);plot(v,u,'.')

end

（2）若以g 为横坐标，迭代稳定结果为纵坐标，得分岔图。syms x0

for n=1:401

g=[0:0.01:4];

t=0;x0=0.8;

while(t<100)

if(t>50)

t=t+1;

x0=g(n)*x0*(1-x0);

u(t,:)=x0;

v(t,:)=g(n);

plot(v(t,:),u(t,:),'.')

hold on

else

t=t+1;

x0=g(n)*x0*(1-x0);

end

end end grid

从上图可以看出g(0~3)一倍解，（3~3.45）二倍解，（3.45~3.55）四倍解，（3.6）八倍解

3.单摆（θ0=0.0，1v 0=0，Q =2，A =12

，w =2

3

）

单摆是在悬挂的细线的另一端连接着一个小球（如图所示）。单摆又称数学摆，是物理学中最简单的模型之一。 可以认为，细线的质量可以忽略，且是刚性的。系统质量集中在可视为质点的小球上。设摆长为l ，小球质量为m ，

相对于平衡的下垂位置的角度为θ，重力加速度为g。则其运动方程为（方程中加上了阻尼，驱动）

d2θ

2+

dθ

+sinθ=Acos(wt)

无阻尼相图：

s=dsolve('D2s+9.8*s=0','s(0)=0.01,Ds(0)=0') s =

cos((7*5^(1/2)*t)/5)/100

v=diff(s,1)

v =

-(7*5^(1/2)*sin((7*5^(1/2)*t)/5))/500

ezplot(s,v)

有阻尼无驱动相图

function dy=L(t,y)

dy=[y(2);-sin(y(1))-0.5*y(2)];

end

[t,Y]=ode45('L',[0,50],[0.01;0]);

plot(Y(:,1),Y(:,2))

有阻尼有驱动相图

function dy=L(t,y)

dy=[y(2);0.5*cos((2/3)*t)-sin(y(1))-0.5*y(2)]; end

[t,Y]=ode45('L',[0,100],[0.01;0]);

>> plot(Y(:,1),Y(:,2))