数学建模实验三--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 f

h x f x f h )()(lim )(0-+='→，即微商h

x f h x f x f )()()(-+≈'。（右端称为差商）从而可在每个结点上用差商来近似替代导数，将微分方程),(y x f y ='转化为代数方程组（此处的代数方程组常称为差分方程）

))(,()()(k k k k x y x f h

x 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 ，则有

⎪⎩⎪⎨⎧+=+==++,),(,,)(1

100k 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、求解微分方程2

2x 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 dx

dy 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 dt

dy 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 的变化关系。