理力-龙格-张新意.doc

龙格-库塔法解洛伦兹方程 结构01 张新意 一．洛伦兹方程 美国科学家Lorenz在研究大气对流模型时得到了以下方程组 X=σY-XY=ρX-Y-XZZ=-βZ+XY 称为Lorenz方程。 现对采用常用具体参数σ=10，β=8/3，ρ=28的Lorenz方程进行求解，描绘出X、Y、Z随时间变化的关系图形和三维图形，并验证该方程的初值敏感性质。 二．求解过程 （1）建立函数 用三维向量dx代表X、Y、Z，程序代码为： function dx=Lorenz(t,x) dx=zeros(3,1); dx(1)=10*(-x(1)+x(2)); dx(2)=28*x(1)-x(2)-x(1)*x(3); dx(3)=x(1)*x(2)-8*x(3)/3; 保存为Lorenz.m （2）用四五阶龙格-库塔法建立X、Y、Z关于时间的函数关系并绘制图形 程序代码： clc clear [t,x]=ode45(lorenz,[0 20],[10,10,10]); figure(1) hold on title(X(t)图像) plot(t,x(:,1)) xlabel(t) ylabel(X) figure(2) hold on title(Y(t)图像) plot(t,x(:,2)) xlabel(t) ylabel(Y) figure(3) hold on title(Z(t)图像) plot(t,x(:,3)) xlabel(t) ylabel(Z) 结果： （3）绘制点（X,Y,Z）的运动轨迹 程序代码： [t,x]=ode45(Lorenz,[0 50],[10,10,10]); hold on title((X,Y,Z)运动轨迹); plot3(x(:,1),x(:,2),x(:,3)); 运行结果： （4）验证Lorenz方程的初值敏感性质，即“蝴蝶效应” 将坐标Y的初始值做微小改变（改变0.01）,绘制出改变前后X(t)、Z(t)曲线，进行比较。 程序代码： clc clear [t1,x1]=ode45(Lorenz,[0 15],[10,10,10]); [t2,x2]=ode45(Lorenz,[0 15],[10,10.01,10]); figure(1) hold on title(X(t)图像) xlabel(t) ylabel(X) plot(t1,x1(:,1),k) plot(t2,x2(:,1),b) figure(2) hold on title(Z(t)图像) xlabel(t) ylabel(Z) plot(t1,x1(:,3),k) plot(t2,x2(:,3),b) 运行结果： 从图中可以明显看出虽然只对Y的初始值进行了微小调整，但对却方程的结果造成了很大影响，“蝴蝶效应”得以验证。