2 微分方程实验1、微分方程稳定性分析绘出下列自治系统相应的轨线，并标出随t 增加的运动方向，确定平衡点，并按稳定的、渐近稳定的、或不稳定的进行分类：,,,+1,(1)(2)(3)(4);2;2;2.dx dx dx dxx x y x dt dt dt dtdy dy dy dy y y x y dt dt dt dt⎧⎧⎧⎧==-==-⎪⎪⎪⎪⎪⎪⎪⎪⎨⎨⎨⎨⎪⎪⎪⎪===-=-⎪⎪⎪⎪⎩⎩⎩⎩解：（1）根据定义，代数方程组的实根即为系统的平衡点，即P(0, 0)，利用直接法判断其稳定性。

在点P(0,0)处，系统的线性近似方程的系数矩阵为1001A ⎡⎤=⎢⎥⎣⎦，解得其特征值λ1=1，λ2=1； p=-(λ1+λ2)=-2<0，q=λ1λ2=1>0；对照稳定性的情况表，可知平衡点(0, 0)是不稳定的。

图形如下：（2）根据定义，代数方程组的实根即为系统的平衡点，即P(0, 0)， 利用直接法判断其稳定性。

解得其特征值λ1=-1，λ2=2；p=-(λ1+λ2)=-1<0，q=λ1λ2=-2<0；易知平衡点(0, 0)是不稳定的。

（3）根据定义，代数方程组的实根即为系统的平衡点，即P(0, 0)，利用直接法判断其稳定性。

解得其特征值λ1=0 + 1.4142i，λ2=0 - 1.4142i；p=-(λ1+λ2)=0，q=λ1λ2=1.4142；易知平衡点(0, 0)是不稳定的。

（4）根据定义，代数方程组的实根即为系统的平衡点，即P(1, 0)，利用直接法判断其稳定性。

解得其特征值λ1=-1，λ2=-2；p=-(λ1+λ2)=3，q=λ1λ2=2；易知平衡点(1, 0)是稳定的。

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

设病菌的数目为N ，单位成员的增长率为r1，则由Malthus 生长律有1dNr N dt=⋅，但是，处于周界表面的那些病菌由于寒冷而受到损伤，它们死亡的数量与N 1/2成比例，其比例系数为r2，求N 满足的微分方程.不用求解，图示其解族.方程是否有平衡解，如果有，是否为稳定的?解：根据题意列出N 满足的微分方程：1212dNr N r N dt =- (1)得到其解为N 1=0, N 2=2221/r r ； 由(1)得：11222121221()()2d N r r N r N r N dt -=-⋅- (2)解得N=2221/4r r画出N （t ）的图形，即微分方程的解族，如下图所示：可以判断出其中N 1=0是不稳定的；N 2=2221/r r 是稳定的。

3、单种群开发模型考虑单种群开发方程：1-x-Ex dx x r dt N=（） 在不求解的情况下，绘出其解族曲线。

(2)用数学表达式证明：在稳定状态下，最优捕捞率为E*= 2r解：由本问题的目标出发，渔场中鱼量达到稳定的平衡状态时的情形，不必知道每一时刻的鱼量变化情况，故不需要解出方程，只需要讨论方程的平衡点并分析其稳定性。

平衡点：满足F(x)=1-x-Ex dx x r dt N=（）= 0 (1) 的点称为方程的平衡点。

解得的两个平衡点为：0(1)Ex N r=-，10x =容易算出两个解E-r 和r-E称平衡点是稳定的是指：对方程(1)的任一个解()x x t =，恒有lim ()*t x t x →∞= (2)判断平衡点x *是否稳定，可根据一阶近似方程：()'(*)(*)dxF x F x x x dt=⋅- (3) 判断。

该方程的一般解为：(*)()*F x t x t C e x =⋅+于是有下述结论：若F'(x*)<0，则x *是稳定平衡点；若F'(x*)>0，则x *不是稳定平衡点。

应用上述近似判别法，所以有当E<r 时， 01F'(x )<0,F'(x )>0⇒ x 0是稳定平衡点，x 1不是；当E>r 时， 01F'(x )>0, F'(x )<0⇒ x 0不是稳定平衡点，x 1是；结果分析：当捕捞适度(即：E<r)时,可使渔场产量稳定在0(1)Ex N r=-， 从而获得持续产量Ex 0，而当捕捞过度（即：E>r ）时，渔场产量将减至x 1=0，破坏性捕捞，从而是不可持续的。

进一步讨论：如何控制捕捞强度E 使得持续产量Ex 0最大：00()(1)E h x Ex N E r==-2(1)02m dh E r N E dx r =-=⇒= 结论：最优捕捞率为*2rE = 。

