湖水污染問題的數(shù)學建模與求解

1、中國傳媒大學2010學年第一學期題目Pristine湖污染問題的建模與求解學生姓名學號班級學生所屬學院任課教師教師所屬學院成績Pristine湖污染問題的建模與求解摘要本文討論了湖水污染濃度變化趨勢的預測問題。通過分析水流輸入輸出湖泊的過程，建立了湖水污染濃度隨時間變化的含參變量的微分方程模型，在河水污染濃度恒定和自然凈化速率呈線性關系的情況下，求得其精確解，帶入具體數(shù)據(jù)得到結論：在PCA聲稱的河水污染濃度下，湖的環(huán)境不會惡化；在工作人員實地測得的河水濃度下，湖的環(huán)境將會惡化。同時建立了計算機模擬模型，帶入具體數(shù)值，運用時間步長法來仿真模擬了在湖水污染濃度穩(wěn)定以前湖水每天的變化情況,輸出自PC

2、A建廠以來每年的湖水污染濃度，得到與微分方程模型相同的結論。在全停產和半停產時，通過前面的兩個模型可以計算湖水污染濃度在自然凈化影響下的恢復到凈化指標所需的年限。并可得到結論：在半停產狀態(tài)下，在選定的自然凈化速率常數(shù)的約束下，只有當河水污染濃度降至原來的%（自然凈化速率呈線性關系），%（自然凈化速率呈指數(shù)關系），才有可能使河水在100年內恢復至l,然后給出整改建議。一、問題重述Pure河是流入Pristine湖的唯一河流。50年前PCA公司在此河旁建起一個生產設施并投入運行。PCA將為處理的湖水排入河中，導致Pristine湖被污染。PCA公司聲稱：已排放的廢水的標準多年從未改變切不會對湖的環(huán)

3、境有影響。現(xiàn)已知：Pristine湖的湖容量為I015L,流入（流出）的水流速度為194L/年。PCA公司聲稱河水污染濃度僅為L,自工廠以來沒有改變過。討論下列問題：（1）建立數(shù)學模型用PCA提供的公開數(shù)據(jù)判斷湖的環(huán)境是否會惡化；（2）以目前湖水污染濃度L，和河水污染濃度L為新數(shù)據(jù)判斷湖的環(huán)境是否會惡化；二、模型的合理假設和符號系統(tǒng)模型的合理假設（1）降水量和增發(fā)量相等；（2）湖中流入量和流出量相等且一直未變；（3）污水量遠小于河水注入量,且污水與河水混合均勻；（4）湖水混合均勻,且流入污水的擴散速度無限大；（5）湖內除Pure河外，無其他污染源；符號系統(tǒng)0：河水污染濃度mol/L；P:湖水污

4、染物濃度mol/L；V：湖泊容量1015L；c:自然凈化速率mol/（L。年）卩：流入（流出）的水流速度1.914L/年；t:從PCA建廠至考察時刻的時間段。三、問題的分析問題分析：對于問題中幾個詞語的理解：1是否會惡化一一湖的環(huán)境惡化即湖水污染濃度大于L,要判斷其是否會惡化，則需計算在某一污染物積累速度（分析影響此速度的因素）下,湖水能達到的最大污染濃度和其變化趨勢，以及湖水經幾年超過L,經過幾年達到最大污染濃度。2. 自然凈化一一自然凈化是獨立的生態(tài)系統(tǒng)進行自我調節(jié)的方式之一，是在空氣，陽光，水和細菌的參與下，進行包括物理沉降，化學反應和生物轉化三大方面的活動，其最終作用是將污染物轉化為無

5、害物質，從而凈化生態(tài)系統(tǒng)。當河水輸入湖泊并均勻混合之后，影響湖水污染濃度的唯一因素便是自然凈化速度。湖水污染問題水流的動態(tài)流程圖：此問題中，我們考察對象是湖水污染濃度的變化趨勢：1在整改之前，其增加的趨勢，超過凈化指標L（即湖水惡化）的可能性和時限；2.在全停產（無污染物輸入）和半停產的情況下，其降低的趨勢，達到凈化指標的時限。3. 在整改之后，其增加的趨勢，未定與凈化指標之下某一水平的時限、在前假設條件的基礎之上，湖水容量不變，出河水外無其他的污染源，故我們可將湖泊作為一個封閉的生態(tài)系統(tǒng)，其簡化的湖水被污染的動態(tài)過程為：受污河輸入湖泊，河水與湖水均勻混合，受污河水進行自我的凈化，湖水數(shù)出湖泊

