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

2、體數值，運用時間步長法來仿真模擬了在湖水污染濃度穩定以前湖水每天的變化情況，輸出自PCA建廠以來每年的湖水污染濃度，得到與微分方程模型相同的結論。 在全停產和半停產時，通過前面的兩個模型可以計算湖水污染濃度在自然凈化影響下的恢復到凈化指標所需的年限。并可得到結論：在半停產狀態下，在選定的自然凈化速率常數的約束下，只有當河水污染濃度降至原來的3.15%（自然凈化速率呈線性關系），4.7%（自然凈化速率呈指數關系），才有可能使河水在100年內恢復至0.001mol/l，然后給出整改建議。1、 問 題 重 述 Pure河是流入Pristine湖的唯一河流。50年前PCA公司在此河旁建起一個生產設施并

3、投入運行。PCA將為處理的湖水排入河中，導致Pristine湖被污染。PCA公司聲稱：已排放的廢水的標準多年從未改變切不會對湖的環境有影響。現已知：Pristine湖的湖容量為L，流入（流出）的水流速度為L/年。PCA公司聲稱河水污染濃度僅為0.001mol/L，自工廠以來沒有改變過。討論下列問題：（1） 建立數學模型用PCA提供的公開數據判斷湖的環境是否會惡化；（2） 以目前湖水污染濃度0.03mol/L，和河水污染濃度0.05mol/L為新數據判斷湖的環境是否會惡化；2、 模型的合理假設和符號系統2.1 模型的合理假設（1）降水量和增發量相等；（2）湖中流入量和流出量相等且一直未變；（3）

4、污水量遠小于河水注入量，且污水與河水混合均勻；（4）湖水混合均勻，且流入污水的擴散速度無限大；（5）湖內除Pure河外，無其他污染源；2.2 符號系統 ：河水污染濃度mol/L； ：湖水污染物濃度mol/L； ：湖泊容量L； ：自然凈化速率mol/（L。年） ：流入（流出）的水流速度L/年； ：從PCA建廠至考察時刻的時間段。 三、問題的分析3.1 問題分析： 對于問題中幾個詞語的理解：1. 是否會惡化湖的環境惡化即湖水污染濃度大于0.001mol/L，要判斷其是否會惡化，則需計算在某一污染物積累速度（分析影響此速度的因素）下，湖水能達到的最大污染濃度和其變化趨勢，以及湖水經幾年超過0.001

5、mol/L，經過幾年達到最大污染濃度。2. 自然凈化自然凈化是獨立的生態系統進行自我調節的方式之一，是在空氣，陽光，水和細菌的參與下，進行包括物理沉降，化學反應和生物轉化三大方面的活動，其最終作用是將污染物轉化為無害物質，從而凈化生態系統。當河水輸入湖泊并均勻混合之后，影響湖水污染濃度的唯一因素便是自然凈化速度。湖水污染問題水流的動態流程圖：此問題中，我們考察對象是湖水污染濃度的變化趨勢：1. 在整改之前，其增加的趨勢，超過凈化指標0.001ml/L（即湖水惡化）的可能性和時限；2. 在全停產（無污染物輸入）和半停產的情況下，其降低的趨勢，達到凈化指標的時限。3. 在整改之后，其增加的趨勢，未

6、定與凈化指標之下某一水平的時限、在前假設條件的基礎之上，湖水容量不變，出河水外無其他的污染源，故我們可將湖泊作為一個封閉的生態系統，其簡化的湖水被污染的動態過程為：受污河輸入湖泊，河水與湖水均勻混合，受污河水進行自我的凈化，湖水數出湖泊。湖中污染物的量直接決定了湖水污染濃度，而污染物的量受到以下兩方面因素的影響：1.河水的污染濃度及其流入速度（根據已知此速度不變），2.湖水的自然凈化速度，前者使其增加，后者使其減少（負增加）。問題一、二、三的實質都是要分析污染濃度的變化趨勢，其去表便在于前一因素的不同。問題一中，河水污染濃度不變，恒為0.001mol/L；問題二中，河水污染濃度可能會變化，受P

7、CA效益的影響而按一定規律波動；問題三中，在全停產或半停產的情況下，和碩污染濃度為0或減為問題二中的一部分。后一因素（自然凈化速度）在三個問題中的作用都是相同的。根據微積分的知識可知，在適當短的時間段之內，通過建立微分方程，可以將連續的過程離散化，從而可以得到湖水污染濃度與時間之間的關系式。利用時間步長發，縮小步長值（從年到月到天），并與微分方程所得的精確解做出比較。四、模型建立與求解問題一：根據PCA的公開申明和所提供數據，可認為：河水污染濃度恒為0.001mol/L。從存在自然凈化和不存在自然凈化兩個方面考慮：（1） 在不考慮自然凈化的情況下：由于假設湖中流入量和流出量相等，而在經過與湖水

