Hat man ein nichtlineares Gleichungssystem
F(x) = 0, F :D−→Rn
durch ein x∗ im Definitionsbereich D ⊂ Rn gelöst, so ist man selten (lange) zufrieden.
Wie die Frau im „Fischer un siene Fru“ wird man die die Modellierung durch das Glei-chungssystem weiter verbessern wollen und z.B. bestimmte Systemparameter λ1, . . . , λk, , die bislang feste Werte λ01, . . . , λ0k , verändern wollen, dass die Lösung x(λ01, . . . , λ0k) die Wirklichkeit noch besser wiedergibt oder bei Ingenieurwissenschaftlicher Umsetzung mehr Geld, Umsatz, Sicherheit, Kundenzufriedenheit, Umweltverträglichkeit etc. etc. generiert.
Gehen die Parameter glatt (am besten differenzierbar) in das System ein, so wird man die dann zu erwartende glatte Änderung des Systems ausnutzen wollen. Der Satz über die implizite Funktion liefert uns im Falle der differenzierbaren Abhängigkeit von x und von λ1, . . . , λk Aussagen über die Fortsetzung der Lösungx(λ01, . . . , λ0k) := x∗ zu Lösungen x(λ1, . . . , λk) für (λ1, . . . , λk) nahe (λ01, . . . , λ0k).
Obwohl die Numerische Mathematik schon länger Algorithmen zur Verfügung stellt, mit denen auch von mehreren Parameteren abhängige Lösungsmannigfaltigkeiten angenähert werden können, gehen wir hier davon aus, dass man bei mehreren Parametern in solchen Parameteränderungstudien zunächst einmal einzelne Parameter getrennt von einander va-riieren wird129.
Wir gelangen dadurch zu einem von einem Zusatzparameter abhängigen System
F(x, λ) = 0, F :D×I −→Rn (200) in dem x in einer (meist offen angenommenen) Teilmenge D des Rn wählbar ist und λ in einem reellen Intervall I.
Dies entspricht genau der bei den Homotopie-Verfahren zur Globalisierung iterativer Me-thoden vorgefundenen Situation, wobei die Parameterabhängigkeit hier nicht künstlich herbeigeführt wird, sondern systemimmanent ist.
Wie bei den Homotopie-Verfahren werden wir auch hier Annahmen über gewissen Regulari-täten von Jacobi-Matrizen machen müssen, wenn wir die von linearen Gleichungssystemen gewohnte Annahme, dass die Lösungsgesamtheit beinGleichungen undn+ 1unbekannten normalerweise eindimensional ist, weiter verwenden wollen.
Ohne solche Annahmen können Lösungsgesamtheiten beliebig scheußlich aussehen, und wenn die Funktionen selbst noch so glatt sind130.
Um uns nicht bei der Analyse von Lösungsmannigfaltigkeiten durch die Beachtung von einschränkenden Definitionsbereichen oder Restriktionen an die Differenzierbarkeoit durch-einander bringen zu lassen. nehmen wir in diesem Unterabschnitt an, dassDundIin (200) der ganz Rn bzw. ganz R sind und dass F überall beliebig häufig differenzierbar ist.
129Alles andere ist recht komplex und nur dann wirklich hilfreich, wenn in diesen Lösungsmannigfaltigkei-ten bestimmt strukturbestimmende Parameterkonstellationen gesucht werden, eine Aufgabe, die für diese Ausarbeitung zu involviert wäre. Siehe jedoch gegebenenfalls: [WCR]
130Nach einem Satz von Whitney gibt es zu jeder abgeschlossenen MengeM imRneine unendlich häufig stetig differenzierbare Funktionf :Rn−→Rderen NullstellenmengeM ist.
Wir werden die Lösungsmenge S :={
(x, λ)∈Rn+1 |F(x, λ) = 0} zunächst nur in der (offenen) „Regularitätsmenge“
R :={
(x, λ)∈Rn+1 |Rang(F′(x, λ)) =n= maximal.}
(201) untersuchen.
Es ist für die mathematische Behandlung der Aufgabe meistens einfacher, wenn man die Sonderrolle von λ aufgibt, und (x, λ) einfach als (n+ 1)-dimensionalen Vektor y ∈ Rn+1 verwendet. Jeder Punkt y0 ∈ S ∩ R ist dann ein Punkt, in dem auf F der Satz über implizite Funktionen angewendet werden kann. y0 ∈ S besagt nämlich, dass F(y0) = 0 ist, und aus y0 ∈ R folgt, dass es unter den n+ 1Spalten der Jacobimatrix F′(y0)∈ R(n,n+1) n linear unabhängige Spalten gibt. Nach den zugehörigen Variablen lässt sich das System F(y) = 0 in einer Umgebung von y0 eindeutig auflösen, wobei die restliche Variable als Parameter herangezogen wird. Es besteht S ∩ R damit aus eindimensionalen „Stücken“, die entweder geschlossen einen Kreis bilden oder bijektiv auf das offene Intervall (0,1) abgebildet werden können. Alle solchen Lösungsabschnitte werden i.a. als „Lösungsäste“
bezeichnet.
−5 −4 −3 −2 −1 0 1 2 3 4 5
−5
−4
−3
−2
−1 0 1 2 3 4 5
λ
x
1 Teilstück ~ Kreis (grün), 11 Teilstück ~ (0,1) (blau), 5 Verzeigungspunkte (rot)
Abbildung 51: Verzweigungsdiagramm
Die rot gezeichneten Schnittpunkte der Kurve in Abbildung 51 zerschneiden diese Kurve in die genannten bijektiv und stetig auf das Intervall (0,1)abbildbaren Kurvenabschnitte.
In den Punkten selbst ist der Rang Jacobimatrix kleiner als n. Andernfalls müsste man in ihnen die Lösungsgesamtheit ja in eindeutiger Weise eindimensional fortsetzen können.
Man nennt diese Punkte Verzweigungspunkte des Diagrammes der Äste und das ganze Diagramm „Verzweigungsdiagramm“. Verzweigungsdiagramme können wesentlich kompli-zierter aussehen, als das aus der obigen Skizze.
Wenn man solche Verzweigungsdiagramme erstellen will, muss man natürlich eigentlich mit solchen Verzweigungspunkten ein wenig umgehen zu können. Einerseits muss man über sie „hinwegkommen“, wenn man bei der Astverfolgung auf sie stößt. Andererseits sind sie besonders interessant, weil das Verhalten der Funktion F in ihnen und in ihrer Umgebung meist strukturgebend für die Diagramme ist131
Leider haben wir nicht die Zeit, auf die Behandlung solcher „singulären Punkte“ einzu-gehen. Um an ihnen zumindest mit den hier zu besprechenden Astverfolgungsmethoden
131Dass Elemente, die bestimmte den Umgang mit ihnen angenehm gestaltende Eigenschaften vermissen lassen, einerseits (glücklicherweise?) nicht so häufig vorkommen und dabei andererseits aber oft ablauf-oder formbestimmend sind, ist schon Laotse im „Dao Dö Djing“ (siehe [LAO]) aufgefallen: „Dreißig Speichen treffen die Nabe - die Leere in der Mitte aber macht das Rad“.
„vorbei zu kommen“ hilft ein Trick aus der Zeit bevor man die Singularitätzen richtig behandeln konnte. Durch eine kleine Störung der Gleichung brechen die „Kreuzungen“
meistens in nichtverbundene glatte Straßen auf. In den folgenden drei Skizzen sehen wir die Lösungsgesamheiten der Gleichung
λx=x3 und ihrer gestörten Versionen λx=x3+ 0.1 und λx=x3−0.1
−2 −1 0 1 2 3 4
−3
−2
−1 0 1 2 3
Pitchfolk−Verzweigung zu λ x = x3
−2 −1 0 1 2 3 4
−3
−2
−1 0 1 2 3
λ x = x3 − 0.01
−2 −1 0 1 2 3 4
−3
−2
−1 0 1 2 3
λ x = x3 − 0.01
Abbildung 52: Reine und gestörte „Pitchfork-Verzweigung“
3.3.1 Prädiktor-Korrektor-Astverfolgungsprobleme.
Die einfachste Zugang zu Lösungsästen ist der in Abbildung 34 schon geschilderte. Man erhöht eine zur Parametrisierung geeignete Variable schrittweise und passt die restlichen n Variablen iterativ an die Änderung an.
Wenn der Pfad bezüglich des gewählten Parameters umkehrt, wird man eine andere Kom-ponente als Parameter wählen müssen. Dass ist programmtechnisch relativ aufwendig. Man muss eventuell häufig das lösbare (n, n)-System neu zusammanstellen.
Anstatt die n Gleichungen mit n + 1 Variablen auf ein (n, n)-System herunterzubrechen (Deflation), ist die Alternative, dass ganze System auf ein(n+1, n+1)-System aufzublasen (Inflation).
Der Schritt von einem Punktyk auf dem Ast zu einem Nachfolger wird man dann in einem Zweistufenprozess durchführen:
1.Stufe: Berechne einen Prädiktor-Punkt y0k durch tangentiale Voraussage.
2.Stufe: Wende Iterationen an auf ein erweitertes System an, bis der Ast (bis auf nahezu Rechnergenauigkeit) wieder erreicht ist.
0.5 1 1.5 2 2.5 3 3.5
0 0.5 1 1.5 2 2.5 3
First idea
Corrector−Steps
Prädiktor
Wir vertrödeln keine Zeit mit Nachvollzug der historischen Entwicklung der für die Reali-sierung der Stufen nötigen Rechenschritte, sondern schreiben einfach einige an.
Stufe 1: Sei der Ast beiyknach der Bogenlängetparametrisiert132, so dassF(y(t)) = 0ist und ∥y′(t)∥2 = 1. Die letzte Bedingung folgt aus der Annahme der Parametrisierung nach der Bogenlänge. Der Ast wird mit der Geschwindigkeit 1 durchlaufen. Dann erhält man den Tangenteneinheitsvektor y′(t) einfach aus der durch Differentiation von F(y(t)) = 0 entstehenden Gleichung
F′(y(t))y′(t) = 0.
Ein MATLAB-Befehl zur Lösung dieser Aufgabe ist
T= null(F′(yk)). (202)
null(A) ist der Befehl, der eine Orthonormalbasis des Nullraumes oder Kerns der Matrix A berechnet. T wäre daher sicher schon die Ableitung der Kurve y(t), wenn wir nur sicher sein könnten, dass T neben der korrekten Richtung und der richti-gen Länge auch noch die Orientierung von y′(t) hätte. Die Berechnung von T aus (202) sichert ja nur die ersten beiden Eigenschaften, nicht aber die Orientierung. Es wäre ohne die korrekte Orientierung nicht ausgeschlossen, dass die T-Vektoren in aufeinanderfolgen Schritten dauernd die Richtung wechselten.
Abhängig von der Größe des behandelten Systems verwendet man unterschiedliche Methoden, eine einheitliche Orientierung zu sichern.
• Für kleine bis moderat große Systeme ist es am einfachsten, die Stetigkeit der Matrix
D(t) =
(F′(y(t)) y′(t)T
)
(203) auszunutzen. Da y′(t) einerseits ungleich Null ist und andererseits senkrecht zu den n linear-unabhängigen Zeilen von F′(y(t)), ergänzt y′(t)T die Matrix F′(y(s))zu einer regulären (n+ 1, n+ 1)-Matrix. Da deren Determinante
D(t) := det
(F′(y(t)) y′(t)T
)
(204) ebenfalls stetig ist und überall ungleich Null, hat D(t) ein festes Vorzeichen.
Analog gibt man T eine einheitliche Orientierung, indem man nach seiner Be-rechnung - gegebenfalls durch Vorzeichenänderung - die Bedingung
det
(F′(yk) T(yk)T
)
>0 (205)
sicherstellt. Dies gibt den TangentialvektorenT(yk)in der Punktfolgey0, y1, y2, ...
entlang y(t) eine einheitliche Orientierung entweder in Richtung vony′(t) oder dagegen.
Fordert man in (205) ein negatives Zeichen, wird die Kurve in Gegenrichtung durchlaufen.
132Bitte nicht die „Realisierung der Bogenlängenparametrisierung über die Inversion der streng mono-tonen Bogenlängenfunktion“ aus den Mathematik-Grundvorlseungen aufarbeiten. Hier geht die Tatsache, dassy(t)nach der Begenlänge allein über die Aussage ein, dass∥y′(t)∥2= 1für allet.
• Für größere Dimensionenn ist es sicher weniger rechenaufwendig nach Vorgabe einer Richtung im ersten Astfolge-Punkt y0 die Tangenten in Folgepunkten so auszurichten, dass ihr inneres Produkt mit ihrer Vorgängerin positiv ist. Al-lerdings setzt dies eine gute Schrittweitensteuerung voraus, die nicht zulässt, dass die Kurve zwischen einem Punkt yk und dem nächsten yk+1 die Richtung vollständig wechselt133.
0 1 2 3 4 5 6 7 8 9 10 11
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
yk+1
yk T(yk+1)
T(yk+1)
Winkel < 180°
Orientierung von T(x
k+1) am Vorgänger T(x
k)
Abbildung 54: Festlegung der Tangentialrichtung
Definition 3.38 ( Einheitstangentenvektorfeld)
Es ist hier ein guter Platz um anzumerken, dass durch die Gleichungen (202) und (205) sowie eine Normierung, also durch
F′(y)T(y) = 0, det
(F′(y) T(y)T
)
>0 und ∥T(y)∥= 1 (206) für y∈ R ein Einheitstangentenvektorfelddefiniert wird, welche auf den Trajektorien y(t) aus F(y(t)) = 0 mit deren Tangentialvektoren übereinstimmen.
0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3
Tangentialschritt und Tangentialfeld
Abbildung 55: Tangentialvektorfeld
133Das ist wirklich ein großes praktisches Problem, und die Literatur ist voll von Steuerungsalgorithmen
Stufe 2: Wir nehmen für die Diskussion dieses Teils des Algorithmus an, dass wir von yk auf der Kurve aus einen Schritt in die richtige Tangentenrichtung getan haben, uns in yk+1p außerhalb des Astes befinden134 und sehnlichst wieder zum Ast zurück wollen.
Eine erste ein gute Idee ist es sicher wieder, die Gleichung F(y) = 0 bei yk+1p zu linearisieren
F(y)≈F(yk+1p ) +F′(yk+1p )(y−yk+1p ) = 0
Wenn yk+1p ∈ R, so dass Rang(F′(ypk+1)) = n, ist die Lösungsgesamtheit dieser Glei-chung eine Gerade, aus der noch ein geeigneter Punkt ausgewählt werden muss (vgl.
Abbildung 56).
0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3
yk
yk+1 p
Lösungslinie der Linearisierung
Lauter mögliche Schritte
"Newton"−Schritte zur Lösung der Linearisierung
Abbildung 56: Lösung der Linearisierung
Für die Auswahl sind verschiedene Zusatzgleichungen vorgeschlagen worden. Von Herber B. Keller stammt die Idee, F(y) = 0 durch die Gleichung eines Kreises um den Startpunktyk+ zu ergänzen, und dieses System
F(y) = 0,
∥y−yk∥22−rk = 0 (207)
mit Newton-Iterationen und Startpunkt yk+1p zu lösen (vgl. Abbildung 57, Keller-Schritt).
Von Schwetlick und anderen stammt der Vorschlag, vom Prädiktionspunkte aus senkrecht auf der bisherigen Tangentialrichtung zu laufen. Das läuft bei Start des Prädiktor-Korrektorschrittes beiyk+1in Abbildung 57) auf die Ausführung der Newton-Iteration für
F(y) = 0,
T(yk+1)T(y−yk+2p ) = 0 (208) mit Startpunktyk+2p hinaus.
134Die Indizes bei ypk+1 deuten an, dass dieser Punkt ein Startwert für die Berechnung des nächsten Astpunktesyk+1 ist, der mit einem Prädiktor gewonnen wurde.
−1 −0.5 0 0.5 1 1.5 2 0
0.5 1 1.5 2 2.5 3
yk+2
yk+1
yk
yk→ yk+1, Keller−Schritt yk+1→ yk+2,
T(yk+1)−Orthogonal−Schritt 2 Wege zurück zum Ast
yk+1 p
T(yk+1) yk+2p
Abbildung 57: Festlegung der Tangentialrichtung
Es waren vermutlich Allgower und Georg (vgl. [AG]), die vorgeschlagen haben, in System (208) die Tangentialrichtung T(yk+1) zu ersetzen durch die im Prädiktions-punkt vorliegende Tangentialrichtung T(yk+2p ), und nach Ausführung eines Newton-Schrittes den Tangentenvektor anzupassen.
−0.5 0 0.5 1 1.5 2 2.5 3 3.5
0 0.5 1 1.5 2 2.5 3
Orthogonalschritt, orthogonal zu T am Prädiktorschritt
yk
yk+1
yk+1 p
T(yk+1 p ) Newton, senkrecht
zu T(y
k+1 p )
Abbildung 58: Schritt senkrecht zur Tangente am aktuellen Ort
In Abbildung 56 hieße dies, dass unter allen Schritten der kürzeste Schritt, senk-recht auf die Linearisierungsgerade gewählt würde. Dafür kann man viele positive Argumente sammeln.
1. Im Zweifel sollte man bei der Verwendung der Linearisierung einer Funktion be-müht sein, keine großen Schritte zu machen, da die Linearisierung die wirkliche Funktion umso besser wiedergibt, je weniger man sich vom Entwicklungspunkt entfernt.
2. Anders als bei einem Schritt senkrecht zur Tangentenrichtung im Startpunkt des Prädiktors wird man nicht planmäßig an der Lösungskurve vorbeigeführt.
3. Stattdessen würde dies Verfahren bei genügender Dämpfung vom Prädiktor-punkt aus orthogonal zum Tangentenvektorfeld T(y), y ∈ R auf die Kurve zu laufen (vgl. Abbildung 59 sowie den Abschnitt ??).
0 0.5 1 1.5 2 2.5 3 3.5 4
0 0.5 1 1.5 2 2.5 3
Kontinuierlicher Korrektur−Pfad vom Prädiktorpunkt senkrecht zu Tangentialfeld.
Abbildung 59: Tangentialvektorfeld 3.3.2 Differentialgleichungsberechnung von Lösungsästen Zur Berechnung eines Lösungsastes eines parameterabhängigen Systems
F(x, λ) = 0, F :Rn×R−→Rn (209) (vgl. 200) lässt sich wie für den Spezialfall (183) die Davidenko-Methode (vgl. Seite 96) heranziehen. Wenn F(x0, λ0) = 0 und detFx(x0, λ0)̸= 0, so kann de Gleichung (209) lokal nachx=x(λ) aufgelöst werden und es gilt die Gleichung:
F(x(λ), λ) = 0. (210)
Wenn man (210) nach λ differenziert, ergibt sich
Fx(x(λ), λ)x′(λ) +Fλ(x(λ), λ) = 0.
Bei Regularität von Fx(x(λ), λ) lässt sich die Gleichung nach x′(λ) auflösen, und x(λ) ist damit lokal durch „Allgemeine Davidenko-Anfangswertaufgabe“ bestimmt:
Fx(x(λ), λ)x′(λ) = −Fx(x(λ), λ)−1Fλ(x(λ), λ), x(λ0) = x0. (211) Leider lässt sich die Lösungs-Kurve mit dieser Metheode nicht um λ-Umkehrpunkte im Regularitätsgebiet (201) herumsteuern.
Hier ist es günstiger, das Einheitstangentenvektorfeld T(y) nach (206) für F(y) = 0, y := (x, λ)
heranzuziehen.
Wenn y(t) die Bogenlängenparametrisierung der Lösungskurve mit y(0) = (x0, λ0) ist, so ist
y′(t) =T(y(t)), y(0) = (x0, λ0) (212)
eine Anfangswertaufgabe, mit der sich der Lösungsast solange verfolgen lässt, wie in in der Regularitätsmenge bleibt135. Obwohl sich dies relativ kompliziert anhört, ist es sehr einfach zu programmieren:
function r e s = i m p l i c i t ( t , y , f u n )
% Annahme : S e i f u n i n C^2(R^n ,R^(n−1)) and Rang ( J a c o b i a n ( f u n ) ) = n−1
% im ( k o n v e x e n und o f f e n e n ) G e b i e t M i n d e r R e g u l a r i t ä t s m e n g e R von f u n .
% Wenn e s e i n y∗ i n M mit f u n ( y∗)=0 g i b t , d e f i n i e r t f u n ( y)= 0 e i n e Kurve
% i n M d u r c h y∗.
%
% 1 . Das v o r l i e g e n d e Programm d e f i n i e r t e i n V e k t o r f e l d , d a s a l s r e c h t e
% S e i t e e i n e r MATLAB−A n f a n g s w e r t a u f g a b e n l ö s e r b e n u t z t werden kann .
%( Die ODE−S u i t e A u f r u f k o n v e n t i o n s i n d e r f ü l l t ) .
% Wird " i m p l i c i t " mit U e b e r g a b e d e s Funktionsnamens a u f " f u n " b e i y0 mit
% f u n ( y0 ) =0 g e s t a r t e t , s o f o l g t d e r ODE−I n t e g r a t o r d e r d u r c h
% f u n ( y)=0 d e f i n i e r t e n Kurve .
% B e i s p i e l
% f u n c t i o n r e s = c i r c l e ( y )
% % d e f i n e s u n i t c i r c l e i n R^2
% r e s = y (1)^2+ y (2)^2−1;
%
% T=10; y0 = [ 1 , 0 ] ’ ; [ t , y ] = ode45 (@( t , y ) i m p l i c i t ( t , y , @ c i r c l e ) , [ 0 ,T ] , y0 )
%
% 2 . I s t f u n ( y0 ) \ ne 0 , s o f o l g t d e r s e l b e A u f r u f
% T=10; y0 = [ 2 , 0 ] ’ ; [ t , y ] = ode45 (@( t , y ) i m p l i c i t ( t , y , @ c i r c l e ) , [ 0 ,T ] , y0 )
% d e r " p a r a l l e l e n Kurve " d e r Punkte , d i e f u n ( y ) = f u n ( y0 ) e r f ü l l e n .
% Im A n w e n d u n g s b e i s p i e l i s t d i e d e r K r e i s um den Ursprung mit Radius 2 .
% two .
[ J , F ] = Jakob ( fun , y ) ; %" J a ko b " b e r e c h n e t d i e J a c o b i m a t r i x "J"
% und den F u n k t i o n s w e r t "F"
n=n u l l( J ) ; % " n u l l " b e r e c h n e t e i n e O r t h o n o r m a l b a s i s
% d e s Kerns von J .
% D i e s i s t d i e T a n g e n t i a l r i c h t u n g wenn Rang ( J)=n−1.
% Man b e a c h t e , d a s s d e r V e k t o r s c h o n e i n h e i t l ä n g e
% h a t .
i f det( [ J ; n ’ ] ) < 0 % Wenn d i e um n^T e r g ä n z t e J a c o b i s c h e k e i n e p o s i t i v e n=−n ; % D e t e r m i n a n t e h a t , d r e h t man d i e R i c h t u n g von n um . end % [ Man s i c h damit a u f e i n e D u r c h l a u f r i c h t u n g f e s t . ]
r e s =n ;
Um unser Programm zu testen, verwenden wir die Gleichung für den Einheitskreis.
function r e s= c i r c l e ( z ) r e s=z (1)^2+ z (2)^2−1;
Der Aufruf kann so aussehen:
x0 = [ 0 ; 1 ] ;
o p t s=o d e s e t ( ’ R e l T o l ’ , 1 e−2 ) ; % D i e s e Option , d i e w i r g l e i c h b e i d e r z w e i−
% t e n I n t e g r a t i o n a n s c h a l t e n , r e d u z i e r t d i e s t a n d a r d m a e s s i g e i n g e s t e l l t e
% I n t e g r a t i o n s g e n a u i g k e i t , s o d a s s w i r g l e i c h d i e Tendenz z u r Abweichung
% von d e r e x a k t e n Loesung bemerken koennen .
[ t , y]=ode45(@( t , y ) i m p l i c i t ( t , y , @ c i r c l e ) , [ 0 , 2 0∗pi] , x0 , o p t s ) ;
plot( y ( : , 1 ) , y ( : , 2 ) ) ; a x i s e q u a l ;
135Dem kann dadurch ein Ende gesetzt werden, dass der Definitionsbereich von F verlassen wird oder
Anders als bei der Davidenko-Gleichung, die eine Parametrisierung des Kreises nach x oder nach y verlangen würde, kann der Integrator in der Kurvenlänge t den Kreis be-liebig häufig umfahren. (Dass der Kreis durch t Bogenlängenparametrisierung erhält, wir einfach dadurch geregelt, dass der in der Differentialgleichung verwendet Richtungsvektor Einheitlänge bekommen hat. Dass bedeutet doch schicht nur, dass die Integralkurve mit Geschwindigkeit 1 durchlaufen wird.
−1 −0.5 0 0.5 1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
Abbildung 60: Ganzer - sauberer - Kreis
Man beachte, dass die letzte Abbildung mit hoher Integrationsgenauigkeit berechnet wurde.
Senken wir diese ab (wie das die Optionen „Ops“ im Programm andeuten), verlassen wir die Lösungskurve:
−1.5 −1 −0.5 0 0.5 1 1.5
−1
−0.5 0 0.5 1
Integration with low integration quality leads to leaving the curve
Abbildung 61: Low quality integration
Das kann man - bei beibehaltener grober Integration - durch einen „Schuss Newton“ ver-hindern. Man addiert zum Tangentialeinheitsvektor, der in Richtung der Kurve zeigt, den-jenigen Vektor δ⊥ unter den Vektoren δ, die die Linearisierung
F(x0) +F′(x0)δ= 0
lösen, der zum Tangentialvektor T(x0) senkrecht steht, der also das System F′(x0)δ⊥(x0) = −F(x0),
T(x0)Tδ⊥(x0) = 0, (213)
löst.
function r e s = i m p l i c i t 2 ( t , y , fun , a l p h a )
% Annahme : S e i f u n i n C^2(R^n ,R^(n−1)) and Rang ( J a c o b i a n ( f u n ) ) = n−1
% im ( k o n v e x e n und o f f e n e n ) G e b i e t M i n d e r R e g u l a r i t ä t s m e n g e R von f u n .
% Wenn e s e i n y∗ i n M mit f u n ( y∗)=0 g i b t , d e f i n i e r t f u n ( y)= 0 e i n e Kurve
% i n M d u r c h y∗.
%
% 1 . Für a l p h a =0 l i e f e r t d a s Programm d a s s e l b e E r g e b n i s w i e i m p l i c i t .m
% A n w e n d u n g s b e i s p i e l
% f u n c t i o n r e s = c i r c l e ( y )
% % d e f i n e s u n i t c i r c l e i n R^2
% r e s = y (1)^2+ y ( 2 ) ^ 2 −1;
%
% T=10; y0 = [ 1 , 0 ] ’ ; [ t , y ] = . . .
% ode45 (@( t , y ) i m p l i c i t 2 ( t , y , @ c i r c l e , 0 ) , [ 0 , T ] , y0 ) ;
%
% 2 . I s t f u n ( y0 ) \ ne 0 and a l p h a =0, s o l i e f e r t
% T=10; y0 = [ 2 , 0 ] ’ ; [ t , y ] = . . .
% ode45 (@( t , y ) i m p l i c i t 2 ( t , y , @ c i r c l e , 0 ) , [ 0 , T ] , y0 ) ;
% w i e d e r d i e P a r a l l e l−Kurve f u n ( y ) = f u n ( y0 ) .
%
%3 . I s t f u n ( y0 ) \ ne 0 a b e r a l p h a >0, s o e r h ä l t d i e I n t e g r a t i o n
% T=10; y0 = [ 2 , 0 ] ’ ; [ t , y ] = ode45 ( @ i m p l i c i t 2 , [ 0 , T] , y0 , ’ ’ , @ c i r c l e , a l p h a )
% e i n e Newton−A u s r i c h t u n g a u f d i e Kurve f u n ( x )=0 hin , w e l c h e mit
% wachsendem a l p h a i n t e n s i v e r w i r d .
%
% Bemerkung : Um zu v e r h i n d e r n , d a s s n u m e r i s c h e F e h l e r dazu f ü h r e n , d i e Kurve
% zu v e r l a s s e n , i s t e s e m p f o h l e n , s t e t s i m p l i c i t 2 .m mit einem k l e i n e n
% p o s i t i v e m a l p h a zu verwenden .
% N e g a t i v e a l p h a−Werte s i n d n i c h t s i n n v o l l , denn s i e t r e i b e n den
% I n t e g r a t o r von d e r Kurve f o r t .
% Zu gro�e p o s i t v e a l p h a−Werte s i n d auch k o n t r a p r o d u k t i v , w e i l s i e d i e
% D i f f e r e n t i a l g l e i c h u n g s t e i f machen . [ J , F ] = Jakob ( fun , y ) ;
N=n u l l( J ) ; JJ =[ J ; N ’ ] ;
i f det( JJ)<0 N=−N;
end
JJ =[ J ; N ’ ] ;
r e s = JJ \[−a l p h a∗F ; 1 ] ;
% Anmerkung
% r e s i s t e i n e Ü b e r l a g e r u n g r e s = r e s 1+a b s ( a l p h a )∗r e s 2 von
% r e s 1=N, d a s d i e Lösung von
% JJ r e s 1 = [ z e r o s ( s i z e (F ) ) ; 1 ]
% i s t , s o w i e von einem Newton−Typ−S c h r i t t , d e r " s e n k r e c h t zum
% T a n g e n t i a l v e k t o r N" z u r L ö s u n g s m a n n i g f a l t i g k e i t von f u n ( y)=0
% z u r ü c k t r e i b t .
%
% JJ r e s 2 = [−F ; 0 ]
%
% Achtung r e s 2 e r f ü l l t
% J r e s 2 = −F ,
% a l s o d i e L i n e a r i s i e r u n g von f u n ( y)=0 i n y und
% N’ r e s 2 = 0 ,
% w e l c h e G l e i c h u n g J zu e i n e r q u a d r a t i s c h e n r e g u l ä r e n M a t r i x macht .
% Die G l e i c h u n g N’ r e s 2 =0 s u c h t u n t e r a l l e n Lösungen von
% J r e s 2 = −F d i e j e n i g e aus , d i e o r t h o g o n a l zum T a n g e n t e n v e k t o r f e l d
Hier ist der Aufruf.
x0 = [ 0 ; 2 ] ;
o p t s=o d e s e t ( ’ R e l T o l ’ , 1 e−2 ) ; %Reduce e x a c t n e s s o f i n t e g r a t i o n [ t , y]=ode45(@( t , y ) i m p l i c i t 2 ( t , y , @ c i r c l e , 0 . 0 3 ) , [ 0 , 4 0∗pi] , x0 , o p t s ) ;
plot( y ( : , 1 ) , y ( : , 2 ) ) ; a x i s e q u a l ;
% Man b e a c h t e , d a s s d e r K r e i s auch mir r e d u z i e r t e r
% I n t e g r a t i o n s g e n a u i g k e i t n i c h t v e r l a s s e n w i r d .
Anmerkungen 3.39 (Newton-Fluss zum Lösungsast)
1. Analog zum Newton-Fluss (187) zur Lösung eines nichtlinearen Systems führt uns die Lösung der Differentialgleichung
x′(t) =δ⊥(x(t)) (214)
mit δ⊥(x(t)) aus (213) „schnellstens“ zurück zum Lösungsast.
In Abbildung 59 ist dies die rote Kurve.
2. Die lineare Abbildung von −F(x0) auf δ⊥(x0) kann auch mit der Pseudoinversen F′(x0)† geschrieben werden als
δ⊥(x0) =−F′(x0)†F(x0). (215) 3.3.3 Berechnung von Umkehrpunkten:
Im Abschnitt 3.2.4 über Homotopien hatten wir schon ausgeführt, dass man in der An-fangszeit der Behandlung parameterabhängiger Gleichungssysteme (200), also
F(x, λ) = 0, F :D×I −→Rn, D ⊂Rn (216) einen Lösungsast durch sukzessive Vergrößerung vonλ mit iterativer Anpassung des zuge-hörigenx-Vektors zu generieren versuchte (vgl. Abbildung 34).
Dort hatten wir auch schon darauf hingewiesen, dass diese Methode Probleme bekommt, wenn der Lösungsast in Bezug auf den Parameterλbei einem Wertλ∗seine Richtung wech-selt (vgl. Abbildung 35). Abgesehen davon, dass solche Umkehrpunkte für die Verfolgung des Astes von Belang schienen136, war natürlich für das dem Modell zugrundeliegende Pro-blem sehr wichtig, dass das Modell voraussagte, dass bei Überschreiten des Parametersλ∗ (sagen wir nach rechts) auch für das beschriebene Problem möglicherweise keine „Lösung“
mehr existierte137. Weil sich ein solches ‘„Aufhören der Existenz“ z.B. in einer explosi-onsartigen Auflösung bestehen könnte138, ist der damals gebräuchliche Name „kritischer Parameter“ für λ∗ relativ einleuchtend.
Um nun herauszubekommen, an welchenλ-Stellen das modellierte System kritisch würde, war es deshalb von großem Interesse, Algorithmen zu bestimmen, mit denen Umkehrpunkte einer implizit definierten Kurve bestimmt werden konnten.
136Hier hatte man die Veränderung vonλauch umzukehren, wenn man dem Ast weiter folgen wollte.
137Jedenfalls keine in der Nähe gelegene, in die das System ohne große und vermutlich plötzliche -Veränderung überwechslen könnte.
138Man hat solche mathematischen Aufgaben damals u.a. im Zusammenhang mit der Beschreibung che-mischer oder auch nuklearer Reaktoren verwendet.
Ende der siebziger Jahre und zu Anfang der 80er Jahre des letzten Jahrhunderts wurde hierfür sehr viele solcher Methoden vorgeschlagen. Diese kann man unterteilen in „ast-gebundene Methoden“ und in „direkte Methoden“.
Die „ast-gebundenen“ Methoden verwendeten eine lokale Parametrisierung des Astes (x(t), λ(t)), t ∈(a, b)
(wobei der Parametertmeist als eine der Komponenten vonxgewählt wurde), und stellten über diese Funktionen Bedingungen an t auf, die den Umkehrpunkt charakterisierten.
λ′(t) = 0 (217)
ist eine solche Bedingung, und Methoden, den Punkt zu gewinnen, erwuchsen aus der Anwendung von Nullstellenverfahren auf das Problem λ′(t) = 0. Nach der Bestimmung eines Parameters t∗, bei demλ′(t∗) = 0 war, ergab sich der Umkehrpunkt selbst als
(x∗, λ∗) := (x(t∗), λ(t∗)).
Da alle Verfahren zur eindimensionalen Nullstellenberechnung nur wirklich gut funktionie-ren, wenn die Ableitung der auf Null zu bringenden Funktion in der Nullstelle von Null verschieden ist, forderte man, dass
λ′′(t∗)̸= 0 (218)
sein möge, und nannte einen Umkehrpunkt, der (218) erfüllte, eineneinfachen Umkehr-punkt. Neben ihrer Bedeutung für die numerische Berechnung einer Nullstelle von λ′(t) sichert sie zugleich auch, dass die Lösungskurve in(x∗, λ∗)auch wirklich umkehrt und nicht nur eine Tangente senkrecht zur λ-Richtung hat:
Beispiel 3.40 Bei
Fa(x, λ) := λ−x3
ist t := x ein geeigneter Parameter um die Lösungskurve in der Nähe des Punktes (0,0) gemäß
(x(t), λ(t) =:= (t, t3)
zu parametrisieren. t∗ := 0 erfüllt sicher (217). Aber (siehe Abbildung 62), die Kurve kehrt bezüglichλ nicht um.
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
λ
x
Umkehrpunkt−Aspirant
Fb(x, λ) :=λ−x2
hat dagegen bei (0,0)einen einfachen λ-Umkehrpunkt.
−1 −0.8 −0.6 −0.4 −0.2 0 0.2
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
λ
x
λ+x2 =0
Einfacher Umkehrpunkt
Abbildung 63: Einfacher Umkehrpunkt
Obwohl das „eindimensionale Rechnen auf einem parametrisierten Ast“ im Prinzip ein einfaches Konzept ist, hat es nicht ganz befriedigt, weil man bei seiner Anwendung die benötigten Werte von (x(t), λ(t)) relativ genau berechnen muss.
Mit den direkten Methoden versuchte man deshalb, ein den Umkehrpunkt (x∗, λ∗)∈ R(n+1) selbst charakterisierendes Gleichungssystem aus n+ 1Gleichungen aufzustellen, zu dessen Lösung man sich nicht auf den Ast zurückziehen muss.
Dazu werden auf die eine oder andere Weise, dien Gleichungen in (216), die den Lösungs-ast festlegen um eine weitere Gleichung ergänzt, die eine charakteristische Qualität des Umkehrpunktes ausdrückt.
Eine der in der Geschichte der Bestimmung von Umkehrpunkten ersten solcher Gleichungen
detFx(x, λ) = 0 (219)
machte von der Information Gebrauch, dass im Umkehrpunkt die Fortsetzbarkeit des Astes in λ-Richtung scheitert, weil die Jacobimatrix bezüglich der Zustandsvariable x singulär wird. Wir demonstrieren das Zusammenwirken von (216) und (219) im folgenden
Beispiel 3.41
Das Gleichungssystem (220) ist eine sehr grobe Diskretisierung des Bratupro-blemes (27) mit nur zwei Diskretisierungspunkten:
2x1 −x2 −λex1 = 0,
−x1 +2x2 −λex2 = 0. (220)
Die Projektion der Lösungsäste auf die (λ, x1)-Ebene ergibt die blauen Linien in der folgenden Abbildung 64.
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0
1 2 3 4 5 6
2D−Bratu
λ
x
Mannigfaltigkeiten singulärer Jacobischer
Umkehrpunkt Verzweigungspunkt
Abbildung 64: Umkehr- und Verzweigungspunkt Die roten Linien sind die Schnitte der Lösungsmannigfaltigkeiten von
det
(2−λex1 −1
−1 2−λex2 )
= 0. (221)
mit der(x1 =x2)-Ebene . Tatsächlich sehen diese Flächen aus, wie in Abbildung 65 dargestellt139
0 2
4
6 0 1 2 3 4 5 6
−1 0 1 2 3 4 5
Y
Vollständige Mannigfaltigkeiten singulärer Jacobischer
X
λ
Abbildung 65: Mannigfaltigkeitem singulärer Jacobischer Diese Flächen schneiden den aus (0,0,0)T kommenden Lösungsast in
(x∗a, λ∗a) = ((1,1), e−1) und (x∗b, λ∗b) = ((3,3),3e−3).
Die Jacobimatrizen (Fx, Fλ) lauten in diesen Punkten Ja =
( 1 −1 −e
−1 1 −e )
bzw. Jb =
(−1 −1 −e3
−1 −1 −e3 )
.
Bei beiden ist der Rang des x-Anteils Fx nur 1, die Determinante ist Null.
Es ist aber Rang(Ja) immer noch zwei, so dass nach dem Satz über implizi-te Funktionen durch (x∗a, λ∗a) = ((1,1), e−1) ein eindimensionaler Lösungsast verläuft.
Der Rang vonJb dagegen bleibt gleich 1, so dass der Satz über Implizite Funk-tionen hier keine eindeutige Fortsetzung garantiert. Tatsächlich ist (wie ja auch die Abbildung 64 zeigt) der Punkt (x∗b, λ∗b) = ((3,3),3e−3) ein Verzweigungs-punkt, durch den (in diesem Fall) zwei Äste laufen140.
Bildet man die Gesamtjacobischen des um (219) erweiterten Systems(216), so ergeben sich in den beiden Punkten die erweiterten quadratischen Jacobi-Matrizen
Ja,plus =
1 −1 −e
−1 1 −e
−1 −1 −12e
bzw. Jb,plus =
−1 −1 −e3
−1 −1 −e3 3 3 2e3
.
Die Tatsache, dass Ja,plus regulär ist, macht uns sicher, dass der Umkehrpunkt (x∗a, λ∗a) durch Anwendung der Newton-Iteration direkt auf das System (216, 219) berechnen kann141.
Diese Beobachtung ist allgemein gültig. Ohne den Beweis anzutreten, halten wir hier fest:
Lemma 3.42
Sei(x∗, λ∗)einfacher Umkehrpunkt der Lösungsgesamtheit von (216). SeiF in einer offenen Umgebung von (x∗, λ∗) zweimal stetig differenzierbar mit Lipschitz-stetiger zweiter Ablei-tung. Dann konvergiert das Newton-Verfahren in Anwendung auf das um (219) erweiterte nichtlineare Gleichungssystem
F(x, λ) = 0,
detFx(x, λ) = 0 (222)
lokal quadratische gegen den Umkehrpunkt.
Es gibt weitere Charakterisierungen der Singularität vonFx(x, λ), und entsprechend hat es auch weitere Methoden gegeben, Umkehrpunkte aufzuspüren. Unter dieast-gebundenen Methoden fällt z.B. die Beobachtung der Eigenwerte von Fx(x(t), λ(t)). Hier verraten sich kritische Punkte durch einen Nulleigenwert. Diese Methode bietet sich dann an, wenn die Eigenwerte bei der Astverfolgung sowieso überwacht werden142.
Als Grundlage für eine direkte Iteration bietet sich das erweiterte System F(x, λ) = 0,
kleinster Eigenwert von(Fx(x, λ)) = 0 (223)
140Wie man diese mathematisch und ingenieurwissenschaftlich ebenfalls sehr interessanten Punkte be-stimmt, wird in diesem Skript nicht erläutert. Es sei aber darauf hingewiesen, dass Ihnen allen aus der Mechanik ein solches Verzweigungsphänomen in Form des Problemes der „Eulerschen Knicklast“ schon begegnet ist.
141Für die Berechnung von Verzweigungspunkten musste man sich etwas Raffinierteres einfallen lassen, dessen Darstellung den Umfang dieser Vorlesung zu groß machen würde.
142Das macht man zum Beispiel dann, wenn man wissen will, ob eine stationäre Lösung der Differential-gleichungx˙ =F(x, λ), also eine Lösung vonF(x, λ) = 0, stabil unter der durch diese Differentialgleichung beschriebenen Zeitentwicklung ist.
nicht an, da die Erweiterungsfunktion noch schwieriger abzuleiten ist als die Determinante.
Ziehlt man statt auf die Existenz eines Nulleigenwertes auf die eines Eigenvektor zu einem Eigenwert Null, so ergibt sich das vielbenutzte System von Moore und Spence.
Lemma 3.43 (Moore-Spence-Erweiterung)
Sei (x∗, λ∗) einfacher Umkehrpunkt der Lösungsgesamtheit von (216), sei Φ∗ der mit einem geeigneten Vektor r ∈ Rn durch rTΦ∗ = 1 normierte Eigenvektor von Fx(x∗, λ∗) zum Eigenwert Null. Dann hat das erweiterte quadratische nichtlineare Gleichungssystem
F(x, λ) = 0, Fx(x, λ)Φ = 0
rTΦ = 1
(224)
inx,Φ∈Rn und λ∈R im Punkt (x∗, λ∗,Φ∗)eine reguläre Ableitungsmatrix, so dass das Newtonverfahren zur Berechnung der Lösung mit lokal quadratischer Konvergenz herange-zogen werden kann, wenn F in einer Umgebung des gewünschten Punktes eine Lipschitz-stetige zweite Ableitung hat.
Anmerkungen 3.44
Einen geeigneten Vektor r erhält man in der Nähe eine Umkehrpunktes z.B. als Differenz zweier aufeinanderfolgender x-Werte von Punkten auf dem Lösungsast, da der Tangential-vektor T(x∗, λ∗) ∈Rn+1 an die Kurve im Umkehrpunkt eine λ-Komponente Null hat. Da dieser Vektor
(Fx(x∗, λ∗), Fλ(x∗, λ∗))T = 0
erfüllt, ist der x-Anteil vonT ein Eigenvektor von Fx(x∗, λ∗)zum Eigenwert Null.
Diese Überlegung führt direkt zur Verwendung des Einheitstangentenvektorfeldes von Seite 118.
Indem wir wie dort die Komponenten x, λ zu einem Gesamtvektor (x, λ) =: y ∈ Rn+1 zusammenfassen, können wir (wieder ohne Beweis) schreiben.
Lemma 3.45 (Erweitertes System für einfache Umkehrpunkte) Sei (y∗) einfacher Umkehrpunkt der Lösungsgesamtheit von
F(y) = 0, F :Rn+1 −→Rn
bezüglich der (n+ 1)-ten Komponente vony(vormals λ). HabeF in einer Umgebung von y∗ ein Lipschitzstetige zweite Ableitung. Dann ist
F(y) = 0,
T(y)n+1 = 0 (225)
ein erweitertes System, dessen lokal eindeutige Lösungy∗ durch das Newton-Verfahren mit lokal quadratischer Konvergenz approximiert werden kann.
Die Beobachtung, dass T(y)n+1 = eTn+1T(y) geschrieben werden kann, führt leicht zur Definition eines L-Umkehrpunktes für eine beliebigen von Null verschiedenen Vektor in Rn+1. Wir werden einen Punkty∗dann einenL-Umkehrpunkt der Lösungsgesamtheit
Beispiel 3.46
So sind die roten Punkte des Kreises in Abbildung 3.49 L1-Umkehrpunkte für L1 :=
(1 1
)
und die blauen L2-Umkehrpunkte für L2 :=
(0 1
) .
−3 −2 −1 0 1 2 3
−2.5
−2
−1.5
−1
−0.5 0 0.5 1 1.5 2 2.5
L1−Umkehrpunkte T
L1
L2 L2
−Umkehrpunkt
Abbildung 66: L-Umkehrpunkte
Definition 3.47 (Einfacher L-Umkehr-Punkt)
1. y∗ in der Regularitätsmenge von F ∈ C1(Rn+1,Rn) ist L-Umkehrpunkt zu L ∈ Rn+1\ {0}, wenn
F(y∗) = 0,
LTT(y∗) = 0 (226)
mit dem Einheitstangentenvektorfeld T(y).
2. Für F ∈C2(Rn+1,Rn) ist mit der Jacobimatix DT(y) von T(y) das Vektorfeld C(y) :=DT(y)T(y)
das sogenannte Krümmungsvektorfeld.
3. Ein L-Umkehrpunkt y∗ heißt einfacher L-Umkehrpunkt, wenn LTC(y∗)̸= 0.
Anmerkungen 3.48
Für eine Trajektorie zum Einheitstangentenvektorfeld T(y) ist C(y) das, was in der Dif-ferentialgeometrie „Hauptnormalenvektor“ oder „Krümmungsvektor“ einer Kurve genannt wird.
Lemma 3.49 ( Berechnung von Umkehrpunkten)
Sei y∗ ein einfacher L-Umkehrpunkt einer Funktion F ∈ C2,1(Rn+1,Rn). Dann ist die Jacobimatrix von
F(y) = 0, LTT(y) = 0
iny∗ regulär, so dass das Newton-Verfahren lokal quadratisch gegeny∗ konvergiert.
Beispiel 3.50
Als Beispiel rechnen wir die oben definierten Größen aus, die bei der Berech-nung der L-Umkehrpunkte (L= (1,1)T) des Einheitskreises
F(x, y) = x2+y2−1
auftreten, dabei wählen wir für den Urbildbereich nicht die Kooerdinatey1 und y2 sondern (x, y).
Es ist
1. F′(x, y) = (2x,2y) und daher 2. T(x, y) = √ 1
x2+y2
(−y x
) .
3. Die Jacobimatrix hiervon ist DT(x, y) = (x2+y2)−3/2
(xy −x2 y2 −xy
)
4. C(x, y) = DT(x, y)T(x, y) =−(x2+y1 2) (x
y )
5. LT(x, y) = (x−y)(x2+y2)−1/2
6. Das erweiterte System x2+y2−1 = 0, (x−y)(x2+y2)−1/2 = 0 7. Die Newton-Gleichung für den Schrittvektor N(x, y):
( 2x 2y xy+y2 −(x2+xy)
)
N(x, y) =−
( x2+y2−1 (x−y)(x2+y2)
)
Zur Veranschaulichung zeigt die nächste Abbildung 67 Trajektoren des
Newton-Flusses (
x′(t) y′(t)
)
=N(x, y).
−1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5
−1
−0.5 0 0.5 1 1.5 2 2.5
Newton−Fluss in die (1,1)−Umkehrpunkte des Einheitskreises
Abbildung 67: Newton-Fluss in (1,1)-Umkehrpunkte
3.3.4 Berandete lineare Systeme
Bei der Behandlung parameterabhängiger Probleme treten als Linearisierungen typischer-weise sogenannte berandete lineare Gleichungssysteme auf. Dies sind Gleichungssteme der
Form (
A b cT d
) (x y
)
= (f
g )
, (227)
wobeiA∈Rn×nregulär ist, b, c, f ∈Rnund d und g ∈Rgegeben undx∈Rn sowiey ∈R zu berechnen sind.
Obwohl man berandete Matrizen auch in anderen Zusammenhängen nützlich findet143, ist ihr Auftreten bei parameterabhängigen Problemen besonders natürlich.
Wir waren ja davon ausgegangen, dass
F(x, λ) = 0, F :Rn×R−→R
für das Studium der Parameterabhängigkeit eines ursprünglich parameterfreien Systems F(x) = 0, F :Rn−→R
entstanden war144.
Wenn die Lösung x0 dazu mit einem Newton-Verfahren berechnet wurde, so war die Jacobimatrix A := F′(x0) sicher regulär. Für einen Newton-Schritte musste ein n × n-Gleichungssystem
Ax=b
gelöst werden. In echten Anwendungsproblemen145 haben Jacobimatrizen fast immer spe-zielle Strukturen, die man bei der Lösung ausnutzt146. Daher kann man annehmen, dass man einen gut ausgearbeiteten Löser für Ax = b hergestellt hat, wenn man sich an das Studium der Parameterabhängigkeit macht. Die Anwendung dieses Lösers schreiben wir abkürzend als Multiplikation mitA−1. Schon in der Linearen Algebra I wurde darauf hin-gewiesen, dass inverse Matrizen tatsächlich nur ganz selten wirklich berechnet werden, und dass das Auftreten vonA−1 in einer Formel meistens zu lesen sei als: „Hier wird ein lineares Gleichungssystem mit der Systemmatrix A gelöst.“ So auch hier.
Die Regularität von A = F′(x0) = Fx(x0, λ0) gewährleistete über den Satz für implizite Funktionen die Fortsetzbarkeit der Lösung x(λ0) = x0 in λ-Richtung. Für die Fortset-zung wurde die Tangente T benötigt, und um diese zu erhalten, brauchte man die ganze Ableitung von F(x, λ).
(Fx(x0, λ0), Fλ(x0, λ0))T = 0.
Dies ist mit b := Fλ(x0, λ0) die erste Block-Zeile des berandeten Systems (227). Indem man Tn+1 gleich 1 setzt147, wird das System zu
Fx(x0, λ0)Tx =−Fλ(x0, λ0)
und man kann den vorhandenen effizienten Löser A−1 direkt einsetzen.
Wenn man sich auf der Tangente ein Stück weit bewegt hat, gilt es zum Lösungsast zurück-zukommen. Nun tritt zum durch die Variableλ„nach rechts ergänzten System“F(x, λ) = 0
143Vgl. rekursive LDLT-Zerlegung etc.
144Genauer war der Parameterλ=λ0hierin noch fest und daher nicht auszeichnungswürdig gewesen.
145Etwa aus den Ingenieurwissenschaften
146Wir werden unten einige typische Fälle ansprechen.
147Was man entfernt von einemλ-Umkehrpunkt tun kann, da der Tangentenvektor einen nicht verschwin-dendenλ-Anteil hat.