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

数学建模实验三

Lorenz 模型与食饵模型

一、 实验目的

1、 学习用Mathematica 求常微分方程的解析解和数值解，并进行定性分析；

2、 学习用MATLAB 求常微分方程的解析解和数值解，并进行定性分析。

二、 实验材料

2.1问题

图3.3.1是著名的洛仑兹(E.N.Lorenz)混沌吸引子，洛仑兹吸引子已成为混沌理论的徽

标，好比行星轨道图代表着哥白尼、开普勒理论一样。洛仑兹是学数学出身的， 1948年起

在美国麻省理工学院(MIT )作动力气象学博士后工作，

1963年他在《大气科学杂志》上

发表的论文《确定性非周期流》是混沌研究史上光辉的著作。以前科学家们不自觉地认为微 分方程的解只有那么几类：1)发散轨道；2)不动点；3)极限环；4)极限环面。除此以外，大 概没有新的运动类型了，

这是人们的一种主观猜测，谁也没有给出证明。事实上这种想法是

非常错误的。1963年美国麻省理工学院气象科学家洛仑兹给出一个具体模型，就是著名的 Lorenz 模型，清楚地展示了一种新型运动体制：混沌运动，轨道既不收敛到极限环上也不 跑掉。而今Lorenz 模型在科学与工程计算中经常运用的问题。例如，数据加密中。我们能 否绘制出洛仑兹吸引子呢？

假设狐狸和兔子共同生活在同一个有限区域内， 有足够多的食物供兔子享用， 而狐狸仅

以兔子为食物.X 为兔子数量，y 表狐狸数量。假定在没有狐狸的情况下，兔子增长率为400%。 如果

没有兔子，狐狸将被饿死，死亡率为 90%。狐狸与兔子相互作用的关系是，狐狸的存

在使兔子受到威胁，且狐狸越多兔子增长受到阻碍越大，设增长的减小与狐狸总数成正比， 比例系数

为0.02。而兔子的存在又为狐狸提供食物， 设狐狸在单位时间的死亡率的减少与兔

子的数量成正比，设比例系数为 0.001。建立数学模型，并说明这个简单的生态系统是如何

变化的。

2.2预备知识

1、求解常微分方程的 Euler 折线法

求初值问题

= f (x,y), • yX) =y

。

(12.1 )

图3.3.1洛仑兹(E.N.Lorenz)混沌吸引子

在区间[X o

,

X n ]上的数值解，并在区间插入了结点 (X o ：：： )

X i

：：：…：:：X n 」(：:：X n

)。由导数

的定义 f (x)二lim 丄乜一h)——，即微商 f (x)：、丄纟一h)——f

(X )

。(右端称为差商)

h

-h

h

从而可在每个结点上用差商来近似替代导数，将微分方程 f (X ,y)转化为代数方程组

(此处的代数方程组常称为差分方程)

—

y （X k ，y （X k ））,k",i, ,n_i h

加上初值条件则可确定一组解。求解这一差分方程即可得到微分方程初值问题的数值解。 变

形上述方程有

y(x k h)二 y(xQ hf (x k , y(xj) ,k=0,1「，n-1

记 x k i = x k h , y(X k ) = y k ，从而 y(X k • h) = y k 1，则有

y o =y(x o ),

<

X2

=X k +h , k =0,1，…,n -1

y k 1 二 y k

hf (X k ,y k ),

这就是求解微分方程初值问题的欧拉

(Euler)折线法。之所以称为欧拉折线法是因为：

就几何

角度而言，所求得的近似解是初值问题精确解的折线逼近， 而且此折线的起点是初值条件所

对应的点。

2、微分方程的Mathematica 求解

(1) 求解命令

有两个命令：DSolve[]与NDSolve 。命令格式分别为 DSolve[方程，y , x]

NDSolve [方程，y , { x , xl , x2门。

其中方程必须为微分方程及相应初始条件， {x , xl , x2 }说明要给出数值解的范围为区

间]x1, x2]o

(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 模型的状态方

程可表示为

{{丫 十 Function [ {K J-

X1（t）二-X（t）X2（t）X3（t）

X2（t）「-；「X2（t）；「X3（t）

X3（t）二—X1（t）X2（t）「X2（t）—X3（t）

若令；「= 10, Q = 28, ：=8/3,且初值为x,。)= x2(0) = 0, x3 (0)=；,；为一个小常

数，假设；=10」°。求微分方程的数值解，并绘制出时间曲线与相空间曲线。

问题(2)是著名的食饵模型，数学模型为

丈=4x-0.02xy

y = —0.9y +0.001xy

2.4练习题

1、求解微分方程、 2xy =xe'的通解。

求解的Mathematica命令为：

DSolve[y'[x]+2*x*y[x]== x*E A(-x A2),y,x] 或者

DSolve[D[y [ x],x]+2*x*y[x]== x*EA(-xA2),y,x]

2、求微分方程xy： y -e x=0在初始条件y = 2e下的特解。应给出的命令为：

DSolve[{x*y'[x]+ y[ x]-EAx==0,y[1]==2E},y,x]

3、求(x2 -1)dy - 2xy-cosx =0在初始条件y(0) = 1下的特解，并画出解的图形。

dx

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

清除要涉及变量的命令为：

Clear[x,y]

求解析解的命令为：

sc=DSolve[{(xA2-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[{(xA2-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、求微分方程组

虫+5x + y = d ,

jdt

—”0

-dt

在初始条件x(0) =1,y(0) =0下的解，并画出解函数y二y(x)的图形。

求解微分方程组的命令为：

Clear[x,y,t] xy=DSolve[{x'[t]+5*x[t]+y[t==EAt,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}}]