1、2 微分方程实验1、微分方程稳定性分析绘出下列自治系统相应的轨线,并标出随t增加的运动方向,确定平衡点,并按稳定的、渐近稳定的、或不稳定的进行分类:解:(1)根据定义,代数方程组的实根即为系统的平衡点,即P(0, 0),利用直接法判断其稳定性。在点P(0,0)处,系统的线性近似方程的系数矩阵为,解得其特征值1=1,2=1;p=-(1+2)=-20;对照稳定性的情况表,可知平衡点(0, 0)是不稳定的。图形如下:(2)根据定义,代数方程组的实根即为系统的平衡点,即P(0, 0),利用直接法判断其稳定性。解得其特征值1=-1,2=2;p=-(1+2)=-10,q=12=-20;易知平衡点(0, 0
2、)是不稳定的。(3)根据定义,代数方程组的实根即为系统的平衡点,即P(0, 0),利用直接法判断其稳定性。解得其特征值1=0 + 1.4142i,2=0 - 1.4142i;p=-(1+2)=0,q=12=1.4142;易知平衡点(0, 0)是不稳定的。(4)根据定义,代数方程组的实根即为系统的平衡点,即P(1, 0),利用直接法判断其稳定性。解得其特征值1=-1,2=-2;p=-(1+2)=3,q=12=2;易知平衡点(1, 0)是稳定的。2、种群增长模型一个片子上的一群病菌趋向于繁殖成一个圆菌落。设病菌的数目为N,单位成员的增长率为r1,则由Malthus生长律有,但是,处于周界表面的那些
3、病菌由于寒冷而受到损伤,它们死亡的数量与N1/2成比例,其比例系数为r2,求N满足的微分方程.不用求解,图示其解族.方程是否有平衡解,如果有,是否为稳定的?解:根据题意列出N满足的微分方程: (1)得到其解为N1=0, N2=;由(1)得: (2)解得N=画出N(t)的图形,即微分方程的解族,如下图所示:可以判断出其中N1=0是不稳定的;N2=是稳定的。3、单种群开发模型考虑单种群开发方程:在不求解的情况下,绘出其解族曲线。(2)用数学表达式证明:在稳定状态下,最优捕捞率为E*= 解:由本问题的目标出发,渔场中鱼量达到稳定的平衡状态时的情形,不必知道每一时刻的鱼量变化情况,故不需要解出方程,只
4、需要讨论方程的平衡点并分析其稳定性。平衡点:满足F(x)= 0 (1)的点称为方程的平衡点。解得的两个平衡点为:,容易算出两个解E-r和r-E称平衡点是稳定的是指:对方程(1)的任一个解,恒有 (2)判断平衡点x*是否稳定,可根据一阶近似方程: (3)判断。该方程的一般解为: 于是有下述结论:若,则x*是稳定平衡点; 若,则x*不是稳定平衡点。应用上述近似判别法,所以有当Er时, x0不是稳定平衡点,x1是; 结果分析:当捕捞适度(即:Er)时,渔场产量将减至x1=0,破坏性捕捞,从而是不可持续的。进一步讨论:如何控制捕捞强度E使得持续产量Ex0最大: 结论:最优捕捞率为 。4、Gompert
5、z模型设渔场鱼量增长服从Gompertz模型:,其中r为固有增长率,N为最大种群数量。若单位时间捕捞量为讨论渔场鱼量的平衡点及其稳定性,求最大持续产量及获得最大产量的捕捞强度和渔场鱼量水平。解:变化规律的数学模型为 记 (1) 令,得 ,则有平衡点为 . 又,推出平衡点是稳定的,而平衡点不稳定.0 (2)最大持续产量的数学模型为:由前面的结果可得 ,令得到最大产量的捕捞强度,从而得到最大持续产量,此时渔场鱼量水平。5、有限资源竞争模型:微分方程是两个物种为了共同的有限资源而竞争的模型,假设c1a1,c2a2。试用微分方程稳定性理论分析:(1)如果,则(2)如果则(3)用图形分析方法来说明上述两
6、种情况解:(1)令得方程的平衡点为P0(0,0),P1(,0),P2(0, ).对平衡点P0(0,0),系数矩阵又c1a1,c2a2则p=-(c1-a1)+(c2-a2) a1,c2a2,则q0,且q0稳定,此时,说明物种1最终要灭亡。(2) 而如果的情况下则方程在P1(,0)稳定,其他点不稳定,此时说明物种2最终会灭亡。6、考虑Lorenz模型其中=10,=28,=8/3,且初值为,x1(0)=x2(0)=0,x3(0)=,为一个小常数,假设=10-10,且0t100。(1)用函数ode45求解,并画出x2x1,x2x3,x3x1的平面图;(2)适当地调整参数,值,和初始值x1(0),x2(
7、0)=0,x3(0),重复一的工作,看有什么现象发生。解:1 .建立自定义函数,在edit中建立“Lorenz.m”的M文件.程序如下:function dy = Lorenz(,y)dy=zeros(3,1);dy(1)=10*(-y(1)+y(2);dy(2)=28*y(1)-y(2)-y(1)*y(3);dy(3)=y(1)*y(2)-8*y(3)/3;end2.在edit中建立“Lzdis.m”的M文件,用来求解和绘图。程序如下:t,y=ode45(Lorenz,0,30,12,2,9);figure(1)plot(t,y(:,1)figure(2)plot(t,y(:,2)figur
8、e(3)plot(t,y(:,3)figure(4)plot3(y(:,1),y(:,2),y(:,3)plot3(y(:,1),y(:,2),y(:,3)3.运行得到如下的结果:Figure(1)是y(1) 即x1 关于t 的变化关系图Figure(2)是y(2) 即x2关于t 的变化关系图Figure(3)是y(3) 即x3关于t 的变化关系图Figure(4)为)x1x2 x3的空间关系图4验证“蝴蝶效应”洛伦兹方程的解对初始值十分敏感,现对x2的初始值稍加修改,将2改为2.01和1.99,让后求解x3的数值解。用edit命令建立“lzsensi.m”的M文件,程序如下:clfholdt
9、,u=ode45(Lorenz,0 15,12,2,9);plot(t,u(:,3),Color,r);t,v=ode45(Lorenz,0 15,12,2.01,9);plot(t,v(:,3),Color,b);t,w=ode45(Lorenz,0 15,12,1.99,9);plot(t,w(:,3),Color,k);运行得到不同初始条件下的x3关于t的图形:黑色线(k)表示初值条件为12,1.99,9时的x3-t 图形绿色线(b)表示初值条件为12,2,9时的x3-t 图形红色线(r)表示初值条件为12,2.01,9时的x3-t图形容易看出:随着时间的推移,三条曲线的吻合程度越来越差
10、,差距越来越大,变化也越来越不明显,成为混沌状态。2.3 加分实验(餐厅废物的堆肥优化问题)一家环保餐厅用微生物将剩余的食物变成肥料。餐厅每天将剩余的食物制成桨状物并与蔬菜下脚及少量纸片混合成原料,加入真菌菌种后放入容器内。真菌消化这此混合原料,变成肥料,由于原料充足,肥料需求旺盛,餐厅希望增加肥料产量。由于无力购置新设备,餐厅希望用增加真菌活力的办法来加速肥料生产.试通过分析以前肥料生产的记录(如表2.1所示),建立反映肥料生成机理的数学模型,提出改善肥料生产的建议。解:根据题意:将食物浆与蔬菜下脚及少量纸片混合成原料,加入真菌菌种,在容器内发酵转化成肥料。为了增加肥料产量,在不购买新设备的
11、条件下,依靠增加真菌活力的方法加速肥料的生产。实验记录给出了食物浆、蔬菜下脚、碎纸的量,并给出了投料日期和产出日期,这样我们可以知道肥料生成的时间长短。并且通过分析温度、湿度及投料比,确定最佳方案生产肥料。于是我们的问题可以描述为:1、在什么温度下生成肥料的速率最快;2、在什么湿度下生成肥料的速率最快;3、在什么样的投料比下生成肥料的速率最快。 为了解决上面提出问题,需要知道肥料生成的的天数,同时计算出对应天数下食物浆 、蔬菜下脚、碎纸之间的比例。除此之外还要建立温度和湿度的图像,通过比较来确立最合适的生成机制。 在解决这个问题的过程中主要运用控制变量法。通过查找资料,将以北方的温度和湿度为模
12、版,建立温度和湿度的图像。首先,进行模型假设: 1 将容器看作封闭的,不考虑质量的损耗。 2 以北方的温度和湿度为标准。 3 真菌的数量相同,初始活力相同。 4容器内生化反应过程中的温度不受人为因素控制,但受外界环境的影响。餐厅没有温度控制方面的投资。 5 反映开始前,真菌和发酵物分别储藏,不发生反应。 6 容器内的真菌分布均匀,且处于发酵的最佳状态。 7 在一定时间内,温度和湿度取平均值。建立数学模型:首先确立食物浆和蔬菜下脚的比例,并以比例为横坐标,肥料生成时间为纵坐标,建立坐标系,做出图像。编号(食物浆/蔬菜下脚)比例碎纸肥料生成天数12.7741902821.4177202733.38
13、09502742.4756102652.821430336 1.9811303678.0666703583.4375004791.86364949100.95000649111.50980749121.36842649食物浆与蔬菜下脚比例和肥料生成天数关系由图可以看出七八月份的产出明显要短,生成速率明显要高,增加碎纸的容器反而分解速率更低。从图中可以看出编号为4的当食物浆与蔬菜下脚比例为2.5左右且无碎纸,时间为7月27日到8月22日时,产出时间最短,生成速率最高,说明此时真菌的活性最大。比较编号为58组和912组可以看出,添加碎纸的组产出时间明显增长,说明碎纸片对真菌的活力起减缓作用,因此,
14、在食物浆与蔬菜下脚比例为2.5左右且无碎纸,投料时间为七月下旬时,真菌活性最大,所需时间最短,速度最快。然后变换图像,使其横坐标按顺序排列:可以看出,在不考虑温度和湿度的情况下,当比例为2.47561时,产出速率最快。通过查找资料,找到北方的平均温度和相对湿度,并作出了图像: 月份平均温度 (摄氏度)相对湿度 百分率%1-0.463206533.47048.473513.47761887721.594823.589921.373101664119.164122.763做出图像: 与我们通常的理解相近,当处于北方夏季时,温度要高(相对湿度也高),用于分解食物的酶活性也高,即此时菌种的活力相对较高
15、。 事实上,如果做出平均气温和湿度的曲线,于上面的曲线比较,也就找出了具有一般性的最佳生成机制。综上分析,要想使增加肥料的产量,必须使菌种在合理的温度和湿度条件下,合理的搭配投料比。由我们的模型知道,当温度在21.523.5 ,相对湿度在73%94%时,投料比为0.40394时肥料的生成速率最快。一阶常微分方程初值问题 数值解法是近似计算中很重要的部分。 常微分方程初值问题的数值解法是求方程的解在点列上的近似值,这里是到的步长,一般略去下标记为。 常微分方程初值问题的数值解法一般分为两大类: (1)单步法:这类方法在计算时,只用到、和,即前一步的值。因此,在有了初值以后就可以逐步往下计算。典型
16、方法如龙格库塔方法。 (2)多步法:这类方法在计算时,除用到、和以外,还要用,即前面步的值。典型方法如Adams方法。 经典的方法是一个四阶的方法,它的计算公式是: 方法的优点是:单步法、精度高,计算过程便于改变步长,缺点是计算量较大,每前进一步需要计算四次函数值。数学模型具有弹性,如果有更多更全面的数据,可以找出参数,建立具有一般性的模型公式。模型的缺点是没有充分利用生物实验的数据来确立模型公式,再有模型对温度和湿度的依赖性很大,也无法计算容器内部热量的产生和扩散。利用四阶龙格-库塔方法求解微分方程的初值问题问题1(1) TestRK4(ode1, 1, 0 -1, 5, inline(-x
17、-1)TestRK4(ode1, 1, 0 -1, 10, inline(-x-1)TestRK4(ode1, 1, 0 -1, 20, inline(-x-1)(2) TestRK4(ode2, 1, 0 1, 5, inline(1./(x+1)TestRK4(ode2, 1, 0 1, 10, inline(1./(x+1)TestRK4(ode2, 1, 0 1, 20, inline(1./(x+1)问题2(1) TestRK4(ode3, 3, 1 0, 5, inline(x.2.*(exp(x)-x)TestRK4(ode3, 3, 1 0, 10, inline(x.2.*(
18、exp(x)-x)TestRK4(ode3, 3, 1 0, 20, inline(x.2.*(exp(x)-x)(2) TestRK4(ode4, 3, 1 -2, 5, inline(2*x./(1-2*x)TestRK4(ode4, 3, 1 -2, 10, inline(2*x./(1-2*x)TestRK4(ode4, 3, 1 -2, 20, inline(2*x./(1-2*x)问题3(1) TestRK4(ode5, 1, 0 1/3, 5, inline(x.2+1/3*exp(-20*x)TestRK4(ode5, 1, 0 1/3, 10, inline(x.2+1/3*
19、exp(-20*x)TestRK4(ode5, 1, 0 1/3, 20, inline(x.2+1/3*exp(-20*x)(2) TestRK4(ode6, 1, 0 1, 5, inline(exp(-20*x)+sin(x)TestRK4(ode6, 1, 0 1, 10, inline(exp(-20*x)+sin(x)TestRK4(ode6, 1, 0 1, 20, inline(exp(-20*x)+sin(x)(3) TestRK4(ode7, 1, 0 0, 5, inline(exp(x).*sin(x)TestRK4(ode7, 1, 0 0, 10, inline(exp(x).*sin(x)TestRK4(ode7, 1, 0 0, 20, inline(exp(x).*sin(x) 考虑到餐厅的实际情况,给出以下建议:1. 增加每批处理的混合物的质量,确定合适的比例。2. 从生物反应方面考虑,增加氧气,充分搅拌,保证容器内混合物充分反应。3. 调节温度,是真菌处于最佳状态。4. 控制碎纸的量。5. 采用循环利用也是比较好的方法。