6、。湖中污染物的量直接決定了湖水污染濃度，而污染物的量受到以下兩方面因素的影響：1.河水的污染濃度及其流入速度（根據(jù)已知此速度不變），2.湖水的自然凈化速度，前者使其增加，后者使其減少（負增加）。問題一、二、三的實質都是要分析污染濃度的變化趨勢，其去表便在于前一因素的不同。問題一中，河水污染濃度不變，恒為L；問題二中，河水污染濃度可能會變化，受PCA效益的影響而按一定規(guī)律波動；問題三中，在全停產或半停產的情況下，和碩污染濃度為0或減為問題二中的一部分。后一因素（自然凈化速度）在三個問題中的作用都是相同的。根據(jù)微積分的知識可知，在適當短的時間段之內，通過建立微分方程，可以將連續(xù)的過程離散化，從而可

7、以得到湖水污染濃度與時間之間的關系式。利用時間步長發(fā)，縮小步長值（從年到月到天），并與微分方程所得的精確解做出比較。四、模型建立與求解問題一：根據(jù)PCA的公開申明和所提供數(shù)據(jù)，可認為：河水污染濃度恒為L。從存在自然凈化和不存在自然凈化兩個方面考慮：（1）在不考慮自然凈化的情況下：由于假設湖中流入量和流出量相等，而在經過與湖水均勻混合后，流出湖水污染濃度明顯減小，故流出污染物的量小于流出污染物的量，污染物將在湖中沉積，從而使湖水污染濃度增加，當其增加至于輸入的河水污染濃度相等時，河水污染濃度達到最大，并穩(wěn)定在這一數(shù)值，不在增加。建立湖水污染濃度隨時間變化的微分方程模型：設在極短時間dt，湖水污染

8、濃度增加d卩，在將湖水被污染這一連續(xù)動態(tài)過程簡化為離散的瞬間靜止狀態(tài)（如問題分析中所述）之后,根據(jù)湖中剩余量=輸入量輸出量，我們可以列微分方程如下：化簡可得：dpdtLl二Vx(p0-p)1）帶入V-1015L,1.9xl。14L/年和p0=0.001mol/L的數(shù)據(jù)，我們可得到p和t的關系式如下.p=0.001xe-o.19xtx（eo.19xt1）通過此關系式我們可知，當t=時，湖水污染濃度將趨近與L，即湖的環(huán)境不會惡化。建立計算機模擬模型：湖水污染濃度的變化時有湖中污染物隨時間的積累而引起的，這個逐步積累的過程我們可以用計算機進行仿真模擬，其實質為完成一個循環(huán)累加的過程，并可改變時間步長

9、，如一年一年的累積，一月一月的累積，一天一天的累積，從而使我們的模擬值逐步精確，可與微分方程求得精確解比較，分析誤差。為提高模擬結果的精確性和運算的效率，我們采用了逐天累加，數(shù)出年污染濃度的方式。模擬程序見附件一，從后面的分析中我們可知：在此河水污染濃度恒定，無自然凈化的最簡單的情況下建立的模型是以后問題的基礎，后面的問題只是改變條件或數(shù)據(jù)，其實質是不變的，故我們在此程序中加入了多個選擇語句，在不同德條件或數(shù)據(jù)下執(zhí)行不同的命令，從而用一個程序解決全部的模擬問題。模擬所得的數(shù)據(jù)如下表1：年份135810年份1520304050根據(jù)模擬數(shù)據(jù)所作的湖水污染濃度變化趨勢圖如下圖1：模擬所得的數(shù)據(jù)顯示：

10、湖水污染濃度將穩(wěn)定于L,這與分析微分方程所得的結果是相符的。（2）在考慮自然凈化的情況下：A. 首先，我們應該了解什么事自然凈化。根據(jù)資料顯示：湖水中的污染物可分為有機污染物和無機污染物兩大類，在多種環(huán)境因素（陽光，空氣，水，水中生物，水中化學物質，重力等）的作用下，通過物理沉降，化學反應和生物轉化一系列復雜的活動，它們的量會發(fā)生變化。有機物在水環(huán)境中的遷移轉化過程如圖2：可以看出整個過程是相當復雜的，不僅過程多，而且在相同的過程中，不同的物質有著不同的結果。此問題中沒有明確給出輸入湖中的污染物種類，也沒有對湖泊環(huán)境做任何描述，無疑給問題的解決增加了極大的難度。為是問題簡化，我們從一般的情況出

