数学建模动态模型.pptx

《数学建模动态模型.pptx》由会员分享，可在线阅读，更多相关《数学建模动态模型.pptx（88页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

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

3、1,0)在坐标轴上第四个解在两条直线 r1-a1x1-b1x2=0 r2-a2x2-b2x1=0 的交点：第9页/共88页如果两条直线不相交，则只存在3个平衡点。在这种情况下，两个种群不能共存。我们希望知道在什么条件下x10且x20.假设aibi,2个种群共存的条件：第10页/共88页平衡态第11页/共88页第五步，回答问题对每种种群存在两类增长限制。第一种是由于与另一种群的竞争，第二种是由于拥挤造成的同一种群内部的竞争。因此，对每一种树，存在着一点，在这一点数木由于拥挤会主动停止增长，且存在另一点，在这一点树木通过竞争阻止另一种群的增长。两种树能够共存的条件是每种树在达到限制自己增长的点之前

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

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

6、,x2=0第17页/共88页第四步，求解模型绘制向量场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)第18页/共88页M=10;%number of samples points x1min=0;x1max=900000;%dom

7、ain 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);%matlab routine axis(x1min x1max x2min x2max);title(Directio

8、n 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)第19页/共88页第20页/共88页四个平衡态解三个为：（0，0），（150000，0），（0，400000）第四个在区域内部，为唯一稳定的平衡态。第21页/共88页第五步，回答问题只要停止捕捞，鲸鱼种群将恢复到原来的水平，生态系统将处于稳定的平衡态。第22页/共88页灵敏性和稳定性对参数a做灵敏性

9、分析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),pretty(x2alpha)第23页/共88页解之得x1=150000*(-1+8000000*a)/Dx2=400000*(-1+1875000*a)/D其中，alpha1=solve(x1alpha);alpha2=solve(x2alpha);for

