湖水污染问题的数学建模与求解

1、中国传媒大学2010学年第 一 学期数学建模与数学实验课程数学建模与数学实验题目Pristine湖污染问题的建模与求解学生姓名学号班级学生所属学院任课教师教师所属学院成绩Pristine湖污染问题的建模与求解摘要本文讨论了湖水污染浓度变化趋势的预测问题。通过分析水流输入输出湖泊的过程，建立了湖水污染浓度随时间变化的含参变量的微分方程模型，在河水污染浓度恒定和自然净化速率呈线性关系的情况下，求得其精确解，带入具体数据得到结论：在 PCA声称的河水污染浓度下，湖的环境不会恶化；在工作人员实地测得的河水浓度下，湖的环境将会恶化。同时建立了计算机模拟模型，带入具体数值，运用时间步长法来仿真模拟了在湖水

2、污染浓度稳定以前湖水每天的变化情况，输出自 PCA建厂以来每年的湖水污染浓度，得到与微分方程模型相同的结论。在全停产和半停产时，通过前面的两个模型可以计算湖水污染浓度在自然净化影响下的恢复到净化指标所需的年限。并可得到结论：在半停产状态下，在选定的自然净化速率常数的约束下，只有当河水污染浓度降至原来的 3.15%（自然净化速率呈线性关系） ，4.7%（自然净化速率呈指数关系） ，才有可能使河水在 100年内恢复至 0.001mol/l ，然后给出整改建议。一、问题重述Pure河是流入 Pristine湖的唯一河流。 50 年前 PCA公司在此河旁建起一个生产设施并投入运行。 PCA将为处理的湖

3、水排入河中，导致Pristine湖被污染。 PCA公司声称：已排放的废水的标准多年从未改变切不会对湖的环境有影响。现已知： Pristine湖的湖容量为 1015 L，流入（流出）的水流速度为1.914 L/ 年。 PCA公司声称河水污染浓度仅为0.001mol/L ，自工厂以来没有改变过。讨论下列问题：（1）建立数学模型用PCA提供的公开数据判断湖的环境是否会恶化；（2）以目前湖水污染浓度 0.03mol/L ，和河水污染浓度 0.05mol/L 为新数据判断湖的环境是否会恶化；二、模型的合理假设和符号系统2.1模型的合理假设（1）降水量和增发量相等；（2）湖中流入量和流出量相等且一直未变；

4、（3）污水量远小于河水注入量，且污水与河水混合均匀；（4）湖水混合均匀，且流入污水的扩散速度无限大；（5）湖内除 Pure 河外，无其他污染源；2.2符号系统0 ：河水污染浓度 mol/L ；：湖水污染物浓度 mol/L ；V ：湖泊容量 1015 L；c ：自然净化速率 mol/ （ L。年）：流入（流出）的水流速度1.914 L/ 年；t ：从 PCA建厂至考察时刻的时间段。三、问题的分析3.1问题分析：对于问题中几个词语的理解：1.是否会恶化湖的环境恶化即湖水污染浓度大于 0.001mol/L ，要判断其是否会恶化，则需计算在某一污染物积累速度（分析影响此速度的因素）下，湖水能达到的最大

5、污染浓度和其变化趋势，以及湖水经几年超过0.001mol/L，经过几年达到最大污染浓度。2. 自然净化自然净化是独立的生态系统进行自我调节的方式之一，是在空气，阳光，水和细菌的参与下，进行包括物理沉降，化学反应和生物转化三大方面的活动，其最终作用是将污染物转化为无害物质，从而净化生态系统。当河水输入湖泊并均匀混合之后，影响湖水污染浓度的唯一因素便是自然净化速度。湖水污染问题水流的动态流程图：此问题中，我们考察对象是湖水污染浓度的变化趋势：1.在整改之前，其增加的趋势，超过净化指标 0.001ml/L （即湖水恶化）的可能性和时限；2.在全停产（无污染物输入）和半停产的情况下，其降低的趋势，达到

