Visual Basic 數值分析 - Least-Squares curve fitting by using an arbitrary order polynomial as approximating function

Least-Squares curve fitting by using an arbitrary order polynomial as approximating function. (線上執行) (原始碼下載)

最小平方近似法 (least-squares approximation) 是用來求出一組離散 (discrete) 數據點的近似函數 (approximating function)，作實驗所得的數據亦常使用最小平方近似法來達成曲線密合 (curve fitting)。以下所介紹的最小平方近似法是使用多項式作為近似函數，除了多項式之外，指數、對數方程式亦可作為近似函數。關於最小平方近似法的計算原理，請參閱市面上的數值分析書籍。

副程式：

Public Sub LeastSquare(d As Long, x() As Single, f() As Single, ReturnCoeff() As Single)
     Dim i As Long, j As Long, no As Long, m() As Single
     Dim row As Long, col As Long, n() As Single
     row = d + 1
     col = d + 2
     no = UBound(x)
     ReDim m(1 To row, 1 To col), n(0 To 2 * d)
     For i = 0 To 2 * d
          For j = 1 To no
               n(i) = n(i) + x(j) ^ i
          Next
     Next
     For i = 1 To row
          For j = 1 To row
               m(i, j) = n(i + j - 2)
          Next
     Next
     For i = 1 To row
          For j = 1 To no
               m(i, col) = m(i, col) + x(j) ^ (i - 1) * f(j)
          Next
     Next
     Call GEMP(m, ReturnCoeff)
End Sub

其中 GEMP 這個副程式，請參閱本站的 Solve a simultaneous linear equations by Gauss Elimination with maximization of pivot elements.。

參數說明：

d：欲 Fit 多項式的羃次(degree)。
x()：自變數之值，陣列基底為 1。
f()：應變數之值，陣列基底為 1。
ReturnCoeff()：傳回多項式之係數，陣列基底為 1。例如：ReturnCoeff(1)代表常數項之係數，ReturnCoeff(2)代表 1 次方之係數，依此類推。

範例：

試利用最小平方近似法，求出以下數據之 2 階多項式近似函數。

x	1.4	2.5	3.1	3.9	5.0	7.1	9.5	11.9	14.1	15.0	16.5	17.2
f(x)	3.8	2.6	4.1	5.2	6.2	6.9	7.2	6.7	5.8	3.8	3.8	2.8

Private Sub Command1_Click()
     Dim i As Long, keyin As String, A() As Single, f() As Single, x() As Single, d As Long
     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, f) Then GoTo WrongValue
     keyin = InputBox("請輸入欲 Fit 多項式的羃次(degree)", "多項式的羃次")
     d = CLng(keyin)
     Call LeastSquare(d, x, f, A)
     For i = 1 To d + 1
          Debug.Print A(i)
     Next
     Exit Sub
WrongValue:
     MsgBox "輸入了無效的數值", vbExclamation, "警告"
End Sub

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

執行以上的程式將可在VB的即時運算視窗顯示以下結果：

1.130045 
 1.284787 
-6.935126E-02

由此可知 2 階多項式為：1.130045 + 1.284787x - 0.06935126x^²

上面的題目用 1、3 階多項式也可做，甚至用 10 階多項式也可求出答案，但是要如何知道以幾階多項式所求的答案才是和原來的數據點誤差最小呢？有一個很笨的方法就是把原來的數據點與求出的近似函數畫在同一張圖形上來比較，這部分可以用 Graph Control 或 Chart Control 辦到，但畢竟肉眼分辨能力有限，其實我們可以用簡單的統計分析來判斷該用幾階多項式以使得誤差最小，這部分將在下次介紹。

LeastSquare 副程式並沒有錯誤捕捉，例如： d 小於 0 的情況，或是 x 和 f 的陣列數目小於 2 或不相等，都會造成執行錯誤，這部分的錯誤捕捉程式碼就留給網友們自行撰寫了。

3/31/1999 補充：

如果數據點很多的時候，使用 LeastSquare 副程式通常會發生很嚴重的捨入誤差 (roundoff error)，解決辦法是將整個專案的 Single 變數改宣告成 Double，請參考以下的範例：

以下的 x 和 f(x) 是擷取某個實驗的數據，並以 6 階多項式來 fit。

Private Sub Command2_Click()
     Dim i As Long, keyin As String, A() As Single, f() As Single, x() As Single
     keyin = "0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 " & _
     " 0.08 0.085 0.09 0.095 0.1 0.105 0.11 0.115 0.12 0.125 0.13 0.135 0.14 0.145 0.15 0.155 0.16 " & _
     " 0.165 0.17 0.175 0.18 0.185 0.19 0.195 0.2  0.205 0.21 0.215 0.22 0.23 0.24 0.25 0.26 0.27 0.28 " & _
     " 0.29 0.3 0.31 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5"
     Call SepStrToNumArray(keyin, " ", 1, x)
     keyin = "0.71 0.732 0.748 0.761 0.772 0.78 0.787 0.793 0.798 0.799 0.797 0.792 0.786 0.78 0.776 " & _
     " 0.765 0.759 0.752 0.743 0.736 0.723 0.716 0.703 0.696 0.683 0.675 0.664 0.657 0.648 0.639 0.633 " & _
     " 0.623 0.618 0.613 0.607 0.6 0.598 0.593 0.589 0.582 0.579 0.577 0.575 0.574 0.573 0.57 0.568 0.567 " & _
     " 0.567 0.567 0.57 0.571 0.573 0.576 0.579 0.582 0.583 0.584 0.588 0.591 0.593 0.5955 0.597 0.597"
     Call SepStrToNumArray(keyin, " ", 1, f)
     Call LeastSquare(6, x, f, A)
     For i = 1 To UBound(A)
          Debug.Print A(i)
     Next
End Sub

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

執行 Command2 之後，即可在 VB 的即時運算視窗顯示以下結果：

接著我們把整個專案的 Single 變數改宣告成 Double (用 VB 的 "取代" 功能)，其他數據維持不變，再執行 Command2 一次，即可在 VB 的即時運算視窗顯示以下結果：

0.706363278353505 
 4.97858806199404 
-85.8230376299251 
 493.616610321853 
-1339.07277833309 
 1759.51765831704 
-904.758357221887

不要懷疑自己的眼睛，這次執行的結果和上一次差異非常大，並不是程式寫錯，而是因為宣告成 Single 變數會造成比宣告成 Double 變數還要嚴重很多的捨入誤差，因此第二次的執行結果才是正確的，也許有網友會問說：如何知道第二次的執行結果正確？因為我使用其他的數值軟體(MATLAB)來作，也是那樣的結果，所以應該是沒問題的。

捨入誤差是在利用電腦做數值運算時一定會遇到的問題，本站很多的程式都是預設宣告成Single，在這建議網友在真正使用時最好改宣告成 Double，雖然這樣會佔用較多的記憶體以及較長的執行時間 (通常感覺不出來)，但是為了獲得更正確的結果，我想這麼作是值得的。如果需要比宣告成 Double 更小的捨入誤差，可以使用更精確的 Decimal 資料型態，但是一般使用 Double 應該就足夠了。

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

Revised : Mar. 31 , 1999 .

Total pageview since 4/6/1999.