10、mat short e double(alpha1),double(alpha2)可见，对任意的a0,x20.第24页/共88页figure ezplot(x1alpha,0 8*10(-7),hold on grid on ezplot(x2alpha,0 8*10(-7),hold on title(Level of coexisting populations vs parameter alpha);第25页/共88页第26页/共88页a=linspace(0,9*10(-7);x1a=subs(x1alpha,alpha,a);x2a=subs(x2alpha,alpha,a);ind

11、=find(x1a0&x2a0);plot(a(ind),x1a(ind),bo);hold on plot(a(ind),x2a(ind),ro);第27页/共88页第28页/共88页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)第29页/共88页稳健性分析稳健性分析是考虑上述模型中f1和f2具有更一般的形式。只要

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

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

14、时间(秒)an=第n次调节后的加速度(米/秒)wn=等待到第n+1次观测前的等待时间(秒)第34页/共88页假设：tn+1=tn+tn=tn+cn+wnvn+1=vn+vn=vn+an-1cn+anwnan=-kvncn0wn=0目标：确定是否有vn0.第35页/共88页第二步，选择建模方法离散时间动力系统理论见教材，p99.第36页/共88页第三步，推导模型公式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-kcx2x2=x1-x2第37页/共88页

15、第四步，求解模型平衡态方程：-kwx1-kcx2=0 x1-x2=0平衡点（0，0）位于上述两直线的交点。下面绘制向量场F=(-kwx1-kcx2,x1-x2)第38页/共88页绘制向量场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);a

16、xis(x1min x1max x2min x2max);title(Direction field(the vectors are rescaled!);hold on xlabel(Current Speed),ylabel(Previous Speed);第39页/共88页第40页/共88页系统演变定义函数function rhs=dockfun(x,c,w,k);rhs=-k*w*x(1)-k*c*x(2);x(1)-x(2);第41页/共88页迭代实现x=8;10;N=12;fprintf(n Current speed Prev.speednn)fprintf(%2.0f%5.2f

17、%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 第42页/共88页第43页/共88页第五步，回答问题假设cw,有x1 -kwx1当kw2时，我们将得到一个稳定的平衡态。当kw1时，系统将逼近平衡态，而不会越过它。据此，可以回答问题如下：只要控制调节不是太剧烈，则控制将起作用。第44页/共88页进一步，两次调节之间的

18、间隔时间越长，调节幅度必须越小。而且，它们之间呈反比。如果两次调节之间的时间间隔增加两倍，调节幅度可以减半。特别地，如果我们以10秒调节一次，则设置加速度为1/10，以免超过目标速度零。第45页/共88页作业：习题4.4 1，11第46页/共88页第五章 动态模型分析连续时间系统的特征值方法离散时间系统的特征值方法相图第47页/共88页连续时间系统特征值方法例5.1 再次考虑例的树木问题。假设硬材树每年增长率为10%，软材树每年增长率为25%。一英亩林地可以提供大约10000吨的硬木或6000吨的软木。竞争的程度还未从数值上确定。两种树能否共存于一个稳定的平衡态？第48页/共88页第一步，提出

19、问题r1r2a1a2第49页/共88页第二步，选择建模方法连续时间动力系统的特征值分析法见教材，p112第50页/共88页第三步，推导模型公式已经确定了r1,r2,a1,a2,我们仍假设biai,不妨取bi=ai/2,则动态系统方程为x=F(x),其中F(x)=(f1,f2).f1(x1,x21-(0.10/10000)(x1)2-(0.05/10000)x1x2f2(x1,x22-(0.25/6000)(x2)2-(0.125/6000)x1x2第51页/共88页第四步，求解模型区域内部有平衡点：x109333,x201333计算偏导数矩阵在平衡点(x10,x20)处的值A。求矩阵A的特征值

20、，得第52页/共88页这两个特征值均具有负实部，所以该平衡点为稳定的。上述过程也可以用计算机求解。第53页/共88页代码实现-求平衡点syms x1 x2 f1=0.1*x1*(1-(1/10000)*x1-0.5*(1/10000)*x2);f2=0.25*x2*(1-(1/6000)*x2-0.5*(1/6000)*x1);x1steady,x2steady=solve(f1,f2);N=length(x1steady);fprintf(The equilibrium points aren)disp(x1steady x2steady)第54页/共88页代码实现-绘制方向场M=15;x1

21、min=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(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,d

22、X1,dX2),axis(x1min x1max x2min x2max);title(Direction field(the vectors are rescaled!);xlabel(Hardwoods),ylabel(Softwoods),hold on 第55页/共88页第56页/共88页代码实现-稳定性分析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,x1num,x2num)la

23、mbda=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 第57页/共88页第58页/共88页第五步，回答问题我们发现硬材树和软材树可以共存于一个平衡态。在一个成熟的稳定的树林中，每英亩大约有9300吨硬材树和1300吨软材树。这个结论基于对两类树种之间竞争

24、程度的近似合理的假设。第59页/共88页灵敏性分析为做灵敏性分析，这里放松假设bi=(1/2)ai,而假设bi=tai.条件:biairi/airj/bj隐含着第60页/共88页定义系统方程syms x1 x2 tf1=0.1*x1*(1-(1/10000)*x1-t*(1/10000)*x2);f2=0.25*x2*(1-(1/6000)*x2-t*(1/6000)*x1);第61页/共88页求平衡态x1steady,x2steady=solve(f1,f2);N=length(x1steady);fprintf(The equilibrium points aren)disp(x1stea

25、dy x2steady)第62页/共88页每一平衡点对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)lambda=eig(A)figurefor tty=subs(lambda,t,tt);plot(tt,y),hold onend end 第63页/共88页平衡态（0，0）第64页/共88页平衡态（150000，0）第65页/共88页平衡态（0，400000）第66页/共88页区域内部的平衡态第67页/共88页5.3

26、相图所谓连续时间动态系统的相图就是连续时间动态系统有代表性的解曲线在状态空间的草图。通常与向量场的勾画一起运用，从而获得动态系统在状态空间动态行为的图形描述。第68页/共88页例5.3 考虑如下图所示的电路图。电路由一个电容，一个电阻和一个电感器构成一个简单闭路。电路中每个元件的作用由在这个回路中的电流和电压之间的关系表示。一个理想的物理模型给出这个关系：第69页/共88页第70页/共88页称函数f(x)为电阻的v-i特征。在古典的RLC电路理论中我们假设f(x)=Rx,其中R表示电阻。基尔霍夫电流定理说明：进入一个节点的电流之和等于流出电流之和。基尔霍夫电压定理说明：闭路上所有电压差之和为零

27、。对情形L=1,C=1/3和f(x)=x3+4x确定这个电路随时间变化的行为。第71页/共88页第一步，提出问题变量：vC,iC,vR,iR,vL,iL假设CdvC/dt=iCvR=f(iR)LdiL/dt=vLiC=iR=iL,vC+vR+vL=0L=1,C=1/3,f(x)=x3+4x目标：确定六个变量随时间变化的性质。第72页/共88页第二步，选择建模方法建立连续时间动态系统模型勾画完整的相图进行分析。第73页/共88页第三步，构造模型设x1=iR=iL=iC,x2=vC,有第74页/共88页代入，得第75页/共88页重新排列得到第76页/共88页第四步，求解模型我们将通过勾画完整的相图

28、分析该动态系统。第77页/共88页第78页/共88页对应的线性系统的相图第79页/共88页第五步，回答问题问题是描述RLC电路的行为。系统的性质可以用两个量描述：通过电阻的电流和电容上的电压降。不论系统的初始状态如何，这两个量最终将趋于零。第80页/共88页灵敏性分析：考察C的灵敏性（原C=1/3)第81页/共88页稳健性考虑RLC电路的v-i特征f(x)(原f(x)=x3+4x).更一般地，假设f(0)=0,f(x)严格单调增加。动态系统方程为：第82页/共88页例考虑非线性RLC电路，L=1,C=1且v-i特征f(x)=x3-x.确定这个电路随时间变化的行为。和前面分析一样，得到如下动态系统：第83页/共88页从向量场难于判断平衡点是否稳定！第84页/共88页相图第85页/共88页第86页/共88页作业第五章习题9，10，11，14，15第87页/共88页感谢您的观看！第88页/共88页