matlab求洛伦兹方程的解,[转载]用Matlab求解洛伦兹方程

1. 洛伦兹方程求解

本文说明用Matlab工具箱求解洛伦兹方程的过程，并给出吸引子的三维动态图象.洛伦兹方程如下：

(1)这是一个自洽的方程组，求解过程如下：

(1) 建立自定义函数

function

dy=Lorenz(t,y) %

y(1)=x,y(2)=y,y(3)=z

dy=zeros(3,1);

dy(1)=10*(-y(1)+y(2));

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

dy(3)=y(1)*y(2)-8*y(3)/3;

(2)用ode45命令求解

[t,y]=ode45(@Lorenz,[0,30],[12,2,9]);

subplot(221);

plot(t,y(:,1));

subplot(222);

plot(t,y(:,2))

subplot(223);

plot(t,y(:,3))

subplot(224);

plot3(y(:,1),y(:,2),y(:,3))

view([20 42]);

(3)求解结果

(4)动态显示吸引子的绘制过程

[t,y]=ode45(@Lorenz,[0 30],[12 2 9]);

clf;

axis([-20

20 -25 25 10 50]);

view([20

42]);

hold

on;

comet3(y(:,1),y(:,2),y(:,3));%显示吸引子的绘制过程

(5)生成动画

[t,]=ode45(@Lorenz,[0

30],[12 2 9]);

m=moviein(100);

axis([-20 20 -25 25 -10

50]);

shading flat;

h=plot3(y(:,1),y(:,2),y(:,3));

for j=1:100

rotate(h,[0 0 1],1.8); %沿Z轴旋转

axis([-20 20 -25 25 -10 50]);

shading flat;

m(:,j)=getframe;

[X,map]=getframe;

imwrite(X,map,['Lz'int2str(j)'.bmp'],'bmp');%写入bmp文件

end

movie(m,5) %循环播放5次

(6)验证蝴蝶效应

clf;

hold on;

[t,u]=ode45('Lorenz',[0 6],[12 2 9]);

plot(t,u(:,3),'color','r');

[t,v]=ode45('Lorenz',[0 6],[12 2.01 9]);

plot(t,v(:,3),'color','b');

[t,w]=ode45('Lorenz',[0 6],[12 1.99 9]);

plot(t,w(:,3),'color','k');