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

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

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

3、的湖水排入河中，导致Pristine 湖被污染。PCA 公司声称：已排放的废水的标准多年从未改变切不会对湖的环境有影响。1514现已知：Pristine湖的湖容量为10 L,流入（流出）的水流速度为9 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 ；15V ：湖泊容量10 L ；c ：自然净化速率mol/ （ L 。年）14：流入（流出）的水流速度 1.9 L/ 年；t ：从 PCA 建厂至考察时刻的时间段。三、问题的分析3.1 问题分析：对于问题中几个词语的理解：1 .是否会恶化湖的环境恶化即湖水污染浓度大于0.001mol/L ，要判断其是否会恶化，则需计算在某一污染物积累速度（分析影响此速度的因素）下，湖水能达到的

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

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

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

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

9、方程如下：dt化简可得：dr V ( 0)(1)带入V 1015L,1.9 1014L/年和0 0.001 mol/L的数据，我们可得到和t的0.19 t 0.19 t关系式如下:0.001 e (e 1)通过此关系式我们可知，当t 时，湖水污染浓度将趋近与 0.001mol/L,即湖的环境不会恶化。建立计算机模拟模型：湖水污染浓度的变化时有湖中污染物随时间的积累而引起 的，这个逐步积累的过程我们可以用计算机进行仿真模拟， 其实质为完成一个循 环累加的过程，并可改变时间步长，如一年一年的累积，一月一月的累积，一天 一天的累积，从而使我们的模拟值逐步精确，可与微分方程求得精确解比较，分 析误差。

10、为提高模拟结果的精确性和运算的效率， 我们采用了逐天累加，数出年 污染浓度的方式。模拟程序见附件一，从后面的分析中我们可知：在此河水污染浓度恒定，无自然 净化的最简单的情况下建立的模型是以后问题的基础，后面的问题只是改变条件 或数据，其实质是不变的，故我们在此程序中加入了多个选择语句， 在不同德条 件或数据下执行不同的命令，从而用一个程序解决全部的模拟问题。模拟所得的数据如下表1：年份13581010 41.734.356.147.828.51年份152030405010 49.439.789.979.9959.999根据模拟数据所作的湖水污染浓度变化趋势图如下图1:模拟所得的数据显示：湖水污

11、染 浓度将稳定于0.001mol/L,这与分析微分方程所得的结果是相符的。(2)在考虑自然净化的情况下：A.首先，我们应该了解什么事自然净化。根据资料显示：湖水中的污染物可分为有机污染物和无机污染物两大类，在多种环境因素(阳光，空气，水，水中生物，水中化学物质，重力等)的作用下，通过物理沉降，化学反应和生物转化一系列复杂的活动，它们的量会发生变化。有机物在水环境中的迁移转化过程如图2:化还at作用可以看出整个过程是相当复杂的， 不仅过程多，而且在相同的过程中，不同的物 质有着不同的结果。此问题中没有明确给出输入湖中的污染物种类， 也没有对湖 泊环境做任何描述，无疑给问题的解决增加了极大的难度。

12、 为是问题简化，我们从一般的情况出发，假设：湖中所进行的反应均为一级反应，有机物的存在不会 对环境参数造成改变，环境固定（自然净化速率恒定），考虑湖泊生态系统中起 主要作用的几种过程作简要分析：1 .物理沉降。不同物质有着不同的沉降速度常数，沉降速率为c ;K 2 .挥发。不同物质有着不同的挥发速率常数v ,在有机物在水体上的大气中的c Kv_分压为0的条件下，挥发速率为C Z （Z为水体深度） 3水解反应。不同物质有着不同的水解速率常数 Kp，在一级反应的条件下，水解速率为c Kp ;4.生物降解反应。不同物质有着不同的降解速率常数 Kb,在一级反应的条件下,生物降解速率c Kb ;B.在最

13、简模拟程序的基础上，从每天的积累量中减去每天的自然净化量K时间1步长，重复循环可得到逐年湖水污染浓度的值，设 K 0.0003天。K越大，湖水污染浓度将越快稳定于一个更小的浓度值。模拟所得的数据如下表2:年份13581010 41.5873.4374.3154.8784.984年份152030405010 45.0875.1035.1065.1065.106根据模拟数据所作的湖水污染浓度变化趋势图如下图3:C.进一步考虑自然净化速率的影响。更贴近于实际的情况是，自然净化速率 c与湖水污染浓度 成指数关系，c随着 的增加而增加，但增加的速率会逐渐减小，B用关系式表达即可为：c A（1 e入A,B

