积分运算？

数值积分法有 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 来判断是否收敛。
例如要求  的积分，则执行以下程式：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 FunctionPublic 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 来判断是否收敛。
例如要求和  的积分，必须先在 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)的积分可能会出现计算结果错误的问题，例如，关于这种问题的解决方式，请参阅市面上的数值分析书籍。

解决方案 »