Genau wie bei linearen Gleichungssystemen kommen in der Praxis nicht nur quadratische Gleichungssysteme vor, also Systeme mit genau so vielen Gleichungen wie Unbekannten.
Abweichungen sind in beiden Richtungen möglich. Sowohl unter- als auch überbestimmte nichtlineare Systeme treten in den Anwendungen auf.
3.4.1 Lösungsmannigfaltigkeiten Bei unterbestimmten Systemen
F(y) = 0
mit F ∈ C1(Rn,Rm) und m < n ist die Lösungsgesamtheit L :=F−1({0}) im Normalfall eine nichtlineare Mannigfaltigkeit der lokalen Dimension n−m.
Den Fall n−m = 1 haben wir oben schon behandelt. Der Schnitt von F−1({0}) mit der Regularitätsmenge
R:={x∈Rn | rang(F′(x)) =m = maximal } (232) war bein−m= 1 homeomorph zum Intervall(0,1)gewesen oder zum Einheitskreis in der Ebene und als Lösungstrajektorie der Differentialgleichungen (212) zum Tangenteinheits-vektorfeld berechenbar.
Das ist bein−m >1alles nicht mehr so einfach, da die lineare Tangentialmannigfaltigkeit Mx0 in einem Lösungspunkt x0 ∈ L eben auch n−m dimensional ist.
Mit einem orthonormalen Tangentensystem T0 ∈ Rn,n−m, dass z.B. durch Gram-Schmidt-Orthonormalisierung aus einem Lösungssatz von
F′(x0)T = 0
gewonnen werden kann, lässt sich Mx0 dann wie folgt schreiben:
{ }
Punkte vonLkann man nahex0 ∈ L nun wieder berechnen, indem man von Startpunkten iny0 ∈Mx0 ausgehend entweder durch Pseudoinvers-Newtonschritte
yn+1 =yn−F′(yn)†F(yn), n= 1,2,3, ... (233) auf die Mannigfaltigkeit zurückstößt150, oder dies - wie in (214) kontinuierlich ausführt.
Wenn die Berechnung der Pseudoinversen in (233) zu aufwendig erscheint, kann man bei Startwerten in Mx0 nahe bei x0 den Newtonschritt auch einfach orthogonal zum Tangen-tensystem T0 ausführen
yn+1 =yn+ ∆(yn).
Dabei löst ∆(yn)das Gleichungssystem
F′(yn)∆(yn) = −F(yn), T0T∆(yn) = 0.
Anmerkungen 3.53 (Begleitende glatte Tangentialbasis)
Innerhalb der Regularitätsmenge hatten wir im Fallen−m= 1zu einerCk-Systemfunktion F(x)ein glattes Tangentialvektorfeld T(x)definiert, welche Ck−1 von x abhing.
Für den Fall n−m >1 haben wir ein solches Vektorfeld bislang nicht brauchen können, da wir ja nicht eine ganze Fläche durch einen Differentialgleichungslöser „abfahren“ lassen können.
Es gibt aber analog zum Satz 3.49 über die Berechnung vonL-Umkehrpunkten einer eines Lösungsastes auch Aufgaben der Berechnung von Umkehrpunkten von Lösungsflächen151. Für ein zu Satz 3.49 analoges Ergebnis benötigt man dann eine Ck−1 Tangentensystem T(x)aus n−m Tangentialvektoren.
Ist T0 = T(x0) das orthonormale Tangentensystem bei x0, so gewinnt man eine Ck−1 -Fortsetzung in eine Umgebung von x0 durch die beiden Schritte
1. Tˆ(x) = [I−F′(x0)T(
F′(x)F′(x0)T)−1
F′(x)]T0
2. Bestimme T(x) aus Tˆ(x)durch Gram-Schmidt-Orthonormalisiereung.
3.4.2 Gauß-Newton-Iterationen
Nun haben wir noch überbestimmte Systeme zu betrachten:
0 =F(x) =
f1(x1, . . . , xn) ...
fm(x1, . . . , xn)
(234)
mit F ∈C1(Rn,Rm) und m > n.
Solche Probleme treten in den Anwendungen relativ häufig auf, wenn die zu lösenden Gleichungen einen Zusammenhang zwischen Größen aus Messwerten erschließen sollen.
In diesem Fall geht man davon aus, dass die Messungen mit Fehlern behaftet sind und führt deshalb - um nicht einer Fehlmessung zu viel Einfluss zu geben - wesentlich mehr Messungn durch als für eine Bestimmung der Daten unbedingt nötig sein würden.
150Wobei alle Techniken zur professionellen Ausführung des Newton-Schrittes bedacht werden sollten, wie etwa die Dämpfung des Schrittes.
151Ohne hier genauer darauf eingehen zu können, merken wir an, dass Verzweigungepunkte von eindimen-sionalemn Lösungsästen im einfachsten Fall interpretiertwerden können als Umkehrpunkte von zugehörigen Lösungsflächen. Vgl. Gegebenenfalls [GrR] und [MH]
Beispiel 3.54 (Ausgleichskreis)
Zur Festlegung der Daten m1, m2, r einer Kreisgleichung (ξ−m1)2+ (η−m2)2 =r2
in der (ξ, η)-Ebene genügt - wie man weiß - die Kenntnis dreier verschiedener Punkte auf dem Kreis. Wenn man sich aber auf drei Punkte verlässt, kann es durch leichte Variation der Daten bei unglücklicher Lage zu starken Veränderungen des Kreises kommen.
−4 −3 −2 −1 0 1 2 3 4 5
−4
−3
−2
−1 0 1 2 3 4
Abbildung 68: Fehlereinfluss bei Kreisberechnung
Bei vielen Anwendungen wird man deshalb mehrfach messen, um aus diesen Daten die wirklichen Werte zu erschließen.
Bei der zementfreien Einpassung der Pfanne einer Hüftgelenk-Endoprothese muss die Pfan-ne äußerst passgenau sitzen. Man muss daher die Ausfräsung dafür entsprechend genau kennen. Da es im Knochenmaterial beim Ausfräsen zu unerwarteten Deformationen kom-men kann, arbeitet man daran, durch Messungen nach diesem Arbeitsgang die geeignete Pfannengröße schnell bestimmen zu können.
Abbildung 69: Kuenstliches Hüftgelenk, Wikimedia Misst man m >3 Punkte
(ξi, ηi), i= 1, . . . , m
eines Kreises, so ergibt das nichtlineare Gleichungssystem (234) in der Form
fi(x) = (ξi−x(1))2+ (ηi−x(2))2−x(3)2, i= 1, . . . , m (235) Die Jacobimatrix J(x) geben wir zur späteren Verwendung auch gleich an:
J(x) =
2(x(1)−ξ1) 2(x(2)−η1) −2x(3) 2(x(1)−ξ2) 2(x(2)−η2) −2x(3)
... ...
.
Es ist klar, dass wir die m Funktionen fi(x1, . . . , xn), i = 1, . . . , m in n < m Variablen i.a. nicht gleichzeitig annullieren können. Wie im Falle des linearen Ausgleichs werden wir sie aber „so null wie möglich“ machen wollen, und das heißt, dass wir wie im quadratischen Fall (m =n) eine TestfunktionT(A, x) :=∥AF(x)∥22 wie in (194) minimieren wollen. Beim Versuch der Lösung eines quadratischen Gleichungssystems gehört die Wahl der regulären metrischen Matrix A zum Aufgabenbereich des numerischen Lösers.
Hier stellen wir uns auf den Standpunkt, dass das Ersetzen des nichtlinearen Gleichungs-systems (234) durch das nichtlineare Ausgleichsproblem
Φ(x) = ∥F(x)∥22 = min! (236) eine Angelegenheit des Aufgabenstellers ist. Er hat zu entscheiden, welchen Einfluss die einzelnen Gleichungen fi(x) = 0 haben sollen, so dass F(x) in (236) gegebenenfalls schon ein Skalierungs-A enthält. Wir werden daher im kommenden Algorithmus eine mögliche Skalierung von F durch den Algorithmus außer Acht lassen.
Wenn wir nun das Minimumproblem
Φ(x) :=F(x)TF(x) = min!
lösen wollen, und Differenzierbarkeit vonF annehmen, finden wir (lokale) Minima offenbar unter den stationären Punkte von Φ, die
∇Φ(x) =F′(x)TF(x) = 0 (237) erfüllen. Wenn die Hessesche, also die zweite Ableitung ∇2Φ(x)von Φexistiert, Lipschitz-stetig vonx abhängt und im stationären Punktx∗ regulär ist, wird das Newton-Verfahren
xk+1 =xk−(
∇2Φ(xk))−1
F′(xk)TF(xk) (238) lokal quadratisch gegebenx∗konvergieren. All unsere Ergebnisse über das Newton-Verfahren und seine Globalisierung könnten wir hier nun anwenden.
Tatsächlich macht man das sehr selten, denn die Hessesche
∇2Φ(x) =F′(x)TF′(x) +
∑m i=1
fi(x)∇2fi(x)
ist nur mit viel Aufwand zu bestimmen. Insbesondere sind die zweiten Ableitungen∇2fi(x), i= 1, . . . , m relativ arbeitsauwändig und bei einer numerischen Approximation zugleich feh-leranfällig.
In der Praxis ersetzt man den Newton-Schritt (238) durch den einfacheren sogenannten Gauss-Newton Schritt
xk+1 =xk−(
F′(xk)TF′(xk))−1
F′(xk)TF(xk). (239) Die Bedingung der Regularität vonF′(xast)wir hier ersetzt durch die Forderung der linea-ren Unabhängigkeit der Spalten von F′(x∗). Dann ist die Matrix (
F′(xk)TF′(xk))
für xk hinreichend nahe bei x∗ invertierbar und der Schritt ausführbar.
Weil in (239) verglichen mit (238) nur den Term mit zweiten Ableitungen in ∇2Φ(xk) fortgelassen hat, kann man dieses Vorgehen als genähertes Newton-Verfahren interpre-tieren. Dies gilt insbesondere wenn im stationären Wert x∗ der Funktionsvektor F(x∗) verschwindet. Denn dann sind die Iterationen in x∗ gleich, und man rechnet außerdem für die Verfahrensfunktion
Ψ(x) = x−(
F′(x)TF′(x))−1
F′(x)TF(x)
leicht nach, dass
Ψ′(x∗) = 0
ist, so dass das Verfahren nach dem Zusatz 3.20 Nr. 2 zum Satz von Ostrowski lokal quadratisch konvergiert.
Für kleine ∥F(x∗)∥-Werte kann man ebenso mit dem Satz von Ostrowski 3.19 die lokal lineare Konvergenz erschließen.
Aufgabe 3.55 Tun Sie’s!
Der Name Gauss-Newton-Verfahren für die Iteration (239) resultiert daraus, dass der Schritt
δk =−(
F′(xk)TF′(xk))−1
F′(xk)TF(xk)
auch verstehbar ist als Lösung des linearen „Gauss’schen“ Ausgleichslösung
∥F(xk) +F′(xk)δ∥2 = min der „Newton’schen“ Linearisierung
0 =F(xk+δ)≈F(xk) +F′(xk)δ.
Das Gegenstück des lokalen Konvergenzsatzes 3.29 für das Newton-Verfahren ist hier Satz 3.56 (Konvergenzverhalten von Gauss-Newton)
Sei F ∈C1,1(M,Rm) mit einer konvexen offenen MengeM ⊂Rn mit m > n. Es sei L die Lipschitzkonstante von F′(x)in M, so dass
∥F′(x)−F′(y)∥ ≤L∥x−y∥für x, y ∈M.
Sei x∗ ∈M ein stationärer Punkt von
Φ(x) :=F(x)TF(x),
so dass F′(x∗)TF(x∗) = 0 ist. Seien weiter die Spalten von F′(x∗) linear unabhängig, so dass die Gauss-Newton-Iteration nahe x∗ ausführbar ist.
Dann gibt es einδ > 0und eine Konstante K >0 so dass für∥xk−x∗∥ ≤δ xk+1 =xk−(
F′(xk)TF′(xk))−1
F′(xk)TF(xk) ausführbar ist und die Abschätzung
∥xk+1−x∗∥ ≤K(
∥xk−x∗∥2+∥F(x∗)∥ · ∥xk−x∗∥) gilt.
Beweis: Sei δ >0 so klein, dass F′(xk)TF′(xk) für ∥xk−x∗∥ < δ stets regulär ist. Dann findet man
xk+1−x∗ = xk−x∗ −(F′(xk)TF′(xk))−1F′(xk)TF(xk)
= (F′(xk)TF′(xk))−1F′(xk)T (
F′(xk)(xk−x∗)−F(xk))
= (F′(xk)TF′(xk))−1F′(xk)T[ (
F′(xk)(xk−x∗)−F(xk) +F(x∗)
| {z }−F(x∗)]
Da wir im Beweis der quadratischen Konvergenz des Newton-Verfahrens nachgewiesen hatten, dass
∥A(xk)∥ ≤ L
2∥c+xk−x∗∥2, und da außerdem
∥F′(xk)−F′(x∗∥ ≤L∥xk−x∗∥
ist. müssen wir nur noch die Koeffizienten dieser beiden Terme im letzten Ergebnis
sam-meln, um zu dem gewünschten Ergebnis zu kommen. 2
Lemma 3.57 (Abstieg von Φ(x) in Gauss-Newton-Richtung) Unter den Voraussetzungen des letzten Satzes ist
δk =−(
F′(xk)TF′(xk))−1
F′(xk)TF(xk) Abstiegsrichtung für Φ in einem nichtstationären xk; denn es ist
d
dtΦ(xk+tδk)|t=0 =−2∥F′(xk)δk∥22 <0 (240)
Beweis: Es ist d
dtΦ(xk+tδk)|t=0 = 2(F′(xk)δk)TF(xk) = 2(δk)T [
F′(xk)TF(xk)]
. (241)
Weil δk die Normalgleichungen
F′(xk)TF′(xk)δk =−F′(xk)TF(xk) erfüllt, kann man[F′(xk)TF(xk)]ganz rechts in (241) ersetzen durch
[F′(xk)TF(xk)] =−F′(xk)TF′(xk)δk,
was gerade (240) ergibt. 2
Die Ergebnis zeigt, dass für das Gauss-Newton-Verfahren eine Armijo-Typ-Dämpfung mög-lich ist, wenn manΦ(x) als Testfunktion benutzt. Das macht auch deshalb besonders viel Sinn, weil man ja tatsächlich den Wert von Φ(·) reduzieren will.
Natürlich kann mit der Gauss-Newton-Richtung auch ein Gauss-Newton-Fluss definiert werden
x′(t) =−(
F′(x)TF′(x))−1
F′(x)TF(x). (242)
Wenn Rechenzeit kein Problem ist, ist dies sicher eine der bequemsten Dämpfungsansätze.
Beispiel 3.58 (Fortsetzung von Beispiel 3.54)
Wir zeigen hier die Anwendung des Gauss-Newton-Flusses zur Berechnung eines Ausgangs-kreises.
% Testprogramm z u r Berechnung e i n e s A u s g l e i c h s k r e i s e s
% mit H i l f e d e r Lösung e i n e s Gauss−Newton−F l u s s e s
% 1 . T e i l E i n g a b e e i n e s K r e i s e s d u r c h d r e i Punkte
% a l s o r i e n t i e r u n g s h i l f e b e i d e r Anschlie�enden E i n g a b e von T e s t d a t e n
a x i s( [−4 , 4 ,−4 , 4 ] ) a x i s e q u a l
hold on z=ginput( 3 )
x=z ( : , 1 ) ; y=z ( : , 2 ) ;
A=[x (2)−x ( 1 ) , y (2)−y ( 1 ) ; x (3)−x ( 1 ) , y (3)−y ( 1 ) ] ;
b=[x (2)^2−x ( 1 ) ^ 2 + y (2)^2−y ( 1 ) ^ 2 ; x (3)^2−x ( 1 ) ^ 2 + y (3)^2−y ( 1 ) ^ 2 ] / 2 ; m=A\b ;
r=sqrt( ( x (1)−m( 1 ) ) ^ 2 + ( y (1)−m( 2 ) ) ^ 2 ) ;
%p l o t ( x , y , ’ or ’ ) ;
p h i=l i n s p a c e( 0 , 2∗pi, 1 0 1 ) ;
xx=m(1)+ r∗cos( p h i ) ; yy=m(2)+ r∗s i n( p h i ) ; plot( xx , yy , ’ b ’ )
% 2 . Nutzung d e s K r e i s e s , um 20 Daten i n s e i n e r Nähe
% e i n z u g e b e n z=ginput( 2 0 )
x=z ( : , 1 ) ; y=z ( : , 2 ) ; f i g u r e
a x i s( [−4 , 4 ,−4 , 4 ] ) a x i s e q u a l
hold on
% e r s t e s P l o t t e n d e r Daten plot( x , y , ’ o r ’ )
% E i n g a b e d e s A n f a n g s k r e i s e s
a0 = input( ’ [ mx, my, r r ] ␣ e i n g e b e n : ’ )
% a0 =[−2 ,−2 ,1] ’;
xx t=a0 (1)+ a0 ( 3 )∗cos( p h i ) ; y y t=a0 (2)+ a0 ( 3 )∗s i n( p h i ) ; plot( xxt , yyt , ’ b ’ )
pause
% 3 . Gauss−Newton−F l u s s zum b e s t e n A u s g l e i c h [ t , a ]= ode45(@( t , a ) r e d u c e ( t , a , x , y ) , [ 0 , 1 0 ] , a0 ) ; M=length( t ) ;
f o r j = 1 : 3 :M
xxn=a ( j , 1 ) + a ( j , 3 )∗cos( p h i ) ; yyn=a ( j , 2 ) + a ( j , 3 )∗s i n( p h i ) ; plot( xxn , yyn , ’ k ’ )
end
%Übermalen d e r Daten .
plot( x , y , ’ o r ’ , ’ MarkerFaceColor ’ , ’ g ’ ) plot( xxn , yyn , ’ r ’ , ’ LineWidth ’ , 2 ) function [ F , J ] = k r e i s ( a , x , y ) ; n = length( x ) ;
f o r k =1:n
F( k ) = ( x ( k)−a ( 1 ) ) ^ 2 + ( y ( k)−a (2))^2−a ( 3 ) ^ 2 ;
J ( k , 1 : 3 ) = [ 2∗( a (1)−x ( k ) ) , 2∗( a (2)−y ( k ) ) , −2∗a ( 3 ) ] ; end
F=F ’ ;
function e r g = r e d u c e ( t , a , x , y ) ; [ F , J ]= k r e i s ( a , x , y ) ;
e r g = −J \F ;
Hinweis zur function „reduce“:Man beachte, dass das Statement für die
Be-unterscheidet. Der Backslash-Operator berechnet im Fall eines regulären Glei-chungssystems die Lösung, im Falle der Überbestimmtheit die Augleichsappro-ximation.
Nachdem in einer ersten Phase 20 Daten für eine Kreis eingegeben werden können, wird ein Startkreis gewählt.
−5 0 5
−4
−3
−2
−1 0 1 2 3 4
Abbildung 70: Daten und Startkreis
Hiernach wird die Gauss-Newton-Fluss-Gleichung von 0 bis 10 integriert. Jeder dritte Kreis der Fluss-Diekretisierung wird geplottet. Der letzte Kreis ist rot gehalten.
−5 0 5
−4
−3
−2
−1 0 1 2 3 4
Ausgleichskreisberechnung mit Gauss−Newton−Fluss
Abbildung 71: Fluss zum Ausgleichskreis Anmerkungen 3.59
Man kann (236) mit einer Matrixfunktion Q(·) ∈ R(m,m), deren Werte Q(x) orthogonal sind, also Q(x)TQ(x) = I erfüllen, äquivalent schreiben als
Φ(x) = ∥Q(x)F(x)∥22 = min!
Wendet man hierauf den Zugang an, so ergibt sich ein vom Gauss-Newton-Schritt abweichender Wert, wennQ(x)nicht konstant ist. Man kann darauf Methoden zur Konvergenzverbesserung aufbauen, cf. [FC].