14、为与环境和污染物种类有关的常数。将此关系式带入微分方程（2）,得到一个新的微分方程，此方程无解。但是我们可以通过计算机模拟求得数值解。A的大小影响着 最终稳定浓度，B的大小影响着 达到最终稳定浓 度的快慢。在前面达到线 性关系的基础上，在 =0/。001mol/L处的c值大小确定A至数量级，并稳定年限定在1A20年间，从而得到一组估计值:A 8 10 4,B 400模拟所得数据如下表3:年份13581010 41.4932.9713.5393.8163.871年份152030405010 43.9033.9063.9063.9063.906根据模拟所得数据所作的湖水污染浓度变化趋势图如下图4:

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

16、1年1。根据我们设定 的K 0.0003天1可计算，湖水最终污染浓度 0 =0.0317mol/L ,超过净化指标 0.001mol/L,故在此条件下，湖的环境将会恶化。B.差分模拟模型：在c与 成线性关系时，得到模拟数据如下表 4:年份13581010 20.821.882.462.883.01年份152030405010 23.133.163.173.173.17在c与 成指数关系时，得到模拟数据如下表 5:年份13581010 20.812.012.813.594.02年份152030405010 24.324.484.564.584.58观察以上数据可知：由于河水污染没得过大，河水污染

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

18、关系时，得到的模拟数据如下表 6:年份15102022354010 40.1850.8011.7927.2729.58152.0497.13年份50829010011214015010 4260.7625.0633.9617.2626.3628.5632.5根据模拟所得数据所作的湖水污染浓度变化趋势如下图5:在c于 成指数关系时，得到的模拟数据如下表 7:年份15102024354510 40.1690.5861.2505.1949.04644.2695.75年份508010012013014015010 4321.7925.5959.6950.6949.6952.5964.8根据模拟所得数据

19、所作的湖水污染浓度变化趋势如下图5:I司题上(2)河水浓度变化, 指数自然洋化0 10.090.080 070 06QQ50 C40 030 020 01.*9M1!1 1*.A*工 子 木 » «* . ! 皋 0* * :J1 « v*1/T , !' *4A1*»1AA* W*a»a11* 9,f /*»圈*vi.i.:*'T! 1 1 , « Hi 1*a*3qa. :14-tiA+9 !T!*d*4i I11 9 Bl¥ JdpIliA口 _ 上1 -iiiiiiii0204080 BO

20、100120140180在第50年时，无论自然净化速率c于湖水污染浓度 成线性关系还是指数关系，可以看到湖水污染浓度都接近工作人员实地测得值，这从检验的角度说明我们的对模型参数的估计是可取的。观察模拟所得数据可知，湖的环境将会恶化。五、模型的评价及评改进5.1 模型的评价(1) .对于首先建立的含参数的微分方程模型，在其有精确解的时候，能给出湖水污染浓度随时间变化的规律，得到很好的预测效果，但微分方程并不是任何条件下都有解，在本问题中，只有当河水污染浓度恒定且自然净化速率呈线性关系时，微分方程可解，故适用的范围受到了限制。(2) .差分模拟法，即用数值方法求解微分方程，它解决了微分方程无精确解

21、的 情况，适用范围很大。使用一个兼顾运算效率和模拟结果精度的时间步长是该方 法的关键。在处理本问题时，采用的时间步长为大。这样处理的精度大约为0.0001,与精确解相当接近。5.2 模型的改进 1） .自然净化速率的改进：资料显示，污染物的很多生物降解过程是通过好氧反应实现的，当污染浓度较大时，需氧量增加，而湖中氧气总量不变，自然净化 速率会有所减小， 所以， 自然净化速率的函数可作为分段函数考虑。 前一段还是用前面的指数模型，后面一段可以考虑为 C e D E（C, D,E 0） 。模型这样改进后，更加符合实际情况。 2） 2） .河水浓度变化规律模型的改进：根据价值规律，一个工厂的发展总是

22、有起有落， 其生产量绕某一中心值上下波动， 故生产河水的浓度的变化应该是具有一定周期性，同时又有所增加。可以考虑在原来的Logistic模型后加上一个傅立叶函数（其周期就是工厂生产的周期）和一个正态随机变量。六、参考文献1数学建模与数学实验高等教育出版社，赵静、但琦主编， 2008年 1 月第三版2 环境化学 ，高等教育出版社出版，戴树桂主编， 1997 年3 月；七、附录附录一 :问题一：（ 1）河水浓度恒为0.001mol/L ，不考虑自然净化，a：采用差分方程模型命令：wenti12（d ， 0.50,0.001， n）b:采用微分方程模型命令： t=1:1:51；P=0.001.*ex

