Weitere Infos hier: http://de.wikipedia.org/wiki/Methode_de ... n_Quadrate
Die hier benutzte Form der Berechnung wird bei höheren Graden (>8)leicht instabiel, was sich dadurch äußert, dass sich das entstandene Gleichungssystem nicht mehr lösen lässt.
Code: Alles auswählen
; Author: alter Mann , adaptiert von Engeln-Müllges/Uhlig "Numerical Algorithms with C"
; Date: 23. 10. 2008
; OS: Windows
; Demo: Ja
;
; Polynomapproximation (Ausgleichpolynom)
; ( wird leicht instabil bei hohen Polynomgraden und großen Stützpunkt-Werten)
#DBL_EPSILON =2.2204460492503131e-016 ; kleinster Wert, so dass 1.0+DBL_EPSILON != 1.0
#EPS_QUAD =(#DBL_EPSILON*#DBL_EPSILON)
EnableExplicit
Procedure.d MathSqrtD (Wert.d)
Protected Wert1.d
If Wert <= 0.0
Wert1 = 0.0
Else
Wert1 = Sqr(Wert)
EndIf
ProcedureReturn Wert1
EndProcedure
Macro SQRT(wert)
MathSqrtD(wert)
EndMacro
Macro SQR2(wert)
(wert)*(wert)
EndMacro
;************************************************************************
;* ein Polynom P in der Darstellung *
;* P(x) = adA[0] + adA[1] * dX + adA[2] * dX^2 + ... + adA[n] * dX^n *
;* nach dem Hornerschema auswerten *
;* *
;* Eingabeparameter: *
;* ================= *
;* lN: Grad des Polynoms *
;* adA: [0..n]-Vektor mit den Koeffizienten des Polynoms *
;* dX: Stelle, an der das Polynom auszuwerten ist *
;* *
;* Funktionswert: *
;* ============== *
;* P(x) *
;************************************************************************
Procedure.d Num_Horner (lN.l, adA.d(1), dX.d) ; Hornerschema zur Polynomauswertung .............*/
Protected ldErg.d = adA(lN)
Protected i.l
For i=lN-1 To 0 Step -1
ldErg = ldErg * dX + adA(i)
Next i
ProcedureReturn ldErg
EndProcedure
;*====================================================================*
;* *
;* Num_ChoCec berechnet die Zerlegungsmatrix einer symmetrischen, *
;* positiv definiten Matrix. Die Zerlegungsmatrix wird auf mat *
;* ueberspeichert. *
;* *
;*====================================================================*
;* *
;* Eingabeparameter: *
;* ================ *
;* n Long n ( n > 0 ) *
;* Dimension von mat, *
;* Anzahl der Komponenten des b-Vektors. *
;* mat Double mat(n,n) *
;* Matrix des Gleichungssystems. *
;* *
;* Ausgabeparameter: *
;* ================ *
;* mat Double mat(n,n) *
;* Dekompositionsmatrix, die die Zerlegung von *
;* mat enthaelt. *
;* *
;* Rueckgabewert: *
;* ============= *
;* = 0 alles ok *
;* = 1 n < 1 gewaehlt *
;* = 2 Matrix ist nicht positiv definit *
;* *
;*====================================================================*
Procedure.l Num_ChoDec (n.l,mat.d(2))
Protected j.l, k.l, i.l
Protected sum.d
If n < 1
ProcedureReturn 1 ; n < 1 Fehler !
EndIf
If mat(0,0) < #EPS_QUAD
ProcedureReturn 2 ; mat ist nicht positiv definit !
EndIf
mat(0,0) = SQRT(mat(0,0))
For j = 1 To n-1 Step 1
mat(j,0) / mat(0,0)
Next j
For i = 1 To n-1 Step 1
sum = mat(i,i)
For j = 0 To i-1 Step 1
sum - SQR2(mat(i,j))
Next j
If sum < #EPS_QUAD
ProcedureReturn 2 ; nicht positiv definit
EndIf
mat(i,i) = SQRT(sum)
For j = i + 1 To n-1 Step 1
sum = mat(j,i)
For k = 0 To i-1 Step 1
sum - mat(i,k) * mat(j,k)
Next k
mat(j,i) = sum / mat(i,i)
Next j
Next i
ProcedureReturn 0
EndProcedure
;*====================================================================*
;* *
;* Num_ChoSol bestimmt die Loesung x des linearen Gleichungssystems *
;* B * x = b mit der unteren Dreieckmatrix B wie sie als Ausgabe *
;* von chodec bestimmt wird. *
;* *
;*====================================================================*
;* *
;* Eingabeparameter: *
;* ================ *
;* n Long n ( n > 0 ) *
;* Dimension von lmat, Anzahl der Komponenten *
;* des b-Vektors u. des Loesungsvektors x. *
;* lmat Double lmat(n,n) *
;* untere Dreieckmatrix, wie sie von chodec als Aus- *
;* gabe geliefert wird. *
;* b Double b(n) *
;* Rechte Seite des Gleichungssystems. *
;* *
;* Ausgabeparameter: *
;* ================ *
;* x Double x(n) *
;* Loesungsvektor des Systems. *
;* *
;* Rueckgabewert: *
;* ============= *
;* = 0 alles ok *
;* = 1 Zerlegungsmatrix unzulaessig oder n < 1 *
;* *
;*====================================================================*
Procedure.l Num_ChoSol (n.l,lmat.d(2),b.d(1),x.d(1))
Protected j.l, k.l
Protected sum.d
If n < 1
ProcedureReturn 1
EndIf
If lmat(0,0) = 0.0
ProcedureReturn 1 ; Unzulaessige Zerlegungsmatrix
EndIf
x(0) = b(0) / lmat(0,0) ; Vorwaertselimination
For k = 1 To n-1 Step 1
sum = 0.0
For j = 0 To k-1 Step 1
sum + lmat(k,j) * x(j)
Next j
If lmat(k,k) = 0.0
ProcedureReturn 1
EndIf
x(k) = (b(k) - sum) / lmat(k,k)
Next k
x(n-1) / lmat(n-1,n-1) ; Rueckwaertselimination
For k = n - 2 To 0 Step -1
sum = 0.0
For j = k + 1 To n-1 Step 1
sum + lmat(j,k) * x(j)
Next j
x(k) = (x(k) - sum) / lmat(k,k)
Next k
ProcedureReturn 0
EndProcedure
;*====================================================================*
;* *
;* Die Funktion cholesky dient zur Loesung eines linearen *
;* Gleichungssystems: mat * x = b *
;* mit der n x n Koeffizientenmatrix mat und der rechten Seite b *
;* nach dem Cholesky-Verfahren. *
;* *
;* Die Eingabematrix muss symmetrisch und positiv definit sein, *
;* d.h. fuer alle n Vektoren y != 0 muss gelten: *
;* T *
;* y * mat * y > 0. *
;* *
;* cholesky arbeitet nur auf der unteren Dreieckmatrix incl. der *
;* Diagonalen, so dass gegenueber dem Gaussverfahren eine erhebliche *
;* Reduktion des Speicherplatzes erreicht wird; es genuegt daher *
;* diesen Teil der Matrix an cholesky zu uebergeben. *
;* *
;*====================================================================*
;* *
;* Anwendung: *
;* ========= *
;* *
;* Lineare Gleichungssysteme mit n x n Koeffizientenmatrix, *
;* die symmetrisch und positiv definit ist. *
;* *
;*====================================================================*
;* *
;* Steuerparameter: *
;* =============== *
;* mod Long mod; *
;* Aufrufart von cholesky: *
;* = 0 Bestimmung der Zerlegungsmatrix und Berechnung *
;* der Loesung des Gleichungssystems. *
;* = 1 Nur Berechnung der Zerlegungsmatrix; wird auf *
;* mat ueberspeichert. *
;* = 2 Nur Loesung des Gleichungssystems; zuvor muss je- *
;* doch die Zerlegungsmatrix bestimmt sein. Diese *
;* Aufrufart wird verwendet, falls bei gleicher *
;* Matrix lediglich die rechte Seite des Systems vari- *
;* iert, z. B. zur Berechnung der Inversen. *
;* *
;* Eingabeparameter: *
;* ================ *
;* n Long n; ( n > 0 ) *
;* Dimension von mat, Anzahl der Komponenten *
;* des b-Vektors u. des Loesungsvektors x. *
;* mat Double mat(n,n) *
;* mod = 0, 1: Matrix des Gleichungssystems. *
;* mod = 2 : Zerlegungsmatrix. *
;* b Double b(n) ( bei mod = 0, 2 ) *
;* Rechte Seite des Gleichungssystems. *
;* *
;* Ausgabeparameter: *
;* ================ *
;* mat Double mat(n,n) ( bei mod = 0, 1 ) *
;* Dekompositionsmatrix, die die Zerlegung von *
;* mat enthaelt. *
;* x Double x(n); ( bei mod = 0, 2 ) *
;* Loesungsvektor des Systems. *
;* *
;* Rueckgabewert: *
;* ============= *
;* = 0 alles ok *
;* = 1 n < 1 gewaehlt oder falsche Eingabeparameter *
;* = 2 Matrix ist nicht positiv definit *
;* = 3 Falsche Aufrufart *
;* *
;*====================================================================*
;* *
;* Benutzte Funktionen: *
;* =================== *
;* Num_ChoDec() : Bestimmt die Dekomposition *
;* Num_ChoSol() : Loest das lineare Gleichungssystem *
;* *
;*====================================================================*
Procedure Num_Cholesky (mod.l,n.l,mat.d(2),b.d(1),x.d(1))
Protected rc.l
If n < 1
ProcedureReturn 1
EndIf
If mod = 0 ; Zerlegung bestimmen und Gleichungssystem loesen
rc = Num_ChoDec (n, mat())
If rc = 0
ProcedureReturn Num_ChoSol (n, mat(), b(), x())
Else
ProcedureReturn rc
EndIf
ElseIf mod = 1 ; Nur Zerlegung berechnen
ProcedureReturn Num_ChoDec (n, mat())
ElseIf mod = 2 ; Nur Loesung bestimmen
ProcedureReturn Num_ChoSol (n, mat(), b(), x())
EndIf
ProcedureReturn 3 ; Falsche Aufrufart
EndProcedure
;************************************************************************
;* die Koeffizienten eines Ausgleichspolynoms n-ten Grades nach der *
;* diskreten Fehlerquadratmethode von Gauss berechnen. *
;* Das Verfahren beruht auf den Gaussschen Normalgleichungen, die mit *
;* dem Choleskyverfahren geloest werden: Die Koeffizienten c[j] des *
;* Vektors *
;* *
;* ( w[0] * (c[0] + c[1] * x[0]^1 + ... + c[n] * x[0]^n) ) *
;* ( w[1] * (c[0] + c[1] * x[1]^1 + ... + c[n] * x[1]^n) ) *
;* ( ... ... ) *
;* ( w[m] * (c[0] + c[1] * x[m]^1 + ... + c[n] * x[m]^n) ) *
;* *
;* sollen so berechnet werden, dass er der rechten Seite *
;* *
;* ( w[0] * y[0] ) *
;* ( w[1] * y[1] ) *
;* ( ... ) *
;* ( w[m] * y[m] ) *
;* *
;* in der Euklidischen Norm moeglichst nahe kommt. *
;* c ergibt sich dann als Loesung des linearen Gleichungssystems *
;* a * c = b (Normalgleichungen) *
;* mit *
;* a[i][j] = w[0] * x[0]^(j+i) + ... + w[m] * x[m]^(j+i), *
;* b[i] = w[0] * y[0] * x[0]^i + ... + w[m] * y[m] * x[m]^i *
;* fuer i,j=0(1)n. *
;* *
;* Eingabeparameter: *
;* ================= *
;* m: Nummer der letzten Stuetzstelle *
;* n: Hoechstgrad des algebraischen Ausgleichspolynoms *
;* x: [0..m]-Vektor mit den Stuetzstellen *
;* y: [0..m]-Vektor mit den Funktionswerten zu den Stuetzstellen *
;* w: [0..m]-Vektor mit den Gewichten *
;* *
;* Ausgabeparameter: *
;* ================= *
;* c: [0..n]-Vektor mit den Koeffizienten des Ausgleichspolynoms *
;* *
;* Funktionswert: *
;* ============== *
;* 0: kein Fehler *
;* 1: Fehler in den Eingabedaten: n zu klein oder m zu klein *
;* 2: Fehler in choly() *
;* *
;************************************************************************
Procedure Num_ApproxPoly (n.l,m.l,x.d(1),y.d(1),w.d(1),c.d(1))
Dim a.d(n,n) ; die [0..n,0..n]-Matrix des Gleichungssystems
Dim b.d(n) ; [0..n]-Vektor mit der rechten Seite des Gleichungssystems
Protected summand.d ; Hilfsvariable zur Berechnung eines Matrixelements
Protected i.l, j.l ; Laufvariablen
If n < 0 Or m < n ; n oder m zu klein?
ProcedureReturn 1
EndIf
; die erste Spalte, die letzte Zeile der Matrix und die rechte Seite des Gleichungssystems berechnen
For i = 0 To n Step 1
a(i,0) = 0.0
a(n,i) = 0.0
b(i) = 0.0
Next i
For j = 0 To m Step 1
summand = w(j)
For i = 0 To n-1 Step 1
a(i,0) + summand
b(i) + summand * y(j)
summand * x(j)
Next i
a(n,0) + summand
b(n) + summand * y(j)
For i = 1 To n Step 1
summand * x(j)
a(n,i) + summand
Next i
Next j
;die Matrix vervollstaendigen, indem immer wieder Zeile i + 1 versetzt in Zeile i kopiert wird
For i = n - 1 To 0 Step -1
For j = 0 To n-1 Step 1
a(i,j + 1) = a(i + 1,j)
Next j
Next i
i = Num_Cholesky(0, n + 1, a(), b(), c()); ; das Gleichungssystem loesen
If i <> 0 ; Fehler beim Loesen des Gleichungsystems?
ProcedureReturn 2
EndIf
ProcedureReturn 0
EndProcedure
;=======================================================================================
;- Demo
Define lib.l = 0
#Image1 = 1
#Window = 2
#ImageGad = 3
#Button = 4
#Option1 = 5
#Option2 = 6
#MaxX = 800
#MaxY = 700
#MaxAnz = 1000
Dim X.d(#MaxAnz)
Dim Y.d(#MaxAnz)
Dim W.d(#MaxAnz)
Dim C.d(#MaxAnz)
Define i.l
Define Anzahl.l = 0
Define Grad.l = 0
Define Event.l
Define Aktion.l = #Option1
Define Input$
Define PosX.l
Define PosY.l
Define Xmin.d
Define Xmax.d
Define XA.l
Define XE.l
Define YA.l
Define YE.l
If MathInit(@lib) = 0
MessageRequester("Achtung!","Kann Lib nicht öffnen")
End
EndIf
If CreateImage(#Image1,#MaxX,#MaxY) = 0
MessageRequester("Achtung!","Kann Image nicht erstellen")
End
EndIf
If OpenWindow(#Window,0,0,#MaxX+10,#MaxY+55,"Ausgleichskurve",#PB_Window_ScreenCentered)
If CreateGadgetList(WindowID(#Window))
ImageGadget (#ImageGad, 5,10, #MaxX,#MaxY,ImageID(#Image1))
ButtonGadget(#Button,5,#MaxY+20,100,25, "Schließen")
OptionGadget(#Option1,110,#MaxY+15,100,20,"Punkte eingeben")
OptionGadget(#Option2,110,#MaxY+35,100,20,"Approximieren")
EndIf
StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
StopDrawing()
SetGadgetState(#ImageGad,ImageID(#Image1))
SetGadgetState(#Option1,1)
Repeat
Event = WaitWindowEvent()
Select Event
Case #PB_Event_Gadget
Select EventGadget()
Case #ImageGad
If Aktion = #Option1 And EventType() = #PB_EventType_LeftClick And Anzahl < #MaxAnz
PosX = WindowMouseX(#Window) - GadgetX(#ImageGad)
PosY = WindowMouseY(#Window) - GadgetY(#ImageGad)
X(Anzahl) = PosX
Y(Anzahl) = PosY
W(Anzahl) = 1.0
Anzahl + 1
StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
For i = 0 To Anzahl-1
Box (Int(X(i)) - 2,Int(Y(i)) - 2,4,4,RGB(250,250,0))
Next
StopDrawing()
SetGadgetState(#ImageGad,ImageID(#Image1))
EndIf
Case #Button
Event = #PB_Event_CloseWindow
Case #Option1
If Aktion <> #Option1
Anzahl = 0
Grad = 0
StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
StopDrawing()
SetGadgetState(#ImageGad,ImageID(#Image1))
Aktion = #Option1
EndIf
Case #Option2
If Aktion <> #Option2
If Anzahl >= 2
Repeat
Grad = Anzahl-1
Input$ = InputRequester("Grad des Ausgleichspolynoms", "Grad (min.:1/max."+Str(Grad)+"):", Str(Grad))
Grad = Val(Input$)
Until Grad>0 And Grad < Anzahl
i = Num_ApproxPoly(Grad,Anzahl-1,X(),Y(),W(),C())
If i > 0
MessageRequester("Achtung!","Fehler (" + Str(i) + ") bei der Approximation!")
StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
StopDrawing()
SetGadgetState(#ImageGad,ImageID(#Image1))
Aktion = #Option1
Anzahl = 0
SetGadgetState(#Option1,1)
Else
StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
Xmin = #MaxX
Xmax = 0.0
For i = 0 To Anzahl-1
Box (Int(X(i)) - 2,Int(Y(i)) - 2,4,4,RGB(250,250,0))
If X(i) < Xmin
Xmin = X(i)
EndIf
If X(i) > Xmax
Xmax = X(i)
EndIf
Next
XA = Int(Xmin)
XE = Int(Xmax)
For i=XA To XE Step 1
YE = Num_Horner(Grad,C(),i)
If i>XA
LineXY(i-1,YA,i,YE,RGB(255,255,255))
EndIf
YA = YE
Next i
StopDrawing()
SetGadgetState(#ImageGad,ImageID(#Image1))
Aktion = #Option2
EndIf
Else
MessageRequester("Achtung!","zu wenig Punkte")
SetGadgetState(#Option1,1)
EndIf
EndIf
EndSelect
EndSelect
Until Event = #PB_Event_CloseWindow
EndIf
CloseWindow(#Window)
[code\]