6、净化指标的时限。3.在整改之后，其增加的趋势，未定与净化指标之下某一水平的时限、在前假设条件的基础之上，湖水容量不变，出河水外无其他的污染源，故我们可将湖泊作为一个封闭的生态系统，其简化的湖水被污染的动态过程为：受污河输入湖泊，河水与湖水均匀混合，受污河水进行自我的净化，湖水数出湖泊。湖中污染物的量直接决定了湖水污染浓度，而污染物的量受到以下两方面因素的影响： 1.河水的污染浓度及其流入速度（根据已知此速度不变） ，2.湖水的自然净化速度，前者使其增加，后者使其减少（负增加）。问题一、二、三的实质都是要分析污染浓度的变化趋势，其去表便在于前一因素的不同。问题一中，河水污染浓度不变，恒为 0.0

7、01mol/L ；问题二中，河水污染浓度可能会变化，受 PCA 效益的影响而按一定规律波动； 问题三中，在全停产或半停产的情况下， 和硕污染浓度为 0 或减为问题二中的一部分。后一因素（自然净化速度）在三个问题中的作用都是相同的。根据微积分的知识可知，在适当短的时间段之内，通过建立微分方程，可以将连续的过程离散化，从而可以得到湖水污染浓度与时间之间的关系式。利用时间步长发，缩小步长值（从年到月到天），并与微分方程所得的精确解做出比较。四、模型建立与求解问题一：根据 PCA 的公开申明和所提供数据，可认为：河水污染浓度恒为0.001mol/L。从存在自然净化和不存在自然净化两个方面考虑：（1）在

8、不考虑自然净化的情况下：由于假设湖中流入量和流出量相等，而在经过与湖水均匀混合后，流出湖水污染浓度明显减小，故流出污染物的量小于流出污染物的量，污染物将在湖中沉积，从而使湖水污染浓度增加，当其增加至于输入的河水污染浓度相等时，河水污染浓度达到最大，并稳定在这一数值，不在增加。建立湖水污染浓度随时间变化的微分方程模型：设在极短时间 dt ，湖水污染浓度增加d，在将湖水被污染这一连续动态过程简化为离散的瞬间静止状态（如问题分析中所述）之后，根据湖中剩余量 =输入量输出量，我们可以列微分方程如下：d( 0)化简可得： dtV······（1

9、）带入 V1015L，1.91014L/ 年和00.001mol/L 的数据，我们可得到和 t 的关系式如下：0.001 e0.19 t(e0.19 t1)通过此关系式我们可知，当t时，湖水污染浓度将趋近与0.001mol/L ，即湖的环境不会恶化。建立计算机模拟模型：湖水污染浓度的变化时有湖中污染物随时间的积累而引起的，这个逐步积累的过程我们可以用计算机进行仿真模拟，其实质为完成一个循环累加的过程，并可改变时间步长，如一年一年的累积，一月一月的累积，一天一天的累积，从而使我们的模拟值逐步精确，可与微分方程求得精确解比较，分析误差。为提高模拟结果的精确性和运算的效率，我们采用了逐天累加，数出年

10、污染浓度的方式。模拟程序见附件一，从后面的分析中我们可知：在此河水污染浓度恒定，无自然净化的最简单的情况下建立的模型是以后问题的基础，后面的问题只是改变条件或数据，其实质是不变的，故我们在此程序中加入了多个选择语句， 在不同德条件或数据下执行不同的命令，从而用一个程序解决全部的模拟问题。模拟所得的数据如下表1：年份1358101.734.356.147.828.51年439.789.979.9959.999根据模拟数据所作的湖水污染浓度变化趋势图如下图1：模拟所得的数据显示： 湖水污染浓度将稳定于 0.001mol/L ，这与分析微分方程所得的结果是相符的。（2）在考

11、虑自然净化的情况下：A. 首先，我们应该了解什么事自然净化。根据资料显示：湖水中的污染物可分为有机污染物和无机污染物两大类，在多种环境因素（阳光，空气，水，水中生物，水中化学物质，重力等）的作用下，通过物理沉降，化学反应和生物转化一系列复杂的活动，它们的量会发生变化。有机物在水环境中的迁移转化过程如图 2：可以看出整个过程是相当复杂的，不仅过程多，而且在相同的过程中，不同的物质有着不同的结果。此问题中没有明确给出输入湖中的污染物种类， 也没有对湖泊环境做任何描述，无疑给问题的解决增加了极大的难度。为是问题简化，我们从一般的情况出发，假设：湖中所进行的反应均为一级反应，有机物的存在不会对环境参数