8、均勻混合后，流出湖水污染濃度明顯減小，故流出污染物的量小于流出污染物的量，污染物將在湖中沉積，從而使湖水污染濃度增加，當其增加至于輸入的河水污染濃度相等時，河水污染濃度達到最大，并穩定在這一數值，不在增加。建立湖水污染濃度隨時間變化的微分方程模型：設在極短時間，湖水污染濃度增加，在將湖水被污染這一連續動態過程簡化為離散的瞬間靜止狀態（如問題分析中所述）之后，根據湖中剩余量=輸入量輸出量，我們可以列微分方程如下：化簡可得： ······（1）帶入L，L/年和mol/L的數據，我們可得到和的關系式如下： 通過此關系式我們可知，當時，湖水污

9、染濃度將趨近與0.001mol/L，即湖的環境不會惡化。建立計算機模擬模型：湖水污染濃度的變化時有湖中污染物隨時間的積累而引起的，這個逐步積累的過程我們可以用計算機進行仿真模擬，其實質為完成一個循環累加的過程，并可改變時間步長，如一年一年的累積，一月一月的累積，一天一天的累積，從而使我們的模擬值逐步精確，可與微分方程求得精確解比較，分析誤差。為提高模擬結果的精確性和運算的效率，我們采用了逐天累加，數出年污染濃度的方式。模擬程序見附件一，從后面的分析中我們可知：在此河水污染濃度恒定，無自然凈化的最簡單的情況下建立的模型是以后問題的基礎，后面的問題只是改變條件或數據，其實質是不變的，故我們在此程序

10、中加入了多個選擇語句，在不同德條件或數據下執行不同的命令，從而用一個程序解決全部的模擬問題。模擬所得的數據如下表1：年份135810 1.734.356.147.828.51年份1520304050 9.439.789.979.9959.999根據模擬數據所作的湖水污染濃度變化趨勢圖如下圖1：模擬所得的數據顯示：湖水污染 濃度將穩定于0.001mol/L，這與分析微分方程所得的結果是相符的。（2） 在考慮自然凈化的情況下：A. 首先，我們應該了解什么事自然凈化。根據資料顯示：湖水中的污染物可分為有機污染物和無機污染物兩大類，在多種環境因素（陽光，空氣，水，水中生物，水中化學物質，重力等）的作用

11、下，通過物理沉降，化學反應和生物轉化一系列復雜的活動，它們的量會發生變化。有機物在水環境中的遷移轉化過程如圖2：可以看出整個過程是相當復雜的，不僅過程多，而且在相同的過程中，不同的物質有著不同的結果。此問題中沒有明確給出輸入湖中的污染物種類，也沒有對湖泊環境做任何描述，無疑給問題的解決增加了極大的難度。為是問題簡化，我們從一般的情況出發，假設：湖中所進行的反應均為一級反應，有機物的存在不會對環境參數造成改變，環境固定（自然凈化速率恒定），考慮湖泊生態系統中起主要作用的幾種過程作簡要分析：1. 物理沉降。不同物質有著不同的沉降速度常數，沉降速率為；2. 揮發。不同物質有著不同的揮發速率常數，在有

12、機物在水體上的大氣中的分壓為0的條件下，揮發速率為(為水體深度)3. 水解反應。不同物質有著不同的水解速率常數，在一級反應的條件下，水解速率為；4. 生物降解反應。不同物質有著不同的降解速率常數，在一級反應的條件下，生物降解速率；B.在最簡模擬程序的基礎上，從每天的積累量中減去每天的自然凈化量時間步長，重復循環可得到逐年湖水污染濃度的值，設。越大，湖水污染濃度將越快穩定于一個更小的濃度值。模擬所得的數據如下表2：年份135810 1.5873.4374.3154.8784.984年份1520304050 5.0875.1035.1065.1065.106根據模擬數據所作的湖水污染濃度變化趨勢圖

13、如下圖3：C.進一步考慮自然凈化速率的影響。更貼近于實際的情況是，自然凈化速率與湖水污染濃度成指數關系，隨著的增加而增加，但增加的速率會逐漸減小，用關系式表達即可為：，A,B為與環境和污染物種類有關的常數。將此關系式帶入微分方程（2），得到一個新的微分方程，此方程無解。但是我們可以通過計算機模擬求得數值解。A的大小影響著最終穩定濃度，B的大小影響著達到最終穩定濃度的快慢。在前面達到線性關系的基礎上，在=0/。001mol/L處的值大小確定A至數量級，并穩定年限定在1020年間，從而得到一組估計值：模擬所得數據如下表3：年份135810 1.4932.9713.5393.8163.871年份15

