求根據已知點擬合曲線函數的代碼(急)


現有一組坐標,X,Y坐標值分別保存在X(),Y()兩個數組中,大約有30組。求用已知的這組點,擬合曲線函數,用來生成曲線,求實現該功能的代碼。我曾試用過很多程序(大多數算法是最小二乘法),但都有問題,會溢出。只要調試成功就給分。

12 个解决方案

#1


試試這個~

'****************************************************************************************************
'** X()------Double 實型一維數組,長度為 n 。存放給定 n 個數據點的 X 坐標。               **
'** Y()------Double 實型一維數組,長度為 n 。存放給定 n 個數據點的 Y 坐標。               **
'** n-------Integer 變量。給定數據點的個數。                                                **

'** a()------Double 實型一維數組,長度為 m 。返回 m-1 次擬合多項式的 m 個系數。            **
'** m-------Integer 變量。擬合多項式的項數,即擬合多項式的最高次數為 m-1。                  **
'**               要求 m<=n 且m<=20。若 m>n 或 m>20 ,則本函數自動按 m=min{n,20} 處理。            **
'**rdblAverageX--Double 變量,返回給定n個數據點的 X 坐標的平均值                                 **
'** dt()------Double 實型一維數組,長度為 3。其中:                                          **
'**               dt(0) 返回擬合多項式與數據點誤差的平方和;                                       **
'**               dt(1) 返回擬合多項式與數據點誤差的絕對值之和;                                   **
'**               dt(2) 返回擬合多項式與數據點誤差絕對值的最大值。                                 **
'****************************************************************************************************

