数学建模实验三--Lorenz模型与食饵模型

数学建模实验三 Lorenz 模型与食饵模型一、实验目的1、学习用Mathematica 求常微分方程的解析解和数值解，并进行定性分析；2、学习用MATLAB 求常微分方程的解析解和数值解，并进行定性分析。

二、实验材料2.1问题图3.3.1是著名的洛仑兹(E.N.Lorenz)混沌吸引子，洛仑兹吸引子已成为混沌理论的徽标，好比行星轨道图代表着哥白尼、开普勒理论一样。

洛仑兹是学数学出身的，1948年起在美国麻省理工学院（MIT ）作动力气象学博士后工作，1963年他在《大气科学杂志》上发表的论文《确定性非周期流》是混沌研究史上光辉的著作。

以前科学家们不自觉地认为微分方程的解只有那么几类：1)发散轨道；2)不动点；3)极限环 ；4)极限环面。

除此以外，大概没有新的运动类型了，这是人们的一种主观猜测，谁也没有给出证明。

事实上这种想法是非常错误的。

1963年美国麻省理工学院气象科学家洛仑兹给出一个具体模型，就是著名的Lorenz 模型，清楚地展示了一种新型运动体制：混沌运动，轨道既不收敛到极限环上也不跑掉。

而今Lorenz 模型在科学与工程计算中经常运用的问题。

例如，数据加密中。

我们能否绘制出洛仑兹吸引子呢？图3.3.1 洛仑兹(E.N.Lorenz)混沌吸引子假设狐狸和兔子共同生活在同一个有限区域内，有足够多的食物供兔子享用，而狐狸仅以兔子为食物.x 为兔子数量,y 表狐狸数量。

假定在没有狐狸的情况下，兔子增长率为400％。

如果没有兔子，狐狸将被饿死，死亡率为90％。

狐狸与兔子相互作用的关系是，狐狸的存在使兔子受到威胁，且狐狸越多兔子增长受到阻碍越大，设增长的减小与狐狸总数成正比，比例系数为0.02。

而兔子的存在又为狐狸提供食物，设狐狸在单位时间的死亡率的减少与兔子的数量成正比，设比例系数为0.001。

建立数学模型，并说明这个简单的生态系统是如何变化的。

2.2预备知识1、求解常微分方程的Euler 折线法求初值问题⎩⎨⎧=='00)(),,(y x y y x f y （12.1）在区间],[0n x x 上的数值解，并在区间插入了结点)()(110n n x x x x <<<<-Λ。

