龙格—库塔法分析lorenz方程课程设计.doc

《数值分析》课程设计实验报告
龙格—库塔法分析Lorenz方程

问题叙述
考虑著名的Lorenz方程
其中s,r,b为变化区域内有一定限制的实参数,该方程形式简单,表面上看并无惊人之处,但由该方程揭示出的许多现象,促使“混沌”成为数学研究的崭新领域,在实际应用中也产生了巨大的影响。
问题分析
Lorenz方程实际上是一个四元一阶常微分方程,用解析法精确求解是不可能的,只能用数值计算,最主要的有欧拉法、亚当法和龙格- 库塔法等。为了得到较高精度的,我们采用经典四阶龙格—库塔方法求解该问题。
实验程序及注释
(1)算法程序
function [T]=Runge_Kutta(f,x0,y0,h,n) %定义算法,其中f为待解方程组,x0是初始自变量,y0是初始函数值,h是步长,n为步数
if nargin<5
n=100; %如果输入参数个数小于5,则步数n=100
end
r=size(y0);r=r(1); %返回初始输出矩阵的行列数,并将值赋给r(1)
s=size(x0);s=s(1); %返回初始输入矩阵的行列数,并将值赋给s(1)
r=r+s;
T=zeros(r,n+1);
T(:,1)=[y0;x0];
for t=2:n+1 %以下是具体的求解过程
k1=feval(f,T(1:r-1,t-1));
k2=feval(f,[k1*(h/2)+T(1:r-1,t-1);x0+h/2]);
k3=feval(f,[k2*(h/2)+T(1:r-1,t-1);x0+h/2]);
k4=feval(f,[k3*h+T(1:r-1,t-1);x0+h]);
x0=x0+h;
T(:,t)=[T(1:r-1,t-1)+(k1+k2*2+k3*2+k4)*(h/6);x0];
end
(2)主程序
function dy=fun(x) %定义函数
s=; %给参数s,r,b赋值
r=;
b=;
dy(1)=s*(x(2)-x(1)); %Lorenz方程表达式
dy(2)=(r*x(1)-x(3)*x(1)-x(2));
dy(3)=x(1)*x(2)-b*x(3);
dy=dy';
(2)运行程序
T=Runge_Kutta('fun',0,[10;10;10],,5000); %调用前面的算法程序
plot3(T(1,:),T(2,:),T(3,:)); %显示三分量的关系图
axis([-20 20 -50 50 0 50]) %定义坐标轴长度
view(3) %设定观察角度
实验数据结果及分析
(1)各初始变量相同时的图像分析
各初始变量取相同的值[10,10,10],运行上述程序后,得到如下图像:
从图中可以看出,各初始变量相同时,曲线总是被吸引回奇怪吸引子附近作来回跳跃。初始变量值取为[-10,-10,-10] ,[20,20,20]时,依然如此。图像如下:
[-10,-10,-10] [20,20,20]
(2)初始值的每个分量变化对图像的影响
y分量:
[0,2,0] [0,5,0]
[0,15,0] [0,20,0]
从上面可以看出,随着初始y值的增大,奇怪吸引子中曲线在其附近