数值分析Lorenz问题与混沌

991790力92王为5/5实验6.1〔Lorenz问题与混沌问题提出：考虑著名的Lorenz方程其中s,r,b为变化区域有一定限制的实参数,该方程形式简单,表面上看并无惊人之处,但由该方程揭示出的许多现象,促使"混沌"成为数学研究的崭新领域,在实际应用中也产生了巨大的影响。实验要求：〔1请读者找出Lorenz方程与上述程序中使用的方程间的关系。〔2对目前取定的参数值SIGMA、RHO和BETA,选取不同的初始值y0〔当前程序中的y0是在坐标原点,运行上述的程序,观察计算的结果有什么特点？解的曲线是否有界？解的曲线是不是周期的或趋于某个固定的点？〔3在问题允许的范围内适当改变其中的参数值SIGMA、RHO和BETA,再选取不同的初始值y0,运行上述的程序,观察并记录计算的结果有什么特点？是否发现什么不同的现象。程序清单：主程序：%BLZPlottheorbitaroundtheLorenzchaoticattractorclfclc%Solvetheordinarydifferentialequationdescribingthe%Lorenzchaoticattractor.Theequationaredifinedin%anM-file,blzeq.m%ThevalueoftheglobalparametersareglobalSIGMARHOBETASIGMA=10.;RHO=10.;BETA=20.;%Thegraphicsaxislimitsaresettovaluesknownto%containthesolution.axis<[050-3030-3030]>view<3>holdontitle<'LorenzAttractor'>y0=[00eps];tfinal=100;[t,y]=ode23<'blzeq',[0tfinal],y0>;plot3<y<:,1>,y<:,2>,y<:,3>>blzeq函数程序%BLZEQEquationoftheLorenzchaoticattractor%ydot=lorenzeq<t,y>%Thedifferentialequationiswritteninalmostlinearform.functionydot=blzeq<t,y>;globalSIGMARHOBETAA=[-BETA0y<2>0-SIGMASIGMA-y<2>RHO-1];ydot=A*y;实验结果及其分析：题中的方程与程序中的方程的关系是变量进行了轮换,x换成了y,y换成了z,z换成了x。原点为原方程的一个奇点,当初始位置稍稍偏离原点如取为[0,eps,0,]〔按原方程中的顺序,下同得到的图像如下：分析：这是一个典型的奇怪吸引子的图像,曲线有界,但他不收敛于某一点也不是周期的,而是在两个位置附近来回的跳跃。取初始位置分别为[10,10,10],[50,50,50]得到的图像如下：分析：个初始变量值相同时,曲线总是被吸引回奇怪吸引子附近作来回跳跃,在初始变量值取道[200,200,200],[-200,-200,-200]时,依然如此。下面分别考察初始值的每个分量变化对图像的影响：y分量：[0,3,0][0,7,0][0,7.69516,0][0,7.69517,0][0,14.23,0][0,14.24,0]分析：从上面可以看见,随着初始x值的增大,奇怪吸引子中曲线在其附近来回跳跃的两个位置中的一个吸引力变弱,另一个吸引力变强,然后在初始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,初始值取[0,eps,0]时,图像如下：对初值进行调整没有发现奇怪吸引子的出现。只调整BETA变量情况如下：SIGMA=10,RHO=28,BETA=9.6/3[0,eps,0]SIGMA=10,RHO=28,BETA=9.7/3[0,eps,0]SIGMA=10,RHO=28,BETA=15/3[0,eps,0]改变SIGMA、RHO的值也有类似的现象。实验6.2〔刚性问题问题提出：考虑下面的初值问题,其中实验内容：刚性比是衡量问题困难程度的重要指标,针对问题合理选择求解刚性问题的方法很重要,Matlab中提供了丰富的函数求解刚性方程,请读者尝试不同方法求解上述方程。实验要求：〔1用Runge-Kutta算法求解方程。〔2分别用刚性问题的算法和一般问题的算法求解方程,与〔1比较他们的计算结果和计算时间,并分析它们的精度。〔3在t=0.1,0.5,1,2,4,8,10等处,计算解的刚性比。〔4尝试编写隐式Runge-Kutta算法和非线性方法的程序,计算方程的解并与前面的计算结果相比较。程序清单：用龙格库塔法求解方程的程序如下：主程序：clfclcaxis<[-0.51.50.82.2,-0.40.4]>view<3>holdony0=[110];tfinal=10;tic;[t,y]=ode23<'func1',[0tfinal],y0>;tocplot3<y<:,1>,y<:,2>,y<:,3>>函数func1.m:functionydot=func1<t,y>;globalSIGMARHOBETAA=[-0.013-1000*y<1>00-2500*y<3>0-0.013-1000*y<1>-2500*y<2>];ydot=A*y;用Matlab中专门处理刚性问题的算法求解方程的程序如下：clfclcaxis<[-0.51.50.82.2,-0.40.4]>view<3>holdony0=[11eps];tfinal=10;tic;[t,y]=ode23s<'func1',[0tfinal],y0>tocplot3<y<:,1>,y<:,2>,y<:,3>>自编的隐式龙格-库塔程序如下：%一阶隐式龙格-库塔法clfclcaxis<[-0.51.50.82.2-0.40.4]>view<3>holdonh=0.0001;y<1,:>=[110];fori=1:1/h;z1=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*<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>>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-<inv<F1>*F>';%迭代求解whilenorm<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>>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-<inv<F1>*F>';end;y<i+1,:>=z3;end;plot3<y<:,1>,y<:,2>,y<:,3>>;tfinal=10;y0=y<1,:>自编的龙格-库塔型非线性法的程序如下：%龙格-库塔型非线性法clfclcaxis<[-0.51.50.82.2-0.40.4]>view<3>holdonh=0.0001;y<1,:>=[110];fori=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*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',[010],y0>tocplot3<y<:,1>,y<:,2>,y<:,3>,'r'>实验结果及其分析：用matlab中的ode23求解的结果如下：占用的CPU时间为：6.9680s〔2用matlab中求解刚性问题的函数ode23s求解的结果如下〔将它与ode23算得的结果画在一起了占用的CPU时间为：0.1870s,这比用ode23要快得多。而从图上可以看出两