数学建模第二章作业答案章绍辉习题2作业讲评1. 继续考虑2.2节的“汽车刹车距离”案例，请问“两秒准则”和“一车长度准则”一样吗？“两秒准则”是否足够安全？对于安全车距，你有没有更好的建议？（“两秒准则”，即后车司机从前车经过某一标志开始，默数2秒之后到达同一标志，而不管车速如何. 刹车距离与车速的经验公式20.750.082678d v v =+，速度单位为m/s ，距离单位为m ）解答（1）“两秒准则”表明前后车距与车速成正比例关系. 引入以下符号：D ～ 前后车距（m ）；v ～ 车速（m/s ）；于是“两秒准则”的数学模型为22D K v v ==. 与“一车长度准则”相比是否一样，依赖于一车长度的选取.比较20.750.082678d v v =+与2D v =，得：()0.082678 1.25d D v v -=-所以当15.12 m/s v <（约合54.43 km/h ）时，有d<D ，即前后车距大于刹车距离的理论值，可认为足够安全；当15.12 m/s v >时，有d>D ，即前后车距小于刹车距离的理论值，不够安全. 也就是说，“两秒准则”适用于车速不算很快的情况.另外，还可以通过绘图直观的解释“两秒准则”够不够安全. 用以下MATLAB 程序把刹车距离实测数据和“两秒准则”都画在同一幅图中（图1）.v=(20:5:80).*0.44704;d2=[18,25,36,47,64,82,105,132,162,196,237,283,334 22,31,45,58,80,103,131,165,202,245,295,353,41820,28,40.5,52.5,72,92.5,118,148.5,182,220.5,266,318,376]; d2=0.3048.*d2;k1=0.75; k2=0.082678; K2=2; d1=[v;v;v].*k1; d=d1+d2;plot([0,40],[0,K2*40],'k') hold onplot(0:40,polyval([k2,k1,0],0:40),':k') plot([v;v;v],d,'ok','MarkerSize',2)title('比较刹车距离实测数据、理论值和两秒准则') legend('两秒准则','刹车距离理论值',...'刹车距离的最小值、平均值和最大值',2) xlabel('车速v （m/s ）') ylabel('距离（m ）') hold off510152025303540020406080100120140160180比较刹车距离实测数据、理论值和两秒准则车速v （m/s ）距离（m ）两秒准则刹车距离理论值刹车距离的最小值、平均值和最大值图1hold onplot([10,35]*0.44704,2*[10,35]*0.44704,'k',... [35,60]*0.44704,3*[35,60]*0.44704,'k',... [60,75]*0.44704,4*[60,75]*0.44704,'k') title('t 秒准则，刹车距离的模型和数据') xlabel('车速v （m/s ）') ylabel('距离（m ）') hold off510152025303540020406080100120140160180车速v （m/s ）距离（m ）t 秒准则，刹车距离的模型和数据t 秒准则刹车距离理论值刹车距离的最小值、平均值和最大值图24. 继续考虑2.3节“生猪出售时机”案例，假设在第t 天的生猪出售的市场价格（元/公斤）为2()(0)p t p gt ht =-+ (1)其中h 为价格的平稳率，取h =0.0002. 其它模型假设和参数取值保持不变.（1） 试比较(1)式与(2.3.1)式，解释新的假设和原来的假设的区别与联系；（2）在新的假设下求解最佳出售时机和多赚的纯利润； （3）作灵敏度分析，分别考虑h 对最佳出售时机和多赚的纯利润的影响；（4）讨论模型关于价格假设的强健性. 解答一（用MATLAB 数值计算）（1）比较(1)式与(2.3.1)式，(1)式表明价格先降后升，(2.3.1)式假设价格匀速下降，(1)式更接近实际（图3）. 两个假设都满足(0)p g '=-，在最佳出售时机附近误差微小（图4）. 绘图的程序p=@(t)12-0.08*t+0.0002*t.^2; figure(1) n=400;plot([0,n],[12,12-0.08*n],'k:',... 0:.1:n,p(0:.1:n),'k') axis([0,400,0,20])title('模型假设(1)式与(2.3.1)式的比较')legend('p(0) - g t (1)式',... 'p(0) - g t + h t^2 (2.3.1)式') xlabel('t （天）')ylabel('p （元/公斤） ') figure(2) n=20;plot([0,n],[12,12-0.08*n],'k:',... 0:.1:n,p(0:.1:n),'k')title('模型假设(1)式与(2.3.1)式的比较')legend('p(0) - g t (1)式',... 'p(0) - g t + h t^2 (2.3.1)式') xlabel('t （天）'), ylabel('p （元/公斤） ')50100150200250300350400024********161820模型假设(1)式与(2.3.1)式的比较t （天）p （元/公斤）p(0) - g t (1)式p(0) - g t + h t 2 (2.3.1)式图3246810121416182010.410.610.81111.211.411.611.812模型假设(1)式与(2.3.1)式的比较t （天）p （元/公斤）p(0) - g t (1)式p(0) - g t + h t 2 (2.3.1)式图4（2）在(1)式和(2.3.1)式组成的假设下，多赚的纯利润为()()23()(0)(0)(0)Q t rp gw c t hw gr t hrt =--+-+保留h ，代入其他具体数值，得()32()900.08 1.6Q t ht h t t =+-+令()2()31800.16 1.60Q t ht h t '=+-+=解得生猪出售时机为()210.161800.1619.230h ht ---=-（舍去负根）多赚的纯利润为()321111900.08 1.6Q ht h t t =+-+.代入h =0.0002，得113.829t =天，110.798Q =元.或者用MATLAB 函数fminbnd 计算，脚本如下： C=@(t)3.2*t; w=@(t)90+t;p=@(t,h)12-0.08*t+h*t.^2;Q=@(t,h)p(t,h).*w(t)-C(t)-90*12; Qh=@(t)-Q(t,0.0002); t1=fminbnd(Qh,0,30) Q1=Q(t1,0.0002)为帮助理解，可用以下脚本绘制图5： figure(2) tp=0:250;plot(tp,Q(tp,0.0002),'k') title('纯利润Q') xlabel('t （天）') ylabel('Q （元） ')050100150200250-600-500-400-300-200-100100纯利润Qt （天）Q （元）图5（3）用以下MATLAB 脚本计算灵敏度(,)t tS t h h h ∆=∆和(,)Q QS Q h h h ∆=∆，将结果列表.结论：h 的微小变化对t 和Q 的影响都很小 Qh=@(t)-Q(t,0.0002*1.01); [tn,Qn]=fminbnd(Qh,0,30); (tn-t1)/t1/0.01 (-Qn-Q1)/Q1/0.01Qh=@(t)-Q(t,0.0002*1.05); [tn,Qn]=fminbnd(Qh,0,30); (tn-t1)/t1/0.05 (-Qn-Q1)/Q1/0.05Qh=@(t)-Q(t,0.0002*1.1); [tn,Qn]=fminbnd(Qh,0,30); (tn-t1)/t1/0.1 (-Qn-Q1)/Q1/0.1表3 数值计算最佳出售时机t 对h 的灵敏度h h +∆h h ∆(%)t t +∆t t ∆(%)(,)t t S t h h h∆=∆0.000202 1 13.886 0.41459 0.41459 0.00021 5 14.121 2.1176 0.42352 0.000221014.431 4.35360.43536表4 数值计算多赚的纯利润Q 对h 的灵敏度h h +∆h h ∆(%)Q Q +∆Q Q ∆(%)(,)Q Q S Q h h h∆=∆ 0.000202 1 10.838 0.36936 0.36936 0.00021 5 11.001 1.8802 0.37604 0.00022 1011.2143.84790.38479（4）市场价格是经常波动的，如果价格下跌，往往会止跌回稳，模型假设(1)式以二次函数来刻画价格止跌回升的变化趋势，如果考虑的时间段长达数月，(1)式比(2.3.1)式更接近实际（见图3），但是本问题的最佳出售时机不超过20天，(1)式与(2.3.1)式在最佳出售时机附近非常近似（见图4），(1)式导致的模型解答可以由(2.3.1)式导致的解答加上灵敏度分析所代替. 所以采用更为简单的(2.3.1)式作为假设更好.具体分析如下：由12()(,)g g t p t h -+∆=，得12(,)1g p t h g gt∆-=-， 代入h =0.0002，t =13.82852279，g =0.08，得0.034571gg∆=-. 由于(,)t g S t g t g∆∆≈，根据课本2.3节，代入(,) 5.5S t g =-，t =10，算得11.901t t +∆=，与t =13.829只相差两天.用于以上分析计算的MATLAB 脚本： dg_g=(12-p(ts,0.0002))/ts/0.08-1 10+dg_g*10*(-5.5)解答二（用MATLAB 的Symbolic Math Toolbox 的MuPAD 软件符号计算）（1）运行以下MuPAD 语句，绘得图6和图7：plot(plot::Function2d(12-0.08*t+0.0002*t^2,t=0..400), plot::Function2d(12-0.08*t,t=0..150, LineStyle=Dashed));plot(plot::Function2d(12-0.08*t+0.0002*t^2,t=0..20), plot::Function2d(12-0.08*t,t=0..20, LineStyle=Dashed),#O);(1)式表明价格先降后升，在实际当中有一定道理. 而 (2.3.1)式假设价格匀速下降. 两个假设都满足(0)p g '=-，在最佳出售时机附近误差微小.图6 假设(2.3.1)式与(1)式的比较图7 假设(2.3.1)式与(1)式的比较(2) 在(1)式和(2.3.1)式组成的假设下，保留h，代入其他具体数值，计算多赚的纯利润. 运行以下MuPAD语句：C:=t->32/10*t:w:=t->90+t:p:=(t,h)->12-8/100*t+h*t^2:Q:=(t,h)-->expand(w(t)*p(t,h)-C(t)-90*12); plot(plot::Function2d(Q(t,0.0002), t=0..290));算得223(2)825,905ht h h t Q t t t =+-+，绘得图8.图8 (,0.0002)Q t 的图像运行以下MuPAD 语句：S:=solve(diff(Q(t,h),t),t) assuming h>0; t1:=S[1];subs(t1,h=0.0002); t2:=S[2];ts:=subs(t2,h=0.0002); Q2:=Q(t2,h);Qs:=subs(Q2,h=0.0002);由方程0Qt∂=∂，解得两根： 2384165625123841656252253240045004450025324004h h h t h h h t -+-+=+-+-=代入h =0.0002，得12192.8381439, 13.82852279t t ==（天）. 2t 符合题意，1t 应该舍去（对应的Q 是负数）. 2t 对应的多赚的纯利润为10.79837809元.（3）接着上一小题，运行以下MuPAD 语句：subs(diff(t2,h)*h/t2, h=0.0002); //t 对h 的灵敏度利用导数算得t 对h 的灵敏度：d (,)0.4124276803d t hS t h h t=⋅=.运行以下MuPAD 语句：subs(diff(Q2,h)*h/Q2,h=0.0002); //Q 对h 的灵敏度，方法一 subs(diff(Q(t,h),h)*h/Q(t,h),t=ts,h=0.0002); //Q 对h 的灵敏度，方法二，更简单用两种方法利用导数算得Q 对h 的灵敏度：d (,)0.367739025d Q hS Q h h Q=⋅=. 结论：h 的微小变化对t 2和Q 2的影响都很小. （4）同解答一5. 继续考虑第2.3节“生猪出售时机”案例，假设在第t 天的生猪体重（公斤）为()000()mt m w w w t w w w e α-=+- (2)其中0(0)90w w ==（公斤），270m w =（公斤），其它模型假设和参数取值保持不变.（1）试比较(2)式与(2.3.2)式，解释新的假设和原来的假设的区别与联系（提示：说明当α (α>0)取何值时，在t =0时可以保持(0)1w r '==；说明当t 增大时，猪的体重会如何变化）.（2）在新的假设下求解最佳出售时机和多赚的纯利润. （3）参数m w 代表猪长成时的最终重量，对m w 做灵敏度分析，分别考虑m w 对最佳出售时机和多赚的纯利润的影响.（4）讨论模型关于生猪体重假设的强健性. 解答一（用MATLAB 数值计算）（1）在(2)式中，为使(0)w r '=，必须00()m m w w w w α-=. 当m w =270，0w =90时，有160α=.新假设(2)式是阻滞增长模型，假设生猪体重的增长率是体重的线性递减函数，于是体重增加的速率先快后慢，时间充分长后，体重趋于m w . 而(2.3.2)式0()w t w rt =+只假设体重匀速增加. 长时间来看，新假设比原假设更符合实际（图9）. 两个假设都满足(0)w r '=，在最佳出售时机附近误差微小（图10）.50100150200250300350400050100150200250300t （天）价格 p （元/公斤）模型假设(2.3.2)式与(2)式的比较p(0) - g t (2.3.2)式p(0) - g t + h 2 (2)式图924681012141618209095100105110115t （天）价格 p （元/公斤）p(0) - g t (2.3.2)式p(0) - g t + h 2 (2)式图10(2) 在(2.3.1)式和(2)式组成的假设下，用MATLAB 函数fminbnd 计算，可以求得生猪出售时机为t =14.434天，多赚的纯利润为Q =12.151元.(3) 编程计算(,)m m m t t S t w w w ∆=∆和(,)m m mQ QS Q w w w ∆=∆，将结果列表.表5 数值计算最佳出售时机t 对m w 的灵敏性m m w w +∆m mw w ∆(%)t t +∆t t ∆(%)(,)m m mt tS t w w w ∆=∆272.7 1 14.977 3.767 3.767 283.5 5 17.057 18.173 3.6345 2971019.46 34.8253.4825表6 数值计算多赚的纯利润Q 对m w 的灵敏性m m w w +∆m m w w ∆(%)Q Q +∆Q Q ∆(%)(,)m m mQ Q S Q w w w ∆=∆272.7 1 13.108 7.872 7.872 283.5 5 17.121 40.897 8.1794 2971022.47584.9638.4963结论：m w 的微小变化对t 和Q 的影响都较小.（4）模型假设(2)式导致的模型解答可以由(2.3.2)式导致的解答加上灵敏度分析所代替，所以实践中采用更为简单的(2.3.2)式作为假设即可. 具体分析过程见解答二之（4）.MATLAB脚本：%% (1) 绘图的程序w=@(t)90*270./(90+180*exp(-t/60));figure(1)n=400;plot([0,n],[90,90+n],'k:',...0:.1:n,w(0:.1:n),'k')axis([0,400,0,300])legend('p(0) - g t (2.3.2)式',... 'p(0) - g t + h^2 (2)式',4) title('模型假设(2.3.2)式与(2)式的比较') xlabel('t（天）')ylabel('价格 p（元/公斤） ')figure(2)n=20;plot([0,n],[90,90+n],'k:',...0:.1:n,w(0:.1:n),'k')legend('p(0) - g t (2.3.2)式',... 'p(0) - g t + h^2 (2)式',2) xlabel('t（天）')ylabel('价格 p（元/公斤） ')%% (2) 最佳出售时机和多赚的纯利润C=@(t)3.2*t;w=@(t,m)90*m./(90+(m-90)*exp(-t/60)); p=@(t)12-0.08*t;Q=@(t,m)p(t).*w(t,m)-C(t)-90*12;Qh=@(t)-Q(t,270);ts=fminbnd(Qh,0,30)Qs=Q(ts,270)%% (3) 灵敏度分析Qh=@(t)-Q(t,270*1.01);[tn,Qn]=fminbnd(Qh,0,30);(tn-ts)/ts/0.01(-Qn-Qs)/Qs/0.01Qh=@(t)-Q(t,270*1.05);[tn,Qn]=fminbnd(Qh,0,30);(tn-ts)/ts/0.05(-Qn-Qs)/Qs/0.05Qh=@(t)-Q(t,270*1.1);[tn,Qn]=fminbnd(Qh,0,30);(tn-ts)/ts/0.1(-Qn-Qs)/Qs/0.1%% (4) 强健性分析dr_r=(w(ts,270)-90)/ts-110+dr_r*10*6.5解答二（用MATLAB 的Symbolic Math Toolbox 的MuPAD 软件符号计算）（1）运行以下MuPAD 语句，算得160α=：solve(subs(diff(90*270/(90+(270-90)*E^(-a*t)),t), t=0)=1, a);运行以下MuPAD 语句，绘得图11：plot(plot::Function2d(90*270/(90+180*E^(-1/60*t)), t=0..400),plot::Function2d(90+t, t=0..180, LineStyle=Dashed), plot::Line2d([0,270],[400,270],LineStyle=Dotted),#O);运行以下MuPAD 语句，绘得图12 ：plot(plot::Function2d(90*270/(90+180*E^(-1/60*t)), t=0..20),plot::Function2d(90+t,t=0..20,LineStyle=Dashed),#O);(2)式()06000()mt m w w w t w w w e -=+-是阻滞增长模型，假设生猪体重的增长率是体重的线性递减函数. 于是，体重w 是时间t 的增函数，体重增加的速率先快后慢，时间充分长后，体重趋于m w . 而(2.3.2)式0()w t w rt =+只假设体重匀速增加. 长时间来看，新假设比原假设更符合实际. 两假设都满足(0)w r '=，在最佳出售时机附近误差微小.图11 假设(2.3.2)式与(2)式的比较图12 假设(2.3.2)式与(2)式的比较w，代入其（2）在由(2)式和(2.3.1)式组成的假设下，保留m他具体数值，计算多赚的纯利润. 运行以下MuPAD语句：C:=t->3.2*t:w:=(t,wm)->90*wm/(90+(wm-90)*E^(-t/60)):p:=t->12-0.08*t:Q:=(t,wm)-->w(t,wm)*p(t)-C(t)-90*12;plot(plot::Function2d(Q(t,270),t=0..30));算得()()6090120.08(,) 3.210809090emm tmw tQ t w tw--=--+-，绘得图13.图13 (,270)Q t的图像运行以下MuPAD语句：T:=solve(diff(Q(t,270),t),t);ts:=T[1];Qs:=Q(ts,270);可解出Q的驻点的数值解14.43357158st=（天），根据函数图像和问题的实际意义，可知这是所求的最佳出售时机，对应的多赚的纯利润为12.15129217s Q =元.（3）接着上一小题，运行以下MuPAD 语句，但是求不出当(,)m Q t w 达到最大值时t 关于m w 的函数解析式：solve(diff(Q(t,wm),t),t);运行以下MuPAD 语句：solve(diff(Q(t,wm),t),wm);可见当(,)m Q t w 达到最大值时m w 关于t 的反函数解析式却有可能求得出，只是MuPAD 给出的表达式很复杂. 其实可以按如下步骤推出m w 关于t 的反函数解析式：g1:=diff(Q(t,wm),t)=0; 算得0Q t∂=∂即： ()()260606030.0812907.2 3.209090902e 90e e m m m m t m t t w t w w w w -----=--⎛⎫++ ⎪⎝⎭观察上式，发现分母大于零，而且去分母之后，合并m w 的同类项，可以表示为m w 的二次方程：g2:=g1*((wm-90)/E^(t/60)+90)^2*25*E^(t/60); //去分母 g2:=collect(g2,wm); //合并wm 的同类项，t 当作参数2606060306060801440016200e 270327038700e e e 648000e 64800012960000e e t m m t t t t t t t w t w ⎛⎫⎛⎫--++-- ⎪ ⎪⎝⎭⎝⎭+--=运行以下MuPAD 语句，由图像（图14）可知在实际问题关心的0<t <30范围内，二次项系数608027030et t -->： plot(plot::Function2d((270-80/E^(t/60)-3*t),t=0..100));图4 二次项系数的符号于是，运行以下MuPAD 语句，解方程：S:=solve(g2,wm);MuPAD 给出解的四种情况，其中第一种是二次项系数非零，正是本问题所要求的解. 但是二次方程有两个根，要检验哪一个根才是当(,)m Q t w 达到最大值时m w 关于t 的反函数解析式. float(subs(S[1][1],t=ts));算得当s t t =时，有0.8519704108m w =-，这是增根，舍去； float(subs(S[1][2],t=ts));算得当s t t =时，有270m w =，这是要找的根；wms:=S[1][2]; //当Q 达到最大值时wm 关于t 的反函数解析式 float(subs(1/(diff(wms,t))*wm/t,t=ts,wm=270));//t 对wm 的灵敏度，利用反函数求导数利用反函数求导数算得t 对m w 的灵敏度：d 1(,) 3.80183985d d d m m m m m w w t S t w w w tt t=⋅=⋅=. Q 对m w 的灵敏度则比较简单，运行以下MuPAD 语句： float(subs(diff(Q(t,wm),wm)*wm/Q(t,wm),t=ts,wm=270)); //Q 对wm 的灵敏度利用导数算得Q 对m w 的灵敏度：d (,)7.786585188d m m m w Q S Q w w Q=⋅=. 结论：m w 的微小变化对t 和Q 存在一定影响，不算厉害.（4）模型假设(2)式以阻滞增长模型来刻画生猪体重的变化趋势，如果考虑的时间段长达数月，(2)式比(2.3.2)式更符合实际，但是本问题的最佳出售时机不超过20天，(2)式与(2.3.2)式在最佳出售时机附近非常近似，(2)式导致的模型解答可以由(2.3.2)式导致的解答加上灵敏度分析所代替. 所以采用更为简单的(2.3.2)式作为假设更好. 具体分析如下：由()90(,)m r r t w t w ++∆=，得(,)90m w t w r r t-∆=-， 代入270m w =，14.43357158s t t ==，r =1，得0.036565352791r r r ∆∆==. 由于(,)t r S t r t r∆∆≈，根据2.3节，代入(,) 6.5S t r =，t =10，r =1，算得12.37674793t t +∆=，与14.43357158s t =只相差两天.以上计算可以用以下MuPAD 语句实现：dr:=float((w(ts,270)-90)/ts-1);10+dr*10*6.5;。