11、發(fā)，假設：湖中所進行的反應均為一級反應，有機物的存在不會對環(huán)境參數(shù)造成改變，環(huán)境固定（自然凈化速率恒定），考慮湖泊生態(tài)系統(tǒng)中起主要作用的幾種過程作簡要分析：1物理沉降。不同物質有著不同的沉降速度常數(shù)九，沉降速率為c二九xp；K2揮發(fā)。不同物質有著不同的揮發(fā)速率常數(shù)v，在有機物在水體上的大氣中的分KxpC二一v壓為0的條件下，揮發(fā)速率為Z（Z為水體深度）3水解反應。不同物質有著不同的水解速率常數(shù)Kp，在一級反應的條件下，水解速率為C二Kxp；4生物降解反應。不同物質有著不同的降解速率常數(shù)Kb，在一級反應的條件下，生物降解速率C=KbXp;B. 在最簡模擬程序的基礎上，從每天的積累量中減去每天的自

12、然凈化量Kx時間步長，重復循環(huán)可得到逐年湖水污染濃度的值，設K=°°°°3天-1。K越大，湖水污染濃度將越快穩(wěn)定于一個更小的濃度值。模擬所得的數(shù)據(jù)如下表2：年份135810年份1520304050根據(jù)模擬數(shù)據(jù)所作的湖水污染濃度變化趨勢圖如下圖3：C. 進一步考慮自然凈化速率的影響。更貼近于實際的情況是，自然凈化速率c與湖水污染濃度p成指數(shù)關系，c隨著p的增加而增加，但增加的速率會逐漸減小，用關系式表達即可為：c=A（1-eBxp），A,B為與環(huán)境和污染物種類有關的常數(shù)。將此關系式帶入微分方程（2），得到一個新的微分方程，此方程無解。但是我們可以通過計算機

13、模擬求得數(shù)值解。A的大小影響著P最終穩(wěn)定濃度，B的大小影響著p達到最終穩(wěn)定濃度的快慢。在前面達到線性關系的基礎上，在p=0/。001mol/L處的c值大小確定A至數(shù)量級，并穩(wěn)定年限定在1020年間，從而得到一組估計值：模擬所得數(shù)據(jù)如下表3：年份135810年份1520304050根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢圖如下圖4：觀察數(shù)據(jù)可知：無論自然凈化速率c與污染濃度P成線性關系還是指數(shù)關系，當河水污染濃度P°=L時，P都將穩(wěn)定于一個小于L的值，也即：湖的環(huán)境不會惡化。問題二：根據(jù)實際情況，在此我們只考慮存在自然凈化的情況。（1）.河水污染濃度恒定為L:與問題一（2）相同，只是

14、Po=L，沿用問題一（2）的微分方程模型和差分模擬模型,我們可以得到以下結果：A.微分方程模型（自然凈化速率c與湖水污染濃度P成線性關系）：p二當時間t=時，卩+KV，即湖水污染濃度穩(wěn)定于一個與K有關的值。當p0=L時，如要p穩(wěn)定于L則K值為9.31年-1。根據(jù)我們設定的K二0.0003天-1可計算，湖水最終污染濃度P0=L,超過凈化指標L，故在此條件下，湖的環(huán)境將會惡化。B.差分模擬模型：在c與P成線性關系時，得到模擬數(shù)據(jù)如下表4：年份135810年份1520304050在c與P成指數(shù)關系時，得到模擬數(shù)據(jù)如下表5：年份135810年份1520304050觀察以上數(shù)據(jù)可知：由于河水污染沒得過大

15、，河水污染沒得從第一年起就超過了凈化指標，并逐年增加，是湖的環(huán)境惡化。（2）.河水污染濃度變化：根據(jù)實際情況，一個工廠的生產量并非恒定不變，每年每月甚至每天也有所不同，從而起其排放污染物的量也將隨生產量的變化而變化的。假設PCA自建立以來的年排污量服從Logistic模型，并考慮到污染物的量受到多方面因素的影響，在每一時刻的梁上加上一個服從正太分布，范圍在此量的10%內的較小量E，則可建立河水污染濃度的模型如下：M1+CerxMxtM,C,r為與工廠效益有關的常數(shù)）設此時工廠處于文本上升的發(fā)展階段，其變化曲線如下圖5：通過計算機模擬產生隨機數(shù)（如圖中星點所示），帶入模擬程序在c于P成線性關系時

