5 Krylov-Typ Rosenbrock-Methoden für differential-algebraische Systeme vom Index 1

26 Pages·2004·0.21 MB·German

Kapitel 5
Krylov-Typ Rosenbrock-Methoden für
differential-algebraische Systeme vom
Index 1
5.1 Steife Anfangswertprobleme
5.1.1 Steifheit
Bei vielen Anwendungen treten sogenannte steife Differentialgleichungen auf. Betrachten wir eine
Anfangswertaufgabe
y0(t) =f(y(t)); y(t) Rny;t [t ;t ]
0 e
2 2 (5.1)
y(t ) =y :
0 0
EigentlichistbereitsderBegriffsteifeDifferentialgleichungungenau,dennSteifheitstellteinPhäno-
men dar, das bei der numerischen Lösung von Anfangswertaufgaben auftritt. So hängt die Steifheit
nebender Differentialgleichungauchvon der gefordertenGenauigkeit in der Lösung,der Längedes
Integrationsintervalls sowieunterUmständenauchvondengegebenenAnfangswertenab.
Steifheit bedeutet, daß in der Lösung neben langsam veränderlichen Komponenten schnell ge-
dämpfteKomponentenauftreten,diefürdasProblemnichtvonInteressesind.SteifeProblemetreten
meistimZusammenhangmitsichsehrschnelleinstellendenGleichgewichtenauf.SozumBeispielin
der Reaktionskinetik,wenn langsameReaktionenzusammen mit schnellenReaktioneneinhergehen,
oderimFallevonAusgleichsprozessenwieDiffusionoderWärmeleitung.
SchnellgedämpfteLösungskomponenten entsprechenEigenwerteninderJacobi-Matrixmitgro-
ßemnegativenRealteil,EigenwertemitgroßemImaginärteilführenzuschnellenOszillationen.Beide
ProblemklassenführenzueinergroßenLipschitz-KonstantenLderrechtenSeitedesSystems.
GrundlegendfürdieÜbertragungderKonvergenzanalyse(Schrittweiteh 0)fürexpliziteRun-
!
ge-Kutta-Verfahren auf endliche Schrittweiten h ist die Annahme, daß hL = (1) gilt, denn unter
O
dieserVoraussetzungsinddieVerfahrenstabil.FürProblememitgroßerLipschitz-KonstantenL,die
abertrotzdemgutartiggestelltsind,sindexpliziteVerfahrennichteffektiv,dahL = (1)dieSchritt-
O
weiten h stark einschränkt. Man greift daher auf spezielle problemangepaßte Verfahren zurück, im
Falle steifer Differentialgleichungen sind dies implizite Runge-Kutta-Verfahren und linear-implizite
Runge-Kutta-Verfahren, sowie die hier nicht näher betrachteten BDF-Methoden als Mehrschrittver-
fahren.AusführlicheDarstellungenundAnalysengeeigneterVerfahrenfindetmanin[24],[46],[99].
103
5.KRYLOV-ROW-METHODEN
5.1.2 StabilitätvonEinschrittverfahren
Bei der Behandlung steifer Differentialgleichungen sind implizite Verfahren trotz des höheren Auf-
wands pro Integrationsschritt den expliziten Verfahren überlegen, da sie wesentlich größere Schritt-
weitenzulassen.DiesliegtandenschlechterenStabilitätseigenschaften expliziterVerfahren,wieman
bereitsandereinfachenlinearenTestgleichung
y0 =(cid:21)y; (cid:21) C (5.2)
2
nachvollziehenkann.BeiAnwendungeinesRunge-Kutta-Verfahrens aufdieGleichung(5.2)mitder
SchrittweitehergibtsicheineRekursionderForm
y =R(h(cid:21))y mit
m+1 m
R(z) =1+bTz(I zA)(cid:0)11l:
(cid:0)
DieFunktionR(z)istrational,fürexpliziteVerfahrenistsieeinPolynom.DielineareGleichung(5.2)
istfürRe(cid:21) 0kontraktiv, d.h.,fürzweiLösungeny (t);y (t)gilt
1 2
(cid:20)
y (t+h) y (t+h) y (t) y (t) fürh 0: (5.3)
1 2 1 2
j (cid:0) j (cid:20) j (cid:0) j (cid:21)
BeilinearenSystemenistdiesäquivalentzu y(t+h) y(t),unddiesfordertmanauchvoneinem
j j (cid:20) j j
Diskretisierungsverfahren, imIdealfallalso R(z) 1fürRez 0.
j j (cid:20) (cid:20)
Die Menge z C : R(z) 1 bezeichnet man als Stabilitätsgebiet. Verfahren, deren Sta-
f 2 j j (cid:20) g
bilitätsgebiet die linke komplexe Halbebene (hier ist gerade die lineare Testgleichung (5.2) stabil)
umfaßt, heißen A-stabil. Bei Stabilität in einem Sektor der Größe 2(cid:11) um die negative reelle Achse
spricht man von A((cid:11))-Stabilität, und liegt Stabilität nur auf der reellen Achse vor, spricht man von
A -Stabilität. A-stabile Verfahren mit R( ) = 0 bezeichnet man als L-stabil. Schärfere Stabilitäts-
0
1
begriffe legen nichtautonome lineare Systeme (AN-Stabilität) und nichtlineare kontraktive Systeme
(B-Stabilität)zugrunde.Ordnungsreduktionsphänomene beisteifenProblemenwerdenmitHilfeder
B-Konvergenz-Ordnungerfaßt.
Bei differential-algebraischen Systemen von höherem Index gestaltet sich die Übertragung der
Stabilitätskonzepte äußerstschwierig,manvergleiche[111],[114].Fürdievonunszubetrachtenden
DAEs vom Index 1 in Hessenbergform beschränken wir uns auf das Konzept der A-Stabilität, um
effizienteVerfahrenauszuwählen.
5.1.3 ImpliziteRunge-Kutta-Verfahren
BeiexplizitenEinschritt-Verfahren istdie StabilitätsfunktionR(z) einPolynom, daherkann einsol-
ches Verfahren nicht A-stabil sein. Implizite Verfahren besitzen rationale Stabilitätsfunktionen, die
Approximationen an die Exponentialfunktion (Lösung von (5.2)) darstellen. Für s-stufige Verfahren
erhältmanimZählerundNennerPolynomevomMaximalgrads.FürdieGauß-bzw.Radau-Verfah-
ren [46] ergeben sich die (s;s)- bzw. (s 1;s)-Padé-Approximationen an die Exponentialfunktion
(cid:0)
alsStabilitätsfunktionen. BeidesindbekanntermaßenA-stabil.Die(s 1;s)-Padé-Approximationist
(cid:0)
zudemauchL-stabil,daderNennergradgrößeralsderZählergradistundsomitR( )= 0gilt.
1
BeiimplizitenRunge-Kutta-Verfahren müssendieVerfahrensgleichungen iterativ gelöstwerden.
EineeinfacheFunktionaliterationverschenktdengewonnenenVorteil,daKontraktivitätgleichbedeu-
tend mit einer Einschränkung hL < C an die Schrittweite ist – aber L ist für steife Systeme eben
sehr groß. Typischerweise wird ein vereinfachtes Newton-Verfahren zur Lösung der Verfahrensglei-
chungenverwendet, wobei die Jacobi-Matrixf (y) zumeistam Startwerty = y approximiert wird
y n
104
5.2Linear-impliziteVerfahren
oder (insofern das Konvergenzverhalten der Iteration zufriedenstellend ist) eine Jacobi-Matrix aus
vorhergehendenSchrittenverwendetwird.
DiesführtinjedemNewton-SchrittaufeinlinearesGleichungssystemderDimensionn s,was
y
(cid:2)
fürgroße Systeme mit einem hohen Aufwand verbunden ist. Es liegt daher nahe, die Verfahrensvor-
schriftzumodifizieren,umdiegutenStabilitätseigenschaftenmiteffizientauflösbarenVerfahrensglei-
chungenzukombinieren.EinersterSchrittindieseRichtungistderÜbergangzudiagonal-impliziten
Verfahren,mansiehe[2],[42],[20],mitdenStufen
i
Y =y +h a f(Y ): (5.4)
mi m ij mj
j=1
X
Die einzelnen Stufengleichungen können nun sukzessive gelöst werden. Ein vereinfachter Newton-
SchrittfürdasInkrementg :=Y y inderi-tenStufeergibtsichmitJ f (y )zu
i mi m y m
(cid:0) (cid:25)
i(cid:0)1
(l+1) (l) (l) (l)
(I ha J)(g g )=g h a f(Y ) ha f(y +g ) (5.5)
(cid:0) (cid:0) ii i (cid:0) i i (cid:0) ij mj (cid:0) ii m i
j=1
X
EssindsomitsukzessivesnichtlineareSystemederDimensionn zulösen.WähltmandieDiagonal-
y
elemente a der Verfahrensmatrix alle gleich, so ist nur eine LU-Zerlegung pro Schritt erforderlich.
ii
DieseVerfahrenheißenSDIRK-Verfahren(singlydiagonallyimplicitRunge-Kutta).
5.2 Linear-implizite Verfahren
BeiderImplementierungimpliziterVerfahrenwerdendieVerfahrensgleichungeniterativmitdemver-
einfachtenNewtonverfahrengelöst,wobeidieIterationzweckmäßigerweiseabgebrochenwird,sobald
die Genauigkeit der Lösung imBereich der gefordertenToleranz(oderum einenbestimmten Faktor
darunter) liegt. Dabei wird das Konvergenzverhalten von der Schrittweite abhängen, mit kleinerer
SchrittweiteverfügtmanaußerdemübereinengenauerenStartwert.
Die Konvergenz der vereinfachten Newton-Iterationkann somit durch die Schrittweiteh gesteu-
ert werden. Man fixiert einfach eine Newton-Iteration pro Stufe — und gelangt zu linear-implizi-
ten Verfahren. Wählt man J = f (y ), so bezeichnet man die Verfahrensklasse als Rosenbrock-
y m
Methoden [85, 107, 46, 99]. Einen allgemeineren Ansatz bieten die von Strehmel und Weiner ent-
wickeltenadaptivenRunge-Kutta-Methoden,diesichdurchguteStabilitätseigenschaftenauszeichnen
[95,97,8,98,99].
5.2.1 Rosenbrock-Methoden
Mit der Jacobi-Matrix J = f (y ) ergibt sich ein vereinfachter Newton-Schritt für die internen
y m
Ableitungenk =f(y +h a k )zu
i m j ij j
P
i(cid:0)1
(1) (0) (0) (0)
(I h(cid:13)J)(k k )= k +f(y +h a k +h(cid:13)k ): (5.6)
(cid:0) i (cid:0) i (cid:0) i m ij j i
j=1
X
105
5.KRYLOV-ROW-METHODEN
(0)
Bei der Wahl des Startwertes k bietet sich eine Kombination der vorhergehenden Stufen an. Wir
i
(1)
erhaltensomitalsStufenlösungk = k nacheinemNewtonschritt
i i
i(cid:0)1
(cid:13)
(0) ij
k = k
i (cid:0) (cid:13) j )
j=1
X
i(cid:0)1 i(cid:0)1 i(cid:0)1
(cid:13) (cid:13)
ij ij
(I h(cid:13)J)(k + k )=f(y +h (a (cid:13) )k )+ k : (5.7)
i j m ij ij j j
(cid:0) (cid:13) (cid:0) (cid:13)
j=1 j=1 j=1
X X X
(0)
DieVerfahrensvorschrift ergibtsichmit(cid:11) = a (cid:13) ,k := k undderüblichenAufdatierung
ij ij (cid:0) ij i (cid:0) i
mitGewichtenb zu
i
i(cid:0)1
Y =y +h (cid:11) k (5.8a)
mi m ij j
j=1
X
(I h(cid:13)J)(k +k ) =f(Y )+k (5.8b)
i i mi i
(cid:0)
s
y =y +h b k : (5.8c)
m+1 m i i
i=1
X
Rosenbrock-Methoden erweisen sich als sehr effiziente Verfahren zur Lösung großer steifer Diffe-
rentialgleichungen. Sie reproduzieren näherungsweise die guten Stabilitätseigenschaften impliziter
Verfahrendurch einesehr einfache Verfahrensvorschrift. Außerdemhabensie denVorteil,leichtim-
plementierbarzu sein. Einige Abstriche müssenin der B-Konvergenzordnung [96] gemacht werden,
da sie Abkömmlinge einfach diagonal-impliziter Verfahren sind, die naturgemäß nicht die verein-
fachende Bedingung C(2) erfüllen können. Im Code ROS4 [46] ist mit der Methode GRK4T von
Kaps/Rentrop[59]einsehreffizientesVerfahrenvierterOrdnungimplementiert.
Rosenbrock-MethodenkönnenaußerdemeffizientzurLösungdifferential-algebraischer Systeme
eingesetztwerden,sofürProblemevomIndex 1in[83],[82],[84]und[3],unddirektfürdieIndex-
3-FormulierungbeimechanischenMehrkörpersystemenin[112],[113],[114],[115].
5.2.2 W-Methoden
In vielen Anwendungen kann man die steifen Komponenten der Lösung lokalisieren. Der Aufwand
zur Lösung der linearen Gleichungssysteme, wie sie bei Rosenbrock-Methoden auftreten, kann ent-
scheidendgesenktwerden, wenn man nichtdie komplette Jacobi-Matrix,sondernnur eine Approxi-
mationT vonniederemRangverwendet,sodaßlineareGleichungssystememitderMatrix(I h(cid:13)T)
(cid:0)
einfach zu lösen sind. Die Matrix J in (5.6) wird dann durch die Matrix T ersetzt. Ordnungsbedin-
gungenfürdiesoerhaltenenW-Methoden[93]werdenunabhängigvonderMatrixT formuliert,man
erhält zu den Ordnungsbedingungenfür Rosenbrock-Methoden zusätzliche Konsistenzbedingungen.
FürdenspeziellenFallT = f (y )+ (h)kannmiteinerkleinerenMengezusätzlicherKonsistenz-
y m
O
bedingungengearbeitetwerden.DieserFallbildetgleichzeitigdietheoretischeGrundlagefürdiebei
der Implementierung von Rosenbrock-Methoden zumeist verwendete Strategie, die gleiche Jacobi-
Matrixfürmehrereaufeinanderfolgende Schrittezuverwenden.
W-MethodenbietenaußerdemeineeinfacheMöglichkeitzurKonstruktionpartitionierterVerfah-
ren[109],mansetztinderJacobi-MatrixeinfachgewisseBlöckezuNull.
106
5.2Linear-impliziteVerfahren
5.2.3 Krylov-ROW-Methoden
Bei vielen Anwendungen ist nur eine gewisse Teilmenge der Lösungskomponenten steif. Um den
Aufwand zur Lösung der linearen Gleichungssysteme bei Rosenbrock-Methoden zu senken, kann
manzumeinenpartitionierteVerfahrenüberW-Methoden,wieobenbeschrieben,konstruieren.Dies
erfordertallerdingsdieexpliziteKenntnisdersteifenundnichtsteifenKomponenten.
EinenAnsatzzurautomatischenPartitionierungbietenKrylov-ROW-Methoden.WiebeidenW-
MethodenwirdJ in (5.8b)zunächstdurch eineApproximationT von niederemRangersetzt.Aller-
dingslassenwirjetztunterschiedlicheMatrizenT = T indenStufeni = 1;:::;szu.DieMatrizen
i
T sindKrylovraum-Approximationen(siehe[104],[86],[101])an J, undsomitgiltimallgemeinen
i
T = J+ (1).DurchgeeigneteApproximationkannjedocherreichtwerden,daßkeinezusätzlichen
i
O
Konsistenzbedingungen imVergleichzumFallT = J (Rosenbrock-Methode)auftreten.
i
AusgangspunktderKonstruktionisteineRosenbrock-Methode.DortistinjederStufeieinlinea-
res Gleichungssystem der Form (I h(cid:13)J)x = w zu lösen. Löst man dieses näherungsweise mit
i i
(cid:0)
FOM (siehe Abschnitt 5.3.2), so ergibt sich die Lösung x im (niedrig-dimensionalen) Krylovraum
i
(J;w ) w durch
K(cid:20)i i 3 i
(I hQ QT J)x =w ; (5.9)
(cid:0) (cid:20)i (cid:20)i i i
wobeidieMatrixQ durcheinenArnoldi-Prozeßgeneriertwird(sieheAbschnitt5.3.1).EineKrylov-
(cid:20)i
ROW-Methodeistsomitgegebendurch
Y =y +h (cid:11) k (5.10a)
mi m ij j
j
X
(I h(cid:13)Q QT J)(k +k )=f(Y )+k (5.10b)
(cid:0) (cid:20)i (cid:20)i i i mi i
y =y +h b k : (5.10c)
m+1 m i i
i
X
DieKoeffizientenergebensichdabeiausderzugrundeliegendenRosenbrock-Methode.Diereduzierte
Jacobi-MatrixbezeichnenwirmitT = Q QT J.
i (cid:20)i (cid:20)i
SchmittundWeiner[91]zeigen,daßuntergeeignetenVoraussetzungendieKrylov-ROW-Metho-
dedieselbeKonvergenzordnungwiediezugrundeliegende Rosenbrock-Methodehat.
Satz5.2.1. Seieny ;k ;Y mitderKrylov-ROW-Methode berechnet, undy~ ;k~ ;Y~ dieent-
m+1 i mi m+1 i mi
sprechendenWertederRosenbrock-Methode. Essei
(Tl Jl)w = hp(cid:0)l ; l = 1;:::;p 1; i = 1;:::;s: (5.11)
i i
(cid:0) O (cid:0)
(cid:16) (cid:17)
Danngilt
k k~ = (hp)
i i
(cid:0) O
y y~ = hp+1 :
m+1 m+1
(cid:0) O
FüreinendetailliertenBeweissowiemöglicheAbsc(cid:0)hwäch(cid:1)ungenderVoraussetzungenverweisen
wir auf [91]. Die Grundidee des Beweises besteht darin, zur Abschätzung der Differenz zwischen
Krylov-undRosenbrock-LösungdieInversenvon(I h(cid:13)J)und(I h(cid:13)T )inNeumann-Reihenzu
i
(cid:0) (cid:0)
entwickelnundunterAusnutzungvon(5.11)derenDifferenzabzuschätzen.
NachobigerDarstellungwerdendieStufengleichungenunabhängigvoneinanderdurchFOMge-
löst. Diese Vorgehensweise ist unbefriedigend, denn die Krylov-Räume der ersten Stufen enthalten
107
5.KRYLOV-ROW-METHODEN
Informationenüber die Eigenraumstrukturder Jacobi-Matrix,die in weiteren Stufen genutztwerden
kann.Für effizienteImplementierungenversuchtman daher, denim erstenSchrittgewonnenen Kry-
lovraumsukzessive mitweiterenStufenzuerweitern.Dabeifordertman:
In Stufe isuchenwir eine NäherungslösungderStufengleichung(I h(cid:13)J)x = w in einem
i i
(cid:15) (cid:0)
(cid:20) -dimensionalenRaum w (nichtnotwendigerweiseeinKrylovraum!).
i K(cid:20)i 3 i
DerRaum werdedurchdieSpaltenderorthogonalenMatrixQ aufgespannt.
(cid:15) K(cid:20)i (cid:20)i
ZurBestimmungderLösungnutzenwireineGalerkin-Bedingung(FOM),d.h.,dasResiduum
(cid:15)
seiorthogonalzu .
K(cid:20)i
Die aufspannenden Matrizen Q können mit einem mehrfachen Arnoldi-Prozeß berechnet werden,
i
mansieheAbschnitt5.3.1.
DieImplementierungdesmehrfachenArnoldi-ProzeßfürKrylov-ROW-MethodendurchWeiner,
Schmitt und Podhaisky im Code ROWMAP[110] resultiertin einem effizienten und robusten Löser
für große Systeme steifer Differentialgleichungen. Weitere Anwendungen von Krylov-Unterraum-
TechnikenfindensichindenArbeitenvonLubichundHochbruckzuexponentiellenIntegratoren,die
speziell für steife und stark oszillierende Probleme konstruiert sind, man siehe [52], [53], [54]. Die
Nutzung iterativer Lösungstechniken für BDF-Verfahren wurde von Brown, Hindmarsh und Petzold
in[7]beschrieben,eineImplementierungfindetsichimCodeVODPK.
5.3 Iterative Lösung linearer Gleichungssysteme
Wir wollen in diesem Abschnitt kurz die für die Krylov-ROW-Methoden benötigten Techniken zur
iterativen LösunglinearerGleichungssystemeumreißen.FürausführlichereDarstellungenverweisen
wiraufdieLiteratur[104],[86],[101].
5.3.1 Der Arnoldi-Prozeß
GesuchtseieineFaktorisierungA = QHQT miteinerHessenbergmatrix H undeinerorthogonalen
Matrix Q. Bei der orthogonalen Strukturierung werden orthogonale Transformationsmatrizenauf A
angewandt, und schrittweise die strukturierte Matrix H erzeugt. Umgekehrt führt eine strukturierte
OrthogonalisierungzurvollständigenArnoldi-Iteration[101].
MitAQ= QH ergibtsichfürdieSpaltenq vonQgerade
k
k+1
Aq = h q ; (5.12)
k jk j
j=1
X
wobeih dieEinträgederHessenbergmatrixH sind.BeimArnoldi-Prozeßnutztman(5.12)alsAn-
jk
satz für q und orthogonalisiertwie bei der modifizierten Gram-Schmidt-Orthogonalisierung [61]
k+1
q gegen q ;:::;q mit anschließender Normierung. Der Arnoldi-Prozeß berechnet auf diese Art
k+1 1 k
eine Gram-Schmidt-Orthogonalisierung der Krylov-Folge q ;Aq ;A2q ;:::. Das Verfahren bricht
1 1 1
ab,wennderKrylovraumdegeneriert,d.h.,wennAkq linearabhängigvonq ;:::;Ak(cid:0)1q ist.
1 1 1
NachnSchrittenergibtsichmit(q ;:::;q )eineorthogonaleBasisdesKrylovraums
1 n
= (A;q ) :=span q ;Aq ;:::An(cid:0)1q : (5.13)
n n 1 1 1 1
K K f g
108
5.3IterativeVerfahren
Bezeichnen wir mit Q die aus den ersten q orthogonalen Vektoren bestehende Matrix, mit H
n n n
2
Rn(cid:2)n bzw.H R(n+1)(cid:2)n linkeobereBlöckevonH,soergibtsich
n+1;n
2
AQ =Q H (5.14)
n n+1 n+1;n
QTAQ =H : (5.15)
n n n
DerOrthoprojektorindenKrylovraumistdurchQ QT gegeben.
n n
Wir erwähnen noch kurz den mehrfachen Arnoldi-Prozeß wie er im Code ROWMAP [110] im-
plementiert ist, für eine ausführliche Darstellung verweisen wir auf die Arbeiten [91], [110]. Beim
Aufbau des Krylovraumes wird dieser nach (cid:20) Schritten nicht durch Aq erweitert, sondern durch
1 (cid:20)1
einen beliebigen Vektor w . Wir orthogonalisieren dann w gegen und erhalten q . Die fol-
2 2 K(cid:20)1 (cid:20)1+1
gendenErweiterungendesRaumesergebensichdanndurchAq ,woraussichnachOrthogonalisie-
(cid:20)1
rungq ergibt, usw., d.h.,Aq wirdgegen q ;:::;q orthogonalisiert,woraussich dann
(cid:20)1+2 (cid:20)1+j 1 (cid:20)1+j+1
q ergibt. Im Prinzip ergibt sich eine Art Shift, Aq wird (nach Orthogonalisierung) nicht di-
(cid:20)1+j+2 j
rektnachq eingefügt,sonderninderi-tenStufedirektnachq (wobeiwirdenallgemeinenFall
j j+i(cid:0)1
voraussetzen,daßdierechtenSeitenw derStufennichtimbereitsgeneriertenAnsatzraumenthalten
i
sind).
In der Matrix QT AQ füllen sich durch den Shift weitere Nebendiagonalen auf, für jede Stufe
(cid:20)i (cid:20)i
i> 1kommtabZeile(cid:20) eineweitereNebendiagonalehinzu.
i
5.3.2 FOM
EinlinearesGleichungssystemAx = b,A Rm(cid:2)m,kannnäherungsweisegelöstwerden,indemesin
2
denKrylovraum (A;b)projiziertwird.Diesbedeutet,daßeineNäherungslösungximKrylovraum
n
K
gesuchtist, so daß die orthogonaleProjektiondes Residuums in den Krylovraum verschwindet.Das
sodefinierteVerfahrenheißtFullyOrthogonalMethod.BezeichnenwirdiemitdemArnoldi-Prozeß
nach n Schritten berechnete Orthonormalbasis des Krylovraumes mit Q , so bestimmt sich die so
n
definierteNäherungslösungx durch
n
QT(Ax b) =0mitx = Q w : (5.16)
n n n n n
(cid:0)
Darausergebensichw undx zu
n n
H w =QTb (5.17)
n n n
x =Q H(cid:0)1QTb: (5.18)
n n n
Dieskannauchsoformuliertwerden:AwirddurchdieniedrigdimensionaleApproximationQ H QT
n n n
ersetzt,undfürdasresultierendesinguläreGleichungssystemwirddieMinimum-Norm-kleinste-Qua-
drate-Lösungverwendet.DieselösthierdasGleichungssystemexakt,dabzumSpaltenraumvonQ
n
gehört.
5.3.3 GMRES
Ein weitereMöglichkeit zur iterativen Lösung eineslinearen GleichungssystemsAx = bstellt GM-
RES (Generalized Minimum Residual) dar. Zum Krylovraum (A;b) bestimmt man die Lösung
n
K
x (A;b)durchdieForderung
n n
2 K
Ax b min: (5.19)
n 2
k (cid:0) k !
109
5.KRYLOV-ROW-METHODEN
MitdemAnsatzx = Q w ergibtsich
n n n
AQ w b = Q H w b
n n 2 n+1 n+1;n n 2
k (cid:0) k k (cid:0) k
= H w QT b ;
k n+1;n n (cid:0) n+1 k2
manerhältw (undsomitdieNäherungx = Q w )alsLösungdesQuadratmittelproblems
n n n n
H w QT b min: (5.20)
k n+1;n n(cid:0) n+1 k2 !
LäßtmandieletzteKomponenteweg,soerhältmandaslösbareGleichungssystem(5.17).
IsteinegrößereAnzahlvonIterationendurchzuführen,soerweistsicheinNeustartdesVerfahrens
nacheinergewissenZahlvonIterationenalsvorteilhaft,dennderAufwandfürdieOrthogonalisierung
wächst von Schritt zu Schritt linear an. Im Falle eines Neustarts wird dann das zuletzt berechnete
ResiduumalsStartwertverwendet.
5.4 Krylov-Typ Rosenbrock-Methoden für DAE’s vom Index 1
Wirbetrachteneindifferential-algebraisches SystemvomIndex1inHessenbergform:
y0 =f(y;z); y Rny; z Rnz; (5.21a)
2 2
0=g(y;z): (5.21b)
IstdieJacobi-Matrixg regulär, so istobiges Systemvom Index 1.DieNebenbedingungkannin ei-
z
ner Umgebung der Lösung nach z = G(y) aufgelöstwerden. Einsetzenin die Differentialgleichung
führtzu einergewöhnlichen Differentialgleichung.Wir setzen voraus, daß die entstehendeDifferen-
tialgleichungy0 = f(y;G(y))steifist.SolcheSystemeentsteheninsbesonderebeiderSemidiskreti-
sierungparabolischerProbleme,und sie sind dannvon hoher Dimension. Oftist ein direkterEinbau
der Randbedingungen nicht zweckmäßig, so daß diese als algebraische Nebenbedingungen an das
System gekoppelt werden. Im letzten Abschnitt dieses Kapitels werden wir mit den Problemen der
ViskoelastizitäteineAnwendungsklassekennenlernen,dieaufgroßeSystemevom Index 1miteiner
großenZahlvonNebenbedingungenführt.
5.4.1 Adaption aufDAEs
Bei der Herleitung von Verfahren für differential-algebraische Systeme bieten sich prinzipiell zwei
Vorgehensweisenan,diealsdirekterundindirekterZugangbezeichnetwerden.
DerindirekteZugang
DerindirekteZugangistaufsemi-expliziteIndex-1-Systemebeschränkt.DieNebenbedingung(5.21b)
wirdzuz = G(y)nachzaufgelöstundindieDifferentialgleichung(5.21a)eingesetzt.Diesegewöhn-
licheDifferentialgleichungwirddannmitdemGrundverfahrengelöst.Füreines-stufigeRosenbrock-
Methode (5.8) ist f(Y ) durch f(Y ;G(Y )) zu ersetzen. Die Jacobi-Matrix der gewöhnlichen
mi mi mi
Differentialgleichungergibtsichmitg(y;G(y)) = 0zuJ = f +f G = f f g(cid:0)1g .Damitwird
y z y y z z y
(cid:0)
110
5.4Krylov-MethodenfürDAEs
(5.8)zu
y = y +h b k ; 0 =g(y ;z ) mit (5.22a)
m+1 m i i m+1 m+1
i
X
i(cid:0)1
Y = y +h (cid:11) k ; 0 =g(Y ;Z ); i = 1;:::;s; (5.22b)
mi m ij j mi mi
j=1
X
I h(cid:13)(f f g(cid:0)1g ) (k +k ) =f(Y ;Z )+k ; i = 1;:::;s: (5.22c)
y z z y i i mi mi i
(cid:0) (cid:0)
(cid:0) (cid:1)
DieBestimmungderUmkehrfunktionz = G(y)wirddurchdieimpliziteRelationindenGleichungen
(5.22b),(5.22a)realisiert.MankanndassokonstruierteVerfahrenaberauchmitderOriginal-Jacobi-
Matrixschreiben.Gleichung(5.22c)mußdanndurchdieäquivalenteGleichung
I 0 f f k +k f(Y ;Z )+k
h(cid:13) y z i i = mi mi i (5.23)
0 0 (cid:0) gy gz li 0
(cid:20)(cid:18) (cid:19) (cid:18) (cid:19)(cid:21)(cid:18) (cid:19) (cid:18) (cid:19)
ersetztwerden.DieVariablenl sindlediglichDummy-Variablen,siewerdennichtweiterverwendet.
i
DerdirekteZugang
BeimdirektenZugangwirddasSystemformalregularisiertdurch
y0 =f(y;z)
(5.24)
"z0 =g(y;z):
Auf das ODE-System (5.24)wird die Verfahrensvorschrift angewandt, und dannder Grenzübergang
für " 0 durchgeführt. Untersuchungen zur Konvergenz so konstruierter Methoden findet man in
!
einerganzenReihevonArbeiten,mansiehe[83],[82],[84]und[3].
Bemerkung 5.4.1. Für implizite Runge-Kutta-Verfahren sind direkter und indirekter Zugang nahe-
zu äquivalent, man vergleiche [43, 46]. Der direkte Zugang führt zu Gleichungen 0 = g(Y ;Z )
mi mi
für die internen Werte Y ;Z , und dies ist de facto eine Auflösung der Nebenbedingungen nach
mi mi
derVariablenz.DieLösungenstimmeniny ;Y ;Z überein,lediglichdiez könnenunter-
m+1 mi mi m+1
schiedlichsein,dasichbeimdirektenZugangz = R( )z + b ! Z ergibt,währendsich
m+1 1 m i i ij mj
beim indirekten Zugang z aus 0 = g(y ;z ) bestimmt. Für steif-genaue Verfahren stimmt
m+1 m+1 m+1
P
aberauchdiesüberein.
DasResultatbeimdirektenZugangergibtsichwiefolgt:
s s
y = y +h b k z =z +h b l : (5.25a)
m+1 m i i m+1 m i i
i=1 i=1
X X
i(cid:0)1 i(cid:0)1
Y = y +h (cid:11) k Z =z +h (cid:11) l (5.25b)
mi m ij j mi m ij j
j=1 j=1
X X
I 0 f f k +k f(Y ;Z )+k
h(cid:13) y z i i = mi mi i (5.25c)
(cid:20)(cid:18)0 0(cid:19)(cid:0) (cid:18)gy gz(cid:19)(cid:21)(cid:18)li+li (cid:19) (cid:18) g(Ymi;Zmi) (cid:19)
Die Größenl sind analog zu den k durch l = i(cid:0)1 (cid:13) =(cid:13)l gegeben. Ein Vorteil des direktenZu-
i i i j=1 ij j
gangsist,daßdieKomponentez nichtdurchAuflösungderNebenbedingungbestimmtwerdenmuß.
P
111
5.KRYLOV-ROW-METHODEN
ImGegensatzzum indirektenZugangerfolgthier eineAufdatierungderKomponentez mit internen
Stufen,dieNebenbedingungg(y;z) = 0istnurinderGrößenordnungderDiskretisierungsfehlerdes
Verfahrenserfüllt.UmdieStabilitätderAufdatierungfürz zugarantieren,wirdintypischenKonver-
genzaussagen R( ) < 1 gefordert [46]. Die Koeffizientenmatrix des linearen Gleichungssystems
j 1 j
(5.25c) für die internen Steigungenk ;l stimmt mit der des Gleichungssystems(5.23) für die inter-
i i
nen Steigungen beim indirekten Zugang überein, lediglich die rechten Seiten unterscheiden sich, da
in(5.25)nochinterneSteigungenfürdiealgebraischenVariablenberechnetwerdenmüssen.
5.4.2 Krylov-ROW-Methodenfür DAEs
WirkönnenzumeineneinfachdieKrylov-ROW-MethodenmitdenbereitserwähntenTechnikenge-
eignet verallgemeinern,um sie auf DAEs anwenden zu können. Zum anderen können wir aber auch
wiedereinenSchrittzurückzudenRosenbrock-Methodengehen,dieseaufDAEserweitern,unddie
linearenGleichungssystemeiterativ lösen.FüreineÜberblicksdarstellung derbeiden Vorgehenswei-
sensieheman[121].
DerindirekteZugang
Unsere ODE hat jetzt die Gestalt y0 = f(y;G(y)), wobei G(y) implizit durch g(y;G(y)) = 0 defi-
niertist. Wir gehen hier davon aus, daß im Verfahren diese Nebenbedingung stetsmit hinreichender
Genauigkeitaufgelöstwird.
MitderJacobi-MatrixJ = f +f G = f f g(cid:0)1g könnenwirdanndasRosenbrock-Verfahren
y z y y z z y
(cid:0)
formulierenundzurKrylov-ROW-Methodeübergehen.EbensokannmandirektindieKrylov-ROW-
MethodefürODEsdieGleichungy0 = f(y;G(y))einsetzen,dasDiagramm
indirekt
ROW(ODE) ROW(DAE)
(cid:0)(cid:0)(cid:0)(cid:0)!
Krylov Krylov (5.26)
? indirekt ?
Krylov?(ODE) Krylov?(DAE)
y (cid:0)(cid:0)(cid:0)(cid:0)! y
kommutiertimFalledesindirektenZugangs.
Die komplette Verfahrensvorschrift ergibt sich aus (5.22), indem in (5.22c) die Jacobi-Matrix
J = f f g(cid:0)1g durchQ QTJ ersetztwird.
y (cid:0) z z y i i
DerdirekteZugangfürdieRosenbrock-Methode
Beim direkten Zugang müssen wir die Reihenfolge der Operationen beachten. Das folgende Dia-
grammverdeutlichtnocheinmaldiezweimöglichenWege:
direkt
ROW("-ODE) ROW(DAE)
(cid:0)(cid:0)(cid:0)(cid:0)!
Krylov Krylov (5.27)
? direkt ?
Krylov(?"-ODE) Krylov?(DAE)
y (cid:0)(cid:0)(cid:0)(cid:0)! y
Zu Beginn wird stets die DAE über die "-Einbettung in eine ODE überführt und die Rosenbrock-
Methodeangewandt.WirdjetztzunächstderGrenzübergang " 0durchgeführt,sogelangtmanzu
!
einerRosenbrock-MethodefürDAEs.DieauftretendenlinearenGleichungssystemekönnendannmit
Krylov-Techniken gelöstwerden–diesistderrechteobereWegin(5.27).AufdemWeglinks-unten
112