12、造成改变，环境固定（自然净化速率恒定），考虑湖泊生态系统中起主要作用的几种过程作简要分析：1. 物理沉降。不同物质有着不同的沉降速度常数，沉降速率为 c；2. 挥发。不同物质有着不同的挥发速率常数Kv ，在有机物在水体上的大气中的分压为 0cK v的条件下，挥发速率为Z( Z 为水体深度)3. 水解反应。不同物质有着不同的水解速率常数 K p ，在一级反应的条件下，水解速率为c K p；4. 生物降解反应。不同物质有着不同的降解速率常数Kb ，在一级反应的条件下，生物降解速率 cKb；B. 在最简模拟程序的基础上，从每天的积累量中减去每天的自然净化量K时间步长，重复循环可得到逐年湖水污染浓度的

13、值，设 K0.0003天 1。 K 越大，湖水污染浓度将越快稳定于一个更小的浓度值。模拟所得的数据如下表2：年份1358101.5873.4374.3154.8784.984年0875.1035.1065.1065.106根据模拟数据所作的湖水污染浓度变化趋势图如下图3：C.进一步考虑自然净化速率的影响。更贴近于实际的情况是，自然净化速率c 与湖水污染浓度成指数关系， c 随着的增加而增加，但增加的速率会逐渐减小，用关系式表达即可为： cA(1eB) ， A,B 为与环境和污染物种类有关的常数。将此关系式带入微分方程（2），得到一个新的微分方程，此方程无解。但是我们可

14、以通过计算机模拟求得数值解。A 的大小影响着最终稳定浓度， B 的大小影响着达到最终稳定浓度的快慢。在前面达到线 性关系的基础上， 在=0/ 。001mol/L 处的 c 值大小确定 A 至数量级，并稳定年限定在 1020 年间，从而得到一组估计值：模拟所得数据如下表 3：年份1358101.4932.9713.5393.8163.871年9033.9063.9063.9063.906根据模拟所得数据所作的湖水污染浓度变化趋势图如下图4：观察数据可知：无论自然净化速率 c与污染浓度成线性关系还是指数关系，当河水污染浓度0 =0.001mol/L 时，都将稳定于一个小于

15、0.001mol/L 的值，也即：湖的环境不会恶化。问题二：根据实际情况，在此我们只考虑存在自然净化的情况。（1）. 河水污染浓度恒定为0.05mol/L:与问题一（ 2）相同，只是0 =0.05mol/L ，沿用问题一（ 2）的微分方程模型和差分模拟模型，我们可以得到以下结果：A. 微分方程模型（自然净化速率c 与湖水污染浓度成线性关系）：0当时间 t时，KV ，即湖水污染浓度稳定于一个与 K 有关的值。当 0 =0.05mol/L时，如要稳定于 0.001mol/L 则 K 值为 9.31 年 1。根据我们设定的 K 0.0003天 1可计算，湖水最终污染浓度0 =0.0317mol/L

16、，超过净化指标 0.001mol/L ，故在此条件下， 湖的环境将会恶化。B. 差分模拟模型：在 c 与成线性关系时，得到模拟数据如下表4：年份1358100.821.882.462.883.01年133.163.173.173.17在 c 与成指数关系时，得到模拟数据如下表5：年份1358100.812.012.813.594.02年324.484.564.584.58观察以上数据可知： 由于河水污染没得过大， 河水污染没得从第一年起就超过了净化指标，并逐年增加，是湖的环境恶化。（2）. 河水污染浓度变化：根据实际情况，一个工厂的生产量并非

17、恒定不变，每年每月甚至每天也有所不同，从而起其排放污染物的量也将随生产量的变化而变化的。假设PCA 自建立以来的年排污量服从Logistic模型，并考虑到污染物的量受到多方面因素的影响，在每一时刻的梁上加上一个服从正太分布，范围在此量的10%内的较小量，则可建立河水污染浓度的模型如下：M0（M,C,r 为与工厂效益有关的常数）1 Cer M t设此时工厂处于文本上升的发展阶段，其变化曲线如下图5：通过计算机模拟产生随机数（如图中星点所示），带入模拟程序。在 c 于成线性关系时，得到的模拟数据如下表6：年份1510202235400.1850.8011.7927.2729.58152.0497.