16、，得到的模擬數(shù)據(jù)如下表6：年份151020223540年份508290100112140150根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢如下圖5在c于P成指數(shù)關系時，得到的模擬數(shù)據(jù)如下表7：年份151020243545年份5080100120130140150根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢如下圖5：在第50年時，無論自然凈化速率c于湖水污染濃度P成線性關系還是指數(shù)關系，可以看到湖水污染濃度都接近工作人員實地測得值，這從檢驗的角度說明我們的對模型參數(shù)的估計是可取的。觀察模擬所得數(shù)據(jù)可知，湖的環(huán)境將會惡化。五、模型的評價及評改進模型的評價（1）.對于首先建立的含參數(shù)的微分方程模型，在其

17、有精確解的時候，能給出湖水污染濃度隨時間變化的規(guī)律，得到很好的預測效果，但微分方程并不是任何條件下都有解，在本問題中，只有當河水污染濃度恒定且自然凈化速率呈線性關系時，微分方程可解，故適用的范圍受到了限制。（2）.差分模擬法，即用數(shù)值方法求解微分方程，它解決了微分方程無精確解的情況，適用范圍很大。使用一個兼顧運算效率和模擬結果精度的時間步長是該方法的關鍵。在處理本問題時，采用的時間步長為天。這樣處理的精度大約為,與精確解相當接近。模型的改進（1）.自然凈化速率的改進：資料顯示，污染物的很多生物降解過程是通過好氧反應實現(xiàn)的，當污染濃度較大時，需氧量增加，而湖中氧氣總量不變，自然凈化速率會有所減小

18、，所以，自然凈化速率的函數(shù)可作為分段函數(shù)考慮。前一段還是用前面的指數(shù)模型，后面一段可以考慮為Cxe-Dxp+E（C,D,E>°）。模型這樣改進后，更加符合實際情況。（2）.河水濃度變化規(guī)律模型的改進：根據(jù)價值規(guī)律，一個工廠的發(fā)展總是有起有落，其生產量繞某一中心值上下波動，故生產河水的濃度的變化應該是具有一定周期性，同時又有所增加?？梢钥紤]在原來的Logistic模型后加上一個傅立葉函數(shù)（其周期就是工廠生產的周期）和一個正態(tài)隨機變量。六、參考文獻1 數(shù)學建模與數(shù)學實驗高等教育出版社，趙靜、但琦主編，2008年1月第三版2 環(huán)境化學，高等教育出版社出版，戴樹桂主編，1997年3月；

19、七、附錄附錄一:問題一：（1）河水濃度恒為L,不考慮自然凈化，a：采用差分方程模型命令：wentil2（'d'，,，'n'）b：采用微分方程模型命令：t=1:1:51；P=.*exp*t）Pplot（t,p）gridon（2）河水濃度恒為L，考慮自然凈化，a：采用差分自然凈化線性自然凈化命令：wentil2('d',0,50,'y','xx')指數(shù)自然凈化命令：wentil2('d',0,50,'y','zs')b采用微分方程模型，線性自然凈化命令：t=1:1:51;

20、p=.*exp.*t).*(-1+1.*exp.*t);pholdonplot(t,p)gridon問題二：(1) 河水濃度恒為L,不考慮自然凈化a：采用差分方程模型命令：wenti12('d','n')b:采用微分方程模型命令：t=1:1:51;p=.*exp.*t).*(-1+1.*exp.*t);pholdonplot(t,p)gridon(2) 河水濃度恒為L，考慮自然凈化，a:采用差分方程模型線性自然凈化命令：wenti12('d','y','xx')指數(shù)自然凈化命令：wneti12('d'

21、;,'y','zs')b采用微分方程模型，線性自然凈化命令t=1:1:51;p=.*exp.*t).*(-19+19.*exp.*t);pholdonplot(t,p)gridon(3) 河水濃度變化，不考慮自然凈化命令：wenti12('d',-1,'n')考慮自然凈化線性自然凈化命令：wenti12('d',0,50,T,'y','xx')指數(shù)自然凈化命令：wenti12('d',0,50,-1,'y','zs')附錄二：主程序：f

22、unctionwenti12(l,chu,n,s,w,b)%按照時間步長法，以每年或每月或每日為時間段，求出濃度的變化，最后輸出每年的濃度，做一個大致的觀察,%格式：wentil2(P1,P2,P3,P4,P5,P6)%Pl:y,m,d,h,min,s表示設定步長值%P2：湖水污染初值%P3:年份%P4：河流污染初值%P5:是否考慮自然降解,n-不考慮,y-考慮%P6:xx采用線性凈化模型，zs采用指數(shù)凈化模型formatcompactformatlongifnargin2%50年,河水污染變化，不考慮自然降解n=50;s=-1;w='n'b='xx'elsei

23、fnargin3%河水污染變化，不考慮自然降解s=-1;w='n'b='xx'endifl'y'q=1elseifl'm'q=12elseifl'd'q=365elseifl'h'q=365*24elseifl'min'q=365*24*60elseifl's'q=365*24*60*60elseerror('error,wrongparameter.tryagain!')endA=zeros(1,n+1);B=zeros(1,n+1);A(1)=ch

24、u;B(1)=chu;var=0;%自然降解varl=O;%河水污染Q，ifw'n'var=0;elseifw'y'var=fen(A(1),l,b);elseerror('error,wrongparameter.tryagain!')endfori=1:1:nifs-1var1=river(i);elsevar1=s;endforj=1:1:qQ，ifw'n'var=0;elseifw'y'var=fen(A(i),l,b);endA(i)=(A(i)*10"15+varl*10"14/q)

25、/(10"15+*10“14/q)-var;B(i)=*var1/q+q)*B(i)-var;endQ，ifw'n'var=0;elseifw'y'var=fen(A(i),l,b);endA(i+1)=(A(i)*10"15+varl*10"14/q)/(10"15+*10“14/q)-var;B(i+1)=*var1/q+q)*B(i)-var;endA;Bholdonplot(1:1:n+1,A,'o')gridonfunctionp=wenti3(l,bi,year,b,a)%格式:wenti3(P

26、l,P2,P3,P4,P5)%Pl:y,m,d,h,min,s表示設定步長值%P2:百分比%P3:可認為的最大凈化年限%P4:xx-線性自然凈化,zs-指數(shù)自然凈化%P5：湖水濃度初始值formatcompactformatlongifl'y'q=1elseifl'm'q=12elseifl'd'q=365;elseifl'h'q=365*24elseifl'min'q=365*24*60elseifl's'q=365*24*60*60elseerror('error,tryagain!&#

27、39;)endA=zeros(100);B=zeros(100);temp=*bi;A(1)=a;fori=1:1:yearforj=1:1:qA(i)=(A(i)*1015+temp*10"14/q)/(10"15+*10"14/q)-fen(A(i),l,b);t=A(i);ifA(i)<endendA(i+1)=(A(i)*1015+temp*10"14/q)/(10"15+*10"14/q)-fen(A(i),l,b);ifA(i)<endendifiyearp=0;elsep=1;end%程序調試所用%sprin

28、tf('停產')sprintf('需要小年%d天',1,)%sprintf('第%d年小天為:.12f',i,jT,t)%sprintf('第%d年小天為:.12f',i,j,A(i)functionlow=wenti3_2(b,year,l,h)%采用二分法求出百分比formatcompactformatlonglow=l;high=h;cen=(low+high)/2cha=high-low;whilecha>10-6%誤差判斷條件m二wenti3('d',cen,yeAR,b,;%調用函數(shù)ifm0hig

29、h=cen;cen=(low+high)/2elselow=cen;cen=(low+high)/2endcha=high-low;endhigh;low;',low,year)sprintf('工廠在生產為原來的f時，湖水才能在%d年內凈化。functionwenti3_3(l,year,bi,b)%格式:wenti3_3(P1,P2,P3,P4)%P1:y,M,d,h,Min,s表示設定步長值%P2:year%P3:百分比%P4:xx-線性自然凈化,zs-指數(shù)自然凈化formatcompactformatlongifl'y'q=1elseifl'm&

30、#39;q=12elseifl'd'q=365;elseifl'h'q=365*24elseifl'min'q=365*24*60elseifl's'q=365*24*60*60elseerror('error,tryagain!')endA=zeros(1,year);A(1)=;N=0;fori=1:1:yeartemp=*bi;forj=1:1:fix(q)A(i)=(A(i)*1015+temp*10"14/q)/(10"15+*10"14/q)-fen(A(i),l,b);e

31、ndA(i+1)=(A(i)*1015+temp*10"14/q)/(10"15+*10"14/q)-fen(A(i),l,b);endN=A(i)被調用的程序：functiony=fen(N,l,b)%自然凈化%格式：fen(P1,P2)%P1：前一時間點的污染值%P2:y,m,d,按年或月或日%線性關系ifb'xx'first=;%first=;ifl'y'y=first*N;elseifl'm'y=first*N/12;elseifl'd'y=first*N/365;elseerror('error,wrongparameter.tryagain!