数学模型实验2

数学建模实验
二．微分方程实验
1. 微分方程稳定性分析
绘出下列自治系统相应的轨线，并标出随t 增加的运动方向，确定平衡点，并按稳定的、渐近稳定、或不稳定的进行分类：
（1）⎪⎩⎪⎨⎧==;,y dt
dy x dt dx （2）⎪⎩⎪⎨⎧=-=;2,y dt dy x dt dx （3）⎪⎩⎪⎨⎧-==;2,x dt dy y dt dx （4）⎪⎩⎪⎨⎧-=-=;
2,
y dt dy x dt dx
（1） 选取平衡点，由,0)(,0)(====y y f x x f 可知为（0,0）
系数矩阵为
,1001⎥⎦
⎤
⎢⎣⎡=A 易得特征值
,
121==λλ则
,01.,02)(2121>==<-=+-=λλλλq p 对照稳定性情况表，平衡点是不稳定的。

（2） 根据（1）题所求方法，取平衡点（0,0），易得特征值
,02,01,2,121<-=<-==-=q p λλ对照稳定性情况表，可知平衡点是不稳定的。

（3） 取平衡点（0,0），易得特征值,04142.1,0,4142.1,4142
.121>==-==q p i i λλ对照稳定性情况表，可知平衡点是不稳定的。

（4） 取平衡点（1,0），易得特征值,022,03,2,121<=>=-=-=q p λλ对照稳定性情
况表，可知平衡点是稳定的。

2. 种群增长模型
一个片子上的一群病菌趋向于繁殖成一个圆菌落。

设病菌的数目为N ，单位成员的增长率为r1，则由Malthus 生长率有
.1N r dt
dN
=但是，处于周界表面的那些病菌由于寒冷而受到损伤，它们死亡的数量与N 1/2成比例，其比例系数为让。

求N 满足的微分方程，不用求解，图示其解族，方程是否有平衡接，如果有，是否为稳定？ 解：
由题可得，N 满足的微分方程为：,)(21
21N r N r dt
dN
N f -==
求取平衡点，令0)(=N f ，可得平衡点为（0，2
122/r r ），
由)).(21(212121212
2N r N r N r r dt
N d --=-，令,022
=dt N d 可求得21224/r r N =， 令2
12221224/,/,0r r N r r N N ===把第一象限划分为三部分，且分别有
.0,0;0,0;0,0222222><<>>>dt
N
d dt dN dt N d dt dN dt N d dt dN 则微分方程的解族图形如下所示，其中，0=N 是不稳定的，2
122/r r N =是稳定的。

解：
（1） 构建方程组⎩
⎨⎧--+-=--+-=)]1([)()]
1([)(2211222222111111x b x b c a x x f x b x b c a x x f ，令,0)(,0)(21==x f x f
易得平衡点)0,(),,0(),0,0(1
11
12222210c b a c p c b a c p p -- 对于P 0，系数矩阵
⎥⎦
⎤⎢⎣⎡--=⎥⎦⎤⎢
⎣⎡--------=221
`1221121222211
222
121111100a c a c x c b x c b a c x c b x c b x c b x c b a c A ， )(2211a c a c p -+--=，已知2211,a c a c >>，所以p<0，该点不稳定。

对于P 1，系数矩阵
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡------=⎥⎦⎤⎢⎣
⎡--------=211122221212211221121222211222
1211111))(()(0)(b a c b b b a c b c c a c a c x c h x c b a c x c b x b x c b x c b a c A λ，
),)
)(().()((,))(()(2
111221221121112212211b a c b b c c a c a c q b a c b b c c a c a c p -----=--+--
-=由题可知，，该点是稳定的。

即),0())(),((2
22
221c b a c t x t x t -→∞→时，，说明物种1最终会灭亡。

对于P 2，系数矩阵
⎥⎥⎥⎥⎦
⎤⎢⎢
⎢⎢⎣⎡
-------=⎥⎦⎤⎢⎣
⎡--------=1211221
11211221121222211222
1211111)(0)()(c c a c a c b a c b a c x c h x c b a c x c b x b x c b x c b a c A λ ),)().((,)()(2
211221112112211c c
a c a c a c q c c a c a c a c p -----=--
-+--= 由题可知，，该点是稳定的。

即)0,(
))(),((1
11
121c b a c t x t x t -→∞→时，，物种2最终要灭亡。

方程组⎩⎨⎧==0)(0)(2
1x f x f 为线性方程组，在平面上匹配两条直线21,l l 将第一象限分为三个
区域。

（1） 当
22
11c a c a >
时，随着时间的增加，物种1将会灭亡，物种2将达到稳定值2
22
2c b a c -。

（2）当
2211c a c a <
时，随着时间的增加，物种1最终能达到稳定值1
11
1c b a c -，物种2最终要灭亡。

4.蝴蝶效应与混沌解
解：
（1）编写Lorenz函数，Matlab程序如下：
Function dx=lorenz(t,x,b,a,c)
dx=[-b*x(1)+x(2)*x(3);
-a*x(2)+a*x(3);
-x(1)*x(2)+c*x(2)-x(3)];
调用ode45函数进行求解，利用plot函数进行绘图，如下所示：
（2） 适当调整参数，将βρα,,值各减1，初始值832110)0(x )0(),0(-=不变，x x ,得到
如下图示：
从图中我们可以看出，图形对参数和初始值的变化敏感很高，随着参数和初始值的轻微变化而改变很大。

解：
（1） 由题可得，微分方程为：⎪⎪⎩
⎪⎪⎨⎧==--=0)0('1)0(222
2x x dt dx h x k mg dt
x d m
其中，去m=1,g=9.8,用Matlab 求的数值解，并作出x(t)的图形： 分别为h=0和h=0.1时，
（2） 由题可知，受附加力后，微分方程为：⎪⎪⎩
⎪⎪⎨⎧==+--=0)0('1)0(sin 222
2x x wt B dt dx h x k mg dt
x d m ，
其中，B=1，w=1,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0，用Matlab 求取数值解并作图： 如下非别为h=0以及h=0.1时的图形，
（3） 如上图所示，在不考虑外力以及摩擦力时，简谐振动会一直进行下去；而当h=0.1
时，振幅逐渐减小并趋近于0，但是周期不变。

在施加一个外力后，振幅发生变化，并且去外力的频率有关系，这就是受迫运动。

解：
（1） 按三段时间利用Matlab 对数据分别拟合，确定增长率r 1,r 2,r 3,图形如下，增长率如图
中所显示：
（2）采用Logistic模型，对所有数据利用Matlab进行拟合，代码如下
t = 0 : 10 : 190;
x = [3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5];
f = inline('p(1)./(1+p(2)*exp(-p(3).*t))','p','t');
p = lsqcurvefit (f, [300, 50, 0.02], t, x);
得出人口的最大容量N m为360.3560（百万），增长率r为0.0234
图形如下所示：
改变参数，t=200,可以得到1990年人口总数为预测为241.7704（百万）。