数值分析报告Lorenz问题与混沌Word文档格式.docx

1、是否发现什么不同的现象。程序清单：主程序：%BLZ Plot the orbit around the Lore nz chaotic attractorelfclc%Solve the ordinary differe ntial equati on describ ing the%Lore nz chaotic attractor.The equati on are difi ned in%an M-file,blzeq.m%The value of the global parameters areglobal SIGMA RHO BETASIGMA=10.;RHO=10.;BETA=2

2、0.;%The graphics axis limits are set to values known to%contain the solutio n.axis（0 50 -30 30 -30 30）view（3）hold ontitle（Lore nz Attractor）y0=0 0 eps;tfin al=100;t,y=ode23（blzeq,0 tfinal,y0）;plot3（y（:,1）,y（:,2）,y（:,3）blzeq函数程序%BLZEQ Equati on of the Lore nz chaotic attractor %ydot=lore nzeq（t,y）%Th

3、e differential equation is written in almost linear form. function ydot=blzeq（t,y）;A=-BETA 0 y（2）0 -SIGMA SIGMA-y（2） RHO -1;ydot=A*y;实验结果及其分析：题中的方程与程序中的方程的关系是变量进行了轮换， x换成了 y, y换成了 z , z换成了 x。原点为原方程的一个奇点，当初始位置稍稍偏离原点如取为 0，eps，0，（按原方程中的顺序，下同）得到的图像如下：取初始位置分别为10，10，10，50，50，50得到的图像如下:分析：个初始变量值相同时，曲线总是被吸引

4、回奇怪吸引子附近作来回跳跃， 在初始变量值取道200, 200, 200，-200，-200, -200时，依然如此。下面分别考察初始值的每个分量变化对图像的影响：y分量： 0，3, 0 0 ，7，0从上面可以看见，随着初始x值的增大，奇怪吸引子中曲线在其附近来回跳跃的两个 位置中的一个吸引力变弱， 另一个吸引力变强， 然后在初始x取某一特定值时，一个位置丧失吸引力，另一位置则将曲线完全吸引过来变成普通吸引子。 初始x继续增大到某一特定值，情况又会变回来。对初始x值单独变化的情况也有类似的现象。 这所明在空 间存在一些区域，当初始位置位于这些区域外时解将出现奇怪吸引子的性质， 而在这些 区域以

5、内解将呈现普通吸引子的性质。从图上可以看出解的曲线为一直线， 这可以从方程的角度来解释。 当x=O,y=O时在方程中dx/dt=O,dy/dt=O ，x，y方向的值不发生变化，仅 z方向的值变化，因此解为一直线。对初值进行调整没有发现奇怪吸引子的出现。只调整BETA变量情况如下：SIGMA=10 RHO=28 BETA=9.6/3 0 ，eps，030、20 r10、0.dO -2030 1SIGMA=10 RHO=28 BETA=15/3 0 , eps, 0Larenz Aitiractar改变SIGMA RHO勺值也有类似的现象。实验6.2 （刚性问题）考虑下面的初值问题 r 0,10,

6、dy = -0.013y1000y1y2* dy = -2500y2y3dy3 = 0.013y1 -1000叶 y2 2500y2y3其中 y（t） =（y1（t）,y2（t）3（t）T,y（0） =（1,1,0）T实验内容：刚性比是衡量问题困难程度的重要指标， 针对问题合理选择求解刚性问题的方法很重要，Matlab中提供了丰富的函数求解刚性方程，请读者尝试不同方法求解上述方 程。（1）用 Runge-Kutta算法求解方程。（2） 分别用刚性问题的算法和一般问题的算法求解方程， 与（1）比较他们的计算结果和计 算时间，并分析它们的精度。（3） 在t=0.1,0.5,1,2,4,8,10 等

7、处，计算解的刚性比。（4） 尝试编写隐式 Runge-Kutta算法和非线性方法的程序，计算方程的解并与前面的计算 结果相比较。用龙格库塔法求解方程的程序如下：axis（-0.5 1.5 0.8 2.2, -0.4 0.4）y0=1 1 0;tfin al=10;tic;fu nc1,0 tfin al,yO）;toc函数 func1.m:function ydot=fu nc1（t,y）;A=-0.013 -1000*y（1） 00 -2500*y（3） 0-0.013 -1000*y（1） -2500*y （2）;用Matlab中专门处理刚性问题的算法求解方程的程序如下 ：clfy0=1

8、1 eps;t,y=ode23s（func1,0 tfinal,y0）自编的隐式龙格-库塔程序如下：9一阶隐式龙格-库塔法axis（-0.5 1.5 0.8 2.2 -0.4 0.4）h=0.0001;y（1,：）=1 1 0;for i=1:1/h;z仁 y（i,:）;%ynz2=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）-h*（0.0065*（z

9、1（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） 00 -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-（i nv（F1）*F）;%迭弋求解while n orm（z3-

10、z2）10;z2=z3;en d;y（i+1,:）=z3;,3）;y0=y（1,:自编的龙格-库塔型非线性法的程序如下：%龙格-库塔型非线性法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

11、,3）-0.5*h*F（3）;F1=-0.013*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）;,0 10,y0）,3）,r占用的CPU时间为：6.9680s（2）用matlab中求解刚性问题的函数 ode23s求解的结果如下（将它与 ode23算得的结果画在一起了）0.1870s，这比用ode23要快得多。而从图上可以看出两者的曲线基本 一致。（2） 解在各时刻的刚性比为:t=0.0001 : 193.7t=0.1,0.5,1,2,4,8,10 : 2.5计算结果表明在t=0到0.01这段0.03这段时间内刚性较大， 解的也变化较大，而在这段时间以外解的变化非常的缓慢。（3） 用自编的一阶隐式龙格-库塔算法，取步长为 0.0001，从t=0算到t=1算得的结果2.016s。用自编的一阶龙格-库塔型非线性算法，取步长为0.0001，从t=0算到t=1算得的结果如下:占用的CPU寸间为:1.2660s.