Visual Basic 數值分析 - Evaluate an integral by using Romberg integration

Evaluate an integral by using Romberg integration. (線上執行) (原始碼下載)

數值積分法有 Trapezoidal rule、Simpson's rule、Romberg integration、Gauss quadrature 等，Trapezoidal 和 Simpson 很適合用於求解等距(equally spaced)數據點的積分，而 Romberg 和 Gauss quadrature 很適合用於求解複雜函數的積分，以下要介紹的是 Romberg integration，關於數值積分的計算原理，請參閱市面上的數值分析書籍。

範例一：多項式的積分

副程式：

Public Function RombergPoly(c() As Single, a As Single, b As Single, tol As Single) As Single
     Dim i As Long, j As Long, T() As Single, L As Long
     Dim x As Single, dx As Single, n As Long, sum As Single
     ReDim T(1 To 3)
     T(1) = (b - a) * (fxP(c, a) + fxP(c, b)) / 2
     T(2) = T(1) / 2 + (b - a) * (fxP(c, (a + b) / 2)) / 2
     T(3) = (4 * T(2) - T(1)) / 3
     j = 3
     Do While Abs(T(UBound(T)) - T(UBound(T) - 2)) > tol
          dx = (b - a) / (2 ^ (j - 1))
          x = a - dx
          n = 2 ^ (j - 2)
          sum = 0
          For i = 1 To n
               x = x + 2 * dx
               sum = sum + fxP(c, x)
          Next
          For i = 2 To UBound(T) Step 2
               T(i - 1) = T(i)
          Next
          T(2) = T(1) / 2 + dx * sum
          ReDim Preserve T(1 To UBound(T) + 2)
          For L = 2 To j
               If L <> j Then
                    T(L * 2) = ((4 ^ (L - 1)) * T(L * 2 - 2) - T(L * 2 - 3)) / ((4 ^ (L - 1)) - 1)
               Else
                    T(UBound(T)) = ((4 ^ (L - 1)) * T(UBound(T) - 1) - _
                    T(UBound(T) - 2)) / ((4 ^ (L - 1)) - 1)
               End If
          Next
          j = j + 1
     Loop
     RombergPoly = T(UBound(T))
End Function

Public Function fxP(c() As Single, x As Single) As Single
     Dim i As Long
     fxP = c(UBound(c))
     For i = UBound(c) - 1 To LBound(c) Step -1
          fxP = c(i) + fxP * x
     Next
End Function

參數說明：

c()：代表多項式的係數，陣列基底為 1，c(1)代表常數項，c(2)代表 1 次方之係數，依此類推。
a：積分下限。
b：積分上限。
tol：容許誤差，使用 absolute convergence criterion 來判斷是否收斂。

例如要求 VB_Numerical-14-1.gif (333 bytes) 的積分，則執行以下程式：

Private Sub Command1_Click()
     Dim Keyin As String, c() As Single, a As Single, b As Single
     Keyin = InputBox("請輸入多項式的係數，低階項的係數先輸入，再輸入高階項的係數" & _
   "，以空白來區隔。" & vbNewLine & "例如：f(x)=10x^3-5x+3，則依序輸入 3  -5  0  10")
     If Not SepStrToNumArray(Keyin, " ", 0, c) Then MsgBox "輸入了非數值": Exit Sub
     a = Val(InputBox("請輸入積分下限"))
     b = Val(InputBox("請輸入積分上限"))
     Debug.Print RombergPoly(c, a, b, 0.000001)
End Sub

其中 SepStrToNumArray 這個副程式，請參閱本站的將字串拆解成陣列。

Keyin 輸入 -1 0 0 2 0 0 0 1，積分下限輸入 -1，積分上限輸入 2，則在 VB 的即時運算視窗顯示以下積分值：

36.375

只要是多項式函數都可經由 RombergPoly 副程式求得積分，那如果要求其它函數的積分該怎麼辦？請參考範例二的介紹。

範例二：任意函數的積分

副程式：

Public Function Romberg(Problem As Long, a As Single, b As Single, tol As Single) As Single
     Dim i As Long, j As Long, T() As Single, L As Long
     Dim x As Single, dx As Single, n As Long, sum As Single
     ReDim T(1 To 3)
     T(1) = (b - a) * (fx(Problem, a) + fx(Problem, b)) / 2
     T(2) = T(1) / 2 + (b - a) * (fx(Problem, (a + b) / 2)) / 2
     T(3) = (4 * T(2) - T(1)) / 3
     j = 3
     Do While Abs(T(UBound(T)) - T(UBound(T) - 2)) > tol
          dx = (b - a) / (2 ^ (j - 1))
          x = a - dx
          n = 2 ^ (j - 2)
          sum = 0
          For i = 1 To n
               x = x + 2 * dx
               sum = sum + fx(Problem, x)
          Next
          For i = 2 To UBound(T) Step 2
               T(i - 1) = T(i)
          Next
          T(2) = T(1) / 2 + dx * sum
          ReDim Preserve T(1 To UBound(T) + 2)
          For L = 2 To j
               If L <> j Then
                    T(L * 2) = ((4 ^ (L - 1)) * T(L * 2 - 2) - T(L * 2 - 3)) / ((4 ^ (L - 1)) - 1)
               Else
                    T(UBound(T)) = ((4 ^ (L - 1)) * T(UBound(T) - 1) - _
                    T(UBound(T) - 2)) / ((4 ^ (L - 1)) - 1)
               End If
          Next
          j = j + 1
     Loop
     Romberg = T(UBound(T))
End Function

Public Function fx(Problem, x As Single) As Single
     Select Case Problem
          Case 1
               fx = Sin(Log(x))
          Case 2
               fx = x ^ 7 + 2 * x ^ 3 - 1
     End Select
End Function

參數說明：

Problem：代表題號。提供這個參數的好處是：不論要計算幾題積分問題，都只要簡單地在 fx 函數中增加屬於該題的程式碼，完全不須更動 Romberg 副程式。
a：積分下限。
b：積分上限。
tol：容許誤差，使用 absolute convergence criterion 來判斷是否收斂。

例如要求 VB_Numerical-14-2.gif (314 bytes) 和 VB_Numerical-14-1.gif (333 bytes) 的積分，必須先在 fx 函數中撰寫相對應的程式碼(如上所述)，然後執行以下程式：

Private Sub Command2_Click()
     Debug.Print Romberg(1, 1, 10, 0.000001)
     Debug.Print Romberg(2, -1, 2, 0.000001)
End Sub

即可在VB的即時運算視窗顯示以下積分值：

7.56091 
 36.375

上述的 Romberg 副程式至少還有二個問題存在，但是對於解決一般的積分問題，仍然可以適用。

無法處理 singularity 的問題。
一些週期性函數(periodic function)的積分可能會出現計算結果錯誤的問題，例如，關於這種問題的解決方式，請參閱市面上的數值分析書籍。

[ 上一個 | 首頁 | 數值分析 | 下一個 ]

Total pageview since 4/6/1999.