牛頓法在求解ＢＯＤ－ＤＯ模型中河流起始點ＢＯＤ的ＶＢ實現

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

[關鍵詞] BOD-DO模型 牛頓法 VB語言

1 BOD-DO模型相關背景

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

(1)

(2)

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

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

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

K ——河流復氧系數，1/d;

t ——河水的流行時間。

其解析式為：(3)

(4)

式中： ——河流起始點的BOD值；

——河流起始點的氧虧值。

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

(5)

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

——飽和溶解氧濃度。

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

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

(6)

(7)

式中： ——臨界點的氧虧值；

tc——由起始點到達臨界點的流行時間。

臨界氧虧發生的時間t 可以由下式計算：

tc (8)

2 實例

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

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

3 牛頓迭代法的流程圖

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

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

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

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

4 牛頓迭代法的VB實現

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

Private Sub C2_Click() ‘計算河流起始點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 = \"方程在區間無解！\": 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 ‘函數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

‘對函數f(x)的求導得f’(x)

……

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

End Function

圖1 河流起始點BOD值的結果

5 討論

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

5.2 可以進一步利用數控技術實現計算機與環境監測儀器之間的通信。

參考文獻

[1] 陸書玉. 環境影響評價[M]. 北京：高等教育出版社， 2001.

[2] 楊克昌. 計算機程序設計典型例題精解[M]. 長沙：國防科技大學出版社，1999.

[3] 譚浩強. C程序設計題解與上機指導[M]. 北京：清華大學出版社，2005.