Also, wenn … das mit der Fehlerquadratur und dem Excel ist schon keine ganz schlechte Idee. Wenn man das Problem auf sowas reduzieren kann und man einen Iterationsknecht sucht, ist der Solver echt nicht schlecht. Wenn man aber selber was verstehen oder machen will, ist es eine schlechte Lösung.
Gut, für eindimensionale Funktionen (also irgendwas mit f(x) = 0) scheinst Du ja die entsprechenden Verfahren zu kennen. Ich werd mal noch ebend kurz etwas resümieren, damit wir beide über das gleiche reden. Und weil mir der Newton am besten gefällt werd ich den erklären.
Erstmal muß ich darauf hinweisen, daß Newton nur Nullstellen sucht. Wenn Du also eine Funktion g(x) = a mit a != 0 hast, mußt Du die umformen zu f(x) = g(x)-a = 0, klar. Was wird denn nun eigentlich genau gemacht? Die Funktion wird in einem Punkt linarisiert, von dieser Linearisierung wird die Nullstelle bestimmt und das ist der neue Punkt für die nächste Linearisierung. Ausgehend von f(x) ist am Punkt x0 eine linarisierte Funktion
phi(x)=f’(x0)*(x - x0) + f(x0).
Die Nullstelle davon ist
phi(x)=0 -> x = x0 - f(x0)/f’(x0).
Das dürfte Dir vielleicht bekanntvorkommen, weil es genau das ist, was man nach NEWTON machen muß.
Und genau das machen wir jetzt für mehrdimensionale Funktionen. Wir gehen von einer Funktion fn(xk) aus. Die hochgestellten Indezes laufen immer von 1…N bzw. 1…K usw. Sind keine Potenzen, sondern Laufindezes, x1 würde dem x in Deinem Beispiel entsprechen, x2 dem y, x3 dem z usw. je nachdem, wieviele Du noch hast. Vorzugsweise sollte N==K sein, also Du solltest genausoviele Gleichungen wie Unbekannte haben. Das ist nicht zwingend, aber notwendige Bedingung für eindeutige Lösbarkeit. In Deinem Fall ist diese Bedingung auch erfüllt. Alles, was wir jetzt noch machen ist die Linearisierung an der Stelle xk0 (jetzt wird vielleicht auch klar, warum hochgestellt und tiefgestellt: die hochgestellten laufen durch Variablen- oder Funktionslisten, die tiefgestellten durch die Iteration):
phin(x)=D(fn)(xk0) * (xk - xk0) + fn(xk0)
Bei näherer Betrachtung ist das fast das gleiche wie bei der eindimensionalen NEWTON-Iteration. Was noch zu klären wäre ist vielleicht das D(fn)(xk0). Das Dingens nennt man Jakobi-Matrix. Und zwar wird für das Element (i)(j) so einer Matrix die erste partielle Ableitung von fi nach xj an der Stelle xj0 gebildet. Wissen Zivis, was eine partielle Ableitung ist? Ach, ich mach später sowieso noch ein Beispiel. Und für das Dingens wird jetzt ganz normal 'ne Nullstelle gesucht. Also
phin(x) = 0 -> xk = xk0 - fn(xk0)*(D(fn)(xk0))-1.
Die ()-1 am Ende um die Jakobi-Matrix bedeutet, daß die Inverse gemeint ist. Das ist auch der Grund, warum genausoviele Gleichungen wie Variablen da sein sollten, weil sonst diese Matrix nicht quadratisch wäre und man sich nur Ärger damit einhandelt. Matrizeninversion? Hm, das sind irgendwie alles nicht gerade Sachen, die man aus dem Abi mitbringt oder während des Zivildienstes lernt. Egal, am Beispiel wird’s vielleicht klar. Das dabei ermittelte xk ist wieder Startpunkt für die neue Linearisierung usw.
Jetzt Dein Beispiel:
0.001 =(x+y+z)(x) /(6-x+y)
0.24 =(x+y+z)/(56-2*z)
4=(x+y+2*z)*(x+z)/(2-3x)
Formen wir erstmal um:
(x + y + z)(x) /(6 - x + y) - 0.001 =0
(x + y + z)/(56 - 2*z) - 0.24 =0
(x + y + 2*z)*(x + z)/(2 - 3x) - 4 =0
Um mal auf die Notation oben zu verweisen wäre jetzt
f1(xk) = f1(x1, x2, x3) = f1(x, y, z) = (x + y + z)(x) /(6 - x + y) - 0.001
f2(xk) = f1(x1, x2, x3) = f2(x, y, z) = (x + y + z)/(56 - 2*z) - 0.24
f3(xk) = f1(x1, x2, x3) = f3(x, y, z) = (x + y + 2*z)*(x + z)/(2 - 3x) - 4
Zuerst mußt Du Dir einen Anfangspunkt suchen. Ich hab keine Ahung, wo der sein könnte, nehmen wir also
xk0 = (x0, y0, z0).
Jetzt die Jakobi-Matrix und die Sache mit den partiellen Ableitungen: man tut einfach so, als sei die Funktion nur von einer Variable abhängig, nämlich von der, nach der man ableitet. Alle anderen sieht man als konstant an. Zum Beispiel bei der … ähm … ih, die enthalten ja alle Brüche … naja, muß ich durch. Nemen wir die dritte und leiten sie nach z ab:
s/w=u*v/w = (x + y + 2*z)*(x + z)/(2 - 3x)
(Die 4 hab ich gleich mal weggelassen, weil die durch die Ableitung eh rausfällt.)
So, wie war das jetzt…
s’ = (u*v)’ = u’*v + u*v’ =
2*(x + z) + (x + y + 2z)*1 =
2x + 2z + x + y + 2z = 3x + y + 4z
(s/w)’ = (s’*w - s*w’)/w² =
((3x + y + 4z)*(2 - 3x) - (x + y + 2*z)*(x + z)*0)/(2 - 3x)² =
(3x + y + 4z)*(2 - 3x)/(2 - 3x)² =
(3x + y + 4z)/(2 - 3x)
Puh, verstanden? Mit den anderen geht das analog, ich will jetzt nicht alles auseinanderklamüsern.
Am Ende erhälst Du eine Jakobi-Matrix, die so aussehen dürfte (Wie stell ich das gleich dar?):
- Zeile:
(12x-x²+2xy+6y+y²+6z+zy)/(6-x+y)² | x*(6-2x-z)/(6-x+y)² | x/(6-x+y)
- Zeile:
1/(56 - 2z) | 1/(56 - 2z) | (28+x+y)/2/(-28+z)²
- Zeile:
(4x-3x²+6z+2y+3yz+6z²)/(-2+3x)² | (x+z)/(2-3x) | (3x+4z+y)/(2-3x)
Bitte nochmal nachrechnen, können sich Fehler eingeschlichen haben. In der ersten Zeile geht’s nur um f1, zweite f2 usw., erste Spalte ist partielle Ableitung nach x, zweite nach y usw. Das läßt sich für beliebig große Systeme entsprechend fortsetzen. In diese Matrix setzt Du jetzt das xk0 ein. Und da das konkrete Zahlenwerte enthält erhälst Du eine Matrix nur mit Zahlen. Jetzt geht’s an die Inversion. Setz 'nen Taschenrechner drauf an oder ein Computerrechensystem. Oder greif Dir ein Buch [1] und find’s selber raus. Wie sowas geht würde den Rahmen deutlich sprengen. Außerdem gibt’s für numerisch bestückte Matrizen ausgefeilte Algorithmen, wie man sowas am günstigsten machen kann, besonders für große und spärlich besetzte Matrizen hat die FEM da gute Methoden entwickelt. Ich schweife ab.
So, dann mußt Du noch in die Funktionen fn das xk0 einsetzen, dann mußt Du bissl wissen, wie man mit Matrizen rechnet, alles zusammenwürfeln und Du kriegst Deinen neuen Punkt für die Linearisierung.
War das jetzt verständlich? Wenn nicht, dann frag nach, am besten per Mail, und wir werden versuchen den Problemen gemeinsam auf die Schliche zu kommen.
Das Zentrum.
[1] Bronstein-Semendjajew: Taschenbuch der Mathematik;