Visual Basic 數值分析 - Solve the integral of discrete evenly and unevenly spaced points by using trapezoidal or Simpson's rule.

Solve the integral of discrete evenly and unevenly spaced points by using trapezoidal or Simpson's rule. (線上執行) (原始碼下載)

離散數據點 (discrete point) 的積分可以用許多方法求得，比較常用的是 trapzoidal rule和 Simpson's rule，這兩種方法都非常簡單，其中 trapzoidal rule 可以用於自變數為等間距和非等間距取樣的積分問題，但是誤差較大；Simpson's rule 則只能用於自變數為等間距取樣的積分問題，但是比較準確。關於其他求離散數據點積分的方法，請參閱市面上的數值分析書籍。trapzoidal rule 就是以梯形面積公式來求積分，我想在這應該不用贅述了。以下就簡單說明 Simpson's rule：

註：關於函數的積分求值問題，請參閱本站的 Evaluate an integral by using Romberg integration.

Simpson's one-third rule：用於當有 N+1 個等間距取樣的離散數據點，且 N 為偶數時，利用以下的公式即可求出積分值。h 為取樣間距，至於公式的推導請參閱市面上的數值分析書籍。

VB_Numerical-16-1.gif (1176 bytes)

Simpson's three-eighths rule ：用於當有 N+1 個等間距取樣的離散數據點，且 N 為 3 的倍數時，利用以下的公式即可求出積分值。h 為取樣間距，至於公式的推導請參閱市面上的數值分析書籍。

VB_Numerical-16-2.gif (1374 bytes)

若 N 為偶數，直接使用 Simpson's one-third rule 即可求出積分；若N為奇數，則 1 ~ (N+1-3) 個數據點可利用 Simpson's one-third rule，(N+1-3 )~ (N+1) 個數據點可利用 Simpson's three-eighths rule，最後再把兩部分的積分值相加即可。

副程式：

Public Function Trapezoid(x() As Single, fx() As Single) As Single
     Dim i As Long, l As Long, u As Long, d As Long, sum As Single
     d = LBound(x): u = UBound(fx): l = LBound(fx)
     If u - l <> UBound(x) - d Then MsgBox " x 和 f(x) 的數據組數不吻合": Exit Function
     If u - l < 1 Then MsgBox "至少要有二個數據點方可進行Trapezoid積分法": Exit Function
     d = d - l
     sum = 0
     For i = l To u - 1
          sum = sum + (fx(i) + fx(i + 1)) * Abs(x(i + d + 1) - x(i + d))
     Next
     Trapezoid = sum / 2
End Function

參數說明：

x()：一維陣列，儲存自變數之值。
fx()：一維陣列，儲存應變數之值。

副程式：

Public Function Simpson(dx As Single, fx() As Single) As Single
     Dim i As Long, j As Long, u As Long, l As Long
     Dim sum As Single, EvenPanel As Boolean
     u = UBound(fx): l = LBound(fx)
     If u - l < 2 Then MsgBox "至少要有三個數據點方可進行Simpson積分法": Exit Function
     EvenPanel = (u - l) Mod 2 = 0
     sum = fx(l)
     For i = 1 To (u - l - IIf(EvenPanel, 0, 3) - 1)
          sum = sum + fx(i + l) * IIf(i Mod 2 = 0, 2, 4)
     Next
     sum = sum + fx(i + l)
     For j = 1 To u - i - l + 1 - IIf(EvenPanel, 1, 0)
          sum = sum + fx(j + i + l - 1) * IIf(j = 1 Or j = 4, 9, 27) / 8
     Next
     Simpson = sum * dx / 3
End Function

參數說明：

dx：自變數的取樣間距。
fx()：一維陣列，儲存應變數之值。

範例一：(等間距取樣)

試求出 f(x) 的積分：

x	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.10	1.1
f(x)	93	87	68	55	42	37	35	39	48	53	51	39

由於是等間距取樣，因此我們可以分別用 trapzoidal rule 和 Simpson's rule 來求積分。

Simpson's rule：

Private Sub Command1_Click()
     Dim keyin As String, fx() As Single, dx As Single
     keyin = InputBox("請輸入應變數f(x)值，以空白來區隔。" & vbNewLine & _
          "例如：10 20 30", "f(x)值")
     If Not SepStrToNumArray(keyin, " ", 1, fx) Then GoTo WrongValue
     keyin = InputBox("請輸入自變數x的取樣間距", "x的取樣間距")
     If Not IsNumeric(keyin) Then GoTo WrongValue
     dx = CSng(keyin)
     Debug.Print Simpson(dx, fx)
     Exit Sub
WrongValue:
     MsgBox "輸入了無效的數值"
End Sub

執行 Command1，f(x) 輸入 93 87 68 55 42 37 35 39 48 53 51 39，在 VB 的即時運算視窗可得積分為 58.39584 。

Trapzoidal rule：

Private Sub Command2_Click()
     Dim keyin As String, x() As Single, fx() As Single, dx As Single
     keyin = InputBox("請輸入自變數x值，以空白來區隔。" & vbNewLine & "例如：1 2 3", "x值")
     If Not SepStrToNumArray(keyin, " ", 1, x) Then GoTo WrongValue
     keyin = InputBox("請輸入應變數f(x)值，以空白來區隔。" & vbNewLine & _
     "例如：10 20 30", "f(x)值")
     If Not SepStrToNumArray(keyin, " ", 1, fx) Then GoTo WrongValue
     Debug.Print Trapezoid(x, fx)
     Exit Sub
WrongValue:
     MsgBox "輸入了無效的數值"
End Sub

執行 Command2，x 輸入0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 ，
f(x) 輸入 93 87 68 55 42 37 35 39 48 53 51 39，在 VB 的即時運算視窗可得積分為 58.1。

範例二：(非等間距取樣)

試求出 f(x) 的積分：

x	0.99	2.1	3.22	4.40	5.7	7.12	8.01	8.37	9.32	9.98
f(x)	4.9	5.7	4.2	7.04	8.31	7.82	5.97	7.01	6.68	4.79

由於這是非等間距取樣，因此我們只能用 trapzoidal rule 來求積分，不能用 Simpson's rule。

執行上述的 Command2，x 輸入 0.99 2.1 3.22 4.4 5.7 7.12 8.01 8.37 9.32 9.98，
f(x) 輸入 4.9 5.7 4.2 7.04 8.31 7.82 5.97 7.01 6.68 4.79，在 VB 的即時運算視窗可得積分為 58.2492。

還有一種方法可以求得非等間距取樣的積分問題，就是利用多項式最小平方法求得最適多項式，再利用解析解 (Analytical Solution) 的方式求出積分值。例如以上的數據利用本站的 Least-Squares curve fitting with statistical analysis by using an arbitrary order polynomial as approximating function. ，在 1~8 階多項式之中可求出最適多項式為 3 階，如下所示：

f(x) = -0.0251x^³ + 0.2971x^² - 0.4911x + 5.0769

積分上限為 9.98，下限為 0.99，對 f(x) 積分後可求得積分值為 57.5259。

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

Total pageview since 4/6/1999.