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

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

2、拟了 在湖水污染浓度稳定以前湖水每天的变化情况，输出自PCA建厂以来每年的湖水 污染浓度，得到与微分方程模型相同的结论。在全停产和半停产时，通过前面的两个模型可以计算湖水污染浓度在自然净 化影响下的恢复到净化指标所需的年限。并可得到结论：在半停产状态下，在选定的自然净化速率常数的约束下，只有当河水污染浓度降至原来的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) 湖中流入量和流出量相等且一直未变；(3)

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

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

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

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

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

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

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

11、下图1：M TO召门右冃址皿扌"O102030405060五模拟所得的数据显示：湖水污染 浓度将稳定于0.001mol/L，这与分析微分方程 所得的结果是相符的。（2）在考虑自然净化的情况下：A.首先，我们应该了解什么事自然净化。根据资料显示：湖水中的污染物可分为有机污染物和无机污染物两大类，在多种环境因素邙阳光，空气，水，水中生物，水中化学物质，重力等）的作用下，通 过物理沉降，化学反应和生物转化一系列复杂的活动，它们的量会发生变化。有机物在水环境中的迁移转化过程如图 2:可以看出整个过程是相当复杂的， 不仅过程多，而且在相同的过程中，不同的物 质有着不同的结果。此问题中没有明确给

12、出输入湖中的污染物种类， 也没有对湖 泊环境做任何描述，无疑给问题的解决增加了极大的难度。 为是问题简化，我们 从一般的情况出发，假设：湖中所进行的反应均为一级反应，有机物的存在不会 对环境参数造成改变，环境固定（自然净化速率恒定），考虑湖泊生态系统中起 主要作用的几种过程作简要分析：1. 物理沉降。不同物质有着不同的沉降速度常数 ，沉降速率为c = ；K2. 挥发。不同物质有着不同的挥发速率常数v，在有机物在水体上的大气中的Kv 'c =分压为0的条件下，挥发速率为Z （Z为水体深度）3. 水解反应。不同物质有着不同的水解速率常数 Kp，在一级反应的条件下，水解速率为c = Kp J

13、；4. 生物降解反应。不同物质有着不同的降解速率常数Kb，在一级反应的条件下，生物降解速率° = Kb ?；B.在最简模拟程序的基础上，从每天的积累量中减去每天的自然净化量 K时间 步长，重复循环可得到逐年湖水污染浓度的值，设 K =0.0003天。K越大，湖 水污染浓度将越快稳定于一个更小的浓度值。模拟所得的数据如下表2：年份135810P 10A1.5873.4374.3154.8784.984年份1520304050p 10A5.0875.1035.1065.1065.106用关系式表达即可为：C二A(1 -e )，A,B为与环境和污染物种类有关的常数。 将此关系式带入微分方程

14、(2)，得到一个新的微分方程，此方程无解。但是我们 七可以通过计算机模拟求得数值解。A的大小影响着'最终稳定浓度，B的大小影响着'达到最终稳定浓 度的快慢。在前面达到线性关系的基础上，在'=0/。001mol/L处的c值大小 确定A至数量级，并稳定年限定在10 20年间，从而得到一组估计值：A =8 10°,B =400模拟所得数据如下表3：年份135810P 10A1.4932.9713.5393.8163.871年份1520304050P 10A3.9033.9063.9063.9063.906根据模拟所得数据所作的湖水污染浓度变化趋势图如下图4：观察数据

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

16、001mol/L则K值为9.31年'。根据我们设定的K =0.0003天1可计算，湖水最终污染浓度：0 =0.0317mol/L，超过净化指标0.001mol/L，故在此条件下，湖的环境将会恶化。B.差分模拟模型：在c与'成线性关系时，得到模拟数据如下表 4：年份135810P E0.821.882.462.883.01年份1520304050P 10"3.133.163.173.173.17在c与"成指数关系时，得到模拟数据如下表 5：年份135810P 10"0.812.012.813.594.02年份1520304050P 10"4

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