18、13年份508290100112140150260.7625.0633.9617.2626.3628.5632.5根据模拟所得数据所作的湖水污染浓度变化趋势如下图5：在 c 于成指数关系时，得到的模拟数据如下表7：年份1510202435450.1690.5861.2505.1949.04644.2695.75年份5080100120130140150321.7925.5959.6950.6949.6952.5964.8根据模拟所得数据所作的湖水污染浓度变化趋势如下图5：在第 50 年时，无论自然净化速率c于湖水污染浓度成线性关系还是指数关系，可以看到湖水污染浓度都接近工作人员实地测得值，这从

19、检验的角度说明我们的对模型参数的估计是可取的。观察模拟所得数据可知，湖的环境将会恶化。五、模型的评价及评改进5.1 模型的评价（1）. 对于首先建立的含参数的微分方程模型，在其有精确解的时候，能给出湖水污染浓度随时间变化的规律，得到很好的预测效果，但微分方程并不是任何条件下都有解，在本问题中，只有当河水污染浓度恒定且自然净化速率呈线性关系时，微分方程可解，故适用的范围受到了限制。（2）. 差分模拟法，即用数值方法求解微分方程，它解决了微分方程无精确解的情况，适用范围很大。使用一个兼顾运算效率和模拟结果精度的时间步长是该方法的关键。在处理本问题时，采用的时间步长为天。这样处理的精度大约为 0.0

20、001, 与精确解相当接近。5.2模型的改进（1）. 自然净化速率的改进：资料显示，污染物的很多生物降解过程是通过好氧反应实现的，当污染浓度较大时，需氧量增加，而湖中氧气总量不变，自然净化速率会有所减小，所以，自然净化速率的函数可作为分段函数考虑。前一段还是用前面的指数模型，后面一段可以考虑为Ce DE (C, D, E0)。模型这样改进后，更加符合实际情况。（2）. 河水浓度变化规律模型的改进：根据价值规律，一个工厂的发展总是有起有落，其生产量绕某一中心值上下波动，故生产河水的浓度的变化应该是具有一定周期性，同时又有所增加。可以考虑在原来的 Logistic 模型后加上一个傅立叶函数（其周期

21、就是工厂生产的周期）和一个正态随机变量。六、参考文献1 数学建模与数学实验 高等教育出版社，赵静、但琦主编， 2008 年 1 月第三版2 环境化学，高等教育出版社出版，戴树桂主编， 1997 年3 月；七、附录附录一 :问题一：（1）河水浓度恒为 0.001mol/L ，不考虑自然净化， a ： 采用差分方程模型命令： wenti12 （d， 0.50,0.001 ，n）b：采用微分方程模型命令： t=1:1:51；Pplot(t,p)grid on（2）河水浓度恒为 0.001mol/L ，考虑自然净化，a：采用差分自然净化线性自然净化命令： wenti12('d',0,5

22、0,0.001,'y','xx')指数自然净化命令： wenti12('d',0,50,0.001,'y','zs')b：采用微分方程模型，线性自然净化命令： t=1:1:51;pholdonplot(t,p)gridon问题二：（1）河水浓度恒为 0.05mol/L, 不考虑自然净化a：采用差分方程模型命令： wenti12('d',0.50,0.05,'n')b: 采用微分方程模型命令：hold onplot(t,p)gridon（2）河水浓度恒为0.05mol/L ，考虑自然净

