湖水污染問題的數(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è)問題。 通過(guò)分析水流輸入輸出湖泊的過(guò)程，建立了湖水污染濃度隨時(shí)間變化的含參變量的微分方程模型，在河水污染濃度恒定和自然凈化速率呈線性關(guān)系的情況下，求得其精確解，帶入具體數(shù)據(jù)得到結(jié)論：在PCA聲稱的河水污染濃度下，湖的環(huán)境不會(huì)惡化；在工作人員實(shí)地測(cè)得的河水濃度下，湖的環(huán)境將會(huì)惡化。 同時(shí)建立了計(jì)算機(jī)模擬模型，帶入具

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

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

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

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

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

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

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

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

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

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

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

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

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

15、于一個(gè)與K有關(guān)的值。當(dāng)=0.05mol/L時(shí)，如要穩(wěn)定于0.001mol/L則K值為。根據(jù)我們?cè)O(shè)定的可計(jì)算，湖水最終污染濃度=0.0317mol/L，超過(guò)凈化指標(biāo)0.001mol/L，故在此條件下，湖的環(huán)境將會(huì)惡化。B.差分模擬模型：在與成線性關(guān)系時(shí)，得到模擬數(shù)據(jù)如下表4：年份135810 0.821.882.462.883.01年份1520304050 3.133.163.173.173.17在與成指數(shù)關(guān)系時(shí)，得到模擬數(shù)據(jù)如下表5：年份135810 0.812.012.813.594.02年份1520304050 4.324.484.564.584.58觀察以上數(shù)據(jù)可知：由于河水污染沒得過(guò)大

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

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

18、凈化速率于湖水污染濃度成線性關(guān)系還是指數(shù)關(guān)系，可以看到湖水污染濃度都接近工作人員實(shí)地測(cè)得值，這從檢驗(yàn)的角度說(shuō)明我們的對(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ù)值方法求解微分方程，它解決了微分方程無(wú)精確解的 情況，適用范圍很大。使用一個(gè)兼顧運(yùn)算效率和模擬

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

20、?？梢钥紤]在原來(lái)的Logistic模型后加上一個(gè)傅立葉 函數(shù)（其周期就是工廠生產(chǎn)的周期）和一個(gè)正態(tài)隨機(jī)變量。6、 參考文獻(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) P plot(t,p) grid on （2） 河水濃度恒為0.001mol/L，考慮自然凈

21、化，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); 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') 指數(shù)自然凈化命令：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') 指數(shù)自然凈化 命令： wenti12('d',0,50,-1,'y',

24、9;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 采用線性凈化模型，zs 采用指數(shù)凈化模型 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 表示設(shè)定步長(zhǎng)值 %P2:百分比 %P3:可認(rèn)為的最大凈化年限 %P4:xx-線性自然凈化,zs-指數(shù)自然凈化 %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 %程序調(diào)試所用 %sprintf('停產(chǎn)') 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); %調(diào)用函數(shù) 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('工廠在生產(chǎn)為原來(lái)的%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 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) 被調(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; else error('error,wrong parameter.try again!') end elseif b 'zs&