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

1、中國(guó)傳媒大學(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é)院成績(jī)Pristine湖污染問題的建模與求解摘要本文討論了湖水污染濃度變化趨勢(shì)的預(yù)測(cè)問題。通過分析水流輸入輸出湖泊的過程，建立了湖水污染濃度隨時(shí)間變化的含參變量的微分方程模型，在河水污染濃度恒定和自然凈化速率呈線性關(guān)系的情況下，求得其精確解，帶入具體數(shù)據(jù)得到結(jié)論：在 PCA聲稱的河水污染濃度下，湖的環(huán)境不會(huì)惡化；在工作人員實(shí)地測(cè)得的河水濃度下，湖的環(huán)境將會(huì)惡化。同時(shí)建立了計(jì)算機(jī)模擬模型，帶入具體數(shù)值，運(yùn)用時(shí)間步長(zhǎng)法來仿真模擬了在湖水

2、污染濃度穩(wěn)定以前湖水每天的變化情況，輸出自 PCA建廠以來每年的湖水污染濃度，得到與微分方程模型相同的結(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將為處理的湖

3、水排入河中，導(dǎo)致Pristine湖被污染。 PCA公司聲稱：已排放的廢水的標(biāo)準(zhǔn)多年從未改變切不會(huì)對(duì)湖的環(huán)境有影響?，F(xiàn)已知： Pristine湖的湖容量為 1015 L，流入（流出）的水流速度為1.914 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）湖中流入量和流出量相等且一直未變；

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

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

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

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

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

9、）帶入 V1015L，1.91014L/ 年和00.001mol/L 的數(shù)據(jù)，我們可得到和 t 的關(guān)系式如下：0.001 e0.19 t(e0.19 t1)通過此關(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í)間步長(zhǎng)，如一年一年的累積，一月一月的累積，一天一天的累積，從而使我們的模擬值逐步精確，可與微分方程求得精確解比較，分析誤差。為提高模擬結(jié)果的精確性和運(yùn)算的效率，我們采用了逐天累加，數(shù)出年

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

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

12、造成改變，環(huán)境固定（自然凈化速率恒定），考慮湖泊生態(tài)系統(tǒng)中起主要作用的幾種過程作簡(jiǎn)要分析：1. 物理沉降。不同物質(zhì)有著不同的沉降速度常數(shù)，沉降速率為 c；2. 揮發(fā)。不同物質(zhì)有著不同的揮發(fā)速率常數(shù)Kv ，在有機(jī)物在水體上的大氣中的分壓為 0cK v的條件下，揮發(fā)速率為Z( Z 為水體深度)3. 水解反應(yīng)。不同物質(zhì)有著不同的水解速率常數(shù) K p ，在一級(jí)反應(yīng)的條件下，水解速率為c K p；4. 生物降解反應(yīng)。不同物質(zhì)有著不同的降解速率常數(shù)Kb ，在一級(jí)反應(yīng)的條件下，生物降解速率 cKb；B. 在最簡(jiǎn)模擬程序的基礎(chǔ)上，從每天的積累量中減去每天的自然凈化量K時(shí)間步長(zhǎng)，重復(fù)循環(huán)可得到逐年湖水污染濃度的

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

14、以通過計(jì)算機(jī)模擬求得數(shù)值解。A 的大小影響著最終穩(wěn)定濃度， B 的大小影響著達(dá)到最終穩(wěn)定濃度的快慢。在前面達(dá)到線 性關(guān)系的基礎(chǔ)上， 在=0/ 。001mol/L 處的 c 值大小確定 A 至數(shù)量級(jí)，并穩(wěn)定年限定在 1020 年間，從而得到一組估計(jì)值：模擬所得數(shù)據(jù)如下表 3：年份1358101.4932.9713.5393.8163.871年9033.9063.9063.9063.906根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢(shì)圖如下圖4：觀察數(shù)據(jù)可知：無論自然凈化速率 c與污染濃度成線性關(guān)系還是指數(shù)關(guān)系，當(dāng)河水污染濃度0 =0.001mol/L 時(shí)，都將穩(wěn)定于一個(gè)小于

15、0.001mol/L 的值，也即：湖的環(huán)境不會(huì)惡化。問題二：根據(jù)實(shí)際情況，在此我們只考慮存在自然凈化的情況。（1）. 河水污染濃度恒定為0.05mol/L:與問題一（ 2）相同，只是0 =0.05mol/L ，沿用問題一（ 2）的微分方程模型和差分模擬模型，我們可以得到以下結(jié)果：A. 微分方程模型（自然凈化速率c 與湖水污染濃度成線性關(guān)系）：0當(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 值為 9.31 年 1。根據(jù)我們?cè)O(shè)定的 K 0.0003天 1可計(jì)算，湖水最終污染濃度0 =0.0317mol/L

16、，超過凈化指標(biāo) 0.001mol/L ，故在此條件下， 湖的環(huán)境將會(huì)惡化。B. 差分模擬模型：在 c 與成線性關(guān)系時(shí)，得到模擬數(shù)據(jù)如下表4：年份1358100.821.882.462.883.01年133.163.173.173.17在 c 與成指數(shù)關(guān)系時(shí)，得到模擬數(shù)據(jù)如下表5：年份1358100.812.012.813.594.02年324.484.564.584.58觀察以上數(shù)據(jù)可知： 由于河水污染沒得過大， 河水污染沒得從第一年起就超過了凈化指標(biāo)，并逐年增加，是湖的環(huán)境惡化。（2）. 河水污染濃度變化：根據(jù)實(shí)際情況，一個(gè)工廠的生產(chǎn)量并非

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

18、13年份508290100112140150260.7625.0633.9617.2626.3628.5632.5根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢(shì)如下圖5：在 c 于成指數(shù)關(guān)系時(shí)，得到的模擬數(shù)據(jù)如下表7：年份1510202435450.1690.5861.2505.1949.04644.2695.75年份5080100120130140150321.7925.5959.6950.6949.6952.5964.8根據(jù)模擬所得數(shù)據(jù)所作的湖水污染濃度變化趨勢(shì)如下圖5：在第 50 年時(shí)，無論自然凈化速率c于湖水污染濃度成線性關(guān)系還是指數(shù)關(guān)系，可以看到湖水污染濃度都接近工作人員實(shí)地測(cè)得值，這從

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

20、001, 與精確解相當(dāng)接近。5.2模型的改進(jìn)（1）. 自然凈化速率的改進(jìn)：資料顯示，污染物的很多生物降解過程是通過好氧反應(yīng)實(shí)現(xiàn)的，當(dāng)污染濃度較大時(shí)，需氧量增加，而湖中氧氣總量不變，自然凈化速率會(huì)有所減小，所以，自然凈化速率的函數(shù)可作為分段函數(shù)考慮。前一段還是用前面的指數(shù)模型，后面一段可以考慮為Ce DE (C, D, E0)。模型這樣改進(jìn)后，更加符合實(shí)際情況。（2）. 河水濃度變化規(guī)律模型的改進(jìn)：根據(jù)價(jià)值規(guī)律，一個(gè)工廠的發(fā)展總是有起有落，其生產(chǎn)量繞某一中心值上下波動(dòng)，故生產(chǎn)河水的濃度的變化應(yīng)該是具有一定周期性，同時(shí)又有所增加?？梢钥紤]在原來的 Logistic 模型后加上一個(gè)傅立葉函數(shù)（其周期

