数值分析报告Lorenz问题与混沌Word文档格式.docx
《数值分析报告Lorenz问题与混沌Word文档格式.docx》由会员分享，可在线阅读，更多相关《数值分析报告Lorenz问题与混沌Word文档格式.docx（14页珍藏版）》请在冰豆网上搜索。
是否发现什么不同的现
象。
程序清单：
主程序：
%BLZPlottheorbitaroundtheLorenzchaoticattractor
elf
clc
%Solvetheordinarydifferentialequationdescribingthe
%Lorenzchaoticattractor.Theequationaredifinedin
%anM-file,blzeq.m
%Thevalueoftheglobalparametersare
globalSIGMARHOBETA
SIGMA=10.;
RHO=10.;
BETA=20.;
%Thegraphicsaxislimitsaresettovaluesknownto
%containthesolution.
axis（[050-3030-3030]）
view（3）
holdon
title（'
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）;
A=[-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]
从上面可以看见，随着初始x值的增大，奇怪吸引子中曲线在其附近来回跳跃的两个位置中的一个吸引力变弱，另一个吸引力变强，然后在初始x取某一特定值时，一个位
置丧失吸引力，另一位置则将曲线完全吸引过来变成普通吸引子。
初始x继续增大到某
一特定值，情况又会变回来。
对初始x值单独变化的情况也有类似的现象。
这所明在空间存在一些区域，当初始位置位于这些区域外时解将出现奇怪吸引子的性质，而在这些区域以内解将呈现普通吸引子的性质。
从图上可以看出解的曲线为一直线，这可以从方程的角度来解释。
当x=O,y=O时在方
程中dx/dt=O,dy/dt=O，x，y方向的值不发生变化，仅z方向的值变化，因此解为一
直线。
对初值进行调整没有发现奇怪吸引子的出现。
只调整BETA变量情况如下：
SIGMA=10RHO=28BETA=9.6/3[0，eps，0]
30、
20r
10、
0.
dO-
-20
■301
SIGMA=10RHO=28BETA=15/3[0,eps,0]
LarenzAitiractar
改变SIGMARHO勺值也有类似的现象。
实验6.2（刚性问题）
考虑下面的初值问题r[0,10],
dy^=-0.013y^1000y1y2
*dy^=-2500y2y3
dy3
~^=—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等处，计算解的刚性比。
（4）尝试编写隐式Runge-Kutta算法和非线性方法的程序，计算方程的解并与前面的计算结果相比较。
用龙格库塔法求解方程的程序如下：
axis（[-0.51.50.82.2,-0.40.4]）
y0=[110];
tfinal=10;
tic;
func1'
[0tfinal],yO）;
toc
函数func1.m:
functionydot=func1（t,y）;
A=[-0.013-1000*y
（1）0
0-2500*y（3）0
-0.013-1000*y
（1）-2500*y
（2）];
用Matlab中专门处理刚性问题的算法求解方程的程序如下：
clf
y0=[11eps];
[t,y]=ode23s（'
func1'
[0tfinal],y0）
自编的隐式龙格-库塔程序如下：
9一阶隐式龙格-库塔法
axis（[-0.51.50.82.2-0.40.4]）
h=0.0001;
y（1,：
）=[110];
fori=1:
1/h;
z仁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）-h*（0.0065*（z1
（1）+z2
（1））+250*（z1
（1）+z2
（1））*（z1
（2）+z2
（2））+625*（z
1
（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）'
;
%迭弋求解
whilenorm（z3-z2）>
10;
z2=z3;
end;
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,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）;
[010],y0）
3）,'
r'
占用的CPU时间为：
6.9680s
（2）用matlab中求解刚性问题的函数ode23s求解的结果如下（将它与ode23算得的结果
画在一起了）
0.1870s，这比用ode23要快得多。
而从图上可以看出两者的曲线基本一致。
（2）解在各时刻的刚性比为:
t=0.0001:
193.7
t=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.
升级会员 