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

数学建模实验三

一、 实验目的

Lorenz模型与食饵模型

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

二、 实验材料

2.1问题

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

标，好比行星轨道图代表着哥白尼、开普勒理论一样。洛仑兹是学数学出身的， 在美国麻省理工学院(MIT )作动力气象学博士后工作，

1948年起

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

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

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

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

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

假设狐狸和兔子共同生活在同一个有限区域内， 没有兔子，狐狸将被饿死，死亡率为 为0.02。而兔子的存在又为狐狸提供食物， 子的数量成正比，设比例系数为 变化的。

有足够多的食物供兔子享用， 而狐狸仅

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

90%。狐狸与兔子相互作用的关系是，狐狸的存

设狐狸在单位时间的死亡率的减少与兔

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

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

2.2预备知识

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

求初值问题

= f (x,y), ? yX) =y。

(12.1 )

1 / 6

在区间[Xo

,

Xn]上的数值解，并在区间插入了结点 ：：：…：:：Xn」(：:：Xn

h

(Xo ：：： )

)。由导数

(X)

Xi

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

。(右端称为差商)

-h

从而可在每个结点上用差商来近似替代导数，将微分方程 (此处的代数方程组常称为差分方程)

h

f (X ,y)转化为代数方程组

y（Xk h） — y（Xk） =f（Xk，y（Xk））,k\

变

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

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

记 xk i = xk h , y(Xk) = yk，从而 y(Xk ? h) = yk 1，则有

yo=y(xo),

<

X2

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

hf (Xk ,yk),

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

就几何

而且此折线的起点是初值条件所

yk 1 二 yk

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

角度而言，所求得的近似解是初值问题精确解的折线逼近， 对应的点。

2、微分方程的Mathematica求解

(1) 求解命令

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

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

其中方程必须为微分方程及相应初始条件， 间]x1, x2]o

(2) 使用的注意事项

{x, xl , x2 }说明要给出数值解的范围为区

① 方程中的函数应写成完整形式 y[x],以表明y是x的函数； ② 方程应写成…==???的形式；

③ 重复使用时，应随时清除要涉及变量的以前定义，方法是

Clear[y];

同时方程中也不含其

④ 使用NDSolve时，所加初始条件的个数应等于微分方程的阶数， 它参数，否则给不出正确结果。

(3) 解的表示形式

Mathematica给出的微分方程的解是以纯函数(或数学中的算子)定义的形式给出的， 例如：DSolve[y'[x]+ 3*y [x]==2x,y,x] 的结果是

{{丫 十 Function [ {KJ-

3、微分方程的MATLAB求解

(1) 求解析解 命令dsolve ;

(2) 求数值解命令 ODE或Simulink。

2.3建立模型

问题(1 )的洛仑兹吸引子可以用下面的微分方程得到，著名的 程可表示为

Lorenz模型的状态方

2 / 6

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*EA(-xA2),y,x]

或者

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

2、 求微分方程xy： y -ex =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下的解，并画出解函数

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

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(x)的图形。

命令为：

y=y/.xy[[1]]; x=x/.xy[[1]];

ParametricPlot[{x[t],y[t]},{t,0,3},PlotRange->{{-10,2},{0,5}}]

3 / 6

注：图中反应出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}]

（减

4 / 6

图3.3.3 相位图

2、用MATLAB求解问题(1 )中Lorenz模型的微分方程。 (1) 打开MATLAB的编辑器；

(2) 在编辑器中用下面的几个语句描述微分方程，并将其保存在 中：

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)];

lorenzeq.m的m文件

(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]);

绘制出时间曲线与相空间曲线，如下图所示。

(曲状态变总的肘间响应團

(b)郴咗间二绯图

图3.3.4时间曲线与相空间曲线

4.2思考问题

1、 运用Mathematica求解Lorenz模型的微分方程组，从而了解系统状态是如何变化的。 2、 求解以下问题(广告的效用)：

某公司生产一种耐用消费品，产品一上市，该公司即开始做广告，一段时期的市场跟踪 调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有购买的百分比 成正比，且估得此比例系数为 (1) (2) (3 80 %?

0.5。

试模拟求解该问题，即购买人口的百分比与(做广告)时间的关系； 建立该问题的数学模型，并求其数值解与模拟结果作以比较； )厂家问：要做多少次广告(设上述单位时间指的是广告次数

X可使市场购买率达到

5 / 6