23、化，a: 采用差分方程模型线性自然净化命令： wenti12('d',0.50,0.05,'y','xx')指数自然净化命令： wneti12('d',0.50,0.05,'y','zs')b：采用微分方程模型，线性自然净化命令t=1:1:51;pholdonplot(t,p)gridon（3）河水浓度变化，不考虑自然净化命令： wenti12('d',0.50,-1,'n')考虑自然净化线性自然净化命令： wenti12('d',0,50,-1,&#

24、39;y','xx')指数自然净化命令： wenti12('d',0,50,-1,'y','zs')附录二：主程序：function wenti12(l,chu,n,s,w,b)%按照时间步长法，以每年或每月或每日为时间段，求出浓度的变化，最后输出每年的浓度，做一个大致的观察 ,%格式： wenti12 （P1,P2,P3,P4,P5,P6)%P1:y,m,d,h,min,s表示设定步长值%P2：湖水污染初值%P3：年份%P4：河流污染初值%P5: 是否考虑自然降解 ,n- 不考虑 ,y- 考虑%P6: xx 采用线性净化

25、模型， zs 采用指数净化模型 format compactformat longif nargin 2 %50 年, 河水污染变化 , 不考虑自然降解 n=50;s=-1;w='n'b='xx'elseif nargin 3 % 河水污染变化 , 不考虑自然降解 s=-1;w='n'b='xx'endif l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365elseif l 'h'q=365*24elseif l 'min

26、9;q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,wrong parameter.try again!')endA=zeros(1,n+1);B=zeros(1,n+1);A(1)=chu;B(1)=chu;var=0; % 自然降解var1=0; % 河水污染if w 'n'var=0;elseif w 'y'var=fen(A(1),l,b);elseerror('error,wrong parameter.try again!')endfor i=

27、1:1:nif s -1var1=river(i);elsevar1=s;endfor j=1:1:qif w 'n'var=0;elseif w 'y'var=fen(A(i),l,b);endA(i)=(A(i)*1015+var1*1.9*1014/q)/(1015+1.9*1014/q)-var;B(i)=0.19*var1/q+(1-0.19/q)*B(i)-var;endif w 'n'var=0;elseif w 'y'var=fen(A(i),l,b);endA(i+1)=(A(i)*1015+var1*1.9*1

28、014/q)/(1015+1.9*1014/q)-var;B(i+1)=0.19*var1/q+(1-0.19/q)*B(i)-var;endA;Bhold onplot(1:1:n+1,A,'o')grid onfunction p=wenti3(l,bi,year,b,a)%格式 : wenti3(P1,P2,P3,P4，P5)%P1:y,m,d,h,min,s表示设定步长值%P2:百分比%P3:可认为的最大净化年限%P4:xx-线性自然净化 ,zs- 指数自然净化%P5：湖水浓度初始值format compactformat longif l 'y'q=1

29、elseif l 'm'q=12elseif l 'd'q=365;elseif l 'h'q=365*24elseif l 'min'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,try again!')endA=zeros(100);B=zeros(100);temp=0.05*bi;A(1)=a;for i=1:1:yearfor j=1:1:qA(i)=(A(i)*1015+temp*1.9*1014/q)/(1015+1.9*1

30、014/q)-fen(A(i),l,b); t=A(i);if A(i)<0.001 break;endendA(i+1)=(A(i)*1015+temp*1.9*1014/q)/(1015+1.9*1014/q)-fen(A(i),l,b);if A(i)<0.001 break;endendif i yearp=0;elsep=1;end%程序调试所用%sprintf('停产 ')sprintf('需要 %d 年%d天 ',i,j)%sprintf('第%d年%d天为 :%.12f',i,j-1,t)%sprintf('第

31、%d年%d天为 :%.12f',i,j,A(i)function low=wenti3_2(b,year,l,h)%采用二分法求出百分比format compactformat longlow=l;high=h;cen=(low+high)/2cha=high-low;while cha>10-6%误差判断条件m=wenti3('d',cen,year,b,0.05); % 调用函数 if m 0high=cen;cen=(low+high)/2elselow=cen;cen=(low+high)/2endcha=high-low;endhigh;low;spri

32、ntf('工厂在生产为原来的 %f时，湖水才能在 %d年内净化。 ',low,year)function wenti3_3(l,year,bi,b)%格式 : wenti3_3(P1,P2,P3,P4)%P1:y,m,d,h,min,s表示设定步长值%P2:year%P3: 百分比%P4:xx-线性自然净化 ,zs- 指数自然净化format compactformat longif l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365;elseif l 'h'q=365*24else

33、if l 'min'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,try again!')endA=zeros(1,year);A(1)=0.03;N=0;for i=1:1:yeartemp=0.05*bi;for j=1:1:fix(q)A(i)=(A(i)*1015+temp*1.9*1014/q)/(1015+1.9*1014/q)-fen(A(i),l,b); endA(i+1)=(A(i)*1015+temp*1.9*1014/q)/(1015+1.9*1014/q)-fen(A(i),l,b );endN=A(i)被调用的程序：function y=fen(N,l,b)%自然净化%格式 :fen(P1,P2)%P1:前一时间点的污染值%P2:y,m,d, 按年或月或日%线性关系if b 'xx'first=0.1095;%first=0.931;if l 'y'y=first*N;elseif l 'm'y=first*N/12;elseif l 'd'y=first*N/365;elseerror('error,wrong parameter.try again!')e