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

1、中國傳媒大學(xué)2010 學(xué)年第一學(xué)期數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)課程數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)題目Pristine湖污染問題的建模與求解學(xué)生姓名學(xué)號(hào)班級(jí)學(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í)間步長法來仿真模擬了在湖水污

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

3、 Pristine湖被污染。PCA公司聲稱：已排放的廢水的標(biāo)準(zhǔn)多年從未改變切不會(huì)對湖的環(huán)境有影響。1514現(xiàn)已知：Pristine湖的湖容量為10 L,流入（流出）的水流速度為9 L/年。PCA公司聲稱河水污染濃度僅為0.001mol/L ,自工廠以來沒有改變過。討論下列問題：（1）建立數(shù)學(xué)模型用PCA提供的公開數(shù)據(jù)判斷湖的環(huán)境是否會(huì)惡化;（2）以目前湖水污染濃度0.03mol/L,和河水污染濃度0.05mol/L為新數(shù)據(jù)判斷湖的環(huán)境是否會(huì)惡化；二、模型的合理假設(shè)和符號(hào)系統(tǒng)2.1 模型的合理假設(shè)（ 1）降水量和增發(fā)量相等；（ 2）湖中流入量和流出量相等且一直未變；（ 3）污水量遠(yuǎn)小于河水注入量

4、，且污水與河水混合均勻；（ 4）湖水混合均勻，且流入污水的擴(kuò)散速度無限大；（ 5）湖內(nèi)除Pure河外，無其他污染源；2.2 符號(hào)系統(tǒng)0 ：河水污染濃度mol/L ；：湖水污染物濃度mol/L ；15V ：湖泊容量10 L；c:自然凈化速率 mol/ （L。年）14：流入（流出）的水流速度1.9 L/ 年；t ：從 PCA 建廠至考察時(shí)刻的時(shí)間段。三、問題的分析3.1 問題分析：對于問題中幾個(gè)詞語的理解：1 .是否會(huì)惡化湖的環(huán)境惡化即湖水污染濃度大于0.001mol/L ,要判斷其是否會(huì)惡化，則需計(jì)算在某一污染物積累速度（分析影響此速度的因素）下，湖水能達(dá)到的最大污染濃度和其變化趨勢，以及湖水經(jīng)

5、幾年超過0.001mol/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)化為無害物質(zhì)， 從而凈化生態(tài)系統(tǒng)。當(dāng)河水輸入湖泊并均勻混合之后，影響湖水污染濃度的唯一因素便是自然凈化速度。湖水污染問題水流的動(dòng)態(tài)流程圖：PCAPure河 輸入 PrisHne湖 輸出此問題中，我們考察對象是湖水污染濃度的變化趨勢：1 .在整改之前，其增加的趨勢，超過凈化指標(biāo)0.001ml/L （即湖水惡化）的可能性和時(shí)限；2 .在全停產(chǎn)（無污染物輸入）和半停產(chǎn)的情

6、況下，具降低的趨勢，達(dá)到凈化指標(biāo) 的時(shí)限。3 .在整改之后，其增加的趨勢，未定與凈化指標(biāo)之下某一水平的時(shí)限、在前假設(shè)條件的基礎(chǔ)之上，湖水容量不變，出河水外無其他的污染源，故我 們可將湖泊作為一個(gè)封閉的生態(tài)系統(tǒng)， 其簡化的湖水被污染的動(dòng)態(tài)過程為：受污 河輸入湖泊，河水與湖水均勻混合，受污河水進(jìn)行自我的凈化，湖水?dāng)?shù)出湖泊。湖中污染物的量直接決定了湖水污染濃度， 而污染物的量受到以下兩方面因素的 影響：1.河水的污染濃度及其流入速度（根據(jù)已知此速度不變），2.湖水的自然凈化速度，前者使其增加，后者使其減少（負(fù)增加）。問題一、二、三的實(shí)質(zhì)都是要分析污染濃度的變化趨勢，其去表便在于前一因素的不同。問題一