18、uot;成线性关系时，得到的模拟数据如下表 6：年份151020223540P0.1850.8011.7927.2729.58152.0497.1310°年份508290100112140150P260.7625.0633.9617.2626.3628.5632.510°根据模拟所得数据所作的湖水污染浓度变化趋势如下图5：在C于'成指数关系时，得到的模拟数据如下表 7:年份151020243545P0.1690.5861.2505.1949.04644.2695.7510°年份5080100120130140150P321.7925.5959.6950.6

19、949.6952.5964.810°根据模拟所得数据所作的湖水污染浓度变化趋势如下图5:十一在第50年时，无论自然净化速率c于湖水污染浓度'成线性关系还是指数关系, 可以看到湖水污染浓度都接近工作人员实地测得值， 这从检验的角度说明我们的 对模型参数的估计是可取的。观察模拟所得数据可知，湖的环境将会恶化。五、模型的评价及评改进5.1模型的评价（1）.对于首先建立的含参数的微分方程模型，在其有精确解的时候，能给出湖 水污染浓度随时间变化的规律，得到很好的预测效果，但微分方程并不是任何条 件下都有解，在本问题中，只有当河水污染浓度恒定且自然净化速率呈线性关系 时，微分方程可解，故

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

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

22、#39; n'b:采用微分方程模型命令： t=1:1:51;P=0.001.*exp(-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.*e

23、xp(0.19.*t);Phold onplot(t,p)grid on问题二：(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 on plot(t,p) grid on(2) 河水浓度恒为0.05mol/L，考虑自然净化，a:采用差分方程模型线性自然净化命令：wenti12('d',0.50,0.05,'y'

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

25、;y','xx')指数自然净化命令：we nti12('d',0,50,-1,'y','zs')附录二：主程序：fun cti on wen ti12(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采用线性净化模型，z

26、s采用指数净化模型format compactformat lo ngif 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 'mi n'q=365*24*60

27、elseif l 's'q=365*24*60*60elseerror('error,wr ong parameter.try agai n!')endA=zeros(1, n+1);B=zeros(1, n+1);A(1)=chu;B(1)=chu;var=0; %自然降解vaM=0; %河水污染if w 'n'var=0;elseif w 'y'var=fe n(A(1),l,b);elseerror('error,wr ong parameter.try agai n!')endfor i=1:1: nif

28、s -1var1=river(i);elsevar1=s;endfor j=1:1:qif w 'n'var=0;elseif w 'y'var=fe n(A(i),l,b);endA(i)=(A(i)*10A15+var1*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=fe n(A(i),l,b);endA(i+1)=(A(i)*10A15+var1*1.9

29、*10A14/q)/(10A15+1.9*10A14/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 onfun cti on p=we nti3(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 lo ngif l 'y&

30、#39;q=1elseif l 'm'q=12elseif l 'd'q=365;elseif l 'h'q=365*24elseif l 'mi n'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,try agai n!')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)/

31、(10A15+1.9*10A14/q)-fe n(A(i),l,b);t=A(i);if A(i)<0.001 break;endendA(i+1)=(A(i)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fe n(A(i),l,b )；if A(i)<0.001 break;endendif i yearp=0；elsep=1；end%程序调试所用%spri ntf('停产')sprintf('需要 %d 年4天',i,j)%spri ntf('第4年4 天为:.12f,i,j-1,t)%spri

32、 ntf('第4年4 天为:.12f,i,j,A(i)fun cti on low=we nti3_2(b,year,l,h)%采用二分法求出百分比format compactformat lo nglow=l;high=h;cen=(low+high)/2 cha=high-low;while cha>10A-6%误差判断条件m=we nti3('d',ce n,year,b,0.05);%调用函数if m 0high=ce n;cen=(low+high)/2elselow=ce n;cen=(low+high)/2endcha=high-low;end hi

33、gh;十七low;sprintf('工厂在生产为原来的f时，湖水才能在d年内净化。：low,year)fun cti on wen ti3_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 lo ngif l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365;elseif l 'h'q=365*24el

34、seif l 'mi n'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,try agai n!')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)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fe n(A(i),l,b);endA(i+1)=(A(i)*10A15+temp*1.9*10A14/q”(10A15+1.9*

35、10A14/q)-fe n(A(i),l,b);endN=A(i)被调用的程序：fun cti on 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;else error('error,wr ong parameter.try aga in!') endelseif b 'zs'%指数