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

数学建模三峡大学理学院俞辉Email:yuhui@精选课件第二局部动态模型动态模型介绍动态模型分析动态模型模拟精选课件第4章动态模型介绍常态分析动力系统离散时间动力系统精选课件4.1常态分析例4.1在一个未被管理的森林，硬材树和软材树竞争可用的土地和水分。越可用的硬材树生长的越慢，但越耐久且提供越有价值的木材。软材树靠生长快，有效消耗水分和土壤养分与硬材树竞争。硬材树靠生长的高度与软材树竞争，它们遮挡了小树的阳光，它们也更耐抗疾病。这两种树能否同时在一片森林中共存，或者一种树会迫使另一种树灭绝？精选课件五步法H和S分别表示硬材树和软材树种群。生物学家习惯使用的计量单位是每英亩上的木材吨数。无限制生长(丰富的空间、阳光、水分、土壤养料等〕：rP种群内的竞争：-aP2〔小种群的增长率线性依赖于种群的大小，即：-aP)种群生长(率)函数：g(P)=rP-aP2，(r为内禀增长率，a<<r是资源限制强度系数)种群间的竞争：-bSH精选课件变量与假设变量：H=硬材树种群〔吨/亩〕S=软材树种群〔吨/亩〕gH=硬材树的生长率〔吨/英亩/年〕gS=软材树的生长率〔吨/英亩/年〕cH=与软材树竞争的损失〔吨/英亩/年〕cs=与硬材树竞争的损失〔吨/英亩/年〕精选课件假设gH=r1H-a1H2gS==r2S-a2S2cH=b1SHcS=b2SHH>=0,S>=0r1,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/a1,0)在坐标轴上第四个解在两条直线r1-a1x1-b1x2=0r2-a2x2-b2x1=0的交点：精选课件如果两条直线不相交，那么只存在3个平衡点。在这种情况下，两个种群不能共存。我们希望知道在什么条件下x1>0且x2>0.假设ai>bi,2个种群共存的条件：精选课件平衡态精选课件第五步，答复以下问题对每种种群存在两类增长限制。第一种是由于与另一种群的竞争，第二种是由于拥挤造成的同一种群内部的竞争。因此，对每一种树，存在着一点，在这一点数木由于拥挤会主动停止增长，且存在另一点，在这一点树木通过竞争阻止另一种群的增长。两种树能够共存的条件是每种树在到达限制自己增长的点之前已经到达它限制另一种树增长的点。精选课件4.2动力系统例4.2蓝鲸和长须鲸是生活在同一海域的相似种群，因此认为他们之间存在竞争。蓝鲸的内禀增长率每年估计为5%，长须鲸为每年8%，环境承载力〔环境能够支付的鲸鱼的最大数量〕估计蓝鲸为150000条，长须鲸为400000条。鲸鱼竞争的程度是未知的。在过去的100年剧烈的捕捞已经使鲸鱼数量减少，蓝鲸大约为5000条，长须鲸大约为70000条。蓝鲸是否会灭绝？精选课件第一步，提出问题变量：B=蓝鲸的数量F=长须鲸的数量gB=蓝鲸种群的增长率〔每年〕gF=长须鲸种群的增长率〔每年〕cB=蓝鲸与长须鲸竞争的影响〔每年的鲸鱼数〕cF=长须鲸与蓝鲸竞争的影响〔每年的鲸鱼数〕精选课件假设gB=0.05B(1-B/150000)gF=0.08F(1-F/400000)cB=cF=aBFB>=0,F>=0,a是正实数目标：确定动力系统是否能够从B=5000,F=70000开始到达稳定的平衡态。精选课件第二步，选择建模方法连续时间动力系统理论见教材p.94.精选课件第三步，构造模型公式令x1=B,x2=F,记x’1=f1(x1,x2),x’2=f2(x1,x2)f1(x1,x2)=0.05x1(1-x1/150000)-ax1x2f2(x1,x2)=0.08x2(1-x2/400000)-ax1x2状态空间为S={(x1,x2):x1>=0,x2>=0}精选课件第四步，求解模型绘制向量场clearall,closeall,clcsymsx1x2alpha=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('Theequilibriumpointsare')disp([x1steadyx2steady])精选课件M=10;%numberofsamplespointsx1min=0;x1max=900000;%domainspecificationx2min=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-componentquiver(X1,X2,dX1,dX2);%matlabroutineaxis([x1minx1maxx2minx2max]);title('Directionfield(thevectorsmayberescaled!)');holdonxlabel('BlueWhales');ylabel('FinWhales');ezplot(f1,[09000000600000]),holdonezplot(f2,[09000000600000])精选课件精选课件四个平衡态解三个为：〔0，0〕，〔150000，0〕，〔0，400000〕第四个在区域内部，为唯一稳定的平衡态。精选课件第五步，答复以下问题只要停止捕捞，鲸鱼种群将恢复到原来的水平，生态系统将处于稳定的平衡态。精选课件灵敏性和稳定性对参数a做灵敏性分析symsalphaf1=.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)精选课件解之得x1=150000*(-1+8000000*a)/Dx2=400000*(-1+1875000*a)/D其中，D=-1+15000000000000*a^2alpha1=solve(x1alpha);alpha2=solve(x2alpha);formatshortedouble(alpha1),double(alpha2)可见，对任意的a<1.25*10-7存在一个稳定的平衡态x1>0,x2>0.精选课件figureezplot(x1alpha,[08*10^(-7)]),holdongridonezplot(x2alpha,[08*10^(-7)]),holdontitle('Levelofcoexistingpopulationsvsparameter\alpha');精选课件精选课件a=linspace(0,9*10^(-7));x1a=subs(x1alpha,alpha,a);x2a=subs(x2alpha,alpha,a);ind=find(x1a>0&x2a>0);plot(a(ind),x1a(ind),'bo');holdonplot(a(ind),x2a(ind),'ro');精选课件精选课件formatbankSx1alpha=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具有更一般的形式。只要向量场具有相同的一般特征，我们的结论仍然是正确的。精选课件4.3离散时间动力系统例4.3宇航员在训练中要求用手动控制做对接演习。作为这个演习的一局部，要求保持一个正在运行的太空船与另一个正在运行的太空船的相对位置。手控制器提供了不同的加速度和减速度，并且在太空船上有一个装置测量这两个飞船的接近速度。建议使用如下的策略进行飞船对接。精选课件首先观察接近速度。如果为零，那么不用再做任何事情。否那么，记住这个接近速度，再看加速度控制器，控制加速度使得它与接近速度相反〔即如果接近速度是正值，那么放慢，如果是负的，那么加快。〕，且正比于这个差值〔即如果发现接近速度到达2倍时，我们将于2倍的速度刹车〕。经过一段时间，再观察接近速度并重复上面的步骤。在什么环境下这个策略才是有效的？精选课件第一步，提出问题设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精选课件变量：tn=第n次观测速度的时间(秒)vn=在tn时刻的速度(米/秒)cn=执行第n次控制调节的时间(秒)an=第n次调节后的加速度(米/秒)wn=等待到第n+1次观测前的等待时间(秒)精选课件假设：tn+1=tn+Δtn=tn+cn+wnvn+1=vn+Δvn=vn+an-1cn+anwnan=-kvncn>0wn>=0目标：确定是否有vn→0.精选课件第二步，选择建模方法离散时间动力系统理论见教材，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=0x1-x2=0平衡点〔0，0〕位于上述两直线的交点。下面绘制向量场F=(-kwx1-kcx2,x1-x2)精选课件绘制向量场clearall,closeall,clcc=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([x1minx1maxx2minx2max]);title('Directionfield(thevectorsarerescaled!)');holdonxlabel('CurrentSpeed'),ylabel('PreviousSpeed');精选课件精选课件系统演变定义函数dockfun.mfunctionrhs=dockfun(x,c,w,k);rhs=[-k*w*x(1)-k*c*x(2);x(1)-x(2)];精选课件迭代实现x=[8;10];N=12;fprintf('nCurrentspeedPrev.speed\n\n')fprintf('%2.0f%5.2f%5.2f\n',0,x(1),x(2))forn=1:Nxnew=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.2f\n',n,x(1),x(2))end精选课件精选课件第五步，答复以下问题假设c<<w,有Δx1≈-kwx1当kw<2时，我们将得到一个稳定的平衡态。当kw<1时，系统将逼近平衡态，而不会越过它。据此，可以答复以下问题如下：只要控制调节不是太剧烈，那么控制将起作用。精选课件进一步，两次调节之间的间隔时间越长，调节幅度必须越小。而且，它们之间呈反比。如果两次调节之间的时间间隔增加两倍，调节幅度可以减半。特别地，如果我们以10秒调节一次，那么设置加速度为1/10，以免超过目标速度零。精选课件作业：习题4.41，11精选课件第五章动态模型分析连续时间系统的特征值方法离散时间系统的特征值方法相图精选课件连续时间系统特征值方法例5.1再次考虑例4.1的树木问题。假设硬材树每年增长率为10%，软材树每年增长率为25%。一英亩林地可以提供大约10000吨的硬木或6000吨的软木。竞争的程度还未从数值上确定。两种树能否共存于一个稳定的平衡态？精选课件第一步，提出问题r1=0.10r2=0.25a1=0.10/10000a2=0.25/6000精选课件第二步，选择建模方法连续时间动力系统的特征值分析法见教材，p112精选课件第三步，推导模型公式已经确定了r1,r2,a1,a2,我们仍假设bi<ai,不妨取bi=ai/2,那么动态系统方程为x’=F(x),其中F(x)=(f1,f2).f1(x1,x2)=0.10x1-(0.10/10000)(x1)2-(0.05/10000)x1x2f2(x1,x2)=0.25x2-(0.25/6000)(x2)2-(0.125/6000)x1x2精选课件第四步，求解模型区域内部有平衡点：x10≈9333,x20≈1333计算偏导数矩阵在平衡点(x10,x20)处的值A。求矩阵A的特征值，得精选课件这两个特征值均具有负实部，所以该平衡点为稳定的。上述过程也可以用计算机求解。精选课件代码实现---求平衡点symsx1x2f1=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('Theequilibriumpointsare\n')disp([x1steadyx2steady])精选课件代码实现---绘制方向场M=15;x1min=0;x1max=11000;x2min=0;x2max=8000;ezplot(f1,[x1minx1maxx2minx2max]),holdonezplot(f2,[x1minx1maxx2minx2max]),holdon[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,dX1,dX2),axis([x1minx1maxx2minx2max]);title('Directionfield(thevectorsarerescaled!)');xlabel('Hardwoods'),ylabel('Softwoods'),holdon精选课件精选课件代码实现---稳定性分析DF=[diff(f1,x1),diff(f1,x2);diff(f2,x1),diff(f2,x2)];fori=1:Nx1num=double(x1steady(i));x2num=double(x2steady(i));A=subs(DF,[x1,x2],[x1num,x2num])lambda=eig(A);fprintf('Theeigenvaluesfortheequilibrium(')fprintf('%1.0f,%1.0f',x1num,x2num);fprintf(')are')fprintf('%1.2f%1.2f\n',lambda(1),lambda(2));plot(x1num,x2num,'ro','MarkerSize',10,'MarkerFaceColor','g');end精选课件精选课件第五步，答复以下问题我们发现硬材树和软材树可以共存于一个平衡态。在一个成熟的稳定的树林中，每英亩大约有9300吨硬材树和1300吨软材树。这个结论基于对两类树种之间竞争程度的近似合理的假设。精选课件灵敏性分析为做灵敏性分析，这里放松假设bi=(1/2)ai,而假设bi=tai.条件:bi<airi/ai<rj/bj隐含着0<t<0.6精选课件定义系统方程symsx1x2tf1=0.1*x1*(1-(1/10000)*x1-t*(1/10000)*x2);f2=0.25*x2*(1-(1/6000)*x2-t*(1/6000)*x1);精选课件求平衡态[x1steady,x2steady]=solve(f1,f2);N=length(x1steady);fprintf('Theequilibriumpointsare\n')disp([x1steadyx2steady])精选课件每一平衡点对t(0<t<0.6)的灵敏性DF=[diff(f1,x1),diff(f1,x2);diff(f2,x1),diff(f2,x2)]fori=1:NA=subs(DF,[x1,x2],[x1steady(i),x2steady(i)])lambda=eig(A)figurefortt=0:0.01:0.6y=subs(lambda,t,tt);plot(tt,y),holdonendend精选课件平衡态〔0，0〕精选课件平衡态〔150000，0〕精选课件平衡态〔0，400000〕精选课件区域内部的平衡态精选课件5.3相图所谓连续时间动态系统的相图就是连续时间动态系统有代表性的解曲线在状态空间的草