Private Sub Iapcir(X() As Double, _
                 Y() As Double, _
                ByVal n As Integer, _
                ByRef a() As Double, _
                ByVal m As Integer, _
                ByRef rdblAverageX As Double, _
                ByRef dt() As Double)
                
   Dim i As Integer, j As Integer, K As Integer
   Dim Z As Double, P As Double, C As Double, G As Double, Q As Double, D1 As Double, D2 As Double
   Dim S(19) As Double, t(19) As Double, B(19) As Double
   
   For i = 0 To m - 1
      a(i) = 0
   Next i
   
   If m > n Then m = n
   If m > 20 Then m = 20
   
   Z = 0#
   
   For i = 0 To n - 1
      rdblAverageX = rdblAverageX + X(i)
      Z = Z + X(i) / (1# * n)
   Next i
   rdblAverageX = rdblAverageX / n
   
   B(0) = 1#
   D1 = 1# * n
   P = 0#
   C = 0#
   
   For i = 0 To n - 1
      P = P + (X(i) - Z)
      C = C + Y(i)
   Next i
   
   C = C / D1
   P = P / D1
   a(0) = C * B(0)
      
   If m > 1 Then
      t(1) = 1#
      t(0) = (-1) * P
      D2 = 0#
      C = 0#
      G = 0#
      For i = 0 To n - 1
         Q = X(i) - Z - P
         D2 = D2 + Q * Q
         C = C + Y(i) * Q
         G = G + (X(i) - Z) * Q * Q
      Next i
      
      C = C / D2
      P = G / D2
      Q = D2 / D1
      D1 = D2
      a(1) = C * t(1)
      a(0) = C * t(0) + a(0)
   End If
   
   For j = 2 To m - 1
      S(j) = t(j - 1)
      S(j - 1) = (-1) * P * t(j - 1) + t(j - 2)
      
      If j >= 3 Then
         For K = j - 2 To 1 Step -1
            S(K) = (-1) * P * t(K) + t(K - 1) - Q * B(K)
         Next K
      End If
      
      S(0) = (-1) * P * t(0) - Q * B(0)
      
      D2 = 0#
      C = 0#
      G = 0#
      
      For i = 0 To n - 1
         Q = S(j)
         
         For K = j - 1 To 0 Step -1
            Q = Q * (X(i) - Z) + S(K)
         Next K
         
         D2 = D2 + Q * Q
         C = C + Y(i) * Q
         G = G + (X(i) - Z) * Q * Q
      Next i
      
      C = C / D2
      P = G / D2
      Q = D2 / D1
      D1 = D2
      a(j) = C * S(j)
      t(j) = S(j)
      
      For K = j - 1 To 0 Step -1
         a(K) = C * S(K) + a(K)
         B(K) = t(K)
         t(K) = S(K)
      Next K
   Next j
   
   dt(0) = 0#
   dt(1) = 0#
   dt(2) = 0#
   
   For i = 0 To n - 1
      Q = a(m - 1)
      
      For K = m - 2 To 0 Step -1
         Q = a(K) + Q * (X(i) - Z)
      Next K
      
      P = Q - Y(i)
      
      If Abs(P) > dt(2) Then
         dt(2) = Abs(P)
      End If
      dt(0) = dt(0) + P * P
      dt(1) = dt(1) + Abs(P)
   Next i
   
End Sub

#2


1.擬合出來的函數不太對啊,比如我用Y=X^2去實驗,擬合出來的函數與這不相符,將原來的X()值代入擬合函數,得出的Y()值也不對。我想用擬合出的函數繪制曲線,現在繪制出來的不太對。請再解釋下吧,謝謝。


#3


最小二乘法,本質上是直線的擬合。二次曲線也罷,指數函數也罷,都是化成一次函數以后,再進行最小二乘法計算。所以,萬能的方法可能沒有(比如一個圓,就頭痛),只能先肉眼判斷一下它像什么曲線,再進行擬合。要萬能一點也就是:多一些擬合的方案,然后選擇一種最好的方案(如:相關系數最大)。

#4


上面的代碼是多項式擬合的
擬合效果不理想

你把會溢出的最小二乘法貼出來看看~

#5


郎格—庫塔法則

#6


以下是我收集的曲線擬合程序正文,使用的時候會有溢出,請各位高手幫我看看。。。。。。謝謝。。。。。。,

引用自http://202.101.71.156/forum/forum_posts.asp?TID=10310&KW=%D7%EE%D0%A1%B6%FE%B3%CB%B7%A8

首先介紹多元線性回歸公式(即最小二乘法計算公式):

回歸公式:f(X1,X2,…,Xk)=Y=A0+A1X1+…+AkXk

要求:{∑[Y-(A0+A1X1+…+AkXk)]^2}min

N為輸入數據個數,輸入數據為:X1,X2,…,Xk,Y

設:

Lij=∑XiXj-∑Xi∑Xj/N

Liy=∑XiY-∑Xi∑Y/N

Lyy=∑Y^2-∑Y∑Y/N

注:∑Xi表示輸入數據中所有Xi的和,∑XiXj表示每次輸入數據中的XiXj的積的和

 

則A0,A1,…,Ak滿足方程組:

A0=(∑Y-A1∑X1-…-Ak∑Xk)/N

┌L11A1+L12A2+…+L1kAk=L1y

│L21A1+L22A2+…+L2kAk=L2y

┤   :       :              :

│   :       :              :

└Lk1A1+Lk2A2+…+LkkAk=Lky

解此方程組,於是A0、A1、……、Ak便為回歸系數。

若把回歸公式改為:f(X)=Y=A0+A1X+A2X^2+…+AkX^k

這樣一來就變成我需要的多次曲線擬合算法。

詳細代碼如下(主函數為GetMINNHvalue):

'=======================================================================

'hRoot()--------雙精度一維數組

'               指定回歸公式,如hRoot(3)=1則表明回歸公式中有三次項

'                               hRoot(2)=0則表明回歸公式中沒有二次項

'hMaxIndex------指定hRoot()的最大下標,最小為0

'hData()--------雙精度二維數組,第二維從0至1

'               點的坐標數據

'LenData--------指定點的個數,即hData()第一維的最大下標

'解出的回歸系數保存在hRoot()中,返回值為相關系數(若為負則出錯)

'-----------------------------------------------------------------------

Public Function GetMINNHvalue(hRoot() As Double, hMaxIndex As Long, hData() As Double, ByVal LenData As Long) As Long
Dim d1 As Double, D2 As Double, d3 As Double
Dim z1 As Long, z2 As Long, z3 As Long
Dim zfx As Long, zfy As Long
Dim Sum() As Double, data() As Double, FC() As Double, Root() As Double
z2 = 0
For z1 = 1 To hMaxIndex
If hRoot(z1) <> 0 Then z2 = z2 + 1
Next
'=============================================================================
zfx = z2: zfy = LenData
ReDim Sum(zfx, zfy - 1) As Double, data(zfy - 1) As Double
'-----------------------------------------------------------------------------
For z1 = 0 To zfy - 1
    data(z1) = hData(z1, 0)
    Sum(zfx, z1) = hData(z1, 1)
Next
z3 = 0
For z2 = 0 To zfx - 1
z3 = z3 + 1
BegIF1: If hRoot(z3) = 0 Then z3 = z3 + 1: GoTo BegIF1
For z1 = 0 To zfy - 1
    d1 = data(z1)
    Sum(z2, z1) = d1 ^ z3
Next
Next
'------------------------------------------------------------------------------
ReDim FC(zfx + 1, zfx) As Double, Root(zfx + 1) As Double
For z2 = 0 To zfx
For z1 = 0 To zfx
    d1 = 0: D2 = 0: d3 = 0
    For z3 = 0 To zfy - 1
        d1 = d1 + Sum(z1, z3) * Sum(z2, z3)
        D2 = D2 + Sum(z1, z3)
        d3 = d3 + Sum(z2, z3)
    Next
    FC(z1, z2) = d1 - D2 * d3 / zfy
Next
Next
'------------------------------------------------------------------------------
Call ExplainEquation(FC(), Root(), zfx - 1)
If Root(zfx) = 0 Then GetMINNHvalue = -1: Exit Function
On Error GoTo ED
d1 = 0
For z1 = 0 To zfy - 1
    d1 = d1 + Sum(zfx, z1)
Next
For z2 = 0 To zfx - 1
    D2 = 0
    For z1 = 0 To zfy - 1
        D2 = D2 + Sum(z2, z1)
    Next
    d1 = d1 - D2 * Root(z2) / Root(zfx)
Next
hRoot(0) = d1 / zfy
z1 = 0
For z2 = 1 To hMaxIndex
If hRoot(z2) <> 0 Then
hRoot(z2) = Root(z1) / Root(zfx)
z1 = z1 + 1
End If
Next
d1 = 0
For z1 = 0 To zfx - 1
    d1 = d1 + Root(z1) * FC(zfx, z1)
Next
GetMINNHvalue = Sqr(d1 / Root(zfx) / FC(zfx, zfx))
On Error GoTo 0
Exit Function
ED:
GetMINNHvalue = -2
End Function

'========================================================================

'行列式求值函數

'------------------------------------------------------------------------

Public Function GetHLValue(n() As Double, ByVal nWide As Long) As Double
Dim nt(127) As Double
Dim z1 As Long, z2 As Long, z3 As Double, z4 As Double, z5 As Double
If nWide > 2 Then
z5 = 0
    For z1 = 0 To nWide - 1
        z4 = n(z1, nWide)
        If z4 = 0 Then GoTo NE
        For z2 = 0 To nWide '- 1
        nt(z2) = n(z1, z2)
        n(z1, z2) = n(nWide, z2)
        Next z2
        z3 = GetHLValue(n(), nWide - 1)
        For z2 = 0 To nWide '- 1
        n(z1, z2) = nt(z2)
        Next z2
        If ((nWide) Mod 2) = 1 Then z5 = z5 + z3 * z4 Else z5 = z5 - z3 * z4
NE:
    Next z1
z3 = GetHLValue(n(), nWide - 1)
If ((nWide) Mod 2) = 0 Then
GetHLValue = z5 + n(nWide, nWide) * z3
Else
GetHLValue = z5 - n(nWide, nWide) * z3
End If
Exit Function
End If
If nWide = 2 Then
For z2 = 0 To nWide
z3 = 1: z4 = 1
    For z1 = 0 To nWide
        z3 = z3 * n((z2 + z1) Mod (nWide + 1), z1)
        z4 = z4 * n((z2 - z1 + nWide + 1) Mod (nWide + 1), z1)
    Next z1
GetHLValue = GetHLValue + z3 - z4
Next z2
ElseIf nWide = 1 Then
GetHLValue = n(0, 0) * n(1, 1) - n(0, 1) * n(1, 0)
ElseIf nWide = 0 Then
GetHLValue = n(0, 0)
End If
End Function

'========================================================================

'求兩數的最大公約數

'------------------------------------------------------------------------

Public Function GetZZXCvalue(ByVal N1 As Long, ByVal N2 As Long) As Long
Dim z1 As Long
N1 = Abs(N1): N2 = Abs(N2)
If N1 < N2 Then z1 = N1 Else z1 = N2
If z1 = 0 Then
GetZZXCvalue = N1 + N2
If GetZZXCvalue = 0 Then GetZZXCvalue = 1
Exit Function
End If
If z1 = 1 Then GetZZXCvalue = 1: Exit Function
Begin:
If N1 < N2 Then
N2 = N2 Mod N1
If N2 = 0 Then GetZZXCvalue = N1: Exit Function
Else
N1 = N1 Mod N2
If N1 = 0 Then GetZZXCvalue = N2: Exit Function
End If
GoTo Begin
End Function

'=======================================================================

'求解多元線性方程組

'-----------------------------------------------------------------------

Public Function ExplainEquation(FC() As Double, Root() As Double, ByVal Height As Long)
Dim z1 As Long, z2 As Long
Root(Height + 1) = GetHLValue(FC(), Height)
For z1 = 0 To Height
    For z2 = 0 To Height
        FC(Height + 2, z2) = FC(z1, z2)
        FC(z1, z2) = FC(Height + 1, z2)
    Next z2
    Root(z1) = GetHLValue(FC(), Height)
    For z2 = 0 To Height
        FC(z1, z2) = FC(Height + 2, z2)
    Next z2
Next z1
For z1 = 0 To Height + 1
    If Root(z1) <> Int(Root(z1)) Then Exit Function
Next
z2 = GetZZXCvalue(Root(0), Root(1))
For z1 = 0 To Height
z2 = GetZZXCvalue(z2, Root(z1 + 1))
Next z1
If Root(Height + 1) < 0 Then z2 = -z2
For z1 = 0 To Height + 1
Root(z1) = Root(z1) / z2
Next z1
End Function

#7


另:samwzhang(分全給我) ,你能說一下什么是郎格—庫塔法則嗎?你說的是不是龍格庫塔法?這個法則和任意曲線的擬合有關系嗎?我急需這方面的資料,有條件的話給我提供下吧,謝謝!
Email:wljin2004@sina.com

#8


龍格—庫塔法則是解決常微分方程初值問題的,用在這里不太合適.
如果你只是想得到一組已確定數據的擬和曲線,其實可以不用在去編那復雜的程序了,已經有很多軟件為我們做好了這一切,比如在matlab中,只要一條指令a=polyfit(x,y,n)就得到了數據x,y對應的n階多項式系數a。
如果你想自己編寫的話,我這也有他人寫的最小二乘法的程序。

#9


我需要比較精確的擬合程序,可以的話發到我的信箱吧,謝謝!

#10


個人覺得:最好不要先肯定是“n階多項式”的形式。如果是一條螺旋線(比如是銀河系的星星組成的點),用“n階多項式”的形式合適嗎?如果能肯定是“n階多項式”的形式,那么用指數曲線擬合的話,也會擬合得很好。

#11


感謝大家,已經按照AprilSong的程序調試通了,擬合出來的函數是Y=A(0)+A(1)(X-rdblAverageX)+…A(I)(X-rdblAverageX)^I,效果還可以,如果還有更精確的,也請大家能發出來,大家一起討論。。。。。。

#12


哎,你會不會微積分啊,數值微分的目的是什么?就是取得系數啊。高次曲線擬合的時候,由於xy實際數量級的差異和誤差被次方放大,所以要轉成微分來處理,計算出系數以后再轉回去啊。
你不會真的這么笨吧?

注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2021 ITdaan.com