Das bekannte Newtonverfahren zur Bestimmung einer Nullstelle, dessen Iterationsformel etwa so geschrieben werden kann:
kann für Systeme von nichtlinearen Gelichungen verallgemeinert werden, indem man sich anstelle der einzelnen Funktionswerte und der Funktion Vektoren vorstellt. Bei der Ableitung wird es kompliziert. Hier wäre eine Matrix einzusetzen, die sich aus allen partiellen Ableitungen aller Funktionen nach allen Variablen ergibt und den Namen "Jacobi-Matrix" bzw. "Jacobische-Determinante" trägt. (Wir haben hier vorausgesetzt, daß es sich um genau n Gleichungen und n Variable handelt).
Da man einen Vektor nicht durch eine Matrix teilen kann, sondern stattdessen mit der Inversen multiplizieren muß, wäre die Iterationsformel so zu schreiben (p() sei der Vektor der x-Werte, F() der Vektor der Funktionen; unten alle mit Pfeil auf dem Kopf):
Damit wären wir schon fertig mit der neuen Iterationsformel.
ABER: Aus Gründen, die eigentlich nicht mehr gelten, nämlich wegen des Rechenaufwandes, mag niemand die immer neuen Inversen Matrizen berechnen, weshalb kluge Leute die nur ab und zu neu berechnen, das geht ja im einfachen Newtonverfahren mit der Ableitung auch! Noch klügere, wie wir, denken noch mal kurz nach. Wenn man der Differenz von palt und pneu einen neuen Namen gibt, etwa "Verbesserungsvektor" oder - praktischer - "u()", wenn man das in der Formel berücksichtigt, alles von links mit J() noch einmal multipliziert, bedenkt, daß J*J-1 die Einheitsmatrix ergibt, welche einen mit ihr multiplizierten Vektor reproduziert, dann erhält man:
was wie ein lineares Gleichungssystem in vektorieller Schreibweise aussieht; glaubt es oder nicht: Es ist eins! Die Unbekannten stecken in u(). Damit ist der Algorithmus programmierbar. Einige lächerliche Kleinigkeiten wie die Dimensionierung etc. lasse ich mal weg, wir sind nicht in der Anfängerabteilung.
| 'Newton-Kantorowitsch runde = 0 Do 'Berechnung der Jacobischen Matrix J => J*x(alt)=x(neu) For i = 1 To n For j = 1 To n B(i, j) = Ableitung(n, i, j, x) Next B(i, n + 1) = - Gleichung(i, x) Next gauss n, B, u 'Auswertung summe = 0 For i = 1 To n x(i) = u(i) + x(i) summe = summe + Abs(u(i)) Next Loop Until summe < eps Or runde = 1000runde = runde + 1 |
Hier sind drei Unterprogrammaufrufe versteckt. Die ersten beiden müssen eigens noch programmiert werden:
| Function Ableitung(n%, Gleinum%, Varnum%, xv() As Double) As Double Dim dx As Double, y1 As Double, y2 As Double, dxv(10) As Double Dim i% For i = 1 To n dxv(i) = xv(i) Next End Function y1 = Gleichung(Gleinum, xv) dx = eps * xv(Varnum) If dx = 0 Then dx = eps dxv(Varnum) = dxv(Varnum) + dx y2 = Gleichung(Gleinum, dxv) Ableitung = (y2 - y1) / dx Function Gleichung(num%, x() As Double) As Double Dim g As Double Select Case num Case 1 g = 3 * x(1) + 4 * x(2) ^ 2 - 6 * x(3) + 5 Case 2 g = x(1) ^ 2 - 3 * x(2) + 5 * x(3) - 27 Case 3 g = -5 * x(1) + x(2) + x(3) ^ 2 - 9 End Select End FunctionGleichung = g 'Lösung dieses Systems: x(1)=1,x(2)=-2,x(3)=4 |
Erforderliche Anmerkung: "Ableitung" ist als numerische Ableitung programmiert, weil man damit unabhängig vom konkreten Funktionensatz ist, dessen Lösung gefunden werden soll. Anders "Gleichung", hier stehen die veränderlichen Gleichungen, die man bei Visual Basic leider nur unter großen Schwierigkeiten zur Laufzeit des Programmes eintragen kann.
Weiterhin wird die Lösung eines linearen Gleichungssystems mittels "gauss ..." aufgerufen. Es ist ja nicht wie bei armen Leuten: Unser Gaußalgorithmus tut hier seinen Dienst.
Mit den Startwerten x(i) = 0.1 für alle drei Parameter und eps = 1E-9 erhält man dieses Ergebnis:
Gute Berechnung!
