数学建模课件 第二部分 动态模型

1、数学建模三峡大学 理学院 俞 辉Email:第二部分 动态模型 动态模型介绍 动态模型分析 动态模型模拟第4章 动态模型介绍 常态分析 动力系统 离散时间动力系统4.1 常态分析 例4.1 在一个未被管理的森林，硬材树和软材树竞争可用的土地和水分。越可用的硬材树生长的越慢，但越耐久且提供越有价值的木材。软材树靠生长快，有效消耗水分和土壤养分与硬材树竞争。硬材树靠生长的高度与软材树竞争，它们遮挡了小树的阳光，它们也更耐抗疾病。这两种树能否同时在一片森林中共存，或者一种树会迫使另一种树灭绝？五步法 H和S分别表示硬材树和软材树种群。生物学家习惯使用的计量单位是每英亩上的木材吨数。 无限制生长(丰富

2、的空间、阳光、水分、土壤养料等）：rP 种群内的竞争：-aP2（小种群的增长率线性依赖于种群的大小 ，即：-aP) 种群生长(率)函数：g(P)= rP -aP2，(r为内禀增长率，a=0,S=0 r1,r2,a1,a2,b1,b2是正实数 目标：是否有H-0或S-0。第二步，选择建模方法 微分方程动力系统理论 见书P.90第三步，构造模型公式 记x1=H,x2=S为两个状态变量，定义在状态空间：(x1,x2):x1=0,x2=0 定常态方程r1x1-a1x12-b1x1x2=0r2x2-a2x22-b2x1x2=0第四步，求解模型 得四个解： 三个解(0，0),(0，r2/a2),(r1/a

3、1,0)在坐标轴上 第四个解在两条直线 r1-a1x1-b1x2=0 r2-a2x2-b2x1=0 的交点：122 11 22 112121 2121 2,rar brb rxxaabbbbaaa 如果两条直线不相交，则只存在3个平衡点。在这种情况下，两个种群不能共存。 我们希望知道在什么条件下x10且x20. 假设aibi,2个种群共存的条件：21122112rrrrabab且平衡态1x2x11rb22rb22ra11ra第五步，回答问题 对每种种群存在两类增长限制。第一种是由于与另一种群的竞争，第二种是由于拥挤造成的同一种群内部的竞争。因此，对每一种树，存在着一点，在这一点数木由于拥挤会主

4、动停止增长，且存在另一点，在这一点树木通过竞争阻止另一种群的增长。两种树能够共存的条件是每种树在达到限制自己增长的点之前已经达到它限制另一种树增长的点。4.2 动力系统 例 4.2 蓝鲸和长须鲸是生活在同一海域的相似种群，因此认为他们之间存在竞争。蓝鲸的内禀增长率每年估计为5%，长须鲸为每年8%，环境承载力（环境能够支付的鲸鱼的最大数量）估计蓝鲸为150000条，长须鲸为400000条。鲸鱼竞争的程度是未知的。在过去的100年剧烈的捕捞已经使鲸鱼数量减少，蓝鲸大约为5000条，长须鲸大约为70000条。蓝鲸是否会灭绝？第一步，提出问题 变量： B=蓝鲸的数量 F=长须鲸的数量 gB=蓝鲸种群的

5、增长率（每年） gF=长须鲸种群的增长率（每年） cB=蓝鲸与长须鲸竞争的影响（每年的鲸鱼数） cF=长须鲸与蓝鲸竞争的影响（每年的鲸鱼数） 假设 gB=0.05B(1-B/150000) gF=0.08F(1-F/400000) cB=cF=aBF B=0,F=0, a是正实数 目标：确定动力系统是否能够从B=5000, F=70000开始达到稳定的平衡态。第二步，选择建模方法 连续时间动力系统理论 见教材p.94.第三步，构造模型公式 令x1=B,x2=F,记x1=f1(x1,x2),x2=f2(x1,x2) f1(x1,x2)=0.05x1(1-x1/150000)-ax1x2 f2(x

6、1,x2)=0.08x2(1-x2/400000)-ax1x2 状态空间为 S=(x1,x2):x1=0,x2=0第四步，求解模型 绘制向量场clear all, close all, clc syms x1 x2 alpha=10(-7); f1=0.05*x1*(1-x1/150000) - alpha*x1*x2;f2=0.08*x2*(1-x2/400000) - alpha*x1*x2; x1steady,x2steady=solve(f1,f2); disp(The equilibrium points are) disp(x1steady x2steady) M=10; % nu

7、mber of samples points x1min=0; x1max=900000; % domain specification x2min=0; x2max=600000; X1,X2=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); dX1=0.05*X1.*(1-X1/150000) - alpha*X1.*X2; % x1-componentdX2=0.08*X2.*(1-X2/400000) - alpha*X1.*X2; % x2-component quiver(X1,X2,dX1,dX2

8、); % matlab routine axis(x1min x1max x2min x2max); title(Direction field (the vectors may be rescaled!); hold on xlabel(Blue Whales); ylabel(Fin Whales); ezplot(f1,0 900000 0 600000), hold on ezplot(f2,0 900000 0 600000) 0123456789x 1050123456x 1052/25 x2 (1-1/400000 x2)-1/10000000 x1 x2 = 0 x1x2四个平

9、衡态解 三个为： （0，0），（150000，0），（0，400000） 第四个在区域内部，为唯一稳定的平衡态。第五步，回答问题 只要停止捕捞，鲸鱼种群将恢复到原来的水平，生态系统将处于稳定的平衡态。灵敏性和稳定性 对参数a做灵敏性分析syms alpha f1=.05*x1*(1-x1/150000)-alpha*x1*x2; f2=.08*x2*(1-x2/400000) - alpha*x1*x2; x1steady,x2steady = solve(f1,f2) x1alpha=x1steady(4); x2alpha=x2steady(4); pretty(x1alpha), pre

10、tty(x2alpha) 解之得x1=150000*(-1+8000000*a)/Dx2=400000*(-1+1875000*a)/D其中，D=-1+15000000000000*a2alpha1 = solve(x1alpha); alpha2 = solve(x2alpha); format short e double(alpha1), double(alpha2) 可见，对任意的a0,x20.figure ezplot(x1alpha,0 8*10(-7), hold on grid on ezplot(x2alpha,0 8*10(-7), hold on title(Level

11、of coexisting populations vs parameter alpha); 012345678x 10-7-6-4-20246810 x 105Level of coexisting populations vs parameter a = linspace(0,9*10(-7); x1a=subs(x1alpha,alpha,a); x2a=subs(x2alpha,alpha,a); ind=find(x1a0 & x2a0); plot(a(ind),x1a(ind),bo); hold on plot(a(ind), x2a(ind),ro); 012345678x

12、10-7-6-4-20246810 x 105Level of coexisting populations vs parameter format bank Sx1alpha = diff(x1alpha, alpha)*(alpha/x1alpha); Sx2alpha = diff(x2alpha, alpha)*(alpha/x2alpha); Sx1a = subs(Sx1alpha, alpha, 10(-7) Sx2a = subs(Sx2alpha, alpha, 10(-7) 稳健性分析 稳健性分析是考虑上述模型中f1和f2具有更一般的形式。 只要向量场具有相同的一般特征，我

13、们的结论仍然是正确的。4.3 离散时间动力系统 例4.3 宇航员在训练中要求用手动控制做对接演习。作为这个演习的一部分，要求保持一个正在运行的太空船与另一个正在运行的太空船的相对位置。手控制器提供了不同的加速度和减速度，并且在太空船上有一个装置测量这两个飞船的接近速度。建议使用如下的策略进行飞船对接。 首先观察接近速度。如果为零，则不用再做任何事情。否则，记住这个接近速度，再看加速度控制器，控制加速度使得它与接近速度相反（即如果接近速度是正值，则放慢，如果是负的，则加快。），且正比于这个差值（即如果发现接近速度达到2倍时，我们将于2倍的速度刹车）。经过一段时间，再观察接近速度并重复上面的步骤。

14、在什么环境下这个策略才是有效的？第一步，提出问题 设vn表示在时间tn观测到的接近速度，tn为第n次观测的时间。 太空船接近速度的改变：vn =vn+1-vn 两次观测之间的时间间隔：tn=tn+1-tn 时间区间被分成两部分： tn=cn+wn cn为调整控制器的时间，wn为下一次观测前的等待时间。 记an为第n次调节后设定的加速度，则vn =an-1cn+anwn 按控制律要求加速度正比于(-vn),因此，an=-kvn变量： tn=第n次观测速度的时间(秒) vn=在tn时刻的速度(米/秒) cn=执行第n次控制调节的时间(秒) an=第n次调节后的加速度(米/秒) wn=等待到第n+1

15、次观测前的等待时间(秒)假设： tn+1=tn+tn=tn+cn+wn vn+1=vn+vn=vn+an-1cn+anwn an=-kvn cn0 wn=0 目标：确定是否有vn0.第二步，选择建模方法 离散时间动力系统理论 见教材，p99.第三步，推导模型公式 vn+1-vn=-kvn-1cn-kvnwn 为简化起见，对所有的n，设cn=c,wn=w 有vn+1-vn=-kwvn-kcvn-1 设x1(n)=vn,x2(n)=vn-1,则 x1 =-kwx1-kcx2 x2=x1-x2第四步，求解模型 平衡态方程： -kwx1-kcx2=0 x1-x2=0 平衡点（0，0）位于上述两直线的交

16、点。 下面绘制向量场F=(-kwx1-kcx2, x1-x2)绘制向量场clear all, close all, clc c = 5,w = 10,k = 0.1; M=10,x1min=-15,x1max=15,x2min=-15,x2max=15;X1,X2=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max);dX1=-k*w*X1-k*c*X2; dX2= X1-X2; quiver(X1,X2,dX1,dX2); axis(x1min x1max x2min x2max); title(Direction

17、 field (the vectors are rescaled!); hold on xlabel(Current Speed), ylabel(Previous Speed); -15-10-5051015-15-10-5051015Direction field (the vectors are rescaled!)Current SpeedPrevious Speed系统演变 定义函数dockfun.mfunction rhs = dockfun(x,c,w,k); rhs = -k*w*x(1)-k*c*x(2); x(1)-x(2);迭代实现x=8;10; N = 12; fpri

18、ntf(n Current speed Prev. speednn) fprintf(%2.0f %5.2f %5.2fn, 0, x(1), x(2)for n=1:N xnew = x + dockfun(x,c,w,k); plot(x(1),xnew(1), x(2),xnew(2),-ro,. MarkerFaceColor,k,MarkerSize,2) x = xnew; fprintf(%2.0f %5.2f %5.2fn, n, x(1), x(2) end -15-10-5051015-15-10-5051015Direction field (the vectors ar

19、e rescaled!)Current SpeedPrevious Speed第五步，回答问题 假设cw,有x1 -kwx1 当kw2时，我们将得到一个稳定的平衡态。 当kw1时，系统将逼近平衡态，而不会越过它。 据此，可以回答问题如下： 只要控制调节不是太剧烈，则控制将起作用。 进一步，两次调节之间的间隔时间越长，调节幅度必须越小。而且，它们之间呈反比。如果两次调节之间的时间间隔增加两倍，调节幅度可以减半。 特别地，如果我们以10秒调节一次，则设置加速度为1/10，以免超过目标速度零。 作业： 习题4.4 1，11第五章 动态模型分析 连续时间系统的特征值方法 离散时间系统的特征值方法 相图

20、连续时间系统特征值方法 例5.1 再次考虑例4.1的树木问题。假设硬材树每年增长率为10%，软材树每年增长率为25%。一英亩林地可以提供大约10000吨的硬木或6000吨的软木。竞争的程度还未从数值上确定。两种树能否共存于一个稳定的平衡态？第一步，提出问题 r1=0.10 r2=0.25 a1=0.10/10000 a2=0.25/6000第二步，选择建模方法 连续时间动力系统的特征值分析法 见教材，p112第三步，推导模型公式 已经确定了r1,r2,a1,a2,我们仍假设biai,不妨取bi=ai/2,则 动态系统方程为x=F(x),其中F(x)=(f1,f2). f1(x1,x2)=0.1

21、0 x1-(0.10/10000)(x1)2-(0.05/10000)x1x2 f2(x1,x2)=0.25x2-(0.25/6000)(x2)2-(0.125/6000)x1x2第四步，求解模型 区域内部有平衡点： x109333,x201333 计算偏导数矩阵在平衡点(x10,x20)处的值A。 求矩阵A的特征值，得671339900 这两个特征值均具有负实部，所以该平衡点为稳定的。 上述过程也可以用计算机求解。代码实现-求平衡点syms x1 x2 f1 = 0.1*x1*(1-(1/10000)*x1- 0.5*(1/10000)*x2); f2 = 0.25*x2*(1-(1/600

22、0)*x2- 0.5*(1/6000)*x1); x1steady,x2steady = solve(f1,f2); N=length(x1steady); fprintf(The equilibrium points aren) disp(x1steady x2steady) 代码实现-绘制方向场M=15; x1min=0; x1max=11000; x2min=0; x2max=8000; ezplot(f1,x1min x1max x2min x2max),hold on ezplot(f2,x1min x1max x2min x2max),hold on X1,X2=meshgrid(

23、x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max); dX1 = 0.1*X1.*(1-(1/10000)*X1-0.5*(1/10000)*X2); dX2 = 0.25*X2.*(1-(1/6000)*X2-0.5*(1/6000)*X1); quiver(X1,X2,dX1,dX2),axis(x1min x1max x2min x2max); title(Direction field (the vectors are rescaled!); xlabel(Hardwoods),ylabel(Softwoods),hold

24、on 01000 20003000 40005000 60007000 80009000 10000 11000010002000300040005000600070008000HardwoodsSoftwoodsDirection field (the vectors are rescaled!)代码实现-稳定性分析DF = diff(f1,x1), diff(f1,x2); diff(f2,x1), diff(f2,x2); for i=1:N x1num = double(x1steady(i); x2num = double(x2steady(i); A=subs(DF,x1,x2,x

25、1num,x2num) lambda = eig(A); fprintf(The eigenvalues for the equilibrium () fprintf(%1.0f, %1.0f,x1num, x2num); fprintf() are ) fprintf(%1.2f %1.2fn,lambda(1), lambda(2); plot(x1num,x2num,ro,MarkerSize,10, MarkerFaceColor,g); end 01000 20003000 40005000 60007000 80009000 10000 1100001000200030004000

26、5000600070008000HardwoodsSoftwoodsDirection field (the vectors are rescaled!)第五步，回答问题 我们发现硬材树和软材树可以共存于一个平衡态。在一个成熟的稳定的树林中，每英亩大约有9300吨硬材树和1300吨软材树。这个结论基于对两类树种之间竞争程度的近似合理的假设。灵敏性分析 为做灵敏性分析，这里放松假设bi=(1/2)ai,而假设bi=tai. 条件: biai ri/airj/bj 隐含着0t0.6定义系统方程syms x1 x2 tf1 = 0.1*x1*(1-(1/10000)*x1- t*(1/10000)*

27、x2); f2 = 0.25*x2*(1-(1/6000)*x2- t*(1/6000)*x1); 求平衡态x1steady,x2steady = solve(f1,f2); N=length(x1steady); fprintf(The equilibrium points aren) disp(x1steady x2steady) 每一平衡点对t(0t0.6)的灵敏性DF = diff(f1,x1), diff(f1,x2); diff(f2,x1), diff(f2,x2) for i=1:N A=subs(DF,x1,x2,x1steady(i) , x2steady(i) lambd

28、a = eig(A) figurefor tt=0:0.01:0.6y=subs(lambda,t,tt);plot(tt,y),hold onend end 平衡态（0，0）00.10.20.30.40.50.60.70.10.150.20.25平衡态（150000，0）00.10.20.30.40.50.60.7-0.1-0.0500.050.10.150.20.25平衡态（0，400000）00.10.20.30.40.50.60.7-0.25-0.2-0.15-0.1-0.0500.050.1区域内部的平衡态00.10.20.30.40.50.60.7-0.25-0.2-0.15-0.

29、1-0.0505.3 相图 所谓连续时间动态系统的相图就是连续时间动态系统有代表性的解曲线在状态空间的草图。 通常与向量场的勾画一起运用，从而获得动态系统在状态空间动态行为的图形描述。例5.3 考虑如下图所示的电路图。电路由一个电容，一个电阻和一个电感器构成一个简单闭路。电路中每个元件的作用由在这个回路中的电流和电压之间的关系表示。一个理想的物理模型给出这个关系：( )CCRRLLdvCidtvf idivdt（电容）（电阻）L（电感）RCL 称函数f(x)为电阻的v-i特征。在古典的RLC电路理论中我们假设f(x)=Rx,其中R表示电阻。基尔霍夫电流定理说明：进入一个节点的电流之和等于流出电流之和。基尔霍夫电压定理说明：闭路上所有电压差之和为零。 对情形L=1,C=1/3和f(x)=x3+4x确定这个电路随时间变化的行为。第一步，提出问题 变量：vC,iC,vR,iR,vL,iL 假设 CdvC/dt=iC vR=f(iR) LdiL/dt=vL iC=iR=iL,vC+vR+vL=0 L=1,C=1/3,f(x)=x3+4x