[摘要] 在BOD-DO模型中，牛頓法通過VB編程實(shí)現(xiàn)，在求解河流起始點(diǎn)BOD值將起到核心作用；文章指出了牛頓法在運(yùn)用過程中可能存在的局限性，并提出進(jìn)一步利用數(shù)控技術(shù)實(shí)現(xiàn)計(jì)算機(jī)與環(huán)境監(jiān)測儀器之間通信的設(shè)想。

[關(guān)鍵詞] BOD-DO模型 牛頓法 VB語言

1 BOD-DO模型相關(guān)背景

溶解氧(DO)、生化需氧量(BOD)是表示水中溶解氧情況的兩個(gè)重要指標(biāo)。前者是反映水體中氧數(shù)量多少的指標(biāo)，若低于4mg/L，則魚類和水中大多數(shù)生物都不能存活；后者是指在20℃下生物氧化5天，微生物消耗溶解氧的數(shù)量。關(guān)于BOD和DO的耦合模型[1]，可以寫作：

(1)

(2)

式中： ——河水中的BOD值，mg/L;

——河水中的氧虧值，mg/L；

K ——河水中BOD衰減（耗氧）系數(shù)，1/d;

K ——河流復(fù)氧系數(shù)，1/d;

t ——河水的流行時(shí)間。

其解析式為：(3)

(4)

式中： ——河流起始點(diǎn)的BOD值；

——河流起始點(diǎn)的氧虧值。

式（4）表示河流的氧虧變化規(guī)律。如果以河流的溶解氧來表示，則

(5)

式中： ——河流中的溶解氧濃度；

——飽和溶解氧濃度。

式（5）稱為S—P氧垂公式。

在溶解氧濃度最低的點(diǎn)——臨界點(diǎn)，河水的氧虧值最大，且變化速率為零，則

(6)

(7)

式中： ——臨界點(diǎn)的氧虧值；

tc——由起始點(diǎn)到達(dá)臨界點(diǎn)的流行時(shí)間。

臨界氧虧發(fā)生的時(shí)間t 可以由下式計(jì)算：

tc (8)

2 實(shí)例

擬建一個(gè)工廠，將廢水經(jīng)過處理后排入附近的一條河流中，在現(xiàn)狀條件下，河流中BOD5的濃度是2mg/L，溶解氧的濃度是8 mg/L，河水水溫是20℃，河流流量是14 m /s ;排放的工業(yè)廢水，BOD5的濃度在處理之前為800 mg/L，水溫為20℃，流量為3.5 m /s ，廢水排放前經(jīng)過處理溶解氧濃度為4 mg/L，假定廢水和河水在排放口附近迅速混合，混合后河道中平均水深達(dá)0.8m， 河寬15m，參數(shù)k1（20℃）為0.23/d，k2為3.0/d， 若河流的溶解氧標(biāo)準(zhǔn)為5.0 mg/L，計(jì)算該工廠排出廢水中允許進(jìn)入河流的最大BOD5濃度。

本例的難點(diǎn)在于將上述的（7）和（8）聯(lián)立求解河流起始點(diǎn)BOD值時(shí)，將得到一個(gè)非手工操作計(jì)算的超越方程，使對工廠排出廢水中允許進(jìn)入河流的最大BOD5濃度的求解陷入僵局。筆者將利用Visual Basic 6.0編程工具實(shí)現(xiàn)對牛頓法的編譯，進(jìn)一步實(shí)現(xiàn)對超越方程的破譯，計(jì)算出工廠排出廢水中允許進(jìn)入河流的最大BOD5濃度。

3 牛頓迭代法的流程圖

牛頓迭代法又稱為牛頓切線法，它是基于下列的思想求根:先求出方程的表達(dá)式f(x)的導(dǎo)數(shù)f′(x)，如果在判定方程的區(qū)間[x1，x2]中無解，即f(x1)與f(x2)同號(標(biāo)記為t=0)，輸出無解提示后結(jié)束；若在判定方程的區(qū)間[x1，x2]中有解，即f(x1)與f(x2)異號(標(biāo)記為t=1)的基礎(chǔ)上，將(x1+x2)/2的初值賦予x1，計(jì)算曲線點(diǎn)(x1，f(x1))的切線與x軸的交點(diǎn)x2: x2=x1-f(x1)/ f′(x1)。依法繼續(xù)，直到x1與x2的差滿足精確要求，即足夠接近真正的x*為止，或作切線達(dá)1000次為止。

從上圖中可以看出: f′(x1)=f(x1)/(x1-x2)

因此， x2=x1-f(x1)/f′(x1)

這就是牛頓迭代法.可以利用它由x1求出x2，然后由x2求出x3，直到求出x*為止.

4 牛頓迭代法的VB實(shí)現(xiàn)

先利用VB計(jì)算出飽和DO值和混合后DO值，進(jìn)一步用牛頓迭代法求解河流BOD初始點(diǎn)值（見圖1）。以下是程序的主體部分：

Private Sub C2_Click() ‘計(jì)算河流起始點(diǎn)BOD

……

For i = a To b Step 0.1[2]

If f(i) * f(b) <= 0 Then x1 = i: x2 = b: t = 1

Next i

If t = 0 Then T1(9).Text = \"方程在區(qū)間無解！\": End

x2 = (x1 + x2) / 2

Do

t = t + 1

x1 = x2

x2 = x1 - f(x1) / f1(x1)

Loop While (Abs(x1 - x2) > 1 / 10 ^ (n + 1) And t < 1000)

T1(9).Text = Int(x2 * 10 ^ n + 0.5) / 10 ^ n

……

End Sub

Private Function f(x As Single) As Single ‘函數(shù)f (x)

Dim c As Single， m As Single， k As Single， q As Single

q = (Val(T1(0).Text) * Val(T1(2).Text) + Val(T1(1).Text) * Val(T1(3).Text)) / (Val(T1(2).Text) + Val(T1(3).Text))

c = Val(T1(6).Text) / (Val(T1(7).Text) - Val(T1(6).Text)) * Log(Val(T1(7).Text) / Val(T1(6).Text))

m = (Val(T1(7).Text) - Val(T1(6).Text)) / Val(T1(6).Text) * (Val(T1(4).Text) - q)

k = Val(T1(7).Text) / Val(T1(6).Text) * (Val(T1(4).Text) - Val(T1(8).Text))

f = Log(k / x) - c * m / x + c

End Function

Private Function f1(x As Single) As Single

‘對函數(shù)f(x)的求導(dǎo)得f’(x)

……

f1 = -1 / x + c * m / (x * x)

End Function

圖1 河流起始點(diǎn)BOD值的結(jié)果

5 討論

5.1 本文用牛頓法所得的河流起始點(diǎn)BOD值為水利工作者提供了方便。但是，牛頓法也有它的局限性，如下左圖和右圖分別表示斜率為0和收斂過程中會(huì)停滯于某一點(diǎn)的情況：

5.2 可以進(jìn)一步利用數(shù)控技術(shù)實(shí)現(xiàn)計(jì)算機(jī)與環(huán)境監(jiān)測儀器之間的通信。

參考文獻(xiàn)

[1] 陸書玉. 環(huán)境影響評價(jià)[M]. 北京：高等教育出版社， 2001.

[2] 楊克昌. 計(jì)算機(jī)程序設(shè)計(jì)典型例題精解[M]. 長沙：國防科技大學(xué)出版社，1999.

[3] 譚浩強(qiáng). C程序設(shè)計(jì)題解與上機(jī)指導(dǎo)[M]. 北京：清華大學(xué)出版社，2005.