23、p(-0.19*t)Pplot(t,p)grid on(2)河水浓度恒为0.001mol/L ，考虑自然净化，a：采用差分自然净化线性自然净化命令： wenti12('d',0,50,0.001,'y','xx')指数自然净化命令：wenti12('d',0,50,0.001,'y','zs')b：采用微分方程模型，线性自然净化命令：t=1:1:51;p=0.001.*exp(-0.19.*t).*(-1+1.*exp(0.19.*t); phold on plot(t,p) grid on问题二：

24、(1)河水浓度恒为0.05mol/L,不考虑自然净化a：采用差分方程模型命令： wenti12('d',0.50,0.05,'n')b：采用微分方程模型命令： t=1:1:51;p=0.05.*exp(-0.19.*t).*(-1+1.*exp(0.19.*t); phold onplot(t,p)grid on( 2)河水浓度恒为0.05mol/L ，考虑自然净化，a: 采用差分方程模型线性自然净化命令：wenti12('d',0.50,0.05,'y','xx')指数自然净化命令：wneti12('d&

25、#39;,0.50,0.05,'y','zs')b:采用微分方程模型，线性自然净化命令t=1:1:51;p=0.00166945.*exp(-0.2995.*t).*(-19+19.*exp(0.2995.*t);phold onplot(t,p)grid on( 3)河水浓度变化，不考虑自然净化命令： wenti12('d',0.50,-1,'n')考虑自然净化线性自然净化命令：wenti12('d',0,50,-1,'y','xx')指数自然净化命令：wenti12('d

26、',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 采用线性净化模型， zs 采用指数净化模型format compactformat longif nargin 2%5

27、0年,河水污染变化,不考虑自然降解n=50;s=-1;w='n'b='xx'elseif nargin 3 %河水污染变化,不考虑自然降解s=-1;w='n'b='xx'endq=1 elseif l 'm'q=12elseif l 'd'q=365elseif l 'h'q=365*24elseif l 'min'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,wrong param

28、eter.try again!')endA=zeros(1,n+1);B=zeros(1,n+1);A(1)=chu;B(1)=chu;var=0; %自然降解var1=0; %河水污染if w 'nvar=0;elseif w 'y'var=fen(A(1),l,b);elseerror('error,wrong parameter.try again!')endfor i=1:1:nif s -1var1=river(i);elsevar1=s;endfor j=1:1:qif w 'n'var=0;elseif w '

29、;y'var=fen(A(i),l,b);endA(i)=(A(i)*10A15+vaU*1.9*10A14/q)/(10A15+1.9*10A14/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)*10A15+vaT1.9*10A14/q)/(10A15+1.9*10A14/q)-var;B(i+1)=0.19*var1/q+(1-0.19/q)*B(i)-var;endA;Bhold on

30、plot(1:1:n+1,A,'o')grid onfunction p=wenti3(l,bi,year,b,a)%格式 : wenti3(P1,P2,P3,P，4 P5)%P1:y,m,d,h,min,s 表示设定步长值%P2:百分比%P3:可认为的最大净化年限%P4:xx线性自然净化，zs-m数自然净化%P5:湖水浓度初始值format long if l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365;elseif l 'h'q=365*24elseif l 'min

31、'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)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fen(A(i),l,b);t=A(i);if A(i)<0.001 break;end endA(i+1)=(A(i)*10A15+temp*1.9*10A14

32、/q)/(10A15+1.9*10A14/q)-fen(A(i),l,b );if A(i)<0.001 break;endend if i yearp=0;elsep=1;end%程序调试所用%sprintf('停产')sprintf('需要 %d 年天',i,j)%sprintf('第 年%dK为:.12f,i,j-1,t)%sprintf('第年认为:.12f,i,j,A(i)function low=wenti3_2(b,year,l,h)%采用二分法求出百分比format compactformat longlow=l;high=

33、h;cen=(low+high)/2cha=high-low;while cha>10A-6%误差判断条件m=wenti3('d',cen,year,b,0.05); %调用函数if m 0high=cen;cen=(low+high)/2elselow=cen;cen=(low+high)/2endcha=high-low;end high;low;sprintf('工厂在生产为原来的f时，湖水才能在d年内净化。,low,year)function wenti3_3(l,year,bi,b)%格式 : wenti3_3(P1,P2,P3,P4)%P1:y,m,d

34、,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 'hq=365*24elseif l 'min'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,try again!')endA=zeros(1,year);A(1)=0.

35、03;N=0;for i=1:1:yeartemp=0.05*bi;for j=1:1:fix(q)A(i)=(A(i)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fen(A(i),l,b);endA(i+1)=(A(i)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fen(A(i),l,b );N=A(i)被调用的程序：function y=fen(N,l,b)%自然净化%格式 :fen(P1,P2)%P1:前一时间点的污染值%P2:y,m,眼年或月或日%线性关系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 paramet