21、就是工廠生產(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；Pplot(t,p)grid on（2）河水濃度恒為 0.001mol/L ，考慮自然凈化，a：采用差分自然凈化線性自然凈化命令： wenti12('d',0,5

22、0,0.001,'y','xx')指數(shù)自然凈化命令： wenti12('d',0,50,0.001,'y','zs')b：采用微分方程模型，線性自然凈化命令： t=1:1:51;pholdonplot(t,p)gridon問題二：（1）河水濃度恒為 0.05mol/L, 不考慮自然凈化a：采用差分方程模型命令： wenti12('d',0.50,0.05,'n')b: 采用微分方程模型命令：hold onplot(t,p)gridon（2）河水濃度恒為0.05mol/L ，考慮自然凈

23、化，a: 采用差分方程模型線性自然凈化命令： wenti12('d',0.50,0.05,'y','xx')指數(shù)自然凈化命令： wneti12('d',0.50,0.05,'y','zs')b：采用微分方程模型，線性自然凈化命令t=1:1:51;pholdonplot(t,p)gridon（3）河水濃度變化，不考慮自然凈化命令： wenti12('d',0.50,-1,'n')考慮自然凈化線性自然凈化命令： wenti12('d',0,50,-1,&#

24、39;y','xx')指數(shù)自然凈化命令： wenti12('d',0,50,-1,'y','zs')附錄二：主程序：function wenti12(l,chu,n,s,w,b)%按照時(shí)間步長(zhǎng)法，以每年或每月或每日為時(shí)間段，求出濃度的變化，最后輸出每年的濃度，做一個(gè)大致的觀察 ,%格式： wenti12 （P1,P2,P3,P4,P5,P6)%P1:y,m,d,h,min,s表示設(shè)定步長(zhǎng)值%P2：湖水污染初值%P3：年份%P4：河流污染初值%P5: 是否考慮自然降解 ,n- 不考慮 ,y- 考慮%P6: xx 采用線性凈化

25、模型， zs 采用指數(shù)凈化模型 format compactformat longif nargin 2 %50 年, 河水污染變化 , 不考慮自然降解 n=50;s=-1;w='n'b='xx'elseif nargin 3 % 河水污染變化 , 不考慮自然降解 s=-1;w='n'b='xx'endif l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365elseif l 'h'q=365*24elseif l 'min

26、9;q=365*24*60elseif l 's'q=365*24*60*60elseerror('error,wrong parameter.try again!')endA=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);elseerror('error,wrong parameter.try again!')endfor i=

27、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)*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;endif w 'n'var=0;elseif w 'y'var=fen(A(i),l,b);endA(i+1)=(A(i)*1015+var1*1.9*1

28、014/q)/(1015+1.9*1014/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=wenti3(l,bi,year,b,a)%格式 : wenti3(P1,P2,P3,P4，P5)%P1:y,m,d,h,min,s表示設(shè)定步長(zhǎng)值%P2:百分比%P3:可認(rèn)為的最大凈化年限%P4:xx-線性自然凈化 ,zs- 指數(shù)自然凈化%P5：湖水濃度初始值format compactformat longif l 'y'q=1

29、elseif 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*60elseerror('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)*1015+temp*1.9*1014/q)/(1015+1.9*1

30、014/q)-fen(A(i),l,b); t=A(i);if A(i)<0.001 break;endendA(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;endendif i yearp=0;elsep=1;end%程序調(diào)試所用%sprintf('停產(chǎn) ')sprintf('需要 %d 年%d天 ',i,j)%sprintf('第%d年%d天為 :%.12f',i,j-1,t)%sprintf('第

31、%d年%d天為 :%.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>10-6%誤差判斷條件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;endhigh;low;spri

32、ntf('工廠在生產(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è)定步長(zhǎng)值%P2:year%P3: 百分比%P4:xx-線性自然凈化 ,zs- 指數(shù)自然凈化format compactformat longif l 'y'q=1elseif l 'm'q=12elseif l 'd'q=365;elseif l 'h'q=365*24else

33、if 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(i)=(A(i)*1015+temp*1.9*1014/q)/(1015+1.9*1014/q)-fen(A(i),l,b); endA(i+1)=(A(i)*1015+temp*1.9*1014/q)/(1015+1.9*1014/q)-fen(A(i),l,b );endN=A(i)被調(diào)用的程序：function y=fen(N,l,b)%自然凈化%格式 :fen(P1,P2)%P1:前一時(shí)間點(diǎn)的污染值%P2:y,m,d, 按年或月或日%線性關(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.try again!')e