数学建模微分方程模型 数学建模微分方程实验

1、数学建模微分方程模型 数学建模微分方程实验 2 微分方程实验 1、微分方程稳定性分析绘出下列自治系统相应的轨线，并标出随t增加的运动方向，确定平衡点，并按稳定的、渐近稳定的、或不稳定的进行分类：dxdxdxdx=x,=-x,=y,=-x+1,dtdtdtdt(1)(2)(3)(4) dydydydy=y;=2y;=-2x;=-2y.dtdtdtdt 解：（1）由f（x）=x=0，f（y）=y=0；可得平衡点为（0,0），1系数矩阵A=00，求得特征值11=1，2=1；p=-(1+2)=-2<0，q=12=1>0；对照稳定性的情况表，可知平衡点（0，0）是不稳定的。图形如下： （2）

2、如上题可求得平衡点为（0,0），特征值1=-1，2=2；p=-(1+2)=-1<0，q=12=-2<0；对照稳定性的情况表，可知平衡点（0，0）是不稳定的。 其图形如下： （3）如上题可求得平衡点为（0,0），特征值1=0 + 1.4142i，2=0 - 1.4142i； p=-(1+2)= 0，q=12=1.4142>0；对照稳定性的情况表，可知平衡点（0，0）是不稳定的。其图形如下： （4）如上题可求得平衡点为（1,0），特征值1=-1，2=-2；p=-(1+2)= 3>0，q=12=2>0；对照稳定性的情况表，可知平衡点（1，0）是稳定的。 其图形如下： 2

3、、种群增长模型一个片子上的一群病菌趋向于繁殖成一个圆菌落.设病菌的数目为N，单位成员的增长率为r1，则由Malthus生长律有dNdt=r1N，但是，处于周界表面的那些病菌由于寒冷而受到损伤，它们死亡的数量与N1/2成比例，其比例系数为r2，求N满足的微分方程.不用求解，图示其解族.方程是否有平衡解，如果有，是否为稳定的?解：由题意很容易列出N满足的微分方程：dNdt1=r1N-r2N2=f(N)令f(N)=0，可求得方程的两个平衡点N1=0,N2=r22/r12 进而求得dNdt22dNdt22=(r1-12r2N-121)(r1N-r2N2)令=0可求得N=r2/4r122则N=N1，N=

4、N2，N=r22/4r12可以把第一象限划为三部分，且从下到上三部分中分别有dNdt0,dNdt220,；dNdt0,dNdt220；dNdt0。 则可以画出N（t）的图形，即微分方程的解族，如下图所示： 由图形也可以看出，对于方程的两个平衡点，其中N1=0是不稳定的；N2=r2/r122是稳定的。 3、有限资源竞争模型1926年Volterra提出了两个物种为共同的、有限的食物来源而竞争的模型 dx1=b1-l1(h1x1+h2x2)x1dtdx2=b-l(hx+hx)x2211222dt假设b1l1b2l2，称bili为物种i对食物不足的敏感度，（1）证明当x1(t0)>0时，物种2

5、最终要灭亡; （2）用图形分析方法来说明物种2最终要灭亡. 解：（1）由上述方程组f（x1）=b1-l1(h1x1+h2x2)x1=0，f（x2）=b2-l2(h1x1+h2x2)x2=0，可得方程的平衡点为P0（0,0），P1（b2b1l1h1,0），P2（0,l2h2）.对平衡点P0（0,0），b1-h1l1x1-h2l1x2系数矩阵A=-h1l2x2b1=b2-h1l2x1-h1l2x20-h2l2x10 b2则p=-（b1+b2）<0，所以该平衡点不稳定。对平衡点P1（b1l1h1,0），-b1-h2l2x1=b2-h1l2x1-h1l2x20b1l2)系数矩阵A=b1-h1l1