7、中，河水污染濃度不變，恒為0.001mol/L;問題二中，河水污染濃度可能會(huì)變化，受 PCA效益的影響而按一定規(guī)律波動(dòng)；問題三中，在全停產(chǎn)或半停產(chǎn)的情況下，和碩污染濃度為0 或減為問題二中的一部分。后一因素 （自然凈化速度）在三個(gè)問題中的作用都是相同的。根據(jù)微積分的知識(shí)可知，在適當(dāng)短的時(shí)間段之內(nèi)，通過建立微分方程，可以將連續(xù)的過程離散化，從而可以得到湖水污染濃度與時(shí)間之間的關(guān)系式。利用時(shí)間步長發(fā)，縮小步長值（從年到月到天），并與微分方程所得的精確解做出比較。四、模型建立與求解問題一：根據(jù) PCA 的公開申明和所提供數(shù)據(jù)，可認(rèn)為：河水污染濃度恒為0.001mol/L。從存在自然凈化和不存在自然凈

8、化兩個(gè)方面考慮：（ 1）在不考慮自然凈化的情況下：由于假設(shè)湖中流入量和流出量相等，而在經(jīng)過與湖水均勻混合后，流出湖水污染濃度明顯減小，故流出污染物的量小于流出污染物的量，污染物將在湖中沉積，從而使湖水污染濃度增加，當(dāng)其增加至于輸入的河水污染濃度相等時(shí)，河水污染濃度達(dá)到最大，并穩(wěn)定在這一數(shù)值，不在增加。建立湖水污染濃度隨時(shí)間變化的微分方程模型：設(shè)在極短時(shí)間dt ，湖水污染濃度增加d ，在將湖水被污染這一連續(xù)動(dòng)態(tài)過程簡化為離散的瞬間靜止?fàn)顟B(tài)（如問題分析中所述）之后，根據(jù)湖中剩余量= 輸入量輸出量，我們可以列微分方程如下：dt化簡可得：dr V ( 0)(1)帶入V 1015L,1.9 1014L/

9、年和0 0.001 mol/L的數(shù)據(jù)，我們可得到和t的0.19 t 0.19 t關(guān)系式如下:0.001 e (e 1)通過此關(guān)系式我們可知，當(dāng)t 時(shí)，湖水污染濃度將趨近與 0.001mol/L,即湖的環(huán)境不會(huì)惡化。建立計(jì)算機(jī)模擬模型：湖水污染濃度的變化時(shí)有湖中污染物隨時(shí)間的積累而引起 的，這個(gè)逐步積累的過程我們可以用計(jì)算機(jī)進(jìn)行仿真模擬， 其實(shí)質(zhì)為完成一個(gè)循 環(huán)累加的過程，并可改變時(shí)間步長，如一年一年的累積，一月一月的累積，一天 一天的累積，從而使我們的模擬值逐步精確，可與微分方程求得精確解比較，分 析誤差。為提高模擬結(jié)果的精確性和運(yùn)算的效率， 我們采用了逐天累加，數(shù)出年 污染濃度的方式。模擬程

10、序見附件一，從后面的分析中我們可知：在此河水污染濃度恒定，無自然 凈化的最簡單的情況下建立的模型是以后問題的基礎(chǔ)，后面的問題只是改變條件 或數(shù)據(jù)，其實(shí)質(zhì)是不變的，故我們在此程序中加入了多個(gè)選擇語句， 在不同德條 件或數(shù)據(jù)下執(zhí)行不同的命令，從而用一個(gè)程序解決全部的模擬問題。模擬所得的數(shù)據(jù)如下表1：年份13581010 41.734.356.147.828.51年份152030405010 49.439.789.979.9959.999根據(jù)模擬數(shù)據(jù)所作的湖水污染濃度變化趨勢圖如下圖1:模擬所得的數(shù)據(jù)顯示：湖水污染 濃度將穩(wěn)定于0.001mol/L,這與分析微分方程所得的結(jié)果是相符的。(2)在考慮自