4、Gompertz 模型设渔场鱼量增长服从Gompertz 模型：xNrx dt dx ln =，其中r 为固有增长率，N 为最大种群数量。

若单位时间捕捞量为Ex h =．讨论渔场鱼量的平衡点及其稳定性，求最大持续产量m h 及获得最大产量的捕捞强度m E 和渔场鱼量水平*0x。

解：()t x 变化规律的数学模型为()Ex xNrx dt t dx -=ln 记 Ex xNrx x F -=ln )( （1） 令()0=x F ，得0ln=-Ex xNrx rENex -=0，01=x ．则有平衡点为1,0x x . 又()E r xNr x F --=ln'，()()∞=<-=1'0',0x F r x F ． 推出平衡点o x 是稳定的，而平衡点1x 不稳定.（2）最大持续产量的数学模型为：⎪⎩⎪⎨⎧≠=-=.0,0ln ..max x Ex x N rx t s Ex h Ex()x f由前面的结果可得 rE ENeh -=r Er Ee r EN Ne dE dh ---=，令.0=dEdh 得到最大产量的捕捞强度r E m =， 从而得到最大持续产量e rN h m /=，此时渔场鱼量水平eNx =*0。

5、有限资源竞争模型：微分方程1111112222221122[(1)][(1)]dx x a c b x b x dt dx x a c b x b x dt⎧=-+--⎪⎪⎨⎪=-+--⎪⎩是两个物种为了共同的有限资源而竞争的模型，假设c1>a1，c2>a2。

试用微分方程稳定性理论分析：（1）如果1212a a c c >，则1()0();x t t →→∞（2）如果1212a a c c <则2()0();x t t →→∞（3）用图形分析方法来说明上述两种情况解：（1）令111111122222221122()[(1)]0()[(1)]0dx f x x a c b x b x dt dx f x x a c b x b x dt ⎧==-+--=⎪⎪⎨⎪==-+--=⎪⎩得方程的平衡点为P 0（0,0），P 1（1111c a c b -,0），P 2（0, 2222c ac b -）. 对平衡点P 0（0,0），系数矩阵112200c a A c a -⎡⎤=⎢⎥-⎣⎦ 又c 1>a 1，c 2>a 2则p=-[(c 1-a 1)+(c 2-a 2)] <0，所以该平衡点不稳定。

以此类推：对平衡点P 1（1111c ac b -,0）：系数矩阵211111112111()()0b c a a c b A c a c a c c -⎡⎤--⎢⎥⎢⎥=-⎢⎥-+-⎢⎥⎣⎦则p=2112111a c c a c a c --+，q= 11211221)()[())]a c c a c a c -----（c ，若1212a a c c >，且假设c1>a1，c2>a2，则q<0不稳定而对于P2（0, 2222c a c b -），有p>0，且q>0稳定，此时1()0();x t t →→∞，说明物种1最终要灭亡。

（2） 而如果1212a a c c <的情况下则方程在P1（1111c a c b -,0）稳定，其他点不稳定，此时2()0();x t t →→∞说明物种2最终会灭亡。

6、考虑Lorenz 模型'1123'223'31223()()()()()()()()()()()()x t x t x t x t x t x t x t x t x t x t x t x t βσσρ⎧=-+⎪=-+⎨⎪=-+-⎩ 其中σ=10，ρ=28，β=8/3，且初值为,x 1（0）=x 2（0）=0，x 3（0）=ε，ε为一个小常数，假设ε=10-10，且0≤t≤100。

（1）用函数ode45求解，并画出x2~x1,x2~x3,x3~x1的平面图； （2）适当地调整参数σ，ρ，β值，和初始值x 1（0），x 2（0）=0，x 3（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)) figure(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) 即x 1 关于t 的变化关系图0510********-20-15-10-505101520Figure(2)是y(2) 即x 2关于t 的变化关系图Figure(3)是y(3) 即x 3关于t 的变化关系图51015202530-25-20-15-10-5051015202505101520253051015202530354045Figure(4)为)x\1\x2 x3的空间关系图4．验证“蝴蝶效应”洛伦兹方程的解对初始值十分敏感，现对x2的初始值稍加修改，将2改为2.01和1.99，让后求解x3的数值解。

用edit命令建立“lzsensi.m”的M文件，程序如下：clfhold[t,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]时的x 3-t 图形 绿色线（b ）表示初值条件为[12,2,9]时的x 3-t 图形 红色线（r ）表示初值条件为[12,2.01,9]时的x 3-t 图形容易看出：随着时间的推移，三条曲线的吻合程度越来越差，差距越来越大，变化也越来越不明显，成为混沌状态。