龙格—库塔法分析Lorenz方程课程设计推荐.docVIP

《数值分析》课程设计实验报告 龙格—库塔法分析Lorenz方程 问题叙述 考虑著名的Lorenz方程 其中s，r，b为变化区域内有一定限制的实参数，该方程形式简单，表面上看并无惊人之处，但由该方程揭示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。 问题分析 Lorenz方程实际上是一个四元一阶常微分方程，用解析法精确求解是不可能的，只能用数值计算，最主要的有欧拉法、亚当法和龙格- 库塔法等。为了得到较高精度的，我们采用经典四阶龙格—库塔方法求解该问题。 实验程序及注释 （1）算法程序 function [T]=Runge_Kutta(f,x0,y0,h,n) ％定义算法，其中f为待解方程组，x0是初始自变量，y0是初始函数值，h是步长，n为步数 if nargin5 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=10.0; ％给参数s，r，b赋值 r=28.0; b=8.0/3; 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],0.01,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值的增大，奇怪吸引子中曲线在其附近来回跳跃的两个位置中的一个吸引力变弱，另一个吸引力变强。初始y继续增大到某一特定值，情况又会变回来。这说明在空间存在一些区域，当初始位置位于这些区域外时解将出现奇怪吸引子的性质，而在这些区域以内解将呈现普通吸引子的性质。 z分量：[0，0，20] 从上图可以看出解的曲线为一直线，这可以从方程的角度来解释。当x=0,y=0时在方程中dx/dt=0,dy/dt=0，x，y 方向的值不发生变化，仅z方向的值变化，因此