6、x1-h2l1x2-h1l2x2h2b1h1 blb2-12l1-则p= b1-b2+b1b2b1l2l1，q= -b1(b2-l1，由题意l1l2，x1(t0)>0，可以得出p>0,q>0,因此该平衡点是稳定的。b1即t时，(x1(t),x2(t)(l1h1,0)，说明物种2最终要灭亡。对平衡点P2（0,b1b2b2l2h2），同理可以得到q<0，在该平衡点不稳定。因此，在l1l2，x1(t0)>0的条件下，物种2最终要灭亡。b1-l1(h1x1+h2x2)=0（2）对于线性方程组b2-l2(h1x1+h2x2)=0在平面上匹配两条直线l1和l2，由题意b1l1

7、b2l2，x1(t0)>0，可将第一象限分为三个区域。在最左边区域，x1,x2都大于0；在中间区域，x1,x2都小于0，在最右边区域，x1,x2分别是大于0和小于0.，由各区域中x1,x2的取值可得到如下图形： 由图也可以看出，随着时间的增加，物种1最终能达到稳定值终要灭亡。4、蝴蝶效应与混沌解 b1l1h1，物种2最考虑Lorenz模型x1(t)=-bx1(t)+x2(t)x3(t)x2(t)=-sx2(t)+sx3(t) x(t)=-x(t)x(t)+rx(t)-x(t)12233其中=10，=28，=8/3，且初值为,x1（0）=x2（0）=0，x3（0）=，为一个小常数，假设=1

8、0-10，且0t100。（1）用函数ode45求解，并画出x2x1,x2x3,x3x1的平面图；（2）适当地调整参数，值，和初始值x1（0），x2（0）=0，x3（0），重复一的工作，看有什么现象发生。解：（1）编写Lorenz函数，function xdot=lorenz1(t,x,b,a,c)xdot=-b*x(1)+x(2)*x(3);-a*x(2)+a*x(3);-x(1)*x(2)+c*x(2)-x(3);对各参数赋值并用ode45函数求解，可得数值解：Columns 1 through 90 0.1250 0.2500 0.3750 0.5000 0.5352 0.5705 0.6

9、057 0.64090 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.00000 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 -0.00000.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 Columns 10 through 180.6761 0.7114 0.7466 0.7818 0.8308 0.8797 0.9286 0.9776 1.01050.0000 0.0000 0.0000 0

10、.0000 0.0000 0.0000 0.0000 0.0000 0.00000.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.00000.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 19 through 271.0434 1.0763 1.1092 1.1421 1.1750 1.2079 1.24091.2797 1.31860.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.00

11、00 0.00000.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.00010.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.0002Columns 28 through 361.3575 1.3964 1.4246 1.4528 1.4810 1.5092 1.5374 1.5656 1.59380.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.00000.0002 0.0003 0.0004 0

12、.0005 0.0008 0.0011 0.0015 0.0021 0.00290.0004 0.0006 0.0008 0.0012 0.0016 0.0023 0.0045 0.0063.Columns 5590 through 559899.5751 99.5921 99.6090 99.6260 99.6462 99.6664 99.7069 99.733816.9457 16.5261 16.2010 15.9854 15.8978 16.0348 17.1933 18.7894-3.3551 -3.7119 -4.1098 -4.5568 -5.1636 -5.8601 -7.54

13、38 -8.8677-5.3476 -5.9274 -6.5941 -7.3519 -8.3766 -9.5370 -12.1944 -14.0725Columns 5599 through 560799.7608 99.7878 99.8148 99.8306 99.8465 99.8623 99.8940 99.909921.2186 24.4709 28.2629 30.5354 32.6533 34.4645 36.7149 37.0528-10.3044 -11.7459 -12.9932 -13.5361 -13.8800 -13.9854 -13.4302 -12.7919-15

14、.7415 -16.8309 -16.9274 -16.3648 -15.3302 -13.8707 -10.1013 -8.1005Columns 5608 through 561399.9257 99.9416 99.9562 99.9708 99.9854 100.0000 36.8957 36.3239 35.5248 34.5435 33.4481 32.2971 -11.9606 -10.9899 -10.0237 -9.0332 -8.0563 -7.1223 -6.2187 -4.5596 -3.2818 -2.2590 -1.4835 -0.9275x2x1,x2x3,x3x

15、1的平面图分别如下：0.0032 99.6867 16.4476 -6.6527 -10.8209 99.8782 35.8466 -13.8351 -12.0805 （2）令参数，值各减1，初始值x1（0），x2（0）不变，x3（0）=10-8 分别得到得到x2x1,x2x3,x3x1的平面图如下： 可以看出，参数，值各减1，初始值x1（0），x2（0）不变，x3（0）数值变为x3（0）=10-8，参数和初始值很小的改变，就会导致最后图形发生很大的变化。 5、用微分方程考察共振现象 设物体沿x轴运动(如图所示)其平衡位置取为原点0，物体的质量为1，在时间t物体的位置为x(t)其所受的恢复为(

16、如弹性力等)与物体所在位置的坐标成正比，即k2x,其中常数k称为恢复系数，运动过程所受的阻力(由于介质及摩擦等)设与速度成正比，即2hdxdt，h>0，称为阻尼系数。（1）根据Newton第二定律，建立相应的微分方程.不妨设初始位置为1，初始速度为0,取k=2, h=0(当h = 0称为简谐振动的方程)和h=0.1，用Matlab软件得到相应的数值解，并在t-x平面上画出x（t）的图形。（2）如果物体还受到附加外力的干扰，且外力是一个依据时间t的函数f(t)(设f (t)=B sinwt)，建立相应的微分方程(该方程称为强迫振动方程).在上述参数不变的情况下，取振幅B=1，分别取w=1,

17、 1.2, 1.4, 1.6, 1.8, 2.0, 2.2,2.4, 2.6, 2.8, 3.0，用Matlab软件得到相应的数值解，并在t-x平面上画出x(t)的图形。（3）分别对上述图形讲行分析，并解释为什么会出现这些现象。 解：（1）根据Newton第二定律，F=ma，可得微分方程：m由题意知边界值：x(0)=1,x(0)=0, 令y1=x，y2=dy1dt=y2dy可将二阶方程转化为：2=mg-k2y1-2hy2dty1(0)=1,y2(0)=0dxdtdxdt22=mg-kx-2h2dxdt ， 其中，m=1；g=9.8；k=2 当h=0时，由matlab解得数值解：Columns

18、1 through 9 0 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000Columns 10 through 18 0.0004 0.0006 0.0009 0.0011 0.0022 0.0032 0.0043 0.0054 0.01081.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0001 1.0001 1.0003Columns 19 through 27 0

19、.0162 0.0216 0.0271 0.0541 0.0812 0.1083 0.1353 0.2155 0.29561.0008 1.0014 1.0021 1.0085 1.0191 1.0339 1.05281.1326 1.2461.Columns 100 through 108 8.2321 8.3593 8.4615 8.5636 8.6658 8.7679 8.87668.9852 9.09393.5019 3.2161 2.9505 2.6641 2.3688 2.0769 1.78341.5214 1.3034Columns 109 through 117 9.2025

20、9.3142 9.4260 9.5377 9.6494 9.7371 9.82479.9124 10.00001.1393 1.0345 1.0001 1.0383 1.1467 1.2773 1.44381.6412 1.8634 h=0时，x（t）图形如下所示： 当h=0.1时，由matlab解得数值解：Columns 1 through 90 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000Columns 10 through 18

21、0.0004 0.0006 0.0009 0.0011 0.0022 0.0054 0.01081.0000 1.0000 1.0000 1.0000 1.0000 1.0001 1.0003Columns 19 through 270.0162 0.0216 0.0271 0.0541 0.0812 0.2156 0.29581.0008 1.0014 1.0021 1.0085 1.01901.1308 1.2417.Columns 109 through 1178.5631 8.6690 8.7749 8.8809 8.98689.3934 9.52892.5857 2.4563 2.3

22、294 2.2106 2.10491.8880 1.8962Columns 118 through 1219.6467 9.7645 9.8822 10.00001.9353 2.0016 2.0909 2.1979 h=0.1时，x（t）图形如下所示：1.0000 1.0000 0.0032 0.0043 1.0000 1.0001 0.1083 0.1353 1.0336 1.0523 9.1223 9.2579 1.9952 1.9213 （2）如果还受到外力干扰，则微分方程变为：2mdx2dxdt2=mg-kx-2hdt+Bsinwt将其化为一阶方程组：dy1dt=y2dy2=mg-k

23、2y1-2hy2+Bsinwt dty1(0)=1,y2(0)=0其中，m=1；g=9.8；k=2；B=1；h=0.1； 用Matlab解得数值解：w=1时Columns 1 through 90 0.0000 0.0000 0.0000 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000.Columns 118 through 1219.8207 9.8804 9.9402 10.00001.8973 1.9210 1.9503 1.9846w=1.2时Columns 1 through 90 0.0000 0.0000 0.0000

24、 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000.Columns 118 through 1219.8281 9.8854 9.9427 10.00001.6999 1.7568 1.8199 1.8886w=1.4时Columns 1 through 90 0.0000 0.0000 0.0000 0.0000 0.0001 1.0000 1.0000 0.0000 0.0001 1.0000 1.0000 0.0000 0.0001 0.0001 1.0000 0.0001 1.0000 0.00010.0002 0.00021

25、.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000.Columns 118 through 1259.6236 9.7309 9.8382 9.9454 9.9591 9.9727 9.9864 10.00002.2401 2.3171 2.4084 2.5103 2.5238 2.5373 2.5510 2.5647w=1.4时Columns 1 through 90 0.0000 0.0000 0.0000 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000.C

26、olumns 118 through 1219.7823 9.8549 9.9274 10.0000 2.1150 2.0676 2.0285 1.9979w=1.6时Columns 1 through 90 0.0000 0.0000 0.0000 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000.Columns 118 through 1219.6933 9.7955 9.8978 10.00000.8062 0.7503 0.7562 0.8231w=1.8时Columns 1 through 90 0.0000 0.0000

27、0.0000 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000.Columns 118 through 1219.7752 9.8501 9.9251 10.0000 0.8629 1.0850 1.3383 1.61740.0000 0.0001 1.0000 1.0000 0.0000 0.0001 1.0000 1.0000 0.0000 0.0001 1.0000 1.0000 0.0001 1.0000 0.0001 1.0000 0.0001 1.0000w=2.0时Columns 1 through 90 0.0000

28、0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.00021.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000.Columns 118 through 1219.7538 9.8359 9.9179 10.00002.3278 2.6023. w=3.0时Columns 1 through 90 0.00000.0002 0.00021.0000 1.00001.0000 1.0000.Columns 118 through 1219.7495 9.83302.2414 2.3303

29、2.8705 3.1243 0.0000 0.0000 0.0000 0.0001 0.0001 1.0000 1.0000 1.0000 1.0000 1.0000 9.9165 10.0000 2.4145 2.4915 如上图所示，w分别为1.0；1.2；.；3.0时x（t）的图形。 当h=0时，即无阻尼振动的情况下，得到t（x）图形如下所示： （3）由上述图形可以看出，不考虑外加作用力时，当没有摩擦等阻力的时候，简谐运动一直进行下去，当有摩擦等阻力的时候，其振幅会逐渐的减小，最后趋近于零，但是周期不会发生变化。在外加一个作用力之后，整个运动情况就发生了很大的变化。首先刚开始的振幅和周期

30、都变化很大，后来随着时间的推移，他们也逐渐开始趋于稳定。但是，稳定之后的振幅和周期都各不相同，与外力的频率有直接的关系，这就是受迫振动，不仅跟系统本身的性质有关，还和外界干扰的情况有很大关系。恢复系数k与w的值越接近时，随着时间增加，物体的振动振幅会一直增大，当w=k=2时，物体的振幅增大的速度最快，可以预见，在t趋于无穷大时，振幅也趋于无穷大，这就是共振现象。也就是当外力与系统处于共振时，会引起振幅无限增大的振动，这在机械和建筑中一般是必须严格避免的。 6、人口模型与曲线拟合问题表2.1列出的是美国1790一1980年的人口统计表.（1）试用Malthus人口模型，按三段时间(1790-18

31、50,1850-1910,1910-1970)分别确定其增长率r。并将数据和不同r值的三段曲线画在同一张图上.（2）利用Logistic模型，重新确定固有增长率r和最大容量Nm，作图，再利用该模型得到的结果预测1990年的美国人口数。解：（1）分段研究，我们先求出增长率，编写命令如下：t = 10 : 10 : 70;x = 3.9 5.3 7.2 9.6 12.9 17.1 23.2;p = polyfit(t, log(x), 1); p得到结果p = 0.0296即17901850年时间段的增长率是0.0296同理可以得到另外两段的增长率，1850-1910年 为0.0226191019

32、70年为0.0129每次画出图形，使用以下命令plot(t, x, bx);hold on在以上每一次求时顺便画出图形，得到最后的总图如下 (2)利用Logistics模型，将以上所有数据进行拟合，编写程序如下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,

33、 0.02, t, x); 得到结果人口最大容量Nm是360.3560（百万），增长率r是0.0234 如下图所示： 输入t_pre = 200;x_pre = p(1) ./ (1 + p(2) * exp(-p(3)*t_pre) );可得最后的结果，即1990年的人口总数x_pre =241.7704万人。如下图： 7、加分实验（餐厅废物的堆肥优化问题） 一家环保餐厅用微生物将剩余的食物变成肥料。餐厅每天将剩余的食物制成桨状物并与蔬菜下脚及少量纸片混合成原料，加入真菌菌种后放入容器内。真菌消化这此混合原料，变成肥料，由于原料充足，肥料需求旺盛，餐厅希望增加肥料产量。由于无力购置新设备，餐

34、厅希望用增加真菌活力的办法来加速肥料生产.试通过分析以前肥料生产的记录(如表2.2所示)，建立反映肥料生成机理的数学模型，提出改善肥料生产的建议。 解答：1.问题条件假设1）每次堆肥的质量相等2）所给的几次堆肥混合物的比例仅由当天的实际情况决定3）所有分离肥仓工作条件相同4）每磅蔬菜可提供的氧气量相等5）细菌消耗的溶解氧气完全由蔬菜叶提供6）每天提供的废物混合物中的化学成分大致相同7）废物混合物在喂给细菌前混合均匀并保持良好的通气环境。2.问题分析堆肥是利用微生物的分解作用将有机废物转化成无害稳定形式的生物化学过程，要提高堆肥产量方法之一就是通过增强细菌的生长繁殖能力以提高分解率。 细菌群体的

35、增长一般要经历延滞期，加速生长期、对数生长期、减数生长期、稳定期、加速死亡期和对数死亡期。另外，对数生长期培养基中所有养分都过剩，细菌可以充分繁殖，其倍增速率恒定，取决于底物浓度、温度、水活度、供氧量。对于当前该餐厅来说底物浓度由每天的剩余食物、蔬菜叶和碎纸屑决定。碎纸屑是吸收水分的调理剂，微菌呼吸所需要的溶解氧由蔬菜也提供，水活度可以通过测定相对湿度来决定，其关系式是：相对湿度：B=P/P0*100%水活度： aw=P/P0 其中P为该溶液蒸汽压，P0为纯水蒸汽压， 在这个堆肥系统中，可供微菌消耗的营养底物和溶解氧都是有限的，他们的消耗会对微菌生长率产生重要影响，一种微菌消耗营养底物的速率和

36、底物浓度之间的关系式为：ds/dt=KmSX/（Ks+S），式中ds/dt表示为底物的有效消耗率，x表示为微菌浓度，Km表示为最大有效系数，在高浓度营养底物中最大的物料消耗率（物质质量/微菌一天的质量）；Ks表示为半速系数（质量/体积），S表示为有限底物的浓度（质量/体积）微菌生长过程是一个生物化学过程，其到经验最佳比例，由于假设各次堆肥后肥料质量相同，因此堆肥时间较短就对应了较好的废料配比，将12组数据按其堆肥日期及完全堆肥时间分成三组，每组模型二：营养底物和氧气是细菌生长的两种底物，物耗公式为：dSi/dt=KmiXSi/(Ksi+Si) （1）i=1时代表营养底物中可降解得有机物。i=2时代表氧气。1）在高浓度底物中，物料的转化过程很迅速，进一步增加底物浓度就不再引起底物转化率的提高，可假设S1>>Ks1,S2>>Ks2，