由导数的定义h x fh x f x f h )()(lim )(0-+='→，即微商hx f h x f x f )()()(-+≈'。

（右端称为差商）从而可在每个结点上用差商来近似替代导数，将微分方程),(y x f y ='转化为代数方程组（此处的代数方程组常称为差分方程）))(,()()(k k k k x y x f hx y h x y =-+,1,,1,0-=n k Λ 加上初值条件则可确定一组解。

求解这一差分方程即可得到微分方程初值问题的数值解。

变形上述方程有))(,()()(k k k k x y x hf x y h x y +=+,1,,1,0-=n k Λ记h x x k k +=+1,k k y x y =)(，从而1)(+=+k k y h x y ，则有⎪⎩⎪⎨⎧+=+==++,),(,,)(1100k k k k k k y x hf y y h x x x y y 1,,1,0-=n k Λ这就是求解微分方程初值问题的欧拉(Euler)折线法。

之所以称为欧拉折线法是因为：就几何角度而言，所求得的近似解是初值问题精确解的折线逼近，而且此折线的起点是初值条件所对应的点。

2、微分方程的Mathematica 求解（1）求解命令有两个命令：DSolve[ ]与NDSolve 。

命令格式分别为DSolve[方程，y ，x]NDSolve ［方程，y ，｛x ，xl ，x2｝］。

其中方程必须为微分方程及相应初始条件，｛x ，xl ，x2｝说明要给出数值解的范围为区间［x1，x2］。

（2）使用的注意事项①方程中的函数应写成完整形式y[x],以表明y 是x 的函数；②方程应写成…==…的形式；③重复使用时，应随时清除要涉及变量的以前定义，方法是Clear[y]；④使用NDSolve 时，所加初始条件的个数应等于微分方程的阶数，同时方程中也不含其它参数，否则给不出正确结果。

（3）解的表示形式Mathematica 给出的微分方程的解是以纯函数（或数学中的算子）定义的形式给出的，例如：DSolve[y'[x]+ 3*y[x]==2x,y,x]的结果是3、微分方程的MATLAB 求解（1）求解析解命令dsolve ；（2）求数值解命令ODE 或 Simulink 。

2.3建立模型问题（1）的洛仑兹吸引子可以用下面的微分方程得到，著名的Lorenz 模型的状态方程可表示为⎪⎩⎪⎨⎧-+-=+-=+-=)()()()()()()()()()()()(322133223211t x t x t x t x t x t x t x t x t x t x t x t x ρσσβ&&&若令，，，3/82810===βρσ 且初值为ε===)0(0)0()0(321x x x ，， 为一个小常数，假设1010-=ε。

求微分方程的数值解，并绘制出时间曲线与相空间曲线。

问题（2）是著名的食饵模型，数学模型为⎩⎨⎧+-='-='xy y y xy x x 001.09.002.04 2.4练习题1、求解微分方程22x xe xy y -=+'的通解。

求解的Mathematica 命令为：DSolve[y'[x]+2*x*y[x]== x*E^(-x^2),y,x] 或者DSolve[D[y[x],x]+2*x*y[x]== x*E^(-x^2),y,x]2、求微分方程0=-+'x e y y x 在初始条件e y x 21==下的特解。

应给出的命令为：DSolve[{x*y'[x]+ y[x]-E^x==0,y[1]==2E},y,x]3、求0cos 2)1(2=-+-x xy dxdy x 在初始条件1)0(=y 下的特解，并画出解的图形。

要求分别求解析解与数值解并作比较。

清除要涉及变量的命令为：Clear[x,y]求解析解的命令为：sc=DSolve[{(x^2-1)y'[x]+2x*y[x]-Cos[x]==0,y[0]==1},y,x]画解析解图像的命令为：y=y/.sc[[1]]g1=Plot[y[x],{x,0,1},PlotStyle->RGBColor[1,0,0]]注：也可将画图范围变为Plot[y[x],{x,0,4}]求数值解的命令为：sn=NDSolve[{(x^2-1)y'[x]+2x*y[x]-Cos[x]==0,y[0]==1}, y,{x,0,1}]画数值解图像的命令为：y=y/.sn[[1]]g2=Plot[y[x],{x,0,1}]比较解析解图像与数值解图像的命令为：Show[g1,g2]4、求微分方程组 ⎪⎩⎪⎨⎧=--=++03,5y x dtdy e y x dt dx t 在初始条件1)0(=x ,0)0(=y 下的解，并画出解函数)(x y y =的图形。

求解微分方程组的命令为：Clear[x,y,t]xy=DSolve[{x'[t]+5*x[t]+y[t==E^t,y'[t]-x[t]-3*y[t]==0,x[0]==1,y[0]==0},{x,y},t]画解的相位图的命令为：y=y/.xy[[1]];x=x/.xy[[1]];ParametricPlot[{x[t],y[t]},{t,0,3},PlotRange->{{-10,2},{0,5}}]注：图中反应出y 随x 的变化关系。

三、实验准备认真阅读实验目的与实验材料后要正确地解读实验，在此基础上制定实验计划（修改、补充或编写程序，提出实验思路，明确实验步骤），为上机实验做好准备。

四、实验思路提示4.1实验步骤1、求解问题（2）中的食饵模型的微分方程组，并画出解的图形和相位图。

（1）以x=800，y=100为初始值，计算x（t）,y（t），当t [0，14]时的数据。

绘出解的图形，并分析捕食者和被捕食者的数量变化规律。

可以先用下面的命令求解析解：Clear[x,y,t]xy=DSolve[{x'[t]==4*x[t]-0.02*x[t]*y[t],y'[t]==-0.9*y[t]+0.001*x[t]*y[t],x[0]==800,y[0]==100},{x,y},t]注：可以发现不能求出解析解。

修改代码如下，可以求数值解：Clear[x,y,t]xy=NDSolve[{x'[t]==4*x[t]-0.02*x[t]*y[t],y'[t]==-0.9*y[t]+0.001*x[t]*y[t],x[0]==800,y[0]==100},{x,y},{t,0,14}]绘出解的图形：y=y/.xy[[1]];x=x/.xy[[1]];Plot[{x[t],y[t]},{t,0,14},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]图3.3.2 捕食者和被捕食者的数量变化（2）以x为横坐标，y为纵坐标绘制相位图。

根据图形分析被捕食者数量增加(减少)对捕食者数量的影响。

绘制相位图的命令：ParametricPlot[{x[t],y[t]},{t,0,14}]图3.3.3 相位图2、用MATLAB求解问题（1）中Lorenz 模型的微分方程。

（1）打开MATLAB的编辑器；（2）在编辑器中用下面的几个语句描述微分方程，并将其保存在lorenzeq.m的m文件中：f unction xdot = lorenzeq(t,x)xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)];（3）新建命令文件：t_final=100; x0=[0;0;1e-10];[t,x]=ode45('lorenzeq',[0,t_final],x0);plot(t,x),figure; plot3(x(:,1),x(:,2),x(:,3)); axis([10 40 -20 20 -20 20]);绘制出时间曲线与相空间曲线，如下图所示。

图3.3.4时间曲线与相空间曲线4.2思考问题1、运用Mathematica求解Lorenz 模型的微分方程组，从而了解系统状态是如何变化的。

2、求解以下问题（广告的效用）：某公司生产一种耐用消费品，产品一上市，该公司即开始做广告，一段时期的市场跟踪调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有购买的百分比成正比，且估得此比例系数为0.5。