数学建模实验二

1、数学建模实验二(微分方程实验)基本实验1. 微分方程稳定性分析 绘出下列自治系统相应的轨线，并标出随t增加的运动方向，确定平衡点，并按稳定的、渐近稳定的、或不稳定的进行分类：解答：（1）由平衡点的定义可得，f(x)=x=0，f(y)=y=0，因此平衡点为(0,0)，系统的线性近似方程的系数矩阵为，解得其特征值.,.对照稳定性的情况表，可知平衡点(0, 0)是不稳定的。自治系统相应轨线为：（2） 由平衡点的定义可得，f(x)=-x=0，f(y)=2y=0，因此平衡点为(0,0)，系统的线性近似方程的系数矩阵为，解得其特征值.,.对照稳定性的情况表，可知平衡点(0, 0)是不稳定的。自治系统相应轨

2、线为：（3） 由平衡点的定义可得，f(x)=y=0，f(y)=-2x=0，因此平衡点为(0,0)，系统的线性近似方程的系数矩阵为，解得其特征值.,.对照稳定性的情况表，可知平衡点(0, 0)是不稳定的。自治系统相应轨线为：（4） 由平衡点的定义可得，f(x)=-x=0，f(y)=-2y=0，因此平衡点为(0,0)，系统的线性近似方程的系数矩阵为，解得其特征值.,.对照稳定性的情况表，可知平衡点(0, 0)是稳定的。自治系统相应轨线为：2. 营养平衡问题 营养以每单位时间R个分子的常速流入一个细胞，并且以其内营养浓度成比例的速度离开，比例常数为K，设N为t时刻的浓度，则上述营养变化速度的数学描述

3、为：即N的变化速度等于营养进入细胞的速度减去它们离开的速度，营养的浓度会达到平衡吗？如果能，平衡解是什么？它是稳定的吗？试用这个方程解的图示解释之。解答：由题意可得N满足的微分方程为：，令，可求得方程的平衡点，当时，；当时，.不难计算出，由题意知,故平衡点是稳定的。由以上分析可做图示分析如下：3. 种群增长模型一个片子上的一群病菌趋向于繁殖成一个圆菌落.设病菌的数目为N，单位成员的增长率为，则由Malthus生长律有，但是，处于周界表面的那些病菌由于寒冷而受到损伤，它们死亡的数量与成比例，其比例系数为，求N满足的微分方程.不用求解，图示其解族.方程是否有平衡解，如果有，是否为稳定的?解答：由题

4、意可得N满足的微分方程为：，令，可求得方程的两个平衡点，分析可得，当时，单调递减，当时，单调递增，当时，单调递减，所以是不稳定的，是稳定的。对求二次微分可得：，令 得，该点即为曲线的拐点。方程的解族如下图所示：由图形也可以看出，是不稳定的，是稳定的。4. 单种群开发模型考虑单种群开发方程用数学表达式证明：在稳定状态下，最优捕捞率为解答：由本问题的目标出发，渔场中鱼量达到稳定的平衡状态时的情形，不必知道每一时刻的鱼量变化情况，故不需要解出方程，只需要讨论方程的平衡点并分析其稳定性。令解得的两个平衡点为：，易得由定理知：若，则是稳定的；若，则不是稳定的。应用上述近似判别法，所以有当E<r时，

5、是稳定平衡点，不是稳定平衡点； 当E>r时，不是稳定平衡点，是稳定平衡点； 所以，当捕捞适度(即：E<r)时,可使渔场产量稳定在，从而获得持续产量，而当捕捞过度（即：E>r）时，渔场产量将减至，从而是不可持续的。令求导可得：所以，最优捕捞率为。5. Compertz 模型 设渔场鱼量自然增长服从Gompertz模型：其中r为固有增长率，N为最大种群数量。若单位时间捕捞量为讨论渔场鱼量的平衡点及其稳定性，求最大持续产量及获得最大产量的捕捞强度和渔场鱼量水平。解答：令：由：，解得该方程的平衡点为，可得平衡点是稳定的，而平衡点不稳定.最大持续产量的数学模型为：由前面的结果可得，令可

6、得最大产量的捕捞强度，从而得到最大持续产量，此时，渔场鱼量水平。6. 有限资源竞争模型 微分方程是两个物种为了共同的有限资源而竞争的模型，假设c1>a1，c2>a2。试用微分方程稳定性理论分析：(1) 如果，则 (2)如果，则 (3)用图形分析方法来说明上述两种情况。解答：（1）令得方程的平衡点为P0（0,0），P1（,0），P2（0, ）.对平衡点P0（0,0），系数矩阵又c1>a1，c2>a2则p=-(c1-a1)+(c2-a2) <0，所以该平衡点不稳定。对平衡点P1（,0），系数矩阵则p=，q=，若，且c1>a1，c2>a2，则q<0不稳

