数值分析 Lorenz问题与混沌

1、991790 力92 王为实验6.1（Lorenz问题与混沌）问题提出：考虑著名的Lorenz方程其中s,r,b为变化区域有一定限制的实参数，该方程形式简单，表面上看并无惊人之处，但由该方程揭示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。实验要求：（1）请读者找出Lorenz方程与上述程序中使用的方程间的关系。（2）对目前取定的参数值SIGMA、RHO和BETA，选取不同的初始值y0（当前程序中的y0是在坐标原点），运行上述的程序，观察计算的结果有什么特点？解的曲线是否有界？解的曲线是不是周期的或趋于某个固定的点？（3）在问题允许的范围内适当改变其中的参数

2、值SIGMA、RHO和BETA，再选取不同的初始值y0，运行上述的程序，观察并记录计算的结果有什么特点？是否发现什么不同的现象。程序清单：主程序：%BLZ Plot the orbit around the Lorenz chaotic attractorclfclc%Solve the ordinary differential equation describing the %Lorenz chaotic attractor.The equation are difined in%an M-file,blzeq.m%The value of the global parameters ar

3、eglobal SIGMA RHO BETASIGMA=10.;RHO=10.;BETA=20.;%The graphics axis limits are set to values known to%contain the solution.axis(0 50 -30 30 -30 30)view(3)hold ontitle('Lorenz Attractor')y0=0 0 eps;tfinal=100;t,y=ode23('blzeq',0 tfinal,y0);plot3(y(:,1),y(:,2),y(:,3)blzeq函数程序%BLZEQ Equ

4、ation of the Lorenz chaotic attractor%ydot=lorenzeq(t,y)%The differential equation is written in almost linear form.function ydot=blzeq(t,y);global SIGMA RHO BETAA=-BETA 0 y(2) 0 -SIGMA SIGMA -y(2) RHO -1;ydot=A*y;实验结果及其分析：题中的方程与程序中的方程的关系是变量进行了轮换，x换成了y，y换成了z，z换成了x。原点为原方程的一个奇点，当初始位置稍稍偏离原点如取为0，eps，0，（

5、按原方程中的顺序，下同）得到的图像如下：分析：这是一个典型的奇怪吸引子的图像，曲线有界，但他不收敛于某一点也不是周期的，而是在两个位置附近来回的跳跃。取初始位置分别为10，10，10，50，50，50得到的图像如下：分析：个初始变量值相同时，曲线总是被吸引回奇怪吸引子附近作来回跳跃，在初始变量值取道200，200，200，-200，-200，-200时，依然如此。下面分别考察初始值的每个分量变化对图像的影响：y分量：0，3，0 0，7，00，7.69516，0 0，7.69517，00，14.23，0 0，14.24，0分析：从上面可以看见，随着初始x值的增大，奇怪吸引子中曲线在其附近来回跳跃

6、的两个位置中的一个吸引力变弱，另一个吸引力变强，然后在初始x取某一特定值时，一个位置丧失吸引力，另一位置则将曲线完全吸引过来变成普通吸引子。初始x继续增大到某一特定值，情况又会变回来。对初始x值单独变化的情况也有类似的现象。这所明在空间存在一些区域，当初始位置位于这些区域外时解将出现奇怪吸引子的性质，而在这些区域以内解将呈现普通吸引子的性质。Z分量：0，0，20分析：从图上可以看出解的曲线为一直线，这可以从方程的角度来解释。当x=0,y=0时在方程中dx/dt=0,dy/dt=0，x，y 方向的值不发生变化，仅z方向的值变化，因此解为一直线。当SIGMA=10，RHO=20，BETA=10，初

7、始值取0，eps，0时，图像如下：对初值进行调整没有发现奇怪吸引子的出现。只调整BETA变量情况如下：SIGMA=10，RHO=28，BETA=9.6/3 0，eps，0SIGMA=10，RHO=28，BETA=9.7/3 0，eps，0SIGMA=10，RHO=28，BETA=15/3 0，eps，0改变SIGMA、RHO的值也有类似的现象。实验6.2（刚性问题）问题提出：考虑下面的初值问题，其中实验内容：刚性比是衡量问题困难程度的重要指标，针对问题合理选择求解刚性问题的方法很重要，Matlab中提供了丰富的函数求解刚性方程，请读者尝试不同方法求解上述方程。实验要求：（1）用Runge-Ku

8、tta算法求解方程。（2）分别用刚性问题的算法和一般问题的算法求解方程，与（1）比较他们的计算结果和计算时间，并分析它们的精度。（3）在t=0.1,0.5,1,2,4,8,10等处，计算解的刚性比。（4）尝试编写隐式Runge-Kutta算法和非线性方法的程序，计算方程的解并与前面的计算结果相比较。程序清单：用龙格库塔法求解方程的程序如下：主程序：clfclcaxis(-0.5 1.5 0.8 2.2, -0.4 0.4)view(3)hold ony0=1 1 0;tfinal=10;tic;t,y=ode23('func1',0 tfinal,y0);tocplot3(y(

9、:,1),y(:,2),y(:,3)函数func1.m:function ydot=func1(t,y);global SIGMA RHO BETAA=-0.013 -1000*y(1) 0 0 -2500*y(3) 0 -0.013 -1000*y(1) -2500*y(2);ydot=A*y;用Matlab中专门处理刚性问题的算法求解方程的程序如下：clfclcaxis(-0.5 1.5 0.8 2.2, -0.4 0.4)view(3)hold ony0=1 1 eps;tfinal=10;tic;t,y=ode23s('func1',0 tfinal,y0)tocplo

10、t3(y(:,1),y(:,2),y(:,3)自编的隐式龙格-库塔程序如下：%一阶隐式龙格-库塔法clfclcaxis(-0.5 1.5 0.8 2.2 -0.4 0.4)view(3)hold onh=0.0001;y(1,:)=1 1 0;for i=1:1/h; z1=y(i,:);%yn z2=z1;%设置z2为yn+1的迭代初值 F=z1(1)-z2(1)-h*(0.0065*z1(1)+0.0065*z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2) z1(2)-z2(2)-625*h*(z1(2)+z2(2)*(z1(3)+z2(3) z1(3)-z2(3)-

11、h*(0.0065*(z1(1)+z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2)+625*(z1(2)+z2(2)*(z1(3)+z2(3); F1=-1-0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1) 0 0 -1-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2) -0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1)-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2)-1; z3=z2-(inv(F1)*F)'%迭代

12、求解 while norm(z3-z2)>10; z2=z3; F=z1(1)-z2(1)-h*(0.0065*z1(1)+0.0065*z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2) z1(2)-z2(2)-625*h*(z1(2)+z2(2)*(z1(3)+z2(3) z1(3)-z2(3)-h*(0.0065*(z1(1)+z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2)+625*(z1(2)+z2(2)*(z1(3)+z2(3); F1=-1-0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1

13、) 0 0 -1-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2) -0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1)-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2)-1; z3=z2-(inv(F1)*F)'end; y(i+1,:)=z3;end;plot3(y(:,1),y(:,2),y(:,3);tfinal=10;y0=y(1,:)自编的龙格-库塔型非线性法的程序如下：%龙格-库塔型非线性法clfclcaxis(-0.5 1.5 0.8 2.2 -0.4 0.4)view(

14、3)hold onh=0.0001;y(1,:)=1 1 0;for i=1:1/h; F=-0.013*y(i,1)-1000*y(i,1)*y(i,2) -2500*y(i,2)*y(i,3) -0.013*y(i,1)-1000*y(i,1)*y(i,2)-2500*y(i,2)*y(i,3);%求导数值x=y(i,1)+0.5*h*y(i,1)*F(1)/(y(i,1)-0.5*h*F(1) y(i,2)+0.5*h*y(i,2)*F(1)/(y(i,1)-0.5*h*F(2) y(i,3)+0.5*h*y(i,3)*F(3)/(y(i,3)-0.5*h*F(3);F1=-0.013*

15、x(1)-1000*x(1)*x(2) -2500*x(2)*x(3) -0.013*x(1)-1000*x(1)*x(2)-2500*x(2)*x(3);y(i+1,1)=y(i,1)+h*F1(1);y(i+1,2)=y(i,2)+h*F1(2);y(i+1,3)=y(i,3)+h*F1(3);end;plot3(y(:,1),y(:,2),y(:,3);y0=y(1,:)tic;t,y=ode23s('func1',0 10,y0)tocplot3(y(:,1),y(:,2),y(:,3),'r')实验结果及其分析：（1） 用matlab中的ode23求解的结果如下：占用的CPU时间为：6.9680s（2）用matlab中求解刚性问题的函数ode23s求解的结果如下（将它与ode23算得的结果画在一起了）占用的CPU时间