中国传媒大学 2010 学年第 一 学期
课程
数学建模与数学实验
题 目
学生姓名 学 号 班 级 学生所属学院 任课教师 教师所属学院 成 绩
Pristine 湖污染问题的建模与求解
摘要
本文讨论了湖水污染浓度变化趋势的预测问题。 通过分析水流输入输出湖泊的过程，建立了湖水污染浓度随时间变化的含参变量的微分方程模型，在河水污染浓度恒定和自然净化速率呈线性关系的情况下，求得其精确解，带入具体数据得到结论：在PCA 声称的河水污染浓度下，湖的环境不会恶化；在工作人员实地测得的河水浓度下，湖的环境将会恶化。 同时建立了计算机模拟模型，带入具体数值，运用时间步长法来仿真模拟了在湖水污染浓度稳定以前湖水每天的变化情况，输出自PCA 建厂以来每年的湖水污染浓度，得到与微分方程模型相同的结论。 在全停产和半停产时，通过前面的两个模型可以计算湖水污染浓度在自然净化影响下的恢复到净化指标所需的年限。并可得到结论：在半停产状态下，在选定的自然净化速率常数的约束下，只有当河水污染浓度降至原来的3.15%（自然净化速率呈线性关系），4.7%（自然净化速率呈指数关系），才有可能使河水在100年内恢复至0.001mol/l，然后给出整改建议。
一、问 题 重 述
Pure 河是流入Pristine 湖的唯一河流。50年前PCA 公司在此河旁建起一个生产设施并投入运行。PCA 将为处理的湖水排入河中，导致Pristine 湖被污染。PCA 公司声称：已排放的废水的标准多年从未改变切不会对湖的环境有影响。
现已知：Pristine 湖的湖容量为10L ，流入（流出）的水流速度为1. 9L/年。PCA 公司声称河水污染浓度仅为0.001mol/L，自工厂以来没有改变过。
讨论下列问题：
（1）建立数学模型用PCA 提供的公开数据判断湖的环境是否会恶化； （2）以目前湖水污染浓度0.03mol/L，和河水污染浓度0.05mol/L为新数据判断湖的环境是否会恶化；
15
14
二、模型的合理假设和符号系统
2.1 模型的合理假设
（1）降水量和增发量相等；
（2）湖中流入量和流出量相等且一直未变；
（3）污水量远小于河水注入量，且污水与河水混合均匀； （4）湖水混合均匀，且流入污水的扩散速度无限大； （5）湖内除Pure 河外，无其他污染源；
二
2.2 符号系统
ρ0：河水污染浓度mol/L；
ρ：湖水污染物浓度mol/L；
15
V ：湖泊容量10L ；
c ：自然净化速率mol/（L 。年）
μ：流入（流出）的水流速度1. 914L/年；
t ：从PCA 建厂至考察时刻的时间段。
三、问题的分析
3.1 问题分析：
对于问题中几个词语的理解：
1. 是否会恶化——湖的环境恶化即湖水污染浓度大于0.001mol/L，要判断其是否会恶化，则需计算在某一污染物积累速度（分析影响此速度的因素）下，湖水能达到的最大污染浓度和其变化趋势，以及湖水经几年超过0.001mol/L，经过几年达到最大污染浓度。
2. 自然净化——自然净化是独立的生态系统进行自我调节的方式之一，是在空气，阳光，水和细菌的参与下，进行包括物理沉降，化学反应和生物转化三大方面的活动，其最终作用是将污染物转化为无害物质，从而净化生态系统。当河水输入湖泊并均匀混合之后，影响湖水污染浓度的唯一因素便是自然净化速度。 湖水污染问题水流的动态流程图：
此问题中，我们考察对象是湖水污染浓度的变化趋势：
1. 在整改之前，其增加的趋势，超过净化指标0.001ml/L（即湖水恶化）的可能性和时限；
2. 在全停产（无污染物输入）和半停产的情况下，其降低的趋势，达到净化指标的时限。
3. 在整改之后，其增加的趋势，未定与净化指标之下某一水平的时限、
三
在前假设条件的基础之上，湖水容量不变，出河水外无其他的污染源，故我们可将湖泊作为一个封闭的生态系统，其简化的湖水被污染的动态过程为：受污河输入湖泊，河水与湖水均匀混合，受污河水进行自我的净化，湖水数出湖泊。湖中污染物的量直接决定了湖水污染浓度，而污染物的量受到以下两方面因素的影响：1. 河水的污染浓度及其流入速度（根据已知此速度不变），2. 湖水的自然净化速度，前者使其增加，后者使其减少（负增加）。问题一、二、三的实质都是要分析污染浓度的变化趋势，其去表便在于前一因素的不同。问题一中，河水污染浓度不变，恒为0.001mol/L；问题二中，河水污染浓度可能会变化，受PCA 效益的影响而按一定规律波动；问题三中，在全停产或半停产的情况下，和硕污染浓度为0或减为问题二中的一部分。后一因素（自然净化速度）在三个问题中的作用都是相同的。
根据微积分的知识可知，在适当短的时间段之内，通过建立微分方程，可以将连续的过程离散化，从而可以得到湖水污染浓度与时间之间的关系式。
利用时间步长发，缩小步长值（从年到月到天），并与微分方程所得的精确解做出比较。
四、模型建立与求解
问题一：
根据PCA 的公开申明和所提供数据，可认为：河水污染浓度恒为0.001mol/L。从存在自然净化和不存在自然净化两个方面考虑： （1）．在不考虑自然净化的情况下：
由于假设湖中流入量和流出量相等，而在经过与湖水均匀混合后，流出湖水污染浓度明显减小，故流出污染物的量小于流出污染物的量，污染物将在湖中沉积，从而使湖水污染浓度增加，当其增加至于输入的河水污染浓度相等时，河水污染浓度达到最大，并稳定在这一数值，不在增加。
建立湖水污染浓度随时间变化的微分方程模型：设在极短时间dt ，湖水污染浓度
增加d ρ，在将湖水被污染这一连续动态过程简化为离散的瞬间静止状态（如问题分析中所述）之后，根据湖中剩余量=输入量—输出量，我们可以列微分方程如下：
ρ0⨯μ⨯dt +ρ⨯V
=ρ+d ρ
V +μ⨯dt
d ρμ
=⨯(ρ0-ρ)
化简可得：dt V ······（1）
1415ρ=0.001mol/L的数据，我们可得到ρ和t 的μ=1.9⨯10V =10带入L ，L/年和0-0. 1⨯t 90⨯t . 19
ρ=0. 00⨯1e ⨯e (- 1) 关系式如下：
四
通过此关系式我们可知，当t =∞时，湖水污染浓度将趋近与0.001mol/L，即湖的环境不会恶化。
建立计算机模拟模型：湖水污染浓度的变化时有湖中污染物随时间的积累而引起的，这个逐步积累的过程我们可以用计算机进行仿真模拟，其实质为完成一个循环累加的过程，并可改变时间步长，如一年一年的累积，一月一月的累积，一天一天的累积，从而使我们的模拟值逐步精确，可与微分方程求得精确解比较，分析误差。为提高模拟结果的精确性和运算的效率，我们采用了逐天累加，数出年污染浓度的方式。
模拟程序见附件一，从后面的分析中我们可知：在此河水污染浓度恒定，无自然净化的最简单的情况下建立的模型是以后问题的基础，后面的问题只是改变条件或数据，其实质是不变的，故我们在此程序中加入了多个选择语句，在不同德条件或数据下执行不同的命令，从而用一个程序解决全部的模拟问题。 模拟所得的数据如下表1：
根据模拟数据所作的湖水污染浓度变化趋势图如下图1：
五
模拟所得的数据显示：湖水污染 浓度将稳定于0.001mol/L，这与分析微分方程所得的结果是相符的。
（2）在考虑自然净化的情况下： A. 首先，我们应该了解什么事自然净化。
根据资料显示：湖水中的污染物可分为有机污染物和无机污染物两大类，在多种环境因素（阳光，空气，水，水中生物，水中化学物质，重力等）的作用下，通过物理沉降，化学反应和生物转化一系列复杂的活动，它们的量会发生变化。 有机物在水环境中的迁移转化过程如图2：
可以看出整个过程是相当复杂的，不仅过程多，而且在相同的过程中，不同的物质有着不同的结果。此问题中没有明确给出输入湖中的污染物种类，也没有对湖泊环境做任何描述，无疑给问题的解决增加了极大的难度。为是问题简化，我们从一般的情况出发，假设：湖中所进行的反应均为一级反应，有机物的存在不会对环境参数造成改变，环境固定（自然净化速率恒定），考虑湖泊生态系统中起主要作用的几种过程作简要分析：
1. 物理沉降。不同物质有着不同的沉降速度常数λ，沉降速率为c =λ⨯ρ； 2. 挥发。不同物质有着不同的挥发速率常数分压为0的条件下，挥发速率为
c =
K v
，在有机物在水体上的大气中的
K v ⨯ρ
Z (Z 为水体深度)
K p
3. 水解反应。不同物质有着不同的水解速率常数
，在一级反应的条件下，水
六
解速率为
c =K p ⨯ρ
；
K b
4. 生物降解反应。不同物质有着不同的降解速率常数生物降解速率
c =K b ⨯ρ
，在一级反应的条件下，
；
B. 在最简模拟程序的基础上，从每天的积累量中减去每天的自然净化量K ⨯时间
-1
K =0.0003天步长，重复循环可得到逐年湖水污染浓度的值，设。K 越大，湖
水污染浓度将越快稳定于一个更小的浓度值。 模拟所得的数据如下表2：
根据模拟数据所作的湖水污染浓度变化趋势图如下图3：
C. 进一步考虑自然净化速率的影响。更贴近于实际的情况是，自然净化速率c 与湖水污染浓度ρ成指数关系，c 随着ρ的增加而增加，但增加的速率会逐渐减小，
B ⨯ρ
c =A (1-e ) ，A,B 为与环境和污染物种类有关的常数。用关系式表达即可为：
将此关系式带入微分方程（2），得到一个新的微分方程，此方程无解。但是我们
七
可以通过计算机模拟求得数值解。
A 的大小影响着ρ最终稳定浓度，B 的大小影响着ρ达到最终稳定浓度的快慢。在前面达到线性关系的基础上，在ρ=0/。001mol/L处的c 值大小
确定A 至数量级，并稳定年限定在10—20年间，从而得到一组估计值：
A =8⨯10-4,B =400
模拟所得数据如下表3：
根据模拟所得数据所作的湖水污染浓度变化趋势图如下图4：
观察数据可知：无论自然净化速率c 与污染浓度ρ成线性关系还是指数关系，当河水污染浓度
ρ0=0.001mol/L时，ρ都将稳定于一个小于0.001mol/L的值，
也即：
湖的环境不会恶化。 问题二：
根据实际情况，在此我们只考虑存在自然净化的情况。
八
（1）. 河水污染浓度恒定为0.05mol/L: 与问题一（2）相同，只是
ρ0=0.05mol/L，沿用问题一（2）的微分方程模型和差
分模拟模型，我们可以得到以下结果：
A. 微分方程模型（自然净化速率c 与湖水污染浓度ρ成线性关系）：
当时间t =∞时，
ρ=
ρ0μ
μ+KV ，即湖水污染浓度稳定于一个与K 有关的值。当
ρ0=0.05mol/L时，如要ρ稳定于0.001mol/L则K 值为9.31年-1。根据我们设定的K =0.0003天可计算，湖水最终污染浓度
-1
ρ0=0.0317mol/L，超过净化指标
0.001mol/L，故在此条件下，湖的环境将会恶化。 B. 差分模拟模型：
在c 与ρ成线性关系时，得到模拟数据如下表4：
在c 与ρ成指数关系时，得到模拟数据如下表5：
观察以上数据可知：由于河水污染没得过大，河水污染没得从第一年起就超过了净化指标，并逐年增加，是湖的环境恶化。 （2）. 河水污染浓度变化：
根据实际情况，一个工厂的生产量并非恒定不变，每年每月甚至每天也有所不同，从而起其排放污染物的量也将随生产量的变化而变化的。假设PCA 自建立以来的年排污量服从Logistic 模型，并考虑到污染物的量受到多方面因素的影响，在
九
每一时刻的梁上加上一个服从正太分布，范围在此量的10%内的较小量ξ，则可建立河水污染浓度的模型如下：
ρ0=
M
+ξr ⨯M ⨯t
1+Ce （M,C,r 为与工厂效益有关的常数）
设此时工厂处于文本上升的发展阶段，其变化曲线如下图5：
通过计算机模拟产生随机数（如图中星点所示），带入模拟程序。 在c 于ρ成线性关系时，得到的模拟数据如下表6：
根据模拟所得数据所作的湖水污染浓度变化趋势如下图5：
十
在c 于ρ成指数关系时，得到的模拟数据如下表7：
根据模拟所得数据所作的湖水污染浓度变化趋势如下图5：
在第50年时，无论自然净化速率c 于湖水污染浓度ρ成线性关系还是指数关系，可以看到湖水污染浓度都接近工作人员实地测得值，这从检验的角度说明我们的对模型参数的估计是可取的。
观察模拟所得数据可知，湖的环境将会恶化。
五、模型的评价及评改进
5.1模型的评价
（1）. 对于首先建立的含参数的微分方程模型，在其有精确解的时候，能给出湖水污染浓度随时间变化的规律，得到很好的预测效果，但微分方程并不是任何条件下都有解，在本问题中，只有当河水污染浓度恒定且自然净化速率呈线性关系时，微分方程可解，故适用的范围受到了限制。 （2）. 差分模拟法，即用数值方法求解微分方程，它解决了微分方程无精确解的 情况，适用范围很大。使用一个兼顾运算效率和模拟结果精度的时间步长是该方 法的关键。在处理本问题时，采用的时间步长为天。这样处理的精度大约为0.0001, 与精确解相当接近。
5.2模型的改进
（1）. 自然净化速率的改进：资料显示，污染物的很多生物降解过程是通过好氧 反应实现的，当污染浓度较大时，需氧量增加，而湖中氧气总量不变，自然净化 速率会有所减小，所以，自然净化速率的函数可作为分段函数考虑。前一段还是
-D ⨯ρ+E
C ⨯e (C , D , E >0) 。模型这样改进用前面的指数模型，后面一段可以考虑为
后，更加符合实际情况。
（2）. 河水浓度变化规律模型的改进：根据价值规律，一个工厂的发展总是有起 有落，其生产量绕某一中心值上下波动，故生产河水的浓度的变化应该是具有一 定周期性，同时又有所增加。可以考虑在原来的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.*exp(-0.19*t) P
plot(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)); p
hold on plot(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)); p
hold on plot(t,p)
grid on
（2）河水浓度恒为 0.05mol/L，考虑自然净化， a:采用差分方程模型 线性自然净化
命令：wenti12('d',0.50,0.05,'y','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 on plot(t,p)
grid on （3）河水浓度变化， 不考虑自然净化
命令：wenti12('d',0.50,-1,'n') 考虑自然净化 线性自然净化
命令： wenti12('d',0,50,-1,'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 采用线性净化模型，zs 采用指数净化模型
format compact format long
if nargin 2 %50年, 河水污染变化, 不考虑自然降解 n=50; s=-1; w='n'; b='xx';
elseif nargin 3 %河水污染变化, 不考虑自然降解 s=-1; w='n'; b='xx'; end if l 'y' q=1 elseif l 'm' q=12
elseif l 'd' q=365 elseif l 'h' q=365*24 elseif l 'min' q=365*24*60 elseif l 's'
q=365*24*60*60 else
error('error,wrong parameter...try again!') end
A=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); else
error('error,wrong parameter..try again!') end
for i=1:1:n if s -1
var1=river(i); else
var1=s; end
for j=1:1:q if w 'n' var=0; elseif w 'y'
var=fen(A(i),l,b); end
A(i)=(A(i)*10^15+var1*1.9*10^14/q)/(10^15+1.9*10^14/q)-var; B(i)=0.19*var1/q+(1-0.19/q)*B(i)-var; end if w 'n'
var=0; elseif w 'y'
var=fen(A(i),l,b); end
A(i+1)=(A(i)*10^15+var1*1.9*10^14/q)/(10^15+1.9*10^14/q)-var; B(i+1)=0.19*var1/q+(1-0.19/q)*B(i)-var; end [A;B] hold on
plot(1:1:n+1,A,'o') grid on
function 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 compact format long if l 'y' q=1 elseif l 'm' q=12 elseif l 'd' q=365; elseif l 'h' q=365*24 elseif l 'min' q=365*24*60 elseif l 's'
q=365*24*60*60 else
error('error,try again!') end
A=zeros(100); B=zeros(100); temp=0.05*bi; A(1)=a;
for i=1:1:year for j=1:1:q
A(i)=(A(i)*10^15+temp*1.9*10^14/q)/(10^15+1.9*10^14/q)-fen(A(i),l,b);
t=A(i);
if A(i)
A(i+1)=(A(i)*10^15+temp*1.9*10^14/q)/(10^15+1.9*10^14/q)-fen(A(i),l,b );
if A(i)
if i year p=0; else p=1; end
%程序调试所用 %sprintf('停产')
sprintf('需要 %d 年%d天~~',i,j)
%sprintf('第%d年%d天为:%.12f',i,j-1,t) %sprintf('第%d年%d天为:%.12f',i,j,A(i))
function low=wenti3_2(b,year,l,h) %采用二分法求出百分比 format compact format long low=l; high=h;
cen=(low+high)/2 cha=high-low;
while cha>10^-6 %误差判断条件
m=wenti3('d',cen,year,b,0.05); %调用函数 if m 0
high=cen;
cen=(low+high)/2 else
low=cen;
cen=(low+high)/2 end
cha=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,h,min,s 表示设定步长值 %P2:year %P3: 百分比
%P4:xx-线性自然净化,zs-指数自然净化 format compact format long if l 'y' q=1 elseif l 'm' q=12 elseif l 'd' q=365; elseif l 'h' q=365*24 elseif l 'min' q=365*24*60 elseif l 's'
q=365*24*60*60 else
error('error,try again!') end
A=zeros(1,year); A(1)=0.03; N=0;
for i=1:1:year temp=0.05*bi; for j=1:1:fix(q)
A(i)=(A(i)*10^15+temp*1.9*10^14/q)/(10^15+1.9*10^14/q)-fen(A(i),l,b); end
A(i+1)=(A(i)*10^15+temp*1.9*10^14/q)/(10^15+1.9*10^14/q)-fen(A(i),l,b ); end N=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; else
error('error,wrong parameter...try again!') end
elseif b 'zs' %指数关系
second=8*10^-4; if l 'y'
y=second*(1-exp(-N*400)); elseif l 'm'
y=second*(1-exp(-N*400))/12; elseif l 'd'
y=second*(1-exp(-N*400))/365; else
error('error,wrong parameter...try again!') end else
error('error,wrong parameter...try again!') end
%河水污染浓度函数的参数方程， function F=fun(x) M=x(1); C=x(2); k=x(3);
F=[eps-M/(1+C);
0.05-M/(1+C*exp(-k*M*50)); 0.1- M/(1+C*exp(-k*M*100))];
function y=river(n) %河水浓度变化
%本函数是求第n 年的河水污染浓度, 返回值是第n 年的河水污染浓度 %分为三个阶段； %一：%河水污染浓度是一个关于时间n 的函数，先求出此函数的三个未知参数M ， C ，k
%二：%在所求出的函数上加上一个小的波动，更能满足实际要求 %三：%计算第n 年的河水污染浓度 %一：求参数
%由于要提供搜索点，所以先在mathmatica 中用solve 命令求出大致的M ，C ，k 值，作为搜索点
%mathematica中的程序为：
% Solve[{0.0001 M/(1 + C),
% 0.05 M/(1 + C*E^(-k*M*50)),
% 0.1 M/(1 + C*E^(-k*M*100))}, {M, C, k}]
%x0=[0.1001;1000;1.37997]; %搜索点 %options=optimset('Display','iter');
%[x,fval]=fsolve('fun',x0,options) %解方程组
x=1000*[0.[1**********]107 1.[1**********]742 0.[1**********]106];
i=1:1:100;
g=x(1)./(1+x(2).*exp(-x(1).*x(3).*i)); %plot(i,g)
%二:加上一个波动,10%的波动范围 h=zeros(1,100); sigma=0.03; for j=1:100
para=sigma*randn; h(j)=g(j)*(1+para); end j=1:100;
%画出波动图, 调试时用 %plot(i,g,j,h,'xr') %计算第n 年的浓度
y=x(1)./(1+x(2).*exp(-x(1).*x(3).*n)); y=y*(1+sigma*randn);