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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

18、，所以，自然凈化速率的函數(shù)可作為分段函數(shù)考慮。前一段還是用前面的指數(shù)模型，后面一段可以考慮為Cxe-Dxp+E（C,D,E>°）。模型這樣改進(jìn)后，更加符合實(shí)際情況。（2）.河水濃度變化規(guī)律模型的改進(jìn)：根據(jù)價(jià)值規(guī)律，一個(gè)工廠的發(fā)展總是有起有落，其生產(chǎn)量繞某一中心值上下波動(dòng)，故生產(chǎn)河水的濃度的變化應(yīng)該是具有一定周期性，同時(shí)又有所增加。可以考慮在原來的Logistic模型后加上一個(gè)傅立葉函數(shù)（其周期就是工廠生產(chǎn)的周期）和一個(gè)正態(tài)隨機(jī)變量。六、參考文獻(xiàn)1 數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)高等教育出版社，趙靜、但琦主編，2008年1月第三版2 環(huán)境化學(xué)，高等教育出版社出版，戴樹桂主編，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)%按照時(shí)間步長法，以每年或每月或每日為時(shí)間段，求出濃度的變化，最后輸出每年的濃度，做一個(gè)大致的觀察,%格式：wentil2(P1,P2,P3,P4,P5,P6)%Pl:y,m,d,h,min,s表示設(shè)定步長值%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表示設(shè)定步長值%P2:百分比%P3:可認(rèn)為的最大凈化年限%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%程序調(diào)試所用%sprin

28、tf('停產(chǎn)')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,;%調(diào)用函數(shù)ifm0hig

29、h=cen;cen=(low+high)/2elselow=cen;cen=(low+high)/2endcha=high-low;endhigh;low;',low,year)sprintf('工廠在生產(chǎn)為原來的f時(shí)，湖水才能在%d年內(nèi)凈化。functionwenti3_3(l,year,bi,b)%格式:wenti3_3(P1,P2,P3,P4)%P1:y,M,d,h,Min,s表示設(shè)定步長值%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)被調(diào)用的程序：functiony=fen(N,l,b)%自然凈化%格式：fen(P1,P2)%P1：前一時(shí)間點(diǎn)的污染值%P2:y,m,d,按年或月或日%線性關(guān)系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!