《数学建模 第二部分 动态模型ppt课件》由会员分享，可在线阅读，更多相关《数学建模 第二部分 动态模型ppt课件（88页珍藏版）》请在金锄头文库上搜索。
1、.,数学建模,三峡大学 理学院 俞 辉 Email:,.,第二部分 动态模型,动态模型介绍 动态模型分析 动态模型模拟,.,第4章 动态模型介绍,常态分析 动力系统 离散时间动力系统,.,4.1 常态分析,例4.1 在一个未被管理的森林，硬材树和软材树竞争可用的土地和水分。越可用的硬材树生长的越慢，但越耐久且提供越有价值的木材。软材树靠生长快，有效消耗水分和土壤养分与硬材树竞争。硬材树靠生长的高度与软材树竞争，它们遮挡了小树的阳光，它们也更耐抗疾病。这两种树能否同时在一片森林中共存，或者一种树会迫使另一种树灭绝？,.,五步法,H和S分别表示硬材树和软材树种群。生物学家习惯使用的计量单位是每英亩
2、上的木材吨数。 无限制生长(丰富的空间、阳光、水分、土壤养料等）：rP 种群内的竞争：-aP2（小种群的增长率线性依赖于种群的大小 ，即：-aP) 种群生长(率)函数：g(P)= rP -aP2，(r为内禀增长率，ar是资源限制强度系数) 种群间的竞争：-bSH,.,变量与假设,变量： H=硬材树种群（吨/亩） S=软材树种群（吨/亩） gH=硬材树的生长率（吨/英亩/年） gS=软材树的生长率（吨/英亩/年） cH=与软材树竞争的损失（吨/英亩/年） cs=与硬材树竞争的损失（吨/英亩/年）,.,假设 gH=r1H-a1H2 gS=r2S-a2S2 cH=b1SH cS=b2SH H=0,S
3、=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=0 r2x2-a2x22-b2x1x2=0,.,第四步，求解模型,得四个解： 三个解(0，0),(0，r2/a2),(r1/a1,0)在坐标轴上 第四个解在两条直线 r1-a1x1-b1x2=0 r2-a2x2-b2x1=0 的交点：,.,如果两条直线不相交，则只存在3个平衡点。在这种情况下，两
4、个种群不能共存。 我们希望知道在什么条件下x10且x20. 假设aibi,2个种群共存的条件：,.,平衡态,.,第五步，回答问题,对每种种群存在两类增长限制。第一种是由于与另一种群的竞争，第二种是由于拥挤造成的同一种群内部的竞争。因此，对每一种树，存在着一点，在这一点数木由于拥挤会主动停止增长，且存在另一点，在这一点树木通过竞争阻止另一种群的增长。两种树能够共存的条件是每种树在达到限制自己增长的点之前已经达到它限制另一种树增长的点。,.,4.2 动力系统,例 4.2 蓝鲸和长须鲸是生活在同一海域的相似种群，因此认为他们之间存在竞争。蓝鲸的内禀增长率每年估计为5%，长须鲸为每年8%，环境承载力（
5、环境能够支付的鲸鱼的最大数量）估计蓝鲸为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=aBF B=0,F=0, a是正实数 目标：确定动力系统
6、是否能够从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(x1,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/1
7、50000) - 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; % number of samples points x1min=0; x1max=900000; % domain specification x2min=0; x2max=600000; X1,X2=meshgrid(x1min:(x1max-x1min)/M:x1max,x2
8、min:(x2max-x2min)/M:x2max); dX1=0.05*X1.*(1-X1/150000) - alpha*X1.*X2; % x1-component dX2=0.08*X2.*(1-X2/400000) - alpha*X1.*X2; % x2-component quiver(X1,X2,dX1,dX2); % matlab routine axis(x1min x1max x2min x2max); title(Direction field (the vectors may be rescaled!); hold on xlabel(Blue Whales); yl
9、abel(Fin Whales); ezplot(f1,0 900000 0 600000), hold on ezplot(f2,0 900000 0 600000),.,.,四个平衡态解,三个为： （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*
10、x2; x1steady,x2steady = solve(f1,f2) x1alpha=x1steady(4); x2alpha=x2steady(4); pretty(x1alpha), pretty(x2alpha),.,解之得 x1=150000*(-1+8000000*a)/D x2=400000*(-1+1875000*a)/D 其中，D=-1+15000000000000*a2 alpha1 = solve(x1alpha); alpha2 = solve(x2alpha); format short e double(alpha1), double(alpha2) 可见，对任意
11、的a0,x20.,.,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);,.,.,a = linspace(0,9*10(-7); x1a=subs(x1alpha,alpha,a); x2a=subs(x2alpha,alpha,a); ind=find(x1a0 ,.,.,format bank Sx1alpha = diff(x1alpha, alpha)
12、*(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 宇航员在训练中要求用手动控制做对接演习。作为这个演习的一部分，要求保持一个正在运行的太空船与另一个正在运行的太空船的相对位置。手控制器提供了不同的加速度和
13、减速度，并且在太空船上有一个装置测量这两个飞船的接近速度。建议使用如下的策略进行飞船对接。,.,首先观察接近速度。如果为零，则不用再做任何事情。否则，记住这个接近速度，再看加速度控制器，控制加速度使得它与接近速度相反（即如果接近速度是正值，则放慢，如果是负的，则加快。），且正比于这个差值（即如果发现接近速度达到2倍时，我们将于2倍的速度刹车）。经过一段时间，再观察接近速度并重复上面的步骤。在什么环境下这个策略才是有效的？,.,第一步，提出问题,设vn表示在时间tn观测到的接近速度，tn为第n次观测的时间。 太空船接近速度的改变：vn =vn+1-vn 两次观测之间的时间间隔：tn=tn+1-t
14、n 时间区间被分成两部分： 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次观测前的等待时间(秒),.,假设：,tn+1=tn+tn=tn+cn+wn vn+1=vn+vn=vn+an-1cn+anwn an=-kvn cn0 wn=0 目标：确定是否有vn0.,.,第
15、二步，选择建模方法,离散时间动力系统理论 见教材，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）位于上述两直线的交点。 下面绘制向量场F=(-kwx1-kcx2, x1-x2),.,绘制向量场,clear all, close all, clc c = 5,w = 10,k = 0.1
16、; 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 field (the vectors are rescaled!); hold on xlabel(Current Speed), ylabel(Previous Speed);,.,.,系统演变,定义函数dockfun.m function rhs = dockfun(x,c,w,k); rhs = -k*w*x(1)-k*c*x(2); x(1)-x(2);,.,迭代实现,x=8;10; N = 12; fprintf(n Current speed Prev. speednn) f