11、然凈化的情況下：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:化還at作用可以看出整個(gè)過程是相當(dāng)復(fù)雜的， 不僅過程多，而且在相同的過程中，不同的物 質(zhì)有著不同的結(jié)果。此問題中沒有明確給出輸入湖中的污染物種類， 也沒有對湖 泊環(huán)境做任何描述，無疑給問題的解決增加了極大的難度。 為是問題簡化，我們從一般的情況出發(fā)，假設(shè)：湖中所進(jìn)行的反應(yīng)均為一級(jí)反應(yīng)，有機(jī)物的存在不

12、會(huì) 對環(huán)境參數(shù)造成改變，環(huán)境固定（自然凈化速率恒定），考慮湖泊生態(tài)系統(tǒng)中起 主要作用的幾種過程作簡要分析：1 .物理沉降。不同物質(zhì)有著不同的沉降速度常數(shù)，沉降速率為c ;K 2 .揮發(fā)。不同物質(zhì)有著不同的揮發(fā)速率常數(shù)v ,在有機(jī)物在水體上的大氣中的c Kv_分壓為0的條件下，揮發(fā)速率為C Z （Z為水體深度） 3水解反應(yīng)。不同物質(zhì)有著不同的水解速率常數(shù) Kp，在一級(jí)反應(yīng)的條件下，水解速率為c Kp4.生物降解反應(yīng)。不同物質(zhì)有著不同的降解速率常數(shù) Kb,在一級(jí)反應(yīng)的條件下,生物降解速率c KbB.在最簡模擬程序的基礎(chǔ)上，從每天的積累量中減去每天的自然凈化量K時(shí)間一一1/ C CCCCJH 1步長

13、，重復(fù)循環(huán)可得到逐年湖水污染濃度的值，設(shè) K 0.0003天。K越大，湖水污染濃度將越快穩(wěn)定于一個(gè)更小的濃度值。模擬所得的數(shù)據(jù)如下表2:年份13581010 41.5873.4374.3154.8784.984年份152030405010 45.0875.1035.1065.1065.106根據(jù)模擬數(shù)據(jù)所作的湖水污染濃度變化趨勢圖如下圖3:C.進(jìn)一步考慮自然凈化速率的影響。更貼近于實(shí)際的情況是，自然凈化速率 c與湖水污染濃度 成指數(shù)關(guān)系，c隨著 的增加而增加，但增加的速率會(huì)逐漸減小，B用關(guān)系式表達(dá)即可為：c A（1 e入A,B為與環(huán)境和污染物種類有關(guān)的常數(shù)。將此關(guān)系式帶入微分方程（2）,得到一

14、個(gè)新的微分方程，此方程無解。但是我們可以通過計(jì)算機(jī)模擬求得數(shù)值解。A的大小影響著 最終穩(wěn)定濃度，B的大小影響著 達(dá)到最終穩(wěn)定濃 度的快慢。在前面達(dá)到線 性關(guān)系的基礎(chǔ)上，在 =0/。001mol/L處的c值大小確定A至數(shù)量級(jí)，并穩(wěn)定年限定在1A20年間，從而得到一組估計(jì)值:A 8 10 4,B 400模擬所得數(shù)據(jù)如下表3:年份13581010 41.4932.9713.5393.8163.871年份152030405010 43.9033.9063.9063.9063.906根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢圖如下圖4:觀察數(shù)據(jù)可知：無論自然凈化速率c與污染濃度 成線性關(guān)系還是指數(shù)關(guān)系，

15、當(dāng)河水污染濃度 0=0.001mol/L時(shí)， 都將穩(wěn)定于一個(gè)小于 0.001mol/L的值，也即: 湖的環(huán)境不會(huì)惡化。問題二：根據(jù)實(shí)際情況，在此我們只考慮存在自然凈化的情況。（1） .河水污染濃度恒定為0.05mol/L:與問題一（2）相同，只是0=0.05mol/L,沿用問題一（2）的微分方程模型和差分模擬模型，我們可以得到以下結(jié)果：A.微分方程模型（自然凈化速率c與湖水污染濃度成線性關(guān)系）：當(dāng)時(shí)間t 時(shí)，KV ,即湖水污染濃度穩(wěn)定于一個(gè)與K有關(guān)的值。當(dāng)0 =0.05mol/L時(shí)，如要 穩(wěn)定于0.001mol/L則K值為分31年1。根據(jù)我們設(shè)定 的K 0.0003天1可計(jì)算，湖水最終污染濃度