7、定。而对于P2（0,），有p>0，且q>0稳定，此时说明物种1最终要灭亡。(2)如果的情况下，方程在P1（,0）稳定，其他点不稳定，此时说明物种2最终会灭亡。(3)对于线性方程组，其中，直线将第一象限分成三个区域。 当时，P2点稳定，通过分析的单调性可得下图： 此时说明物种1最终要灭亡。 当时，P1点稳定，通过分析的单调性可得下图： 此时说明物种2最终要灭亡。7. 蝴蝶效应与混沌解考虑Lorenz模型其中，且初值为,，为一个小常数，假设，且。（1）用函数ode45求解，并画出x2x1,x2x3,x3x1的平面图；（2）适当地调整参数，值，和初始值x1(0)，x2(0)=0，x3(0

8、)，重复一的工作，看有什么现象发生。解答： （1）创建plot_yuanzhe.m文件，在plot_yuanzhe.m文件中编写下面的语句：在plot_yuanzhe.m文件中，编写下面的语句：f=(t,x)-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3);t_final=100;x0=0;0;1e-10;t,x=ode45(f,0,t_final,x0);subplot(2,2,1),plot(x(:,2),x(:,1)xlabel('x_2');ylabel('x_1');subplot

9、(2,2,2),plot(x(:,2),x(:,3)xlabel('x_2');ylabel('x_3');subplot(2,2,3),plot(x(:,3),x(:,1)xlabel('x_3');ylabel('x_1');运行程序，可得下面的图形：（2）修改此题的参数，令，且初值不变，保持不变，运行上面的程序，可得下面的图形：修改此题的参数，令，且初值不变，保持不变，运行上面的程序，可得下面的图形：可以发现，修改参数和初始值，图形会发生很大变化。8. 药物动力学模型解答：题目中没有，这里我拟将作为肠道中的初始药物浓度。（1

10、）根据课件中的“药物分部中的房室模型”，该题为二室模型，其满足的微分方程组如下：（2）求解微分方程：创建D8_yuanzhe.m文件，编写下面的语句：x,y=dsolve('Dx=-a*x','Dy=a*x-b*y','x(0)=c','y(0)=0')可以得到下面的结果：x = (a2*c)/(exp(a*t)*(a - b) - (a*b*c)/(exp(a*t)*(a - b)/a y = (a*c)/(exp(b*t)*(a - b) - (a*c)/(exp(a*t)*(a - b)最小二乘曲线拟合的方法求解参数：创建D

11、8_2_yuanzhe.m文件，编写下面的语句：t=1 2 3 4 6 8 10 12 16;y=0.7 1.2 1.4 1.4 1.1 0.8 0.6 0.5 0.3;f=(a,t)(-a(1)*-a(3)/exp(-a(2)*t).*(a(1)-a(2)-(a(1)*a(3)/exp(a(1)*t)*(a(1)-a(2);a=1;2;3;a,b,c=lsqcurvefit(f,a,t,y)可以得到下面的结果：a = 0.1830 0.4347 5.9981b = 0.0356c =0.1080 -0.0039 -0.0652 -0.0691 0.0333 0.0741 0.0432 -0.

12、0384 -0.0707所以：参数的估计值为：.加分实验（化学动力学模型）解答：（1）根据题意建立模型：（2）创建D9_yuanzhe.m文件，编写下面的语句：x,y,z=dsolve('Dx=-a*x-d*x','Dy=a*x-b*y-c*y','Dz=b*y+d*x','x(0)=1','y(0)=0','z(0)=0')可以得到下面的结果：x = -(a*b*(a*b - a*d + b*d + c*d - d2)/(a + d)*(a - b - c + d) - (d2*(a*b - a*

13、d + b*d + c*d - d2)/(a + d)*(a - b - c + d) - (a2*(a*b - a*d + b*d + c*d - d2)/(a + d)*(a - b - c + d) + (a*c*(a*b - a*d + b*d + c*d - d2)/(a + d)*(a - b - c + d) - (2*a*d*(a*b - a*d + b*d + c*d - d2)/(a + d)*(a - b - c + d) + (b*d*(a*b - a*d + b*d + c*d - d2)/(a + d)*(a - b - c + d) + (c*d*(a*b - a

14、*d + b*d + c*d - d2)/(a + d)*(a - b - c + d)/(a*b*exp(a*t + d*t) - d2*exp(a*t + d*t) - a*d*exp(a*t + d*t) + b*d*exp(a*t + d*t) + c*d*exp(a*t + d*t) y = -(a*b2*d2*exp(a*t + d*t)/(b + c)*(a - b - c + d) - (a2*b2*c*exp(a*t + d*t)/(b + c)*(a - b - c + d) - (a2*b3*exp(a*t + d*t)/(b + c)*(a - b - c + d) +

15、 (a2*b2*d*exp(a*t + d*t)/(b + c)*(a - b - c + d) + (a2*b*exp(b*t + c*t)*(a*b - a*d + b*d + c*d - d2)/(a + d)*(a - b - c + d) - (a*b3*d*exp(a*t + d*t)/(b + c)*(a - b - c + d) + (a*b*c*d2*exp(a*t + d*t)/(b + c)*(a - b - c + d) - (a*b*c2*d*exp(a*t + d*t)/(b + c)*(a - b - c + d) - (2*a*b2*c*d*exp(a*t + d*t)/(b + c)*(a - b - c + d) + (a2*b*c*d*exp(a*t + d*t)/(b + c)*(a - b - c + d) + (a*b*d*exp(b*t + c*t)*(a*b - a*d + b*d + c*d - d2)/(a + d)*(a - b - c + d)/(a*b2*exp(a*t + d*t)*exp(b*t + c*t) - b*d2*exp(a*t + d*t)*exp(b*t + c*t) + b2*d*exp(a*t + d*t