微分动态模型.doc

第二部分 动力系统模型第四章 动力系统的解的定性分析4.1 连续动力系统的解的定性分析例4.1 设森林中生长着硬树木H(吨)0及软树木S(吨)0，假定树木量的年增长量（吨/英亩/年）为dH/dt=(r1-a1*H-b1*S)*H, dS/dt=(r2-a2*S-b2*H)*S, 其中r1，r20分别是硬木及软木的内禀增长率，r1-a1*H-b1*S 是硬木的实际增长率，0a1*H0体现了软木对硬木增长的影响. 对于软木的情况也有对应的说明.解: 首先来求出平衡点，即硬木和软木量与时间无关的解(H,S)，显然(H,S)=(0,0)是解(零解)，但这不是我们感兴趣的定常解，不予考虑，同样，还有两个定常解(H,S)=(0,r2/a2),(H,S)=(r1/a1,0), 由于只有一种树木存在，也不予考虑. 我们看是否有两种树种都存在的定常解(H=(r1*a2-r2*b1)/(a1*a2-b1*b2), S=(r2*a1-r1*b2)/(a1*a2-b1*b2); 假设种群内部的竞争比种群之间的竞争强，则a1b1, a2b2, 从而a1*a2-b1*b20; 所以两种树木共存的条件是r1*a2-r2*b10; r2*a1-r1*b20; 即a1/b2r1/r2b1/a2; 以下我们要讨论的是，一般情况下两种树种的量是否会趋于这个共存的平衡点. 用动力系统的语言来讲就是这个定常解是否是全局吸引的. 首先来证明这方程的解的正半轨最终将位于一个有界的闭域：0=H=2*r1/a1; 0=S=2*r1/a1时，dH/dt=2*r2/a2时, dS/dt= -r2*S. 因此终究会位于该闭域内. 因此任何轨道的W极限集存在并有界，二维动力系统的有界的W极限集是极限环或是平衡点与连接平衡点的非闭轨. 首先来看在这平衡点附近微分方程的线性部分，在平衡点方向场的Jacobi为-a1*H，-b1*H; -b2*S, -a2*S; 特征方程为 x2+(a1*H+a2*S)*x+(a1*a2-b1*b2)*H*S=0;根据Hurwitz定理，特征值都具有负的实部，所以平衡点是渐近稳定的，即平衡点至少是局部吸引的，但是否是全局吸引的还需进一步研究。这时，我们还要考察其他平衡点的吸引性. 在零解附近，方向场的Jacobi为r1,0; 0, r2, 故零解是反向吸引的，因此非零解不会趋于零解, 对于定常解(H,S)=(0,r2/a2)，Jacobi为r1-b1*S, 0; -b2*S, r2-2*a2*S=r1-b1*r2/a2,0; -b2*r2/a2, -r2; 特征方程为x2+(r2-r1+b1*r2/a2)*x-r1*r2+b1*r22/a2=0; 特征值为(r1*a2-r2*b1)/a2, -r2, 故平衡点(H,S)=(0,r2/a2)是鞍点；现在的问题是对于一般的共存解(H,S), 是否会趋于这个鞍点呢？在S轴上有两条轨道，0Sr2/a2, r2/a2=S=r2/a2), dH/dt=0, dS/dt=(r2-a2*S)*S=0, , 在r1-a1*H-b1*S=0, 上的一段， (0=H=(r1*a2-r2*b1)/(a1*a2-b1*b2), dH/dt=0,dS/dt=(r2-a2*S-b2*H)*S=(r2-2*(r1-a1*H)/b1-b2*H)*S=-(r1*a2-r2*b1)-(a1*a2-b1*b2)*H)/b1*S=0, 在r2-a2*S-b2*H=0上的一段(0=H=0 即在这三角形区域内的轨线都不会走出这三角形区域，而包围共存解(H,S)的闭轨必走出这三角形区域. 所以系统不存在闭轨. 从而一般的共存解(H,S)必定趋向共存的平衡点. 就是说，在自然状态下，两个树种最终会达到共存的平衡态.例5 单摆方程y = dsolve(D2y =-sin(y),y(0)=a,Dy(0)=0)得y = RootOf(Int(1/(2*cos(_f)-2*cos(a)(1/2),_f=a._Z)+t) RootOf(Int(-1/(2*cos(_f)-2*cos(a)(1/2),_f=a._Z)+t)这个结果不是令人满意的，需要作人工的处理，方法如下：积分一次得(Dy)2= 2(cos(y)-cos(a)=4(sin(a/2)2-sin(y/2)2),t=int(1/2/sqrt(sin(a/2)2 sin(y/2)2), 0, y),令sin(y/2)=sin(a/2)*sin(z), 化为t=int(1/sqrt(1- sin(a/2)2sin(z)2), 0, z)，是第一类不完全椭圆积分，记作t=EllipticF(sin(z), k) , k= sin(a/2),在MATLAB中的MAPLE符号运算中有这个函数，调用方法如maple(EllipticF(0.5,0.9), 注：要调用MAPLE函数的方法就是把MATLAB中的命令用单引号括起写在maple()中的括号中即可. EllipticF的反函数是椭圆函数JacobiAM，两者之间的关系为EllipticF(sin(z), sin(a/2) = InverseJacobiAM(z, sin(a/2)，当sin(z)=1时的第一类不完全椭圆积分就是第一类完全椭圆积分EllipticK(sin(a/2), 当a=0时的值等于pi/2，当a=pi时等于inf.这就说明了单摆的周期是振幅的函数，随着振幅的增大，周期增大到无穷大.最后单摆方程的解可以写成y=2*asin(sin(JacobiAM(t,sin(a/2)*sin(a/2)= 2*asin(JacobiSN(t,sin(a/2)*sin(a/2).这个结果要比前一结果更具体了， 并且可以用MATLAB中的MAPLE功能来求值。但如果不需要精确求值，可以使用MATLAB中的Jacobi 椭圆函数SN,CN,DN=ellipj(u,m), 其中u=int(1/sqrt(1- m*(sin (x)2), 0, phi)，SN=sin(phi), CN=cos(phi), DN= sqrt(1- m*(sin (phi)2), am(u)=phi, m= k2=sin(a/2)2. 可见单摆方程的解用MATLAB函数可表示为SN,CN,DN=ellipj(t, sin(a/2)2)，y=2*asin(sin(a/2)*SN).MATLAB还有完全椭圆积分K，E=ellipke(m), m=k2, K是第一类椭圆积分的值，E是第二类椭圆积分的值. 当m= sin(a/2)2，K，E=ellipke(m)，则4*K就是单摆的周期.4.2 动力系统解的数值计算(定量分析)在动力系统的定性分析中，我们不知道解与共存定常解接近的动态情况，比方说在何时接近到何程度。这时的定量分析就需要进行数值计算，要进行数值计算就需要给出参数的值.若在例4.1中取，r1=0.1; r2=0.25; a1=0.05; a2=0.1; b1=a1/2;b2=a2/2; 这些参数满足例4.1的条件，共存平衡解为(H,S)=(1, 2);画出方向场的箭图r1=1.e-1; r2=2.5e-1; a1=5.e-2; a2=1e-1; b1=a1/2;b2=a2/2;H=1.2*(r1*a2-r2*b1)/(a1*a2-b1*b2);S=1.2*(r2*a1-r1*b2)/(a1*a2-b1*b2);A=linspace(0,H,20); B=linspace(0,S,20);x,y=meshgrid(A, B);DY1=(r1-a1*x-b1*y).*x;DY2=(r2-a2*y-b2*x).*y;quiver(x,y, DY1,DY2),从箭图上可见共存平衡解(1,2)是全局吸引的.下面我们设这树林被人类砍伐后两种树的量为(0.5,1)(吨)，问要经过150年后树木恢复的情况. 由于方程无解析解，用MATLAB来数值求解. 首先写一个ODE函数function DYDT=ODEFUN41(T, Y)r1=0.1; r2=0.25; a1=0.5e-1; a2=1e-1; b1=a1/2;b2=a2/2;DYDT=(r1-a1*Y(1)-b1*Y(2)*Y(1); (r2-a2*Y(2)-b2*Y(1)*Y(2);设在T=0时，初始值为(0.5, 1), 求两种树的增长情况.数值求解并作图的程序为Tspan=0,150;Y0=0.5, 1;Sol=ode45(ODEFUN41, Tspan, Y0);T=linspace(0,150,151);S=deval(Sol,T);plot(T,S(1,:),T,S(2,:); %两种树的生长曲线figure;plot(S(1,:),S(2,:);%相轨线在本题中，积分的时间较长(150), 由此产生的误差会较大，我们可以改变时间单位，使得积分的长度变小. 对于本题，我们可以作自变量的变换T=40*T1，方程化为DYDT1=(4-2*Y(1)-Y(2)*Y(1); (1-4*Y(2)-2*Y(1)*Y(2);这样，时间区间就化为0，150/40=0，3.75另外，要注意，由于动力系统的非定常解不可能在有限时间内到达平衡点. 所以对于本模型，树木不可能在有限时间内恢复到平衡态.练习：找自变量和因变量的线性变换，使得例4.1的微分方程化为一个三参数的无量纲的微分方程:dx/dT=(1-x-B1*y)*x; dy/dT=(R2-R2*y-B1*x)*y 4.3 离散动力系统例4.3 在两艘太空船对接时用手工控制两太空船的相对速度，设第n个观测时刻t(n)测得速度为v(n),在宇航员的反应时间c(n)后，宇航员调整速度控制器, 再经过等待时间w(n)后再次测量速度，因此下一次测量的时刻是t(n+1)=t(n)+c(n)+w(n)，设第n次调解后加速度为a(n),则t(n+1)时刻的速度为v(n+1)=v(n)+a(n-1)*c(n)+a(n)*w(n),控制规则是a(n)=-k*v(n),k是正常数. 问题是，是否当n趋向无穷大时，v(n)趋于零？解：设c(n)=c,w(n)=w是常数，则该问题是一个2阶常系数齐次线性差分方程v(n+1)+(w*k-1)v(n)+c*k*v(n-1)=0初值条件是v(0)=v00, v(1)=v1. 特征方程为x2+(w*k-1)*x+c*k=0, 特征根为x1,x2=(1-w*ksqrt(w*k-1)2-4*c*k)/2; 当两特征值不相等时方程的通解为v(n)=C1*x1n+C2*x2n, 当特征值是一对共轭虚数 r*exp(i*beta), beta0时，v(n)=rn*(C1*cos(n*theta)+C2* sin(n*theta)；当特征值是相同实根x时，v(n)=xn*(C1+C2*n)；可见当且仅当两特征值的绝对值均小于1时，v(n)趋于零.6.1 离散非线性动力系统例6.1 一支R单位(师)的部队(R方)和一支B单位(师)的部队(B方)作战，设直接交火的伤亡率与对方人数成正比，由炮火产生的伤亡率与对方的部队数和己方的人数的乘积成正比。因此，可列出差分方程如下：R(n+1)=R(n)-(a1+b1*R(n)*B(n);B(n+1)=B(n)-(a2+b2*B(n)*R(n);R(0)=R0, B(0)=B0其中 a1,b1,a2,b20是常数，a2表示R方步兵的武器装备的强度，b2表示R方炮火的强度.假定R0B0，而B方想要加强武器装备来取胜，问B方武器要比R方强几倍才能取胜？即取a1，b1至少多大时可使在第n=N个战斗小时时，R(N)=0而B(N)0.如果可以分别地任意取a1，b1，则一个简单的方法可确定a1，b1的下确界，把方程改写为(R(n+1)/R0)=(R(n)/R(0)-(a1*B0/R0+b1*B0*(R(n)/R0)*(B(n)/B0);(B(n+1)/B0)=(B(n)/B0)-(a2*R0/B0+b2*R0*(B(n)/B0)*(R(n)/R0) ;则可见只要加强炮火的强度，取b1到下确界b1=R0/B0*b2, a1=(R0/B0)2*a2时, 差分方程的解为R(n)/R0=B(n)/B0，即R方和B方同时为零. 即B方加强炮火的强度至少为部队人数的比的倒数倍, 而步兵武器强度为部队人数的比的倒数的平方倍. 但现在的问题是：设 a1=k*a2, b1=k*b2, 问k至少多大时，B方可取胜？ 对于这个问题，没有解析解，只能用数值计算，所以要给出参数的值，设R0=5(师)，B0=2(师)，a2=0.05, b2=0.005; 按照前面的分析，对于a1至少要是a2的6.25倍，对于b1至少应是b2的2.5倍，所以现在k是比2.5大而比6.25小的数.分别取k=3，4，5，6，利用计算机来计算，看k为何值时，B方可取胜.解：首先确定N=168(小时)，（14个战斗日，每日战斗12小时），编MATLAB程序如下function R,B,N=zhandou(K)R=5; B=2; a2=0.05; b2=0.005;for N=1:168RT=R-K*(a2+b2*R)*B;BT=B-(a2+b2*B)*R;R=RT;B=BT;if R=0|Bt, 并仍用t记新的自变量，方程化为dx1/dt= x2 - f(x1); dx2/dt = - x1*L/C; 这是Lienard方程dx1/dt= x2 - f(x1); dx2/dt = - g(x1);的特殊形式, 此系统的惟一平衡点是(0,0), 要求其稳定性，方程的线性部分的系数矩阵为-f(0)/L, 1/L; -1/C,0; 特征方程为 x2+f(0)/L*x+1/(L*C)=0; f(0)=4, 故根据Hurwitz定理，由于特征方程的系数都是正的, 特征值都具有负的实部，所以平衡点是渐近稳定的. 因此此平衡点至少是局部吸引的, 取V函数V=L/C*x12+x22, 全导数 dV/dt = -2*x1*f(x1)/C0, G(x)=int(g, 0, x), G(Inf)= Inf, 2. 存在 c10, 使得c1x0; 0xc2时 f(x)0;3. 存在k2max(c2, -c1), 使得当xM时 F(x)=k1; x= -M时， F(x)0, 使得abs(x)0, 和任意xA， 存在yA及t=0，使得|y-x|delta(这个性质称为动力系统phi(t，x)对初值的敏感性). 2 phi 有拓扑传递性: 对于任意两个开集B，C，存在t0，使得phi(t,B) phi(t,C)(看似随机的行为)3. phi 的周期点在A中稠密则称A是动力系统phi(t，x)的混沌集，若A又是测度大于零的集合，则称动力系统phi(t，x)是混沌的动力系统.分形最早的例子是Weierstrass的处处不可微的连续函数及Cantor的三分集, 后来在1975年由美国IBM公司的数学家Benoid B. Mandelbrot 首次提出fractal这个新术语，在1982年出版“自然界中的分形几何”，于1985年荣获Barnard奖. 著名的分形的例子还有Siepinski垫子、Siepinski毯子、Siepinski海绵、Vicssek 图形、Julia集，Mandelbrot集、Koch曲线，分形没有精确的数学定义，以下是一个常见的分形定义: 一个集合A称为分形, 如1. A具有精细的结构,即有任意小比例的细节;2. A的整体和局部都不能用传统的几何语言来描述;3. A通常有某种自相似的形式,可能是近似的或是统计的.4. A的某种分形维数大于其拓扑维数;5. 很多情况下A以非常简单的方式定义或叠代产生.例6.5 设有一个连续的线性系统 dx/dt = - c*x，c0, x(0)=x0 0; 它的解为x(t)=x0*exp(- t); 再考虑对应的离散系统，设t(n)=n*h, h0是步长，x(t)=y(n), 若用向前差商代替导数，则得差分方程y(n+1)=(1-c*h)*y(n)，从而y(n)=x0*(1-c*h)n, x(t)=x0*(1-c*h)(t/h), 其中t=n*h， 可见，当h趋于0时，解趋于x0*exp(-c*t), 而当h增大时，离散系统与连续系统的解的差别就越来越大. 原因是导数用差商来替代一般都是近似的. 当0c*h1; 而当1c*h2时，y(n)振荡地趋于无穷大. 如果认为用中心差商代替导数会得到更好的近似效果，但这是不对的. 因为这时得到一个二阶差分方程y(n+1)-y(n-1)=-2*h*c*y(n)，特征方程有两个实根，其乘积为-1，通解为y(n)=C1*(sqrt(1+(h*c)2)-h*c)n+C2*(-(sqrt(1+(h*c)2)+h*c)n, 其中前一个解对应于一个小于1的正特征值，解是衰减的，后一个解对应于小于-1的特征值，解是振荡趋于无穷的，由于连续模型的解是衰减的，故C2=0，解应为y(n)=x0*(sqrt(1+(h*c)2)-h*c)n， 但在数值计算时，由于误差的作用，多次叠代后解将被后一个解“淹没”而振荡地趋于无穷！例6.5说明了把一个连续系统离散化时, 时间步长h不能取得太大. 不然的话，得到的离散系统得解的性质可能会与连续系统的解的性质完全不同. 更有趣的是，对于非线性方程还会有本质上的变化：如轨道的极限集的变化，周期解的产生，混沌与分形现象！例6.5 设有一个连续的非线性系统, dx/dt=a*x(1-b*x); a0;b0, x(0)=x0, 0x01/b; 这模型称为logistic模型，用于描述比如人口等增长率不大的生物的增长. 很容易求得它的解为x(t)=1/(b+(1/x0-b)*exp(-a*t)这解的图像是一条S形曲线，当t趋于无穷大时，x趋于稳定平衡点1/b ，现离散这模型如下：记y(n)=x(t(n), t(n)=n*h, 用差商代替导数，得y(n+1)=y(n)+a*h*y(n)*(1-b*y(n)=(1+a*h)*y(n)*(1 a*b*h/(1+a*h)*y(n); 令z(n)=a*b*h/(1+a*h)*y(n); 则得离散系统z(n)=c*z(n)*(1-z(n); c=1+a*h，因为解z(n)满足 0 z(n)0，使得c=4，这是对h的第一个要求. 即使这样还不够， 经过以下的讨论可知，还要求0a*h2. 因为这个非线性离散模型称为虫口模型，可以描述昆虫、细菌等增长率很大的离散的增长行为，这时内禀增长率a会很大. 当1c3时,系统只有2个不动点：不稳定的不动点0及稳定的不动点1- 1/c，分别对应于离散系统的y=0及1/b，也分别对应于原连续系统的不稳定平衡点x=0及稳定平衡点x= 1/b，但当3c1+sqrt(6)= 3.449489742783178.时不动点0仍然是不稳定的;的稳定不动点1- 1/c变成了不稳定不动点，且从这个不动点分离出稳定的周期为2的离散极限环，由二点 (1+csqrt(1+c)*(3+c)/(2*c)组成, 这是连续系统没有的现象. 当1+sqrt(6)c3.544.到 3.5699456 时，极限环的个数及周期不断地倍增至无穷，当c进一步增大时，对于有些c值会出现混沌和分形现象，解对初值具有极大的敏感性！ 由于计算机的一般浮点运算有误差，因此，用计算机来验证时，只能得到稳定的极限集.对于离散的动力系统，可以使用计算机代数系统，如MATLAB的符号运算，来得到精确的结果. 以下是一段MATLAB的符号运算的程序，用于虫口模型的研究，请解读一下各语句的作用：syms x cf1=c*x*(1-x);F1=f1-x;F1=simple(F1);F1=factor(F1),S1=solve(F1),D1=diff(f1);D11=subs(D1,x,S1(1),D12=sub