16、 0 =0.0317mol/L ,超過凈化指標(biāo) 0.001mol/L,故在此條件下，湖的環(huán)境將會(huì)惡化。B.差分模擬模型：在c與 成線性關(guān)系時(shí)，得到模擬數(shù)據(jù)如下表 4:年份13581010 20.821.882.462.883.01年份152030405010 23.133.163.173.173.17在c與 成指數(shù)關(guān)系時(shí)，得到模擬數(shù)據(jù)如下表 5:年份13581010 20.812.012.813.594.02年份152030405010 24.324.484.564.584.58觀察以上數(shù)據(jù)可知：由于河水污染沒得過大，河水污染沒得從第一年起就超過了凈化指標(biāo)，并逐年增加，是湖的環(huán)境惡化。（2）

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

18、40.1850.8011.7927.2729.58152.0497.13年份50829010011214015010 4260.7625.0633.9617.2626.3628.5632.5根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢如下圖5:在c于 成指數(shù)關(guān)系時(shí)，得到的模擬數(shù)據(jù)如下表 7:年份15102024354510 40.1690.5861.2505.1949.04644.2695.75年份508010012013014015010 4321.7925.5959.6950.6949.6952.5964.8根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢如下圖5:I司題上(2)河水濃度變化, 指

19、數(shù)自然洋化0 10.090.080 070 06QQ50 C40 030 020 01.*9M1!1 1*.A*工 子 木 » «* . ! 皋 0* * :J1 « v*1/T , !' *4A1*»1AA* W*a»a11* 9,f /*»圈*vi.i.:*'T! 1 1 , « Hi 1*a*3qa. :14-tiA+9 !T!*d*4i I11 9 Bl¥ JdpIliA口 _ 上1 -iiiiiiii0204080 BO 100120140180在第50年時(shí)，無論自然凈化速率c于湖水污染濃

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

21、是該方 法的關(guān)鍵。在處理本問題時(shí)，采用的時(shí)間步長為大。這樣處理的精度大約為0.0001,與精確解相當(dāng)接近。5.2 模型的改進(jìn)(1) .自然凈化速率的改進(jìn)：資料顯示，污染物的很多生物降解過程是通過好氧反應(yīng)實(shí)現(xiàn)的，當(dāng)污染濃度較大時(shí)，需氧量增加，而湖中氧氣總量不變，自然凈化速率會(huì)有所減小，所以， 自然凈化速率的函數(shù)可作為分段函數(shù)考慮。前一段還是用前面的指數(shù)模型，后面一段可以考慮為C e D E（C, D,E 0）。模型這樣改進(jìn)后，更加符合實(shí)際情況。（ 2） .河水濃度變化規(guī)律模型的改進(jìn)：根據(jù)價(jià)值規(guī)律，一個(gè)工廠的發(fā)展總是有起有落， 其生產(chǎn)量繞某一中心值上下波動(dòng)，故生產(chǎn)河水的濃度的變化應(yīng)該是具有一定周期

22、性，同時(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月；七、附錄附錄一 :問題一：（1）河水濃度恒為0.001mol/L,不考慮自然凈化，a：采用差分方程模型命令：wenti12（ d ， 0.50,0.001， n）b:采用微分方程模型命令：t=1:1:51；P=0.001.*exp(-0.19*t)Pplot(t,p)grid on(2)河水濃度恒為0.001mol/L,考

23、慮自然凈化，a：采用差分自然凈化線性自然凈化命令： wenti12('d',0,50,0.001,'y','xx')指數(shù)自然凈化命令：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); phold on plot(t,p) grid on問題二：(1)河水濃度恒為0.05mol/L,不考慮自然凈化a：采用差分方程模型命令： wenti12(&

24、#39;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 onplot(t,p)grid on( 2)河水濃度恒為0.05mol/L ，考慮自然凈化，a: 采用差分方程模型線性自然凈化命令：wenti12('d',0.50,0.05,'y','xx')指數(shù)自然凈化命令：wneti12('d',0.50,0.05,'y','zs')b:采用微分方程模型

25、，線性自然凈化命令t=1:1:51;p=0.00166945.*exp(-0.2995.*t).*(-19+19.*exp(0.2995.*t);phold onplot(t,p)grid on( 3)河水濃度變化，不考慮自然凈化命令： wenti12('d',0.50,-1,'n')考慮自然凈化線性自然凈化命令：wenti12('d',0,50,-1,'y','xx')指數(shù)自然凈化命令：wenti12('d',0,50,-1,'y','zs')附錄二:主程序：fun

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

27、'elseif nargin 3 %河水污染變化,不考慮自然降解s=-1;w='n'b='xx'endq=1 elseif l 'm'q=12elseif l 'd'q=365elseif l 'h'q=365*24elseif l 'min'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,wrong parameter.try again!')endA=zeros(1,n+1);B=zeros(1,

28、n+1);A(1)=chu;B(1)=chu;var=0; %自然降解var1=0; %河水污染if w 'nvar=0;elseif w 'y'var=fen(A(1),l,b);elseerror('error,wrong parameter.try again!')endfor i=1:1:nif s -1var1=river(i);elsevar1=s;endfor j=1:1:qif w 'n'var=0;elseif w 'y'var=fen(A(i),l,b);endA(i)=(A(i)*10A15+vaU*

29、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=fen(A(i),l,b);endA(i+1)=(A(i)*10A15+vaT1.9*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 onfunction p=went

30、i3(l,bi,year,b,a)%格式: wenti3(P1,P2,P3,P，4 P5)%P1:y,m,d,h,min,s 表示設(shè)定步長值%P2:百分比%P3:可認(rèn)為的最大凈化年限%P4:xx線性自然凈化，zs-m數(shù)自然凈化%P5:湖水濃度初始值format long if l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365;elseif l 'h'q=365*24elseif l 'min'q=365*24*60elseif l 's'q=365*24*60*60

31、elseerror('error,try again!')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)/(10A15+1.9*10A14/q)-fen(A(i),l,b);t=A(i);if A(i)<0.001 break;end endA(i+1)=(A(i)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fen(A(i),l,b );if A(i)<

32、0.001 break;endend if i yearp=0;elsep=1;end%程序調(diào)試所用%sprintf('停產(chǎn)')sprintf('需要 %d 年天'浦)%sprintf('第 年%dK為:.12f,i,j-1,t)%sprintf('第年認(rèn)為:.12f,i,j,A(i)function low=wenti3_2(b,year,l,h)%采用二分法求出百分比format compactformat longlow=l;high=h;cen=(low+high)/2cha=high-low;while cha>10A-6%誤差判

33、斷條件m=wenti3('d',cen,year,b,0.05); %調(diào)用函數(shù)if m 0high=cen;cen=(low+high)/2elselow=cen;cen=(low+high)/2endcha=high-low;end high;low;sprintf('工廠在生產(chǎn)為原來的f時(shí)，湖水才能在d年內(nèi)凈化。,low,year)function wenti3_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-)t數(shù)自然凈

34、化format compactformat longif l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365;elseif l 'hq=365*24elseif l 'min'q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,try again!')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(

35、i)=(A(i)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fen(A(i),l,b);endA(i+1)=(A(i)*10A15+temp*1.9*10A14/q)/(10A15+1.9*10A14/q)-fen(A(i),l,b );N=A(i)被調(diào)用的程序：function y=fen(N,l,b)%自然凈化%格式:fen(P1,P2)%P1:前一時(shí)間點(diǎn)的污染值%P2:y,m,眼年或月或日%線性關(guān)系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;elseerror('error,wrong parameter.tr