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

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

一、实验目的

1、学习用Mathematica 求常微分方程的解析解与数值解,并进行定性分析;

2、学习用MA TLAB 求常微分方程的解析解与数值解,并进行定性分析。

二、实验材料

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、微分方程的MA TLAB 求解

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