14、20304050 3.9033.9063.9063.9063.906根據模擬所得數據所作的湖水污染濃度變化趨勢圖如下圖4：觀察數據可知：無論自然凈化速率與污染濃度成線性關系還是指數關系，當河水污染濃度=0.001mol/L時，都將穩定于一個小于0.001mol/L的值，也即：湖的環境不會惡化。問題二：根據實際情況，在此我們只考慮存在自然凈化的情況。（1） .河水污染濃度恒定為0.05mol/L:與問題一（2）相同，只是=0.05mol/L，沿用問題一（2）的微分方程模型和差分模擬模型，我們可以得到以下結果：A. 微分方程模型（自然凈化速率與湖水污染濃度成線性關系）：當時間時，即湖水污染濃度穩定

15、于一個與K有關的值。當=0.05mol/L時，如要穩定于0.001mol/L則K值為。根據我們設定的可計算，湖水最終污染濃度=0.0317mol/L，超過凈化指標0.001mol/L，故在此條件下，湖的環境將會惡化。B.差分模擬模型：在與成線性關系時，得到模擬數據如下表4：年份135810 0.821.882.462.883.01年份1520304050 3.133.163.173.173.17在與成指數關系時，得到模擬數據如下表5：年份135810 0.812.012.813.594.02年份1520304050 4.324.484.564.584.58觀察以上數據可知：由于河水污染沒得過大

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

17、表6：年份151020223540 0.1850.8011.7927.2729.58152.0497.13年份508290100112140150 260.7625.0633.9617.2626.3628.5632.5根據模擬所得數據所作的湖水污染濃度變化趨勢如下圖5：在于成指數關系時，得到的模擬數據如下表7：年份151020243545 0.1690.5861.2505.1949.04644.2695.75年份5080100120130140150 321.7925.5959.6950.6949.6952.5964.8根據模擬所得數據所作的湖水污染濃度變化趨勢如下圖5：在第50年時，無論自然

18、凈化速率于湖水污染濃度成線性關系還是指數關系，可以看到湖水污染濃度都接近工作人員實地測得值，這從檢驗的角度說明我們的對模型參數的估計是可取的。觀察模擬所得數據可知，湖的環境將會惡化。 五、模型的評價及評改進5.1模型的評價（1） .對于首先建立的含參數的微分方程模型，在其有精確解的時候，能給出湖水污染濃度隨時間變化的規律，得到很好的預測效果，但微分方程并不是任何條件下都有解，在本問題中，只有當河水污染濃度恒定且自然凈化速率呈線性關系時，微分方程可解，故適用的范圍受到了限制。（2） .差分模擬法，即用數值方法求解微分方程，它解決了微分方程無精確解的 情況，適用范圍很大。使用一個兼顧運算效率和模擬

19、結果精度的時間步長是該方 法的關鍵。在處理本問題時，采用的時間步長為天。這樣處理的精度大約為0.0001, 與精確解相當接近。5.2模型的改進（1）.自然凈化速率的改進：資料顯示，污染物的很多生物降解過程是通過好氧 反應實現的，當污染濃度較大時，需氧量增加，而湖中氧氣總量不變，自然凈化 速率會有所減小，所以，自然凈化速率的函數可作為分段函數考慮。前一段還是用前面的指數模型，后面一段可以考慮為。模型這樣改進后，更加符合實際情況。（2）.河水濃度變化規律模型的改進：根據價值規律，一個工廠的發展總是有起 有落，其生產量繞某一中心值上下波動，故生產河水的濃度的變化應該是具有一 定周期性，同時又有所增加

20、。可以考慮在原來的Logistic模型后加上一個傅立葉 函數（其周期就是工廠生產的周期）和一個正態隨機變量。6、 參考文獻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，考慮自然凈

21、化，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：采用差分方程模型 命令：wenti

22、12('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'

23、)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',

24、9;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年,河水污染變化,不考慮自然

25、降解 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

26、('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; e

27、nd for j=1:1:q if w 'n' var=0; elseif w 'y' var=fen(A(i),l,b); end A(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; end if w 'n' var=0; elseif w 'y' var=fen(A(i),l,b); end A(i+1)=(A(i)*1015+var1*1.9*1014/q)/(1015+1.9*1014/q

28、)-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

29、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)*1015+temp*1.9*1014/q)

30、/(1015+1.9*1014/q)-fen(A(i),l,b); t=A(i); if A(i)<0.001 break; end end A(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; end end if i year p=0; else p=1; end %程序調試所用 %sprintf('停產') sprintf('需要 %d 年%d天',i,j) %sprintf('第%d年%d天為:%.12f

31、9;,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=(lo

32、w+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 els

33、eif 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)*1015+temp*1.9*1014/q)/(1015+1.9*1014/q)-f

34、en(A(i),l,b); end A(i+1)=(A(i)*1015+temp*1.9*1014/q)/(1015+1.9*1014/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&