实验2--差分方程和数值微分

实验2差分方程和数值微分题目3某湖泊每天有4310m的河水流入，河水中污物浓度为30.02/gm，经渠道排水后湖泊容积保持4320010m不变，现测定湖泊中污物浓度为30.2/gm，建立差分方程计算湖泊中1年内逐月（按30天计算）下降的污物浓度，问要多长时间才能达到环保要求的浓度30.04/gm？为了把这个时间缩短为一年，应将河水中污物浓度降低到多少？【模型及其求解】记湖泊总水量0V，每天流入湖泊的河水量inV，湖泊最初污物浓度0x，目前河流的污物浓度r0p，假定外部环境没有任何其他的污染物进入湖泊，并且湖泊本身不会制造污染物，流入湖泊的河水的污物浓度不变，不考虑每天任意时刻河水与湖水的混合过程，另外，假定湖很大，河流入水口与湖泊出水口相距很远（认为当天从湖泊流出的水的浓度是不变的），同时，湖水的总量任意时刻都保持不变（不会出现先蓄积再排出的现象）。基于以上假定，则第k天湖水的污物浓度kx只和第k-1天湖水的污物浓度1kx及流过的河水污物浓度rp有关，可以得到如下关系式：00100ininrkkVVVpxxVV记00inVVaV，00inrVpbV，则：1kkxaxb代入题目中给定的数据，在MATLAB中计算。源程序为：%第二章第3题源程序clearallV0=200*10^4;%参数初始值设定Vin=10^4;K=0;pr0=0.02;x0=0.2;a=(V0-Vin)/V0;P=[00.010.020.040.060.080.10.140.160.20.24]%河水不同的污物浓度forj=1:11pr=P(j);b=Vin*pr/V0;x(1)=x0;x(2)=a*x(1)+b;fork=2:720x(k)=a*x(k-1)+b;%差分关系式C(j,k)=x(k);%湖水污物浓度变化矩阵endC(j,1)=x0;if(j==3)%题目给出的河水污染物浓度下的情况m=1:12;B(m+1)=x(30*m);B(1)=x0;endendfigure,plot(0:12,B,'b*',0:12,B,'r-');%一年内按月份污物浓度变化(pr=0.02)title('图1.一年内湖泊污染物浓度(按月)')xlabel('月份')ylabel('污物浓度(g/m^{3})')gtext('pr=0.02');figure,plot(1:720,C(1,:),'m',1:720,C(2,:),'b',...%pr取不同值污染浓度曲线比较1:720,C(3,:),'y',1:720,C(4,:),'r',1:720,C(5,:),'g',...1:720,C(6,:),'c',1:720,C(7,:),'k',1:720,C(8,:),'m',...1:720,C(9,:),'g',1:720,C(10,:),'r',1:720,C(11,:),'k'),...axis([072000.25]);title('图2.河水污物浓度不同时湖泊污染物浓度变化情况(按天)');xlabel('天')ylabel('污物浓度(g/m^{3})'),gtext('pr=0');gtext('pr=0.01');gtext('pr=0.02');gtext('pr=0.04');gtext('pr=0.06');gtext('pr=0.08');gtext('pr=0.1');gtext('pr=0.14');gtext('pr=0.16');gtext('pr=0.2');gtext('pr=0.21');a=0.9950;k=360;xk=0.04;x1=0.2;%如果一年内污染物浓度达要求求解prXb=(1-a)*(xk-x1*a^(k-1))/(1-a^(k-1));prX=b*V0/Vin;问题解决Step1：首先，得到一年内1～12月份的湖水污物浓度数据，如表1。月份123456污物浓度0.17560.15390.13520.11910.10530.0934月份789101112污物浓度0.08310.07430.06670.06020.05460.0498表1.一年内湖水污物浓度（月份）（浓度单位：g/m3）MATLAB绘制出曲线如图1：0246810120.040.060.080.10.120.140.160.180.2图1.一年内湖泊污染物浓度(按月)月份污物浓度(g/m3)pr=0.02由此曲线可以清楚地看到湖水的污染物浓度不断降低，但是一年时间仍然不能使得湖水污物浓度达到环保要求。利用Matlab进一步计算可以得到如表2所示结果：k12…439440…x(k)0.20.1991…0.04000.03990.0399…表2.河水污物浓度0.02g/m3时，湖水污物浓度达标时间（浓度单位：g/m3）湖泊污物浓度将在第439天达到30.04/gmStep2：由差分方程1kkxaxb逐步递推得到关系式：11221111(1)1kkkkkaxaxbaaaaxba从而得到b的表达式：111(1)()1kkkaxaxba按照题目要求，如果要求把湖水污染物浓度降低到环保标准（即30.04/gm）所用时间减少到一年，则可以通过上式求出。参数为00inVVaV=0.9950，k=360（1年），kx30.04/gm，1x30.2/gm代入式中，求得：54.147910b再由关系式：0inrVpbV计算得到：30.0083/rpgm，也就是说当河水的污染物浓度为30.0083/rpgm时，按照模型计算，可以在一年时间使得湖水达到环保要求。【结果分析与讨论】在建立本题目的模型过程中，首先作了一些假定，归结起来主要有：1．外部环境没有任何其他的污染物进入湖泊，并且湖泊本身不会制造污染物；2．流入湖泊的河水的污物浓度不变；3．不考虑每天任意时刻河水与湖水的混合过程；4．假定湖很大，河流入水口与湖泊出水口相距很远（认为当天从湖泊流出的水的浓度是不变的）；5．湖水的总量任意时刻都保持不变（不会出现先蓄积再排出的现象）。在这五条假定的前提下，建立了形如：1kkxaxb的差分方程数学模型，并对问题相应的求解。利用MATLAB画出了不同的河水污染物浓度pr的情况下，湖水污染物浓度随时间的变化曲线，如图2：010020030040050060070000.050.10.150.20.25图2.河水污物浓度不同时湖泊污染物浓度变化情况(按天)天污物浓度(g/m3)pr=0pr=0.01pr=0.02pr=0.04pr=0.06pr=0.08pr=0.1pr=0.14pr=0.16pr=0.2pr=0.21从图上曲线的变化趋势可以很清楚的看到，河水的污染物浓度是影响湖水质量的重要因素。在低于湖水的污染物浓度的情况下（0.2rp），河水的污染物浓度越低，湖水的污染物浓度下降得越快，湖水质量提高的越快；当达到0.2时，湖水的浓度会保持不变；超过0.2时，湖水的污染物浓度会越来越高。以上这三点实验结论都是符合客观事实的。在现实中，河水是不断流动的，而湖水是不流动的，所以湖泊的污染情况直接取决于河水的污染情况，因此该模型能够明确的给出一种湖泊类不流动水系污染治理的思路，即应将治理工作的重点集中在治理同湖泊相连的河流的治理。在实际应用中，该模型还应该考虑到湖泊原本积累的污染物会导致新的污染，从而加剧湖水的污染状况，因此要在差分方程中引入一个阻滞项，这样上面的曲线斜率就会降低，即湖泊污染程度改善速度会降低。题目5据报道，某种山猫在较好、中等及较差的自然环境下，年平均增长率分别为1.68%，0.55%和-4.50%，假定开始时有100只山猫，按以下情况讨论山猫数量逐年变化过程及趋势：（1）3种自然环境下25年的变化过程（作图）；（2）如果每年捕获3只，会发生什么情况？山猫会灭绝吗？如果每年捕获一只呢？（3）在较差的自然环境下，如果使山猫数量稳定在60只左右，每年要人工繁殖多少只？【第(1)问模型建立及求解】记第k年该山猫的数量为，自然环境下的年平均增长率为r，且a=1+r，则第k+1年山猫的数量为=a，a=1+r，k=0,1,2,…(1)在较好、中等及较差的自然环境下，以r=0.0168，0.0055，-0.0450以及=100，b=0代入，用MATLAB计算并作图，程序如下：%------------------------------作业题2_5函数M文件源程序---------------------------functionx=Exf2_5(x0,n,r,b)%建立名为Exf2_5的函数M文件，x0,n,r,b可调节a=1+r;x=x0;%赋初值fork=1:nx(k+1)=a*x(k)+b;%按照（1）式进行迭代计算end%------------------------------作业题2_5脚本M文件源程序---------------------------clearallx0=100;n=25;b=0;%x0:山猫初始数量；n：年数；b:人工繁殖数（若捕获则为负数）k=(0:n)';y1=Exf2_5(x0,n,0.0168,b);%给定x0,n,r,b,调用Exf2_5函数计算y2=Exf2_5(x0,n,0.0055,b);y3=Exf2_5(x0,n,-0.0450,b);round([k,y1',y2',y3']),%对结果四舍五入取整plot(k,y1,k,y2,':',k,y3,'--'),%将3条线画在一个图上xlabel('年数')ylabel('山猫数量')title('山猫数量在3种自然环境下25年的变化过程（b=0）')%加入X轴标记，Y轴标记和标题gtext('r=0.0168'),gtext('r=0.0055'),gtext('r=-0.0450'),%在图上标记各图线的增长率得到的结果如下面的表1和图1所示：表1：山猫数量在3种自然环境下25年的变化过程（b=0）k（r=1.68%，b=0）（r=0.55%，b=0）（r=-4.50%，b=0）01001001001102101962103101913105102874107102835109103796111103767112104728114104699116105661011810663111201066012122107581312410755141261085215128109501613110948171331104618135110441913711142201401124021142112382214411336231471133524149114332515211532图1：山猫数量在3种自然环境下25年的变化过程（b=0）【第(1)问结果分析及理论计算】由Matlab计算出来的数据表格和变化趋势图可见，无人工干预时，在较好（r=1.68%）和中等（r=0.55%）自然环境下，该山猫数量逐年增多（但在较好环境下增长速度更快）；而在较差（r=-4.50%）自然环境下，山猫数量逐年减少，趋于灭绝。理论分析如下：在自然环境下（1）式的解的形式为显然当a1（即r0时），∞，山猫数量逐年增多；而a1（即r0）时，，表明在较差的自然环境下山猫将濒于灭绝。可见，理论计算分析结果与Matlab计算结果是一致的。【第(2)问模型建立及求解】如果每年捕获3只和每年捕获1只，在三种自然环境下分别以b=-3和b=-1代入，用MATLAB计算得到的结果如下面的表2和图2、图3所示：表2：每年分别捕获3只和一只山猫时，其数量在3种自然环境下的变化过程k（r=1.68%，b=-3）（r=0.55%，b=-3）（r=-4.50%，b=-3）（r=1.68%，b=-1）（r=0.55%，b=-1）（r=-4.50%，b=-1）01001001001001001001999893101100952979585101998939693781029984495907210398795