Mathematische Methoden in der Dynamik

Modul 1337
Baudynamik und Erdbebeningenieurwesen
Vorlesung
Mathematische Methoden in der Dynamik
Vorlesungsskript: Prof. Thomas Apel
Institut für Mathematik und Bauinformatik
Fakultät für Bauingenieurwesen und Umweltwissenschaften
Februar 2016
Vorwort
Das Modul 1623, Baudynamik und Erdbebeningenieurwesen, wird für Studierende des Master-Studiengangs Bauingenieurwesen und Umweltwissenschaften
und des Master-Studiengangs Mathematical Engineering gelesen. Die Vorlesung
Mathematische Methoden in der Dynamik bildet einen Modulbestandteil im Umfang von zwei Trimesterwochenstunden. Das Skript wurde im Frühjahrstrimester
2010 neu erstellt und seither jedes Jahr überarbeitet. Für ihre Hilfe bei der
Texterfassung und beim Rechnen von Beispielen sei an dieser Stelle Leonid
Blank und Bernhard Gehrmann gedankt.
Das Skript ist kein Lehrbuch. Es ist an vielen Stellen knapper gehalten als es
Lehrbücher sind, deshalb wird die gleichzeitige Verwendung von Lehrbüchern
empfohlen; Tipps gibt die Literaturliste auf Seite 4.
Das Skript gibt die Stoffauswahl des Vorlesenden wieder, so dass klar wird,
welcher Stoff prüfungsrelevant ist. Außerdem erspart das Vorliegen des Skripts,
dass während der Vorlesungen alles angeschrieben werden muss. In den Vorlesungen werden Herleitungen und Beispiele meist ausführlicher besprochen;
Sie sollten das Skript selbständig um Bemerkungen ergänzen, die Ihnen zum
Verständnis des Stoffes notwendig erscheinen.
Wichtige Webseiten
• https://dokumente.unibw.de/pub/bscw.cgi/6070897
Dokumente zur Vorlesung
2
24. März 2016
Inhaltsverzeichnis
Literatur
4
1. Einführende Beispiele und Erklärungen
5
1.1. Der Feder-Masse-Schwinger . . . . . . . . . . . . . . . . . . . . . 5
1.2. Schwingungen mehrerer Massen . . . . . . . . . . . . . . . . . . . 9
1.3. Systeme erster Ordnung . . . . . . . . . . . . . . . . . . . . . . . 15
2. Numerische Verfahren für Differentialgleichungen erster
2.1. Allgemeine Herangehensweise . . . . . . . . . . . . .
2.2. Einschrittverfahren . . . . . . . . . . . . . . . . . . .
2.3. A-Stabilität, absolute Stabilität . . . . . . . . . . . .
2.4. A-Stabilität für die Schwingungsgleichung . . . . . .
2.5. Mehrschrittverfahren . . . . . . . . . . . . . . . . . .
Ordnung
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.
.
.
.
.
16
16
17
21
25
28
3. Numerische Verfahren für Schwingungsprobleme
3.1. Allgemeine Herangehensweise . . . . . . . . .
3.2. Newmark-Verfahren . . . . . . . . . . . . . .
3.3. Houbolt-Verfahren . . . . . . . . . . . . . . .
3.4. Kollokationsverfahren . . . . . . . . . . . . .
3.5. Hilber-Hughes-Taylor-Verfahren . . . . . . . .
3.6. Diskussion der Verfahren . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
31
31
31
37
38
40
42
.
.
.
.
.
.
.
50
50
53
54
60
64
66
68
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4. Numerische Lösung von Eigenwertproblemen
4.1. Matrizen mit spezieller Struktur . . . . . . . . . .
4.2. Eigenwertberechnung und Nullstellensuche . . . . .
4.3. Der QR-Algorithmus . . . . . . . . . . . . . . . . .
4.4. Die Potenzmethode und verwandte Verfahren . . .
4.5. Das Lanczos-Verfahren . . . . . . . . . . . . . . . .
4.6. Symmetrische verallgemeinerte Eigenwertprobleme
4.7. Zusammenfassung . . . . . . . . . . . . . . . . . .
A. Grundlagen Eigenwertprobleme
24. März 2016
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
69
3
Literatur
[Bat02] Bathe, K.-J.: Finite-Elemente-Methoden. Springer, Berlin, 1986, 1990,
1995, 2002.
[GvL07] Golub, G. H. und Ch. F. van Loan: Matrix Computations. The
Johns Hopkins University Press, Baltimore, Maryland, 1983, 1989,
1996, 2007.
[Hug00] Hughes, T. J. R.: The finite element method. Linear static and
dynamic finite element analysis. Dover Publications., Mineola, NY,
2000. Reprint of the 1987 original ed.
[HW96] Hairer, E. und G. Wanner: Solving ordinary differential equations
II. Stiff and differential-algebraic problems. Springer, Berlin, 1996.
[MV03] Meyberg, K. und P. Vachenauer: Höhere Mathematik, Band 1.
Springer, Berlin, . . . , 1999, 2001, 2003.
[QSS02] Quarteroni, A., R. Sacco und F. Saleri: Numerische Mathematik
2. Springer, Berlin, 2002.
[SK91]
Schwetlick, H. und H. Kretzschmar: Numerische Verfahren
für Naturwissenschaftler und Ingenieure – eine computerorientierte
Einführung. Fachbuchverlag, Leipzig, 1991.
[TB00]
Trefethen, L. N. und D. Bau: Numerical linear algebra. SIAM,
Philadelphia, 2000.
[Wat10] Watkins, D. S.: Fundamentals of matrix computations. Wiley, Hoboken, NJ, 1991, 2002, 2010.
[ZTZ05] Zienkiewicz, O. C., R. L. Taylor und J.Z. Zhu: The Finite Element
Method: Its Basis and Fundamentals. Elsevier, Oxford, 2005.
4
24. März 2016
1. Einführende Beispiele und Erklärungen
1.1. Der Feder-Masse-Schwinger
Bsp 1.1 (Feder-Masse-Schwinger)
X
XXX
XX
XX
XXX
X
X
XXX
elastische Feder
Gleichgewichtslage
u(t)
u̇(t)
ü(t)
...
...
...
Auslenkung
Geschwindigkeit
Beschleunigung
~ ?u(t)
HH
Pendelmasse
Dämpfungsvorrichtung
(zähe Flüssigkeit)
Zur mathematischen Modellierung dieses Systems betrachten wir die Kräfte, die
auf die Masse m wirken. Nach dem Hookeschen Gesetz gilt für die Rückstellkraft
der Feder F1 = −ku mit der Federkonstanten k. Die Reibung ist proportional
zur Geschwindigkeit, also wirkt eine Kraft F2 = −cu̇. Zusätzlich kann noch eine
äußere Kraft F3 = F (t) wirken. Nach dem Newtonschen Gesetz gilt folglich
mü = F1 + F2 + F3 = −ku − cu̇ + F (t),
also die gewöhnliche Differentialgleichung 2. Ordnung
mü + cu̇ + ku = F (t).
(1.1)
Dabei sind die Sonderfälle
• freie Schwingung:
F (t) = 0
• ungedämpfte Schwingung:
b=0
eingeschlossen. Zur vollständigen Modellierung gehören noch die Anfangsbedingungen
u(0) = u0 , u̇(0) = v0 ,
(1.2)
also die Vorgabe einer Anfangsauslenkung und einer Anfangsgeschwindigkeit.
E 1.2 In Beispiel 1.1 haben wir die Schwingungsgleichung (1.1) hergeleitet.
Manchmal ist es sinnvoll, die Gleichung so zu skalieren, dass der Koeffizient vor
der zweiten Ableitung Eins ist. Wir dividieren also durch m und erhalten mit
2δ = c/m ≥ 0, ω02 = k/m und f (t) = F (t)/m die Gleichung
ü + 2δ u̇ + ω02 u = f (t).
Dann
ist ω0 die Eigenfrequenz des ungedämpften Systems, siehe E 1.3, und
p
2
ω0 − δ 2 die Eigenfrequenz des schwach gedämpften Systems, siehe E 1.4.
24. März 2016
5
E 1.3 (Freie ungedämpfte Schwingung) Wir betrachten den Fall δ = 0, f (t) ≡
0, das heißt die Schwingungsgleichung
ü(t) + ω02 u(t) = 0
mit den Anfangsbedingungen (1.2). Durch Einsetzen der allgemeinen Lösung
u(t) = C1 cos(ω0 t) + C2 sin(ω0 t)
in die Anfangsbedingungen werden die freien Parameter C1 und C2 berechnet.
Unter Berücksichtigung von u̇(t) = ω0 (−C1 sin(ω0 t) + C2 cos(ω0 t)) erhält man
!
u(0) = C1 = u0 ,
!
u̇(0) = C2 ω0 = v0 ,
also C1 = u0 und C2 = v0 /ω0 , so dass die Lösung
u(t) = u0 cos(ω0 t) +
v0
sin(ω0 t)
ω0
(1.3)
lautet.
E 1.4 (Freie gedämpfte Schwingung: schwache Dämpfung) Wir wollen nun
δ > 0 voraussetzen. Die zur homogenen Gleichung gehörende charakteristische
Gleichung lautet
λ2 + 2δλ + ω02 = 0.
Ist δ < ω0 , dann erhalten wir mit Hilfe des Parameters ωd :=
p
ω02 − δ 2 > 0
λ1,2 = −δ ± iωd .
Die allgemeine Lösung der homogenen Gleichung (freie Schwingung) lautet also
u(t) = e−δt [C1 cos(ωd t) + C2 sin(ωd t)].
Gegenüber der ungedämpften Schwingung verkleinert sich die Kreisfrequenz,
denn es ist ωd < ω0 . Außerdem verkleinert sich die Amplitude der Schwingung
im Laufe der Zeit.
Die Konstanten C1 und C2 können aus den Anfangsbedingungen (1.2) bestimmt
werden. Man erhält C1 = u0 und C2 = (v0 + δu0 )/ωd , also
v0 + δu0
u(t) = e−δt u0 cos(ωd t) +
sin(ωd t) .
ωd
(1.4)
Die rote Kurve in Abbildung 1 stellt die Lösung im Fall ω0 = 1, δ = 0.2, u0 = 1,
v0 = 0 dar.
6
24. März 2016
E 1.5 (Freie gedämpfte Schwingung: Kriechfall) Bei starker Dämpfung, δ >
ω0 , erhalten wir zwei negative reelle Lösungen der charakteristischen Gleichung,
q
λ1,2 = −δ ± δ 2 − ω02 < 0.
Die allgemeine Lösung der homogenen Gleichung ist dann eine Linerkombination
zweier monoton fallender Exponentialfunktionen,
u0 (t) = C1 eλ1 t + C2 eλ2 t .
Das mechanische System ist infolge zu starker Dämpfung (Reibung) zu keiner
echten Schwingung mehr fähig.
Die Konstanten C1 und C2 können wieder aus den Anfangsbedingungen (1.2)
bestimmt werden. Man erhält als Lösung
u(t) =
(v0 − u0 λ2 )eλ1 t − (v0 − u0 λ1 )eλ2 t
.
λ1 − λ2
(1.5)
Die grüne Kurve in Abbildung 1 stellt die Lösung im Fall ω0 = 1, δ = 1.8,
u0 = 1, v0 = 0 dar.
1,0
w
0,8
= 1,
0
x
0
= 1,
v
0
= 0
0,6
d = 1.8
0,4
0,2
d=1
0
2
4
6
8
10
12
t
K
0,2
d = 0.2
Abbildung 1: Freie gedämpfte Schwingungen
E 1.6 (Freie gedämpfte Schwingung: aperiodischer Grenzfall) Im Grenzfall
δ = ω0 erhalten wir zwei gleiche negative reelle Lösungen der charakteristischen
Gleichung,
λ1,2 = −δ < 0.
Die allgemeine Lösung der homogenen Gleichung besitzt somit die Form
u0 (t) = (C1 t + C2 )e−δt .
Das mechanische System verhält sich ähnlich zum Kriechfall, wie die blaue Kurve
in Abbildung 1 illustriert. Bei stärkerer Dämpfung, δ > ω0 , klingt die Lösung
langsamer ab, bei geringerer Dämpfung, δ < ω0 , kommt es zum Durchschwingen.
Bestimmt man die Konstanten C1 und C2 aus den Anfangsbedingungen (1.2),
erhält man die Lösung
u(t) = (u0 + v0 t + δu0 t)e−δt .
24. März 2016
(1.6)
7
E 1.7 (Erzwungene Schwingung) Erzwungene Schwingungen entstehen durch
eine mit einer Erregerfrequenz ω periodisch wirkende Kraft, die in der rechten
Seite der Differentialgleichung modelliert wird,
ü + 2δ u̇ + ω02 u = f (t) := K0 sin(ωt).
Die allgemeine Lösung
u(t) = uhom (t) + uinh (t)
der inhomogenen Differentialgleichung ergibt sich als Summe der allgemeinen
Lösung uhom (t) der homogenen Differentialgleichung (siehe (1.4), (1.5) oder
(1.6)) und einer speziellen Lösung uinh (t) der inhomogenen Gleichung. Letztere
hat die Form einer ungedämpften Schwingung mit der Erregerfrequenz ω,
uinh (t) = A sin(ωt − φ),
wobei die Amplitude A zu
K0
A= p 2
(ω0 − ω 2 )2 + 4δ 2 ω 2
berechnet werden kann (Übung).
Im Fall δ > 0 klingt der Lösungsanteil uhom (t) mit der Zeit ab, so dass das
mechanische System nach einer Einschwingphase ungedämpft schwingt. Dabei
wird die von der Erregerfrequenz ω abhängige Amplitude A besonders groß (d. h.
maximal), wenn ω gleich der Resonanzkreisfrequenz
q
(1.7)
ωr = ω02 − 2δ 2
ist (Übung). Diesen Sonderfall bezeichnet man als Resonanzfall, die Schwingungsamplitude A beträgt dann
A=
K
p 0
.
2δ ω02 − δ 2
Bei kleiner Dämpfung, δ 1, kommt es bei einer Erregung mit der Resonanzfrequenz des Systems, ω = ωr ≈ ω0 , zu riesigen Amplituden, der
Resonanzkatastrophe. Diese kann zur Zerstörung des Systems führen, siehe zum Beispiel die Zerstörung verschiedener Brücken (Tacoma: http://
www.youtube.com/watch?v=3mclp9QmCGs, Ibáñez: http://www.youtube.com/
watch?v=t9kR9T9wZsM oder Wolgograd: http://www.youtube.com/watch?v=
WEQrt_w7gN4).
8
24. März 2016
ungedaempft
gedaempft
Kriechfall
1
0.5
0
−0.5
−1
−1
−0.5
0
0.5
1
Abbildung 2: Phasendiagramm für verschiedene Schwingungen
E 1.8 (Phasendiagramm) Man kann die Funktionen u(t) und u̇(t)/ω als Kurven mit dem Kurvenparameter t in einem Phasendiagramm darstellen. Bei
einer freien, ungedämpften Schwingung erhält man einen Kreis. Bei freien,
(schwach) gedämpften Schwingungen ergibt sich eine Spirale. Abbildung 2 zeigt
das Phasendiagramm für u0 = 1, v0 = 0, ω0 = 1, t ∈ [0, 5π] in den Fällen
δ = 0 (ungedämpfte Schwingung), δ = 0,1 (schwache Dämpfung) und δ = 1,1
(Kriechfall).
1.2. Schwingungen mehrerer Massen
Bsp 1.9 Wir betrachten ein System aus N Massepunkten, die mit Federn/Dämpfern
verbunden sind. Gesucht sind die Auslenkungen


xi,1 (t)
xi (t) =  xi,2 (t)  ,
i = 1, . . . , N,
xi,3 (t)
dieser Massepunkte aus der Ruhelage, die man zu einem Vektor


x1


x = x(t) =  ...  ∈ R3N
xN
zusammenfassen kann. Das Kräftegleichgewicht lautet
M ẍ + C ẋ + Kx = f (t).
(1.8)
Dabei bedeuten M eine Massenmatrix (Diagonalmatrix), C eine Dämpfungsmatrix, K eine Steifigkeitsmatrix und f (t) eine Erregerfunktion.
24. März 2016
9
E 1.10 (Wellengleichung) Wir betrachten die longitudinale Bewegung eines
elastischen Stabes.
0
-
x
L
Bezeichne
L
[m]
die Länge des Stabes,
A
[m2 ]
die konstante Querschnittsfläche des Stabes,
ρ
[kg/m3 ]
die Dichte,
[m]
die Verschiebung des Querschnittselements, das
sich ursprünglich an der Stelle x ∈ (0, L) befand
= Auslenkung aus dem Ruhezustand,
1
die Dehnung,
[N/m2 ]
den Elastizitätsmodul,
[N/m2 ]
die Spannung nach dem Hookeschen Gesetz,
[ mkg2 s ]
die Impulsdichte,
u(x, t)
=
∂u
∂x
E
σ = Eε
∂u
ρ
∂t
dann lautet die Newtonsche Bewegungsgleichung (das Kräftegleichgewicht) in
einem beliebigen Intervall (a, b) ⊂ (0, L)
Z
Z b
∂u
∂σ
d b
ρA
dx = [σ(b) − σ(a)] A =
A dx.
|
{z
}
dt a
∂t
a ∂x
|
{z
}
Spannungen
|
{z
}
Gesamtimpuls
|
{z
}
Kraft
Kraft
Da das Integrationsgebiet (a, b) beliebig ist, gilt folglich
∂
∂u
∂σ
ρA
=
A
∂t
∂t
∂x
∂2u
∂σ
∂
∂u
ρ 2 =
=
E
∂t
∂x
∂x
∂x
Ist E 6= E(x), dann gilt die Wellengleichung
2
∂2u
2∂ u
−
a
= 0,
∂t2
∂x2
a2 =
E
.
ρ
Zur vollständigen Beschreibung benötigt man zwei Anfangsbedingungen
u(x, 0) = u0 (x)
ut (x, 0) = v0 (x)
(Anfangsauslenkung)
(Anfangsgeschwindigkeit)
und an jedem Rand eine Randbedingung. Möglichkeiten sind:
oder
oder E
10
u(0, t) = 0
∂u
(0, t) = 0
∂x
∂u
(0, t) + ku(0, t) = 0
∂x
(feste Einspannung)
(freies Ende, entspricht σ(0, t) = 0)
(elast. Aufhängung, σ(0, t) = −ku(0, t))
24. März 2016
Bsp 1.11 (Schwingende Saite) Wir suchen die Auslenkung u(x, t) bei einer
schwingenden Saite.
Saite
u
x
0
1
Sie wird durch die dieselbe Anfangs-Randwertaufgabe
2
∂2u
2∂ u
−
a
= 0,
∂t2
∂x2
u(x, 0) = u0 (x)
ut (x, 0) = v0 (x)
u(0, t) = u(1, t) = 0
(1.9)
Anfangsauslenkung,
Anfangsgeschwindigkeit,
Randbedingung,
beschrieben, wobei der Parameter a2 das Verhältnis aus Steifigkeit (E-Modul,
Vorspannung) und Masse (Dichte) beschreibt.
E 1.12 (Diskretisierung der Wellengleichung) Wir diskretisieren diese partielle Differentialgleichung in der Ortsvariablen mit Hilfe der Methode der finiten
Differenzen (Finite Elemente können analog verwendet werden). Dazu definieren
wir ein Ortsgitter mit der Schrittweite h = N1 ,
xj = jh,
j = 0, . . . , N,
und approximieren die Lösung u durch die Vektorfunktion u(t) mit den Komponenten uj (t) ≈ u(xj , t). Die zweite partielle Ableitung im Ort wird durch eine
zweite Differenz ersetzt,
∂2u
1
≈ 2 (uj+1 − 2uj + uj−1 ),
2
∂x
h
so dass die Differentialgleichung (1.9) im Punkt (xj , t) durch die Differenzengleichung
a2
üj − 2 (uj+1 − 2uj + uj−1 ) = 0, j = 1, . . . , N − 1,
h
approximiert wird. Dieses System gewöhnlicher Differentialgleichungen kann in
der Form
M ü + Ku = 0.
(1.10)
geschrieben werden, wobei die Massematrix M hier die Einheitsmatrix ist, aber
bei Verwendung der Finite-Elemente-Methode keine Diagonalmatrix ist. Man
beachte, dass das System dieselbe Form hat, wie das Schwingungsproblem bei
mehreren Punktmassen.
Ähnliche Modelle erhält man auch in der Akustik, bei Schwingungen elastischer Membranen und elastischer dreidimensionaler Körper, sowie bei elektromagnetischen Schwingungen.
24. März 2016
11
E 1.13 (Zusammenhang mit Eigenwertproblemen) Zur Lösung der Saitenschwingungsgleichung (1.9) wählen wir den Ansatz
u(x, t) = u(x)e iωt ,
wobei u(x) die Amplitude und ω = 2π
T die Kreisfrequenz zur Periodendauer T
1
beschreibt . Einsetzen in (1.9) liefert
u(x)( iω)2 e iωt − a2 u00 (x)e iωt = 0,
was mit ω 2 = λ zur Rand-Eigenwertaufgabe
−a2 u00 (x) = λu(x),
u(0) = u(1) = 0,
(1.11)
führt. Die Problemstellung lautet also: Finde Zahlen λ ∈ C und Funktionen
u(x) 6≡ 0, so dass (1.11) erfüllt ist.
In diesem einfachen Fall kann man die Lösungen exakt angeben,
λ k = a2 k 2 π 2 ,
uk (x) = sin(kπx),
k = 1, 2, 3, . . .
Die Lösung der Ausgangsaufgabe (1.9) kann man durch die Reihe
u(x, t) =
∞
X
Ck sin(kπx)e iakπt
k=−∞
darstellen, wobei die komplexen Koeffizienten Ck so ermittelt werden, dass sie
die Anfangsbedingungen u(x, 0) = u0 (x), ut (x, 0) = v0 (x) erfüllen.
Aufg 1.14 Wiederholen Sie Ihre Kenntnisse zu Matrx-Eigenwertproblemen aus
dem 1. Trimester. Als Hilfe kann der Anhang dieses Skripts dienen, siehe Seiten 69 ff.
E 1.15 (Diskretisierung des Rand-Eigenwertproblems) Wenn wir das Differenzenverfahren zur näherungsweisen Lösung der Rand-Eigenwertaufgabe (1.11)
verwenden, erhalten wir das Matrix-Eigenwertproblem


 2


− h12
u1
u1
h2
2
 u2 
 − 12
  u2 
− h12
h2



 h

1
2
1

 u3 



− h2
− h2


  u3 
h2
2
a 
  ..  = λ  ..  ,
..
..
..






.
.
.

 . 
 . 
2
1
1

 uN −2 
− h2
− h2   uN −2 
h2
1
2
uN −1
uN −1
− h2
h2
|
{z
} | {z }
u
K
kurz
Ku = λM u.
1
Dieser Ansatz ist schön kompakt aufschreibbar. Man kann auch u(x, t) = u1 (x) cos(ωt) +
u2 (x) sin(ωt) ansetzen und erhält, dass u1 und u2 Lösungen der gleichen RandEigenwertaufgabe (1.11) sind.
12
24. März 2016
Eigenfunktionen, h=1/40
0.25
0.2
0.15
0.1
0.05
0
−0.05
−0.1
−0.15
−0.2
−0.25
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Abbildung 3: Näherungen der ersten drei Eigenfunktionen
Lösen wir dieses, erhalten wir N − 1 Näherungen für die unendlich vielen
Eigenwerte des kontinuierlichen Eigenwertproblems (1.11). Dabei werden die
kleinsten Eigenwerte recht gut angenähert, die großen nicht. Für den Fehler in
der Eigenwertnäherung kann man
|λ − λh | ≤ Ch2 λ2
zeigen.
Bsp 1.16 (Anwendung des Differenzenverfahrens) Sei a = 1. In folgender Tabelle sind die Fehler bei der Approximation der drei kleinsten Eigenwerte λ1 = π 2 ,
λ2 = 4π 2 und λ3 = 9π 2 für verschiedene Schrittweiten angegeben. Welche Konvergenzordnung lesen Sie ab? Wie beeinflusst die Größe des Eigenwerts den
Approximationsfehler?
N = h−1
|λ1 − λ1,h |
|λ2 − λ2,h |
|λ3 − λ3,h |
10
20
40
8,0908 e−2
2,0277 e−2
5,0723 e−3
1,2818 e+0
3,2363 e−1
8,1108 e−2
6,3835 e+0
1,6317 e+0
4,1018 e−1
Die Näherungen der ersten drei Eigenfunktionen sind in Abbildung 3 dargestellt.
E 1.17 (Entkopplung des Systems) Wir haben gesehen, dass sich die Bewegung von Mehrkörpersystemen durch das Differentialgleichungssystem
M ü + C u̇ + Ku = F (t).
darstellen lassen. Diskretisiert man die Wellengleichung nur im Ort ergibt sich
ein Differentialgleichungssystem mit derselben Struktur.
24. März 2016
13
Im Fall von Rayleigh-Dämpfung ,
C = aM + bK,
erhalten wir
M ü + (aM + bK)u̇ + Ku = F (t).
(1.12)
Seien nun λi = ωi2 > 0 die Eigenwerte und die Spalten der Matrix V die
zugehörigen Eigenvektoren des verallgemeinerten Eigenwertproblems
Kv = λM v,
dann erhält man mit der Diagonalmatrix Λ, die auf der Hauptdiagonale die
Eigenwerte enthält, die Beziehung
KV = M V Λ.
Außerdem benutzen wir die Matrix V für eine Transformation der gesuchten
Vektorfunktion u auf die neue Unbekannte y durch die Gleichung
u = V y.
Eingesetzt in das Differentialgleichungssystem (1.12), erhält man
M V ÿ + (aM + bK)V ẏ + KV y = F (t)
M V ÿ + (aM V + bM V Λ)ẏ + M V Λy = F (t)
ÿ + (aI + bΛ)ẏ + Λy = V −1 M −1 F (t),
welches ein System von nicht gekoppelten Differentialgleichungen darstellt. Seien
nun fi (t) die Komponenten der Vektorfunktion V −1 M −1 F (t), dann kann man
dieses System in der Ein-Freiheitsgrad-Form“
”
ÿi + (a + bλi )ẏi + λi yi = fi
schreiben.
Mit λi = ωi2 und δi = 12 (a + bλi ) ergibt sich die aus Erklärung 1.2 bekannte
Form
ÿi + 2δi ẏi + ωi2 yi = fi .
Es ist aber auch die Form
ÿi + 2ξi ωi ẏi + ωi2 yi = fi .
mit dem Dämpfungsverhältnis ξi = 12 (aωi−1 + bωi ) gebräuchlich. Später werden
wir auch wieder u statt y für die gesuchte Auslenkung schreiben.
Man kann also ein numerisches Verfahren zur Lösung des Differentialgleichungssystems (1.12) anhand der skalaren Schwingungsgleichung analysieren. Man
beachte dabei, dass für den Fall, dass das System (1.12) durch eine Diskretisierung der räumlichen Variablen entstanden ist, und dass die Frequenzen ωi die
Eigenfrequenzen des diskretisierten Systems sind. Da nur die niedrigen Frequenzen auch etwas mit dem kontinuierlichen System zu tun haben, siehe Beispiel
1.16, sollten die durch die hohen Frequenzen entstehenden Lösungsanteile vom
Verfahren gedämpft werden. Die zeitliche Veränderung der niederfrequenten
Anteile sollte aber gut approximiert werden.
Bem 1.18 Verfahren zur numerischen Berechnung von Eigenwerten von Matrizen lernen wir in Kapitel 4 kennen.
14
24. März 2016
1.3. Systeme erster Ordnung
E 1.19 (Überführung in ein System 1. Ordnung) Wir betrachten das Modell
des Einmassenschwingers beziehungsweise die Schwingungsgleichung für eine
einzelne Schwingform,
ü + 2δ u̇ + ω 2 u = f,
u(0) = u0 ,
u̇(0) = v0 .
Die Differentialgleichung 2. Ordnung kann mit Hilfe der neuen Variablen v = u̇
in ein System 1. Ordnung umgeformt werden,
u̇ = v,
u(0) = u0 ,
v̇ = −2δv − ω 2 u + f, v(0) = v0 ,
beziehungsweise in vektorieller Schreibweise
u̇
0
1
u
0
=
+
,
v̇
−ω 2 −2δ
v
f
u(0)
v(0)
(1.13)
=
u0
v0
.
E 1.20 Es gibt zwei Zugänge zur numerischen Lösung von Anfangswertaufgaben.
• Man kann jede Anfangswertaufgabe als Spezialfall des allgemeinen Anfangswertproblems
ẏ(t) = f (t, y(t)),
y(0) = y0 ,
auffassen und Allzweck-Methoden mit sorgfältiger Fehlerkontrolle verwenden. Es stehen sehr leistungsfähige Software-Pakete zur Verfügung. Nicht
beachtet wird dabei die Struktur der Aufgabe, z. B., dass ein Teil der Unbekannten Positionen/Auslenkungen und ein anderer Teil Geschwindigkeiten
darstellen.
• Man kann spezielle Eigenschaften des betrachteten Problems identifizieren
(z. B., dass hohe Frequenzen durch die Diskretisierung ohnehin nicht
vernünftig abgebildet werden und deshalb herausgedämpft werden sollten)
und spezielle numerische Methoden entwickeln, die diese Eigenschaften
ebenfalls besitzen bzw. im Diskreten nachbilden.
Wir werden beide Zugänge diskutieren.
24. März 2016
15
2. Numerische Verfahren für Differentialgleichungen
erster Ordnung
2.1. Allgemeine Herangehensweise
E 2.1 (Häufige Vorgehensweise bei der numerischen Lösung)
der Anfangswertaufgabe
ẏ = f (t, y),
Zur Lösung
y(0) = y0 ,
diskretisiert man die Zeit t ∈ R+ durch die Einführung eines Gitters
ωτ = {t0 , t1 , t2 , . . .}
(Menge von Stützstellen) mit der lokalen Schrittweite τk = tk+1 − tk . Die
numerische Lösung sucht man in Form einer Gitterfunktion uτ : ωτ → Rn ,
wobei
uk := uτ (tk ) ≈ y(tk ), k = 0, 1, 2, . . . ,
gelten soll. Wir bezeichnen mit
τ := max τk
k
die maximal verwendete Schrittweite.
E 2.2 (Diskretisierungsfehler) In der Regel werden in jedem Zeitschritt lokale
Fehler erzeugt, die sich in komplizierter Weise übertragen und den globalen Diskretisierungsfehler
eτ (tk ) = y(tk ) − uτ (tk )
bilden. Gilt für τ ≤ τ∗ und für alle tk ∈ ωτ
|eτ (tk )| ≤ Cτ p ,
dann heißt das entsprechende Diskretisierungsverfahren konvergent mit der
Ordnung p.
E 2.3 (Wahl der Schrittweite) Die Schrittweite muss so gewählt werden, dass
der Fehler innerhalb einer vorgegebenen Toleranz liegt. Dazu verwendet man in
jedem Zeitschritt einen Fehlerschätzer. Liegt der geschätzte Fehler innerhalb der
Toleranz, wird der Zeitschritt akzeptiert, andernfalls wird der Schritt mit einer
kleineren Schrittweite wiederholt. Liegt der geschätzte Fehler deutlich unterhalb
der Toleranz, kann im nächsten Zeitschritt eine größere Schrittweite gewählt
werden. Dieses Vorgehen wird als adaptive Schrittweitensteuerung bezeichnet.
Wir werden darauf jedoch im Rahmen dieser Vorlesung nicht weiter eingehen.
Der Einfachheit halber beschränken wir uns hier auf eine konstante Schrittweite τ ,
d. h., es gilt tk = kτ , k = 1, 2, . . .
16
24. März 2016
2.2. Einschrittverfahren
Def 2.4 Explizite Einschrittverfahren sind Verfahren der Form
u0 = y0
uk+1 = uk + τ Φ(tk , uk , τ, f ),
k = 0, 1, 2, . . .
zur Lösung von Anfangswertaufgaben. Die Funktion Φ heißt Inkrement- bzw.
Zuwachsfunktion.
Bsp 2.5 (Eulerschen Polygonzugverfahren) Beim Eulerschen Polygonzugverfahren (auch explizites Euler-Verfahren) approximiert man die Differentialgleichung
ẏ(t) = f (t, y(t))
im Punkt t = tk durch die Differenzengleichung
uk+1 − uk
= f (tk , uk ).
τ
Äquivalent umgeformt bedeutet das
u0 = y0
uk+1 = uk + τ f (tk , uk ),
k = 0, 1, 2, . . .
Das Eulersche Polygonzugverfahren ist das einfachste Beispiel eines expliziten
Einschrittverfahrens und konvergiert mit erster Ordnung.
Bsp 2.6 Wir lösen die Anfangswertaufgabe
ÿ(t) + ω 2 y(t) = 0,
y(0) = 1,
ẏ(0) = 0,
mit dem expliziten Euler-Verfahren. Die Abbildung 4 zeigt das Phasendiagramm
und das t, y-Diagramm. Man erkennt, dass nur mit sehr kleinen Schrittweiten
halbwegs brauchbare Ergebnisse erzielt werden.
Bsp 2.7 (Verfahren zweiter Ordnung) Beispiele für explizite Einschrittverfahren zweiter Ordnung sind das zweistufige Verfahren von Heun (1900),
uk+1 = uk + τ 21 f1,k + 12 f2,k ,
f1,k = f (tk , uk ),
f2,k = f (tk + τ, uk + τ f1,k ),
und das Halbschrittverfahren von Collatz (1969),
uk+1 = uk + τ f2,k ,
f1,k = f (tk , uk ),
f2,k = f tk + 12 τ, uk + 12 τ f1,k .
Bsp 2.8 Das klassische Runge-Kutta-Verfahren
uk+1 = uk + τ
1
6 f1,k
+ 31 f2,k + 13 f3,k + 16 f4,k ,
f1,k = f (tk , uk ),
f2,k = f tk + 12 τ, uk + 12 τ f1,k
f3,k = f tk + 21 τ, uk + 12 τ f2,k
f4,k = f (tk + τ, uk + τ f3,k )
hat die Ordnung 4.
24. März 2016
17
Phasendiagramm, t ∈ [0, 2π]
2
1.5
1
v/ω
0.5
exact, ω beliebig
Euler explizit, τ=2π/40
Euler explizit, τ=2π/20
0
−0.5
−1
−1.5
−2
−2
−1.5
−1
−0.5
0
u
0.5
1
1.5
2
Auslenkung, t ∈ [0, 10π]
2
1
u
0
exact, omega=1
Euler explizit, ω=1, τ=2π/40
Euler explizit, ω=1, τ=2π/20
−1
−2
0
5
10
15
20
25
30
t
Abbildung 4: Phasendiagramm und t, u-Diagramm für das explizite EulerVerfahren
Def 2.9 (Verfahren vom Runge-Kutta-Typ) Ein r-stufiges Einschrittverfahren
vom Runge-Kutta-Typ ist ein Verfahren der Form
uk+1 = uk + τ
r
X
βi fi,k
i=1
Man unterscheidet explizite, semi-implizite und implizite Verfahren. Für ein
explizites Verfahren gilt
fi,k = f tk + γi τ, uk + τ
i−1
X
βi,s fs,k
βi,s fs,k
βi,s fs,k
s=1
Für ein semi-implizites Verfahren gilt
fi,k = f tk + γi τ, uk + τ
i
X
s=1
Für ein implizites Verfahren gilt
fi,k = f tk + γi τ, uk + τ
r
X
s=1
Ü 2.10 Was bedeuten diese kleinen Unterschiede für den Rechenaufwand?
(Beachten Sie, dass die Funktion f vektorwertig sein kann.)
18
24. März 2016
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
exact, ω beliebig
Euler implizit, τ=2π/40
Euler implizit, τ=2π/20
Euler implizit, τ=2π/10
Euler implizit, τ=2π/3
v/ω
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
exact, ω=1
Euler implizit, ω=1, τ=2π/40
Euler implizit, ω=1, τ=2π/20
Euler implizit, ω=1, τ=2π/10
Euler implizit, ω=1, τ=2π/3
0
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 5: Phasendiagramm und t, u-Diagramm für das implizite EulerVerfahren
Bsp 2.11 (Implizites Eulerverfahren) Mit
uk+1 = uk + τ f1,k
f1,k = f (tk + τ, uk + τ f1,k ) = f (tk+1 , uk+1 )
erhält man das implizite Eulerverfahren (Euler rückwärts). Es kann auch durch
die Vorschrift
u0 = y0 ,
uk+1 = uk + τ f (tk+1 , uk+1 ),
k = 0, 1, 2, . . .
beschrieben werden.
Bsp 2.12 Wir betrachten wie in Beispiel 2.6 die Anfangswertaufgabe
ü(t) + ω 2 u(t) = 0,
u(0) = 1,
u̇(0) = 0,
lösen dieses Mal aber mit dem impliziten Euler-Verfahren. Die Abbildung 5 zeigt
das Phasendiagramm und das t, u-Diagramm. Man erkennt, dass die eigentlich
ungedämpfte Schwingung stark gedämpft wird.
Bsp 2.13 Durch die größere Zahl von Koeffizienten ist es möglich, ein einstufiges
Verfahren zweiter Ordnung zu konstruieren.
uk+1 = uk + τ f1,k
1
1
f1,k = f (tk + τ, uk + τ f1,k )
2
2
Zur Berechnung von f1,k ist die Lösung einer im Allgemeinen nichtlinearen
Gleichung nötig.
24. März 2016
19
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
v/ω
0.2
exact, ω beliebig
Trapezregel, ω=1, τ=2π/20
Trapezregel, ω=1, τ=2π/10
Trapezregel, ω=1, τ=2π/3
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
0.5
u
0
exact, ω=1
Trapezregel, ω=1, τ=2π/20
Trapezregel, ω=1, τ=2π/10
Trapezregel, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 6: Phasendiagramm und t, u-Diagramm für die Trapezregel
E 2.14 (Trapezregel) Das semi-implizite Einschrittverfahren mit γ1 = β1,1 =
β1,2 = 0, γ2 = 1, β2,1 = β2,2 = β1 = β2 = 12 ,
uk+1 = uk +
τ
[f (tk , uk ) + f (tk+1 , uk+1 )]
2
ist von 2. Ordnung. Es wird als Trapezregel bezeichnet. (Warum?)
Bsp 2.15 Wir betrachten wie in Beispiel 2.6 die Anfangswertaufgabe
ü(t) + ω 2 u(t) = 0,
u(0) = 1,
u̇(0) = 0,
also
u̇
v̇
=
0
1
−ω 2 0
u
v
,
u(0)
v(0)
=
u0
v0
,
und lösen dieses Mal aber mit der Trapezregel,
uk+1
vk+1
=
uk
vk
τ
+
2
0
1
−ω 2 0
uk
vk
+
uk+1
vk+1
Im Phasendiagramm in Abbildung 6 erkennt man auf den ersten Blick gar
keine Fehler, das liegt daran, dass das Verfahren unabhängig von ω frei von
numerischer Dämpfung ist. Fehler treten trotzdem auf, sie äußern sich in einer
geringfügigen Phasenverschiebung, wie im t, u-Diagramm gut zu sehen ist.
20
24. März 2016
2.3. A-Stabilität, absolute Stabilität
E 2.16 Wir betrachten die Dahlquistsche Testgleichung
ẏ = λy,
y(0) = 1,
(2.14)
mit einem Parameter λ ∈ C. Für die exakte Lösung y(t) = eλt gilt
0
für Re λ < 0,
lim |y(t)| =
t→∞
∞
für Re λ > 0.
Für Re λ = 0 gilt |y(t)| = 1.
Eine wünschenswerte Eigenschaft ist also, dass für Re λ ≤ 0
|uk | ≤ 1
∀k
(2.15)
gilt. Man beweist die Eigenschaft (2.15) im Allgemeinen über die Beziehung
|uk+1 | ≤ |uk |.
Def 2.17 Ein numerisches Verfahren zur Approximation von (2.14) heißt für
z := τ λ absolut stabil (A-stabil), wenn (2.15) gilt.
Der Bereich der absoluten Stabilität ist die Menge
A = {z = τ λ ∈ C : (2.15) ist erfüllt},
d. h., A ist die Menge der Werte des Produktes τ λ, für die die numerische
Methode Lösungen liefert, die für tn → ∞ beschränkt bleiben.
Ein numerisches Verfahren zur Approximation von (2.14) heißt unbedingt absolut
stabil, wenn (2.15) für beliebige Re λ ≤ 0 und beliebige Schrittweiten τ > 0 gilt.
Folgerung 2.18 Eine Methode ist unbedingt absolut stabil, wenn C− ⊂ A gilt,
wobei C− = {z ∈ C : Re z ≤ 0}.
Bsp 2.19 Wir untersuchen das explizite Eulerverfahren. Es gilt
uk+1 = uk + τ f (tk , uk )
= uk + τ λuk
= uk (1 + τ λ)
Das Verfahren ist nur dann absolut stabil, wenn |1 + τ λ| ≤ 1 bzw., nach kurzer
Rechnung, wenn die Stabilitätsbedingung für das explizite Eulerverfahren
0<τ ≤−
2 Re λ
|λ|2
erfüllt ist. In Abbildung 7 ist der Stabilitätsbereich illustriert. Man sieht, dass die
Stabilitätsbedingung besonders ungünstig ist, wenn Im λ Re λ oder |λ| 1.
24. März 2016
21
Im
A
−1
1
Re
Abbildung 7: Stabilitätsbereich des expliziten Euler-Verfahrens
Euler vorwaerts, t ∈ [0, 10]
3
τ=0.41
τ=0.40
τ=0.39
τ=0.15
exakt
2
1
u
0
−1
−2
−3
0
1
2
3
4
5
t
6
7
8
9
10
Abbildung 8: t, u-Diagramm für Beispiel 2.20
Bsp 2.20 (Numerischer Test) Sei λ = −5, T = 10. Dann lautet die Stabilitätsbedingung τ < 0,4. Im numerischen Test beobachtet man folgendes Verhalten:
τ
τ
τ
τ
= 0,41
= 0,40
= 0,39
= 0,15
⇒
⇒
⇒
⇒
Lösung schaukelt sich auf,
Lösung oszilliert unverändert,
Lösung fällt betragsmäßig,
monoton fallende Lösung,
siehe Abbildung 8.
Bsp 2.21 (Explizite Einschrittverfahren) Die Grafik in Abbildung 9 stellt die
Ränder der Stabilitätsbereiche für explizite Einschrittverfahren 1. (grün), 2. (türkis), 3. (blau) und 4. Ordnung (rot) dar. Von unbedingter Stabilität ist man
weit entfernt, hier helfen nur implizite Verfahren.
22
24. März 2016
3
2
1
0
−1
−2
−3
−3
−2.5
−2
−1.5
−1
−0.5
0
0.5
Abbildung 9: Stabilitätsbereiche für explizite Einschrittverfahren 1. (grün), 2.
(türkis), 3. (blau) und 4. Ordnung (rot)
Im
A
−1
1
Re
Abbildung 10: Stabilitätsbereich das implizite Euler-Verfahren
Bsp 2.22 (Implizites Euler-Verfahren) Zur Analyse der absoluten Stabilität
rechnen wir
uk+1 = uk + τ λuk+1
(1 − τ λ)uk+1 = uk
uk+1 =
⇒
uk =
uk
1 − τλ
1
u0
=
k
(1 − τ λ)
(1 − τ λ)k
für u0 = 1.
Damit ist (2.15) äquivalent zu |1 − τ λ| ≥ 1. Der Bereich ist in Abbildung 10
illustriert. Das Verfahren ist unbedingt absolut stabil.
24. März 2016
23
Euler implizit, t ∈ [0, 0.5]
15
τ=0.5
τ=0.25
τ=0.15
τ=0.05
exakt
10
u
5
0
−5
0
0.05
0.1
0.15
0.2
0.25
t
0.3
0.35
0.4
0.45
0.5
Abbildung 11: t, u-Diagramm für das Beispiel in Bemerkung 2.23
E 2.23 Für ẏ = 5y, y(0) = 1 wächst die Lösung y = e5t für t → ∞. Für τ > 0,4
berechnet das implizite Eulerverfahren eine betragsmäßig fallende Lösung. Man
sieht, dass A-Stabilität eine sinnvolle Eigenschaft für Aufgaben mit Re λ < 0
ist, aber nicht für Aufgaben mit Re λ > 0.
Das Diagramm in Abbildung 11 zeigt die numerische Lösung für die Schrittweiten
τ = 0,05 (sinnvolle Lösung), τ = 0,15 (monoton wachsende Lösung mit großem
Fehler), τ = 0,25 (betragsmäßig exponentiell wachsende, aber alternierende
Lösung) und τ = 0,5 (betragsmäßig monoton fallende Lösung).
E 2.24 Explizite Verfahren benötigen wegen des kleinen Bereichs absoluter
Stabilität typischerweise sehr kleine Schrittweiten. Sie sind im Allgemeinen nur
konkurrenzfähig, wenn das bei impliziten Verfahren zu lösende Gleichungssystem
großdimensioniert und nichtlinear ist. So etwas kommt zum Beispiel bei der
Crashsimulation vor.
Bsp 2.25 Für das einstufige Verfahren zweiter Ordnung aus Beispiel 2.13 gilt
R(z) = 1 +
z
1−
z
2
=
1 + z/2
.
1 − z/2
(2.16)
Damit gilt |R(z)| ≤ 1 genau dann wenn
z z 1 + ≤ 1 − .
2
2
Mit z = x + iy erhält man die äquivalenten Ungleichungen
1 + z 2 ≤ 1 − z 2
2
2
2
2
2
2
1 + x2 + y2
≤ 1 − x2 + y2
x ≤ −x,
d. h., |R(z)| ≤ 1 gilt für alle z mit Re z ≤ 0. Das Verfahren ist unbedingt absolut
stabil.
Ü 2.26 Analysieren Sie die absolute Stabilität der Trapezregel (Beispiel 2.14).
24
24. März 2016
2.4. A-Stabilität für die Schwingungsgleichung
E 2.27 Wir betrachten die Schwingungsgleichung
ü + 2δ u̇ + ω 2 u = f,
u(0) = u0 ,
u̇(0) = v0 .
In Erklärung 1.19 haben wir diese mit Hilfe der neuen Variablen v = u̇ in ein
System 1. Ordnung umgeformt,
u̇
0
1
u
0
u(0)
u0
=
+
,
=
.
v̇
−ω 2 −2δ
v
f
v(0)
v0
Dieses können wir nun mit einem Verfahren für Systeme erster Ordnung numerisch lösen. Was bedeutet A-Stabilität in diesem Zusammenhang?
• Wir können das System erster Ordnung entkoppeln und so die Werte des
Parameters λ in der Dahlquistschen Testgleichung herausfinden, die zur
angebenenen Gleichung gehören.
• Wir können die Analyse der A-Stabilität auf den Fall eines Systems aus
zwei Differentialgleichungen 1. Ordnung anpassen.
Wir werden nun beide Zugänge diskutieren.
E 2.28 (Diagonalisierung des Systems erster Ordnung) Wir diagonalisieren
die im
√ System auftretende Matrix. Die Eigenwerte λ1,2 = −δ ± iωd mit
ωd = ω 2 − δ 2 sind echt komplex und die zugehörigen Eigenvektoren
1
1
und
λ1
λ2
sind nicht orthogonal. Mit den Bezeichnungen
λ1
1 1
Λ=
, V =
λ1 λ2
λ2
z1
z2
=V
−1
u
v
erhält man
V
u̇
v̇
= V ΛV
−1
u
v
+
0
f
,
0
u
= ΛV
+ V −1
f
v
ż1
z1
0
=Λ
,
+ V −1
ż2
z2
f
−1
u̇
v̇
−1
,
also zwei entkoppelte Gleichungen. Im homogenen Fall, f = 0, sind das die
Dahlquistschen Gleichungen
żi = λi zi .
Wir merken uns, dass uns die Dahlquistsche Testgleichung für λ = −δ ± iωd
interessiert.
Der Vorteil der Diagonalisierung besteht darin, dass wir das Verhalten von
Näherungsverfahren zur Lösung der Anfangswertaufgabe nur noch für eine
skalare Differentialgleichung (allerdings komplexwertig) analysieren müssen. Der
Nachteil besteht drin, dass die unterschiedliche physikalische Bedeutung der
Variablen u und v nicht mehr gesehen wird, denn die zi sind Linearkombinationen
aus u und v.
24. März 2016
25
Satz 2.29 (A-Stabilität nach Hughes) Ein numerisches Verfahren der Form
uj+1
vj+1
=A
uj
vj
+ cj
ist stabil, wenn die Eigenwerte der Matrix A betragsmäßig kleiner oder gleich
Eins sind. Dabei reicht es aus, den skalaren Fall uj , vj ∈ R1 zu betrachten. Falls
A einen doppelten Eigenwert2 hat, so muss dieser echt kleiner als Eins sein.
[Hug00, S. 497]
E 2.30 Die Erklärung von Hughes ist etwa die Folgende. Wir bezeichnen
yj :=
uj
vj
,
so dass das Verfahren in der Form
yj+1 = Ayj + cj
geschrieben werden kann. Nun betrachten wir eine Störung in den Anfangsdaten,
d. h., statt mit dem exakten Anfangswert y0 wird mit dem gestörten Anfangswert
ỹ0 = y0 + ε gerechnet. Wir erhalten
ỹ1 = Aỹ0 + c0 = A(y0 + ε) + c0 = y1 + Aε,
ỹ2 = Aỹ1 + c1 = A(y1 + Aε) + c1 = y2 + A2 ε,
..
.
ỹn = yn + An ε.
Sei nun A diagonalisierbar, A = V ΛV −1 , dann ist An = V Λn V −1 , und es gilt
wegen kΛn k = max{|λ1 |n , |λ2 |n } die Abschätzung
kAn εk = kV Λn V −1 εk ≤ kV kkΛn kkV −1 εk = C(V, ε) max{|λ1 |n , |λ2 |n },
wobei
(
≤C
wenn max{|λ1 |, |λ2 |} ≤ 1,
C max{|λ1 | , |λ2 | }
→ ∞ sonst.
n
n
Damit kleine Störungen in den Anfangsdaten nur kleine Änderungen in der
Lösung verursachen, fordern wir also, dass die Eigenwerte der Matrix A betragsmäßig kleiner oder gleich Eins sind.
Das gleiche Argument kann verwendet werden, wenn man den lokalen Diskretisierungsfehler analysiert. Im Fall, dass die Matrix A nicht diagonalisierbar ist,
muss für das gleiche Argument vorausgesetzt werden, dass die Eigenwerte von
A echt kleiner als Eins sind.
2
Eigentlich muss nur ein doppelter defektiver Eigenwert echt kleiner als Eins sein, siehe
die Begründungen in [Hug00, S. 497]. Defektive Eigenwerte sind solche, bei denen die
geometrische Vielfachheit kleiner als die algebraische Vielfachheit ist.
26
24. März 2016
E 2.31 Die Erklärung kann aber auch anders interpretiert werden. Wir wollen,
dass (insbesondere hochfrequente) Lösungsanteile, die durch die allgemeine
Lösung der homogenen Differentialgleichung in der Lösung auftreten, nicht
durch das numerische Verfahren verstärkt werden, sondern eher gedämpft werden.
Solche Lösungsanteile erfüllen in der numerischen Lösung die Beziehung
yj+1 = Ayj .
Ist A diagonalisierbar, A = V ΛV −1 , dann gilt
V −1 yj+1 = ΛV −1 yj
und somit die Abschätzung
kV −1 yj+1 k = kΛV −1 yj k ≤ kΛkkV −1 yj k = max{|λ1 |, |λ2 |}kV −1 yj k.
Stabilität im Sinne von Hughes bedeutet also, dass die Norm von V −1 yj für alle
j beschränkt bleibt.
Bsp 2.32 Wir diskretisieren die Schwingungsgleichung
u̇
0
1
u
=
v̇
−ω 2 −2δ
v
mit dem impliziten Euler-Verfahren,
uk+1
=
vk+1
1
−τ
uk+1
=
ω 2 τ 1 + 2δτ
vk+1
uk
vk
uk
vk
+τ
0
1
−ω 2 −2δ
uk+1
vk+1
,
und untersuchen die Stabilität. Das Verfahren ist stabil, wenn die Eigenwerte
der Matrix auf der linken Seite betragsmäßig größer oder gleich Eins sind. Es
gilt für δ ∈ [0, ω]
1−λ
−τ
0 = 2
ω τ 1 − λ + 2δτ = (1 − λ)(1 − λ + 2δτ ) + ω 2 τ 2
= (1 − λ)2 + 2δτ (1 − λ) + ω 2 τ 2
p
1 − λ± = −δτ ± i ω 2 τ 2 − δ 2 τ 2
p
λ± = 1 + δτ ∓ i ω 2 τ 2 − δ 2 τ 2
|λ± |2 = (1 + δτ )2 + ω 2 τ 2 − δ 2 τ 2
= 1 + 2δτ + ω 2 τ 2 ≥ 1.
Für δ ≥ ω gilt
λ± = 1 + δτ ±
p
p
δ 2 τ 2 − ω 2 τ 2 = 1 + τ δ ± δ 2 − ω 2 ≥ 1.
Das Verfahren ist also für alle δ ≥ 0 unbedingt stabil.
24. März 2016
27
2.5. Mehrschrittverfahren
E 2.33 Einschrittverfahren haben gewisse Nachteile:
• hohe Anzahl von Funktionswertauswertungen, die man für eine hohe
Konsistenzordnung benötigt,
• Funktionswertauswertungen aus einem früheren Schritt werden nicht mehrmals benutzt.
Abhilfe bringen Mehrschrittverfahren. Wir besprechen hier einige Klassen solcher
Verfahren, wobei wir der Einfachheit halber wieder voraussetzen wollen, dass
die Schrittweite konstant gleich τ ist.
Def 2.34 Bei einem r-Schritt-Verfahren werden zur Berechnung von uk+1 die
bereits berechneten Größen uk , uk−1 , . . . , uk−r+1 verwendet. Sie haben die allgemeine Form:
uk+1 + αr−1 uk + . . . + α0 uk−r+1 =
= τ [βr f (tk+1 , uk+1 ) + βr−1 f (tk , uk ) + . . . + β0 f (tk−r+1 , uk−r+1 )]
bzw.
uk+1 +
r−1
X
i=0
αr−1−i uk−i = τ
r−1
X
βr−1−i f (tk−i , uk−i ).
i=−1
mit |α0 | + |β0 | > 0. Für βr = 0 ist das Verfahren explizit, sonst implizit.
Bem 2.35 (Anlaufrechnung) Für k = 0, . . . , r − 2 muss man sich zunächst
Startwerte verschaffen, z. B. mit einem Einschrittverfahren.
Satz 2.36 (2. Dahlquist-Barriere) Es gilt:
• Kein explizites lineares Mehrschrittverfahren ist unbedingt A-stabil.
• Kein implizites lineares Mehrschrittverfahren ist unbedingt A-stabil, wenn
seine Konsistenzordnung größer als 2 ist.
Bsp 2.37 Ein explizites Zweischrittverfahren zweiter Ordnung ist die Mittelpunktregel
uk+1 = uk−1 + 2τ f (tk , uk ).
Wir testen das Verfahren am Beispiel aus 2.15. Das Phasendiagramm in Abbildung 12 zeigt Instabilitäten für im Vergleich zur Periode großen Schrittweiten.
Außerdem ist die Periodenverkürzung gut zu erkennen.
Bsp 2.38 Ein implizites Zweischrittverfahren vierter Ordnung ist
uk+1 = uk−1 +
τ
[f (tk−1 , uk−1 ) + 4f (tk , uk ) + f (tk+1 , uk+1 )] .
3
Haben Sie einen Wiedererkennungserfolg? Wir testen das Verfahren ebenfalls am
Beispiel aus 2.15. Man sieht in Abbildung 13 die hohe Genauigkeit im Fall, dass
die Schrittweite die Stabilitätsbedingung einhält. Für zu große Schrittweiten
wird das Verfahren jedoch instabil.
28
24. März 2016
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
v/ω
0.2
exact, ω beliebig
Mittelpunktregel, ω=1, τ=2π/20
Mittelpunktregel, ω=1, τ=2π/10
Mittelpunktregel, ω=1, τ=2π/3
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 2π]
1
u
0.5
0
exact, ω=1
Mittelpunktregel, ω=1, τ=2π/20
Mittelpunktregel, ω=1, τ=2π/10
Mittelpunktregel, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 12: Phasendiagramm und t, u-Diagramm zur Mittelpunktregel
Phasendiagramm, t ∈ [0, 2π]
1
v/ω
0.5
exact, omega beliebig
Zweischrittverfahren 4. Ordnung, ω=1, τ=2π/20
Zweischrittverfahren 4. Ordnung, ω=1, τ=2π/10
Zweischrittverfahren 4. Ordnung, ω=1, τ=2π/3
0
−0.5
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
0
exact, omega beliebig
Zweischrittverfahren 4. Ordnung, ω=1, τ=2π/20
Zweischrittverfahren 4. Ordnung, ω=1, τ=2π/10
Zweischrittverfahren 4. Ordnung, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 13: Phasendiagramm und t, u-Diagramm zum Zweischrittverfahren
vierter Ordnung
24. März 2016
29
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
v/ω
0.2
exact, ω beliebig
BDF2, ω=1, τ=2π/20
BDF2, ω=1, τ=2π/10
BDF2, ω=1, τ=2π/3
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 2π]
1
0.5
u
0
exact, ω=1
BDF2, ω=1, τ=2π/20
BDF2, ω=1, τ=2π/10
BDF2, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 14: Phasendiagramm zum BDF2-Verfahren
Bsp 2.39 (BDF-Verfahren) Man approximiert die linke Seite der Gleichung
ẏ(tk+1 ) = f (tk+1 , y(tk+1 ))
durch einen Differenzenquotient möglichst hoher Konsistenzordnung auf der
Basis der Stützstellen tk+1 , tk , tk−1 , . . . , tk−r+1 . Es ergeben sich die impliziten
BDF (backward differentiation formulae), die besonders für steife Differentialgleichungen geeignet sind.
Für r = 1 ergibt sich das implizite Euler-Verfahren (erster Ordnung)
1
[uk+1 − uk ] = f (tk+1 , uk+1 ).
τ
Für r = 2 ergibt sich das BDF2-Verfahren (zweiter Ordnung)
3
1
uk+1 − 2uk + uk−1 = τ f (tk+1 , uk+1 ).
2
2
Wir testen das Verfahren am Beispiel aus 2.15. Das Phasendiagramm in Abbildung 14 zeigt deutlich die Dämpfung bei großer Schrittweite.
30
24. März 2016
3. Numerische Verfahren für Schwingungsprobleme
Das Kapitel basiert auf [Hug00, Kapitel 9].
3.1. Allgemeine Herangehensweise
E 3.1 (Aufgabenstellung) Zu lösen ist das System gewöhnlicher Differentialgleichungen 2. Ordnung
M ẍ + C ẋ + Kx = f
mit den Anfangsbedingungen
x(0) = x0 ,
ẋ(0) = v0 ,
wobei x(t) die Auslenkung des schwingenden Systems, M die Massenmatrix, C
die Dämpfungsmatrix, K die Steifigkeitsmatrix und f den Vektor der äußeren
Kräfte (Kraftvektor) bedeuten. Die Massenmatrix ist positiv definit, die Dämpfungsmatrix und die Steifigkeitsmatrix sind zumindest positiv semidefinit.
E 3.2 (Allgemeines zur Diskretisierung) Wir diskretisieren die Zeitachse der
Einfachheit halber uniform mit der Schrittweite τ ,
tj = jτ,
j = 0, 1, 2, . . .
und bestimmen für diese Zeitpunkte Näherungslösungen
uj ≈ x(tj ),
vj ≈ ẋ(tj ),
aj ≈ ẍ(tj ).
3.2. Newmark-Verfahren
E 3.3 (Herleitung des Verfahrens) Wir betrachten zum Zeitpunkt tj+1 = tj +τ
das Gleichungssystem
M ẍ(tj+1 ) + C ẋ(tj+1 ) + Kx(tj+1 ) = f (tj+1 ),
x(tj+1 ) = x(tj ) + τ ẋ(tj ) +
τ2
ẍ(ξ1 ),
2
ẋ(tj+1 ) = ẋ(tj ) + τ ẍ(ξ2 ),
mit ξ1 , ξ2 ∈ (tj , tj+1 ) und diskretisieren es durch
M aj+1 + Cvj+1 + Kuj+1 = fj+1 := f (tj+1 ),
uj+1 = uj + τ vj + τ 2 21 − β aj + βaj+1 ,
vj+1 = vj + τ [(1 − γ)aj + γaj+1 ].
(3.17)
(3.18)
(3.19)
Die Parameter β und γ bestimmen die Genauigkeit und Stabilität des Verfahrens.
24. März 2016
31
E 3.4 (Implementierung) Für die Implementierung stellt man das Gleichungssystem nach einer der drei Variablen uj+1 , vj+1 oder aj+1 um und löst dieses.
Bathe [Bat02] stellt mit Hilfe von (3.18) den Vektor aj+1 als Funktion von uj+1
dar, mit dessen Hilfe folgt aus (3.19) eine Darstellung von vj+1 als Funktion
von uj+1 und aus diesen beiden Beziehungen und (3.17) ein Gleichungssystem
für uj+1 , das dann zu lösen ist.
Im Gegensatz dazu setzt Hughes (3.18) und (3.19) in (3.17) ein und bezieht alles
auf aj+1 . Zu lösen ist dann in jedem Zeitschritt das lineare Gleichungssystem
(M + γτ C + βτ 2 K) aj+1 = f˜j+1
(3.20)
τ2
= fj+1 − C(vj + (1 − γ)τ aj ) − K uj + τ vj + (1 − 2β)aj .
2
(3.21)
mit
f˜j+1
E 3.5 (Stabilität des Newmark-Verfahrens) Die Matrix A (vgl. Satz 2.29) ist
für das Newmark-Verfahren nicht leicht herzuleiten; sie hängt von der höchsten
auftretenden Frequenz ω des ungedämpften Systems ab. Es ergeben sich aber
folgende Resultate [Hug00, 9.1.2]:
• Für
2β ≥ γ ≥
1
2
(3.22)
ist das Verfahren unbedingt stabil.
• Für γ ≥ 12 und 2β < γ ist das Verfahren stabil unter einer Beschränkung
der Schrittweite, τ ≤ τ0 (ω, β, γ).
E 3.6 (Konsistenzfehler) Der Konsistenzfehler ist von der Ordnung 2 für γ =
und von der Ordnung 1 für γ 6= 21 .
1
2
Werte γ =
6 12 haben nur deshalb eine Bedeutung, weil man damit hohe FehlerFrequenzen, die aus der Diskretisierung im Ort entstehen, dämpfen kann. Um
stets eine Dämpfung großer Frequenzen zu erzielen, fordert man eine strengere
Stabilitätsbedingung als (3.22), nämlich
γ ≥ 12 ,
β≥
1
4
γ+
1 2
.
2
Man sieht aus
β ≥ 14 (γ + 12 )2 =
1
4
(γ − 12 )2 + 2γ ≥ 21 γ,
dass diese Bedingung tatsächlich strenger ist.
E 3.7 (Konvergenz) Konvergenz folgt aus Konsistenz und Stabilität.
32
24. März 2016
E 3.8 (Methode der konstanten mittleren Beschleunigung) Ein
Spezialfall entsteht durch die Wahl γ = 12 und β = 14 :
wichtiger
M aj+1 + Cvj+1 + Kuj+1 = fj+1 ,
uj+1 = uj + τ vj + 14 τ 2 (aj + aj+1 ),
vj+1 = vj + 21 τ (aj + aj+1 ).
Dieses Verfahren ist implizit, unbedingt stabil und konvergent von 2. Ordnung.
Das Verfahren entsteht auch, wenn man annimmt, dass im Intervall (tj , tj+1 )
die Beschleunigung konstant gleich 12 (aj + aj+1 ) ist. Deshalb wird es auch als
Methode der konstanten mittleren Beschleunigung bezeichnet.
E 3.9 (Trapezregel) Das Newmark-Verfahren mit γ = 12 und β = 14 wird auch
als Trapezregel bezeichnet. Wir wollen im Folgenden die Beziehung zu dem in
Erklärung 2.14 besprochenen impliziten Einschrittverfahren aufzeigen.
Wir betrachten die Differentialgleichung
M ẍ + C ẋ + Kx = f
als System erster Ordnung für die Variable y = [x, ẋ]T ,
ẏ =
ẋ
ẍ
=
ẋ
−1
M (f − C ẋ − Kx)
und approximieren dieses mit der Trapezregel aus 2.14,
uj+1
uj
=
+
vj+1
vj
τ
vj
vj+1
+
+
M −1 (fj − Cvj − Kuj )
M −1 (fj+1 − Cvj+1 − Kuj+1 )
2
(3.23)
Dann definieren wir
aj = M −1 (fj − Cvj − Kuj ),
das ist die erste Formel des Newmark-Verfahrens, und vereinfachen (3.23) zu
uj+1
vj+1
=
uj
vj
+
τ
2
vj
aj
+
vj+1
aj+1
.
Die zweite Zeile ist die dritte Formel des Newmark-Verfahrens. Durch Einsetzen
der zweiten Zeile in die erste,
uj+1 = uj + 12 τ vj + 12 τ vj+1
= uj + 12 τ vj + 12 τ [vj + 12 τ (aj + aj+1 )]
= uj + τ vj + 41 τ 2 (aj + aj+1 ),
erhält man die zweite Formel des Newmark-Verfahrens.
24. März 2016
33
E 3.10 (Störmer-Verlet-Verfahren, Zentrales Differenzenverfahren)
Ein zweiter wichtiger Spezialfall des Newmark-Verfahrens entsteht durch die
Wahl von γ = 12 und β = 0:
M aj+1 + Cvj+1 + Kuj+1 = fj+1 ,
uj+1 = uj + τ vj + 21 τ 2 aj ,
vj+1 = vj +
1
2 τ (aj
+ aj+1 ).
(3.24)
(3.25)
Löst man dieses System nach aj+1 auf, vergleiche (3.20), dann erhält man
(M + 1 τ C) aj+1 = f˜j+1 .
2
Damit ist das Störmer-Verlet-Verfahren3 ein explizites Verfahren, wenn M und
C Diagonalmatrizen sind. Das Verfahren ist nur bedingt stabil. Man kann zeigen,
dass die Stabilitätsbedingung
2
τ≤
ω
lautet, wobei ω die höchste auftretende Frequenz ist.
E 3.11 (Zentrales Differenzenverfahren) Der Name zentrales Differenzenverfahren kommt daher, dass man aus (3.24) und (3.25) auch
uj+1 − uj−1
uj+1 − 2uj + uj−1
,
aj =
vj =
2τ
τ2
schlussfolgern kann [Hug00, S. 495]. (Übung)
E 3.12 (Einschritt- oder Mehrschrittverfahren) Das Newmark-Verfahren, wie
wir es hier eingeführt haben, entspricht einem Einschrittverfahren. Man kann es
aber auch als Zweischrittverfahren
(M + γτ C + βτ 2 K)uj+1 + −2M + (1 − 2γ)τ C + ( 21 − 2β + γ)τ 2 K uj +
+ M − (1 − γ)τ C + ( 1 + β − γ)τ 2 K uj−1 = τ 2 f¯j
2
mit der einzigen Variablen u formulieren.
Man beachte aber, dass das Newmarkverfahren nicht für beliebige Differentialgleichungen 1. Ordnung angewendet werden kann. Es werden nicht alle
Komponenten der Verktorfunktion y = [x, ẋ]T gleich behandelt. Deshalb passt
das Verfahren auch nicht in die Theorie von Dahlquist, vergleiche Satz 2.36.
Ü 3.13 Man bestimme f¯j als Funktion von fj−1 , fj und fj+1 [Hug00, S. 527].
E 3.14 Man kann das Newmark-Verfahren auch auf 3. Ordnung verallgemeinern
[ZTZ05, 17.3.3]:
M aj+1 + Cvj+1 + Kuj+1 = fj+1
1
uj+1 = uj + τ vj + 12 τ 2 aj + τ 3 [β3 bj+1 + (1 − β3 )bj ]
6
1 2
vj+1 = vj + τ aj + τ [β2 bj+1 + (1 − β2 )bj ]
2
aj+1 = aj + τ (β1 bj+1 + (1 − β1 )bj )
Erkennen sie das Konstruktionsprinzip?
3
nach Carl Størmer (1874–1957, norwegischer Mathematiker und Geophysiker) und Loup
Verlet (geb. 1931, französischer Physiker), erste Beschreibung des Verfahrens aber schon
1687 durch Isaac Newton (1642–1726, englischer Naturforscher und Philosoph)
34
24. März 2016
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
exact, ω beliebig
Newmark, γ=0.5, β=0.25, ω=1, τ=2π/20
Newmark, γ=0.5, β=0.25, ω=1, τ=2π/10
Newmark, γ=0.5, β=0.25, ω=1, τ=2π/3
v/ω
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
0
exact, ω=1
Newmark, γ=0.5, β=0.25, ω=1, τ=2π/20
Newmark, γ=0.5, β=0.25, ω=1, τ=2π/10
Newmark, γ=0.5, β=0.25, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 15: Phasendiagramm und t, u-Diagramm zum Newmark-Verfahren
mit γ = 0.5, β = 0.25 (Trapezregel)
Bsp 3.15 Wir testen das Newmark-Verfahren mit verschiedenen Parametern
am Beispiel aus 2.15,
ü(t) + ω 2 u(t) = 0,
u(0) = 1,
u̇(0) = 0,
mit ω = 1.
• Für die Trapezregel (γ = 0.5, β = 0.25) erkennen wir in Abbildung 15 die
nicht vorhandene Dämpfung, die unbedingte Stabilität und eine Periodenverlängerung.
• Für das Störmer-Verlet-Verfahren (γ = 0.5, β = 0) erhalten wir ebenfalls
geringe Dämpfung, aber fehlende Stabilität für zu große Schrittweite (τ =
2π/3 ist schon instabil) und eine Periodenverkürzung, siehe Abbildung 16.
• Beim gedämpften Newmark-Verfahren (γ = 0.6) unter Einhaltung der
Stabilitätsbedingung (β = 0.3025) sehen wir in Abbildung 17 deutlich die
Dämpfung, aber auch die unbedingte Stabilität.
24. März 2016
35
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
exact, omega beliebig
Newmark, γ=0.5, β=0.0, ω=1, h=2π/20
Newmark, γ=0.5, β=0.0, ω=1, h=2π/10
Newmark, γ=0.5, β=0.0, ω=1, h=2π/3
v/ω
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
0.5
u
0
exact, omega=1 beliebig
Newmark, γ=0.5, β=0.0, ω=1, h=2π/20
Newmark, γ=0.5, β=0.0, ω=1, h=2π/10
Newmark, γ=0.5, β=0.0, ω=1, h=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 16: Phasendiagramm und t, u-Diagramm zum Newmark-Verfahren
mit γ = 0.5, β = 0.0 (Störmer-Verlet-Verfahren)
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
exact, omega beliebig
Newmark, γ=0.6, β=0.3025, ω=1, τ=2π/20
Newmark, γ=0.6, β=0.3025, ω=1, τ=2π/10
Newmark, γ=0.6, β=0.3025, ω=1, τ=2π/3
v/ω
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
0.5
u
0
exact, omega=1 beliebig
Newmark, γ=0.6, β=0.3025, ω=1, τ=2π/20
Newmark, γ=0.6, β=0.3025, ω=1, τ=2π/10
Newmark, γ=0.6, β=0.3025, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 17: Phasendiagramm und t, u-Diagramm zum Newmark-Verfahren
mit γ = 0.6, β = 0.3025
36
24. März 2016
3.3. Houbolt-Verfahren
E 3.16 (Houbolt-Verfahren) Das implizite Dreischrittverfahren
M aj+1 + Cvj+1 + Kuj+1 = fj+1
1
vj+1 =
[11uj+1 − 18uj + 9uj−1 − 2uj−2 ]
6τ
1
aj+1 = 2 [2uj+1 − 5uj + 4uj−1 − uj−2 ]
τ
ist unbedingt stabil und konvergent von 2. Ordnung, [Hug00, 9.3.3].
Bsp 3.17 Wir testen das Houbolt-Verfahren an unserem Standardbeispiel 2.15
mit ω = 1. In Abbildung 18 erkennen wir starke Dämpfung, die unbedingte
Stabilität und eine Periodenverlängerung.
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
v/ω
0.2
exact, ω beliebig
Houbolt, τ=2π/20
Houbolt, τ=2π/10
Houbolt, τ=2π/3
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
0
exact, ω=1
Houbolt, τ=2π/20
Houbolt, τ=2π/10
Houbolt, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 18: Phasendiagramm und t, u-Diagramm zum Houbolt-Verfahren
Ü 3.18 Wir rechnen die Konvergenz zweiter Ordnung nach. Es gilt
...
x(t − τ ) = x(t) − τ ẋ(t) + 12 τ 2 ẍ(t) + 16 τ 3 x (t) + O(τ 4 ),
...
x(t − 2τ ) = x(t) − 2τ ẋ(t) + 2τ 2 ẍ(t) + 34 τ 3 x (t) + O(τ 4 ),
3 ...
x (t) + O(τ 4 ),
x(t − 3τ ) = x(t) − 3τ ẋ(t) + 29 τ 2 ẍ(t) + 27
6 τ
also
1
[11x(t) − 18x(t − τ ) + 9x(t − 2τ ) − 2x(t − 3τ )] = ẋ(t) + O(τ 2 ),
6τ
1
[2x(t) − 5x(t − τ ) + 4x(t − 2τ ) − x(t − 3τ )] = ẍ(t) + O(τ 2 ).
τ2
24. März 2016
37
3.4. Kollokationsverfahren
E 3.19 (Rechenvorschrift) [Hug00, S. 530] Das Gleichungssystem
M aj+θ + Cvj+θ + Kuj+θ = fj+θ
fj+θ = (1 − θ)fj + θfj+1
uj+θ = uj + θτ vj + (θτ )2 [( 12 − β)aj + βaj+θ ]
vj+θ = vj + θτ [(1 − γ)aj + γaj+θ ]
definiert zunächst aj+θ und mit
aj+θ = (1 − θ)aj + θaj+1
die Variable aj+1 . Die Werte für uj+1 und vj+1 erhält man durch die NewmarkFormeln (3.18) und (3.19),
uj+1 = uj + τ vj + τ 2 12 − β aj + βaj+1 ,
vj+1 = vj + τ [(1 − γ)aj + γaj+1 ].
E 3.20 (Kollokations-Parameter) Der Parameter θ heißt Kollokations-Parameter. Typisch sind Werte im Bereich 1,3 bis 1,5.
• Für θ = 1 entsteht das Newmark-Verfahren.
• Für γ =
1
2
und β =
1
6
entsteht die Wilsonsche θ-Methode .
E 3.21 (Eigenschaften des Kollokationsverfahrens) Genau für γ =
Verfahren konvergent mit der Ordnung 2. Gilt zusätzlich θ ≥ 1 und
1
2
ist das
θ
2θ2 − 1
≥β≥
2(θ + 1)
4(2θ3 − 1)
dann ist das Verfahren unbedingt stabil.
Ü 3.22 (Stabilität der Wilsonschen θ-Methode) Für welche Werte von θ ist
die Wilsonsche θ-Methode unbedingt stabil?
E 3.23 Die besten Verfahren erhält man für β ∈ [ 16 , 14 ] und ein zugehöriges,
in Abhängigkeit von β zu wählendes θ = θ(β), siehe [Hug00, Tabelle 9.3.1].
Beispiele:
β
1
6
1
5
1
4
θ(β)
1,42
1,16
1,00
Bsp 3.24 Wir testen das Kollokationsverfahren mit verschiedenen Parametern
am Beispiel aus 2.15,
ü(t) + ω 2 u(t) = 0,
u(0) = 1,
u̇(0) = 0,
mit ω = 1.
• Für die Parameter γ = 12 , β = 15 , θ = 1,16 erkennen wir in Abbildung 19 geringe Dämpfung, die unbedingte Stabilität und eine Periodenverlängerung.
• Für die Wilsonsche θ-Methode (γ = 12 , β = 16 ) mit θ = 1,42 sehen wir
in Abbildung 20 die etwas stärkere Dämpfung und ebenfalls unbedingte
Stabilität sowie Periodenverlängerung.
38
24. März 2016
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
exact, ω beliebig
Kollokation, θ=1.16, ω=1, τ=2π/20
Kollokation, θ=1.16, ω=1, τ=2π/10
Kollokation, θ=1.16, ω=1, τ=2π/3
v/ω
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
0
exact, ω=1
Kollokation, θ=1.16, ω=1, τ=2π/20
Kollokation, θ=1.16, ω=1, τ=2π/10
Kollokation, θ=1.16, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 19: Phasendiagramm und t, u-Diagramm zum Kollokationsverfahren
mit γ = 0,5, β = 0,2 und θ = 1,16
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
exact, ω beliebig
Kollokation, θ=1.42, ω=1, τ=2π/20
Kollokation, θ=1.42, ω=1, τ=2π/10
Kollokation, θ=1.42, ω=1, τ=2π/3
v/ω
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
0
exact, ω=1
Kollokation, θ=1.42, ω=1, τ=2π/20
Kollokation, θ=1.42, ω=1, τ=2π/10
Kollokation, θ=1.42, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 20: Phasendiagramm und t, u-Diagramm zum Wilsonschen θVerfahren mit θ = 1,42
24. März 2016
39
3.5. Hilber-Hughes-Taylor-Verfahren (HHT-Verfahren, α-Verfahren)
E 3.25 Es gibt verschiedene Möglichkeiten, das Verfahren darzustellen. Wir
verwenden hier, anders als in [Hug00], die Variante aus der Originalarbeit mit
α ∈ [0, 1), wobei man für α = 0 das Newmark-Verfahren erhält,
M aj+1 + (1 − α)Cvj+1 + αCvj + (1 − α)Kuj+1 + αKuj
= f ((1 − α)tj+1 + αtj ).
Vervollständigt wird der Formelsatz durch die Newmark-Formeln (3.18) und
(3.19),
uj+1 = uj + τ vj + τ 2 12 − β aj + βaj+1 ,
vj+1 = vj + τ [(1 − γ)aj + γaj+1 ].
Wählt man α ∈ 0, 13 , γ = 12 + α und β = 14 (1 + α)2 , so ergibt sich ein unbedingt
stabiles Verfahren mit der Konvergenzordnung 2. Je größer α ist, desto größer
ist die numerische Dissipation (Dämpfung der Schwingung).
Ü 3.26 (Alte Klausuraufgabe) Berechnen Sie für die Anfangswertaufgabe
ẍ(t) + 0,1 ẋ(t) + x(t) = 0,
x(0) = 1,
ẋ(0) = 0
einen Iterationsschritt mit dem Hilber-Hughes-Taylor-Verfahren. Verwenden Sie
die Parameter α = 0,1, β = 0,3025, γ = 0,6 sowie die Zeitschrittweite τ = 0,1.
Bsp 3.27 Wir testen das HHT-Verfahren mit kleinem (α = 0,05) und großem
(α = 0,3) Parameter sowie γ = 12 + α und β = 14 (1 + α)2 an unserem StandardBeispiel aus 2.15 mit ω = 1.
• Für α = 0,05 sehen wir in Abbildung 21 Dämpfung, unbedingte Stabilität
und Periodenverlängerung.
• Für α = 0,3 erkennen wir in Abbildung 22 die etwas stärkere Dämpfung
und ebenfalls unbedingte Stabilität sowie Periodenverlängerung.
Bem 3.28 (Andere Schreibweise) Analog zu den Kollokationsverfahren kann
man das α-Verfahren auch in der Form
M aj+1 + Cvj+α + Kuj+α = fj+α ,
uj+α = (1 − α)uj+1 + αuj ,
vj+α = (1 − α)vj+1 + αvj ,
fj+α = f ((1 − α)tj+1 + αtj ),
schreiben, wobei zur Vervollständigung des Systems noch die Newmark-Formeln
(3.18) und (3.19) hinzugefügt werden müssen.
40
24. März 2016
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
v/ω
0.2
exact, ω beliebig
HHT, α=0.05, ω=1, τ=2π/20
HHT, α=0.05, ω=1, τ=2π/10
HHT, α=0.05, ω=1, τ=2π/3
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
0
exact, ω=1
HHT, α=0.05, ω=1, τ=2π/20
HHT, α=0.05, ω=1, τ=2π/10
HHT, α=0.05, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 21: Phasendiagramm und t, u-Diagramm zum HHT-Verfahren mit
α = 0,05
Phasendiagramm, t ∈ [0, 2π]
1
0.8
0.6
0.4
v/ω
0.2
exact, ω beliebig
HHT, α=0.3, ω=1, τ=2π/20
HHT, α=0.3, ω=1, τ=2π/10
HHT, α=0.3, ω=1, τ=2π/3
0
−0.2
−0.4
−0.6
−0.8
−1
−1
−0.5
0
u
0.5
1
Auslenkung, t ∈ [0, 10π]
1
u
0.5
0
exact, ω=1
HHT, α=0.3, ω=1, τ=2π/20
HHT, α=0.3, ω=1, τ=2π/10
HHT, α=0.3, ω=1, τ=2π/3
−0.5
−1
0
5
10
15
20
25
30
t
Abbildung 22: Phasendiagramm und t, u-Diagramm zum HHT-Verfahren mit
α = 0,3
24. März 2016
41
E 3.29 (Verallgemeinerte α-Verfahren) Die Herausforderung bei der Konstruktion von Verfahren für Schwingungsprobleme besteht darin, dass man hohe Frequenzen im Allgemeinen dämpfen möchte, aber niedrige Frequenzen nicht. Man
kann die α-Verfahren als Verallgemeinerung der Newmark-Verfahren auffassen,
wobei der neue Parameter α eingeführt wurde, um das Dämpfungsverhalten zu
steuern, ohne die Approximation zweiter Ordnung zu verlieren.
Bei den Verallgemeinerten α-Verfahren nach Chung und Hubert (1993) wurde
dazu noch ein weiterer Parameter αM eingeführt:
M aj+αM + Cvj+α + Kuj+α = fj+α ,
uj+α = (1 − α)uj+1 + αuj ,
vj+α = (1 − α)vj+1 + αvj ,
aj+αM = (1 − αM )aj+1 + αM aj ,
fj+α = f ((1 − α)tj+1 + αtj ),
wobei zur Vervollständigung des Systems ebenfalls die Newmark-Formeln (3.18)
und (3.19) hinzugefügt werden müssen. Für αM = 0 entsteht die HHT-Methode.
3.6. Diskussion der Verfahren
E 3.30 Ein konkurrenzfähiges Verfahren sollte folgende Eigenschaften besitzen
[Hug00, S. 535 ff.]:
1. Unbedingte Stabilität, wenn es auf lineare Probleme angewandt wird. Bedingt stabile Verfahren benötigen oft kleine Zeitschritte, die von hohen
Frequenzen bestimmt werden, obwohl oft nur niedrige Frequenzen aufgelöst werden sollen. Unbedingte Stabilität erzwingt allerdings implizite
Verfahren.
2. Pro Iteration sollte nur ein lineares Gleichungssystem für eine der Variablen u, v, a, . . . gelöst werden müssen und nicht ein gekoppeltes System aus
diesen Variablen. Es muss ohnehin schon genug gespeichert und gerechnet
werden. Es sind aber Verfahren mit größeren linearen Gleichungssystemen
entwickelt worden.
3. Konvergenz 2. Ordnung: Dahlquist hat bewiesen, dass es kein unbedingt
A-stabiles lineares Mehrschrittverfahren mit Konvergenzordnung größer
als 2 gibt (zweite Dahlquist-Barriere). Die kleinste Fehlerkonstante unter
allen unbedingt stabilen Verfahren 2. Ordnung hat die Trapezregel.
4. Steuerbare numerische (algorithmische) Dissipation ist wünschenswert, um
unphysikalische, hohe Frequenzen zu dämpfen. Steuerbar ist die numerische
Dissipation bei den Kollokationsmethoden und bei den α-Verfahren. Die
Trapezregel besitzt hingegen gar keine numerische Dissipation.
5. Keine Anlaufrechnung. Diese benötigt mehr Implementations- und Analyseaufwand, Mehrschrittverfahren benötigen auch mehr Speicheraufwand,
[Hug00].
Der beste Algorithmus für Schwingungsprobleme in der Strukturdynamik ist
folglich ein unbedingt stabiles Einschrittverfahren 2. Ordnung, das optimal
zwischen effektiver numerischer Dissipation und hoher Genauigkeit (kleiner
Fehlerkonstante) balanciert.
42
24. März 2016
E 3.31 Alle in 3.1 bis 3.4 behandelten Verfahren sind (ggf. nach Einstellung der
richtigen Parameter)
•
•
•
•
implizit,
unbedingt stabil,
konvergent von 2. Ordnung,
dissipativ (außer Trapezregel)
und benötigen etwa den gleichen numerischen Aufwand.
E 3.32 (Dämpfung hoher Frequenzen) Der Spektralradius ρ (das ist der Betrag des betragsgrößten Eigenwertes) der Transfer-Matrix ist ein Maß für die
Dämpfung; je kleiner dieser ist, desto stärker wirkt die numerische Dämpfung
des Verfahrens: für ρ = 1 ist keine Dämpfung vorhanden. In Abbildung 23 ist
der Spektralradius gegen τ /T aufgetragen, wobei T = 2π/ω die Periode der
exakten Lösung x(t) der Schwingungsdifferentialgleichung
ẍ(t) + ω 2 x(t) = 0
(ohne Dämpfung, c = 0) ist. Mit einer gewissen Schrittweite τ kann man niedrige
Frequenzen auflösen (τ /T klein), aber hohe Frequenzen (τ /T groß) nicht.
Abbildung 23: Spektralradien der Transfer-Matrix im single-frequency-mode
für verschiedene Verfahren [Hug00, S. 533], τ := ∆t, α := −α
Im Diagramm 23 ist erkennbar, dass der Unterschied der vorgestellten Verfahren
im Niederfrequenzbereich nur gering ist, sich die einzelnen Verfahren aber im
hochfrequenten Bereich deutlich unterscheiden:
• Es wirkt keine Dämpfung für die Trapezregel, d. h. für das NewmarkVerfahren mit γ = 21 , β = 41 .
• Starke Dämpfung hoher Frequenzen erhält man für das Houbolt-Verfahren,
das Kollokations-Verfahren und für die α-Methode mit großem α. (Man
beachte, dass α ∈ [0, 13 ) gewählt werden sollte.)
24. März 2016
43
E 3.33 (Dämpfung niederer Frequenzen) Das Ziel der Rechnung ist, niedere
Frequenzen möglichst exakt zu berechnen. Daher ist im niederfrequenten Bereich
eine starke Dämpfung nachteilig.
Sei T die Periode der exakten Lösung x(t) der Schwingungsdifferentialgleichung ohne Dämpfung (c = 0), uτ (t) die mit der Schrittweite τ berechnete
Näherungslösung und Tτ die numerisch bestimmte Periodendauer (was in der
Praxis nicht einfach ist), siehe auch die Illustration in Abbildung 24.
1
2 pi xi
0
−1
0
T/2
T
T_tau
t
Abbildung 24: Illustration des Periodenfehlers
In der Abbildung 25 ist das algorithmische Dämpfungsverhältnis ξ := [x(T ) −
uτ (Tτ )]/2π gegen τ /T aufgetragen.
Aus Abbildung 25 erkennt man:
• Die Trapezregel ist optimal, sie dämpft überhaupt nicht.
• Die starke Dämpfung kleiner Frequenzen ist der Hauptnachteil gedämpfter
Newmark-Verfahren (γ > 12 ).
• HHT-Verfahren waren schon bei hohen Frequenzen besser als die Kollokationsverfahren. Das gleiche trifft auch bei niedrigen Frequenzen zu.
• Bei den HHT-Verfahren ist kleines α besser als großes, bei großen Frequenzen war es genau umgekehrt.
• Das Houbolt-Verfahren ist nicht konkurrenzfähig.
E 3.34 (Relativer Periodenfehler) Wie schon in Erklärung 3.33 dargestellt, ist
die numerisch bestimmte Periode Tτ im Allgemeinen von der exakten Periode T
verschieden. Die Größe (Tτ − T )/T bezeichnet man als relativen Periodenfehler.
Dieser ist in den Abbildungen 26 und 27 für verschiedene Verfahren gegen die
relative Schrittweite τ /T aufgetragen.
Aus Abbildung 26 kann man schlussfolgern, dass die Trapezregel oder eine
α-Methode mit kleinem α geeignete Verfahren sind, wenn der relative Periodenfehler von Relevanz ist. Diese dämpfen allerdings hohe Frequenzen schlecht.
In Abbildung 27 wird der relative Periodenfehler für verschiedene ungedämpfte
Newmark-Verfahren (γ = 12 ) verglichen. (Zur Erinnerung: Newmark mit γ = 12
besitzt keine algorithmische Dämpfung. Für γ > 12 wird zwar eine Dämpfung
erreicht, jedoch reduziert sich die Konvergenzordnung von 2 auf 1.) Man sieht,
dass kleineres β zu kleinerem Periodenfehler führt. Allerdings sind NewmarkVerfahren mit β < 14 nur bedingt stabil. Bei zu großer Schrittweite explodiert
der Periodenfehler.
44
24. März 2016
0.14
Houbolt
Trapezregel
alpha Methode alpha=0.05
Kollokation theta= 1.4, beta= 1/6
Newmark beta=0.3025, gamma=0.6
alpha Methode alpha=0.3
0.12
0.1
Zeta
0.08
0.06
0.04
0.02
0
0
0.1
0.2
0.3
0.4
0.5
Zeitschrittweite
0.6
0.7
0.8
0.9
1
Abbildung 25: Grad der algorithmischen Dämpfung für α-Methoden,
Kollokations-Verfahren und die Houbolt-, Park-, und WilsonMethode. Oben: Kopie aus [Hug00, S. 534] (τ := ∆t, α := −α),
unten: eigene Rechnung (Blank 2011)
24. März 2016
45
0.9
Houbolt
Trapezregel
alpha Methode alpha=−0.05
Kollokation theta= 1.4, beta= 1/6
Kollokation theta= 1.215798, beta= 0.19
alpha Methode alpha=0.3
Kollokation theta= 1.420815, beta= 1/6
0.8
0.7
0.6
(Ttau −T)/T
0.5
0.4
0.3
0.2
0.1
0
−0.1
0
0.05
0.1
0.15
0.2
Zeitschrittweite
0.25
0.3
0.35
0.4
Abbildung 26: Relative Periodenfehler für verschiedene Verfahren. Oben: Kopie
aus [Hug00, S. 535], (α := −α), unten: eigene Rechnung (Blank
2011)
46
24. März 2016
1.6
Newmark beta=1/12
Newmark beta=1/6
Newmark beta=1/5
Newmark beta=1/4
1.4
1.2
1
(Ttau −T)/T
0.8
0.6
0.4
0.2
0
−0.2
−0.4
0
0.1
0.2
0.3
0.4
0.5
Zeitschrittweite
0.6
0.7
0.8
0.9
1
Abbildung 27: Periodenfehler für ungedämpfte Newmark-Verfahren (γ = 12 ) im
Vergleich mit dem Wilson- und Houbolt-Verfahren. Oben: Kopie
aus [Hug00, S. 536] , unten: eigene Rechnung (Blank 2011)
24. März 2016
47
E 3.35 (Überschwingen) Wir betrachten hohe Frequenzbereiche, die durch
numerische Verfahren mit
τ = 10 T
nicht aufgelöst werden können. Diese sollten die numerische Lösung möglichst
nicht instabil werden lassen. Jeder einzelne Anteil einer Schwingung mit hoher
Frequenz sollte sich nicht im Laufe der Zeit “aufschaukeln”. In Abbildung 28
werden verschiedene Verfahren daraufhin untersucht.
Abbildung 28: Vergleich des Überschwingens veschiedener Verfahren [Hug00, S.
539]
48
24. März 2016
Folgende Schlussfolgerungen können gezogen werden:
• Bei der Trapez-Regel ist kein Überschwingen erkennbar, weder in u noch
in v.
• Beim gedämpften Newmark- und beim α-Verfahren ist kein Überschwingen
in u, jedoch mildes Überschwingen in v erkennbar.
1
100
0.5
50
Auslenkung
Auslenkung
• Bei den Kollokations-Verfahren und beim Wilson-Verfahren ist ein mildes
Überschwingen in v erkennbar, aber für u werden unbrauchbare Ergebnisse
berechnet.
0
−0.5
−1
0
−50
0
50
100
150
200
Zeit
250
300
350
−100
400
0
50
100
Trapezregel
Geschwindigkeit
−0.5
Auslenkung
0
50
100
150
200
Zeit
250
300
350
350
400
0
−1
1
400
0.5
200
0
−0.5
0
50
100
150
200
Zeit
250
300
350
400
0
50
100
150
200
Zeit
250
300
350
400
0
−200
0
50
100
150
200
Zeit
250
300
350
−400
400
Kollokation beta=0.20, theta=1.159772
Newmark gamma=0.6, beta=0.3025
6
Geschwindigkeit
2
Geschwindigkeit
300
1
−2
400
Auslenkung
Geschwindigkeit
0
1
0
−1
−2
250
2
0.5
−1
200
Zeit
Kollokation beta=0.24, theta=1.021712
1
−1
150
0
50
100
150
200
Zeit
250
300
350
4
2
0
−2
−4
400
0
50
100
150
200
Zeit
250
300
350
400
0
50
100
150
200
Zeit
250
300
350
400
350
400
600
1
400
Auslenkung
Auslenkung
0.5
0
200
0
−200
−0.5
−400
−1
−600
0
50
100
150
200
Zeit
250
300
350
400
Wilson = Kollokation beta=1/6, theta=1.4
Hughes alpha=0.1
6
Geschwindigkeit
Geschwindigkeit
2
1
0
−1
−2
4
2
0
−2
−4
0
50
100
150
200
Zeit
250
300
350
400
0
50
100
150
200
Zeit
250
300
Abbildung 29: Vergleich des Überschwingens veschiedener Verfahren, eigene
Rechnung (Blank 2011)
24. März 2016
49
4. Numerische Lösung von Eigenwertproblemen
4.1. Matrizen mit spezieller Struktur
Def 4.1 Eine reelle n-reihige Matrix A heißt orthogonal, wenn A−1 = AT bzw.
AAT = AT A = I gilt.
Satz 4.2 (Eigenschaften orthogonaler Matrizen) Für jede orthogonale Matrix
gilt:
1.
2.
3.
4.
Die Spalten bilden ein Orthonormalsystem.
Die Zeilen bilden ein Orthonormalsystem.
Der Betrag der Eigenwerte ist gleich 1.
Der Betrag der Determinante ist gleich 1.
Beweis Die ersten beiden Eigenschaften folgen aus AAT = AT A = I. Die dritte
Eigenschaft soll im Rahmen von Aufgabe 4.12 selbst bewiesen werden. Die vierte
Eigenschaft folgt aus der dritten, da die Determinante gleich dem Produkt der
Eigenwerte ist.
q.e.d.
E 4.3 Günstig für die numerische Praxis sind Ähnlichkeitstransformationen mit
orthogonalen Matrizen Y , da dann Y −1 = Y T gilt. Man spricht dann von einer
orthogonalen Ähnlichkeitstransformation B = Y T AY .
Oft werden orthogonale Matrizen mit Q bezeichnet, so dass man meist
B = QT AQ
liest, wenn es um eine orthogonale Ähnlichkeitstransformation geht.
Satz 4.4 (QR-Faktorisierung) Jede reelle Matrix A kann in das Produkt aus
einer orthogonalen Matrix Q und einer rechten oberen Dreiecksmatrix R zerlegt
werden,
A = Q · R.
Die QR-Faktorisierung hat folgende Eigenschaften:
1. Die Überführung von A in R mit Hilfe von orthogonalen Transformationen
kann stabiler ausgeführt werden als bei der LU-Faktorisierung.
2. Das Problem ist äquivalent zur Überführung einer Basis a1 , . . . , an in eine
Orthonormalbasis q 1 , . . . , q n .
A = a1 a2 · · · an ,
Q = q1 q2 · · · qn
3. Die Faktorisierung ist auch für Rechteckmatrizen möglich:
A
=
=
50
Q
R
A
bzw.
=
Q
R
=
24. März 2016
E 4.5 (Gram-Schmidt-Orthogonalisierung) Die Spalten der gegebenen Matrix
A und der gesuchten Matrix Q werden wie folgt bezeichnet:
Q = q1 q2 · · · qn .
A = a1 a2 · · · an ,
Der Gram-Schmidt-Algorithmus ist vermutlich aus der linearen Algebra bekannt:
1: q = ka1 k a1
1
1
2: for k = 1(1)n − 1 do
P
3:
q̃ k+1 = ak+1 − kj=1 (q j , ak+1 )q j
q k+1 =
4:
1
kq̃ k+1 k q̃ k+1
end for
5:
Die Matrix R erhält man aus (q j , ak+1 ) und Normbildung oder aus R = QT A.
Weil Rundungsfehler zum Verlust der Orthogonalität führen, ist es besser, den
Schritt in Zeile 3 wie folgt auszuführen:
3.1:
3.2:
3.3:
3.4:
3.5:
(0)
q k+1 = ak+1
for j = 1(1)k do
(j)
(j−1)
(j−1)
q k+1 = q k+1 − (q j , q k+1 )q j
end for
(k)
q̃ k+1 = q k+1
Ü 4.6 Für Stabilitätsbetrachtungen studiere man [TB00, Sect. 9, Ex. 3]. Stabilere Algorithmen zur Berechnung der QR-Faktorisierung erhält man mit Hilfe
der Givens-Rotation oder Householder-Spiegelung, siehe Literatur.
Def 4.7 Eine reelle n-reihige Matrix heißt symmetrisch, wenn A = AT gilt.
Satz 4.8 (Eigenschaften symmetrischer Matrizen)
trische Matrix gilt:
Für jede reelle, symme-
1. Alle Eigenwerte sind reell.
2. Eigenvektoren zu verschiedenen Eigenwerten sind orthogonal.
3. Algebraische und geometrische Vielfachheit jedes Eigenwerts sind gleich.
Beweis Sei λ ∈ C ein Eigenwert einer reellen, symmetrischen Matrix A und
x der zugehörige Eigenvektor. Wir bezeichnen mit x den Vektor, der aus x
entsteht, wenn man bei jedem Element zur konjugiert komplexen Zahl übergeht.
Folglich gilt
Ax = Ax = λx = λ x
und damit
λxT x = (Ax)T x = xT AT x = xT Ax = λxT x.
Da der Eigenvektor x nicht der Nullvektor ist, gilt
X
X
xT x =
xi xi =
|xi |2 > 0.
i
i
Also können wir die vorhergehende Gleichung durch xT x teilen und erhalten
λ = λ, d. h., λ ist reell.
24. März 2016
51
Seien λ und µ zwei verschiedene Eigenwerte einer reellen, symmetrischen Matrix
A und seien x und y die zugehörigen Eigenvektoren, d. h. Ax = λx und Ay = µy.
Folglich gilt
λxT y = (Ax)T y = xT AT y = xT Ay = µxT y,
also (λ − µ)xT y = 0, was wegen λ 6= µ die Orthogonalität xT y = 0 nach sich
zieht.
Für den Beweis der dritten Eigenschaft verweisen wir auf die Literatur, z. B.
[MV03, §6.7, Satz 7.1].
q.e.d.
Satz 4.9 Zu einem k-fachen Eigenwert einer reellen, symmetrischen Matrix
können k orthogonale Eigenvektoren angegeben werden.
Beweis Der Beweis kann konstruktiv über das Schmidtsche Orthogonalisierungsverfahren geführt werden.
q.e.d.
Folgerung 4.10
1. Werden die (orthogonalen) Eigenvektoren x1 , . . . , xn einer reellen, symmetrischen n × n-Matrix auf die Länge Eins normiert,
qi =
xi
,
|xi |
und zur Modalmatrix
Q=
q1 . . . qn
zusammengefasst, so ist Q eine orthogonale Matrix.
2. Eine symmetrische Matrix kann mit Hilfe einer orthogonalen Ähnlichkeitstransformation diagonalisiert werden,
QT AQ = Λ = diag(λi ).
Man bezeichnet diese auch als Hauptachsentransformation.
E 4.11 Wenn man von einer Matrix weiß, dass sie symmetrisch ist, kann man
gleich mehrere angenehme Eigenschaften benutzen:
beliebige Matrix A
symmetrische Matrix A
Die Eigenwerte sind i. Allg. komplex.
Die Eigenwerte sind reell.
Die geometrische Vielfachheit kann
kleiner als die algebraische Vielfachheit sein.
Geometrische und algebraische Vielfachheit stimmen überein.
Die Eigenvektoren zu verschiedenen Eigenwerten sind zwar linear unabhängig,
aber nicht notwendig orthogonal.
Die Eigenvektoren zu verschiedenen
Eigenwerten sind orthogonal, d. h., es
existiert ein Orthonormalsystem aus
Eigenvektoren.
Die Ähnlichkeitstransformation hat die
Form A = XΛX −1 .
Die Ähnlichkeitstransformation hat die
Form A = QΛQT mit einer orthogonalen Matrix Q.
52
24. März 2016
Ü 4.12 Eine reelle n-reihige Matrix heißt schiefsymmetrisch, wenn
A = −AT
gilt. Man zeige analog zum Beweis von Satz 4.8, dass die Eigenwerte einer
schiefsymmetrischen Matrix sämtlich rein imaginär sind.
Mit der gleichen Beweisidee kann man auch zeigen, dass die Eigenwerte einer
orthogonalen Matrix auf dem Einheitskreis in der komplexen Ebene liegen, d. h.,
sie haben den Betrag 1.
Ü 4.13 Versuchen Sie aus dem Plot von Spektren die Eigenschaften der Matrizen
zu erkennen. Woran erkennen Sie, dass es sich um eine (a) komplexe, (b) reelle,
(c) symmetrische (hermitesche), (d) schiefsymmetrische (schiefhermitesche) bzw.
(e) orthogonale (unitäre) Matrix handelt?
4.2. Eigenwertberechnung und Nullstellensuche
E 4.14 (Betrachtungen zur Kondition)
• Kleine Änderungen von A implizieren kleine Änderungen der λj .
Wielandt-Hoffmann-Theorem: Sei A symmetrisch mit Eigenwerten λ1 ≤
· · · ≤ λn . Die gestörte symmetrische Matrix B = A + δA habe die Eigenwerte µ1 ≤ · · · ≤ µn . Dann gilt:
n
X
j=1
(µj − λj )2 ≤ kB − Ak2F :=
n
X
(bij − aij )2
i,j=1
Die Voraussetzung der Symmetrie ist wesentlich für diesen Satz. Im Allgemeinen gelten viel schwächere Resultate, siehe Übung 4.15.
• Bei richtiger Orientierung ändert sich vj nur wenig mit A, wenn λj einfach
und genügend von den übrigen Eigenwerten λi getrennt ist.
Gibt es dicht benachbarte Eigenwerte, dann reagiert vj empfindlich auf
Störungen in A.
• Kleine Änderungen der Koeffizienten an eines Polynoms können zu großen
Änderungen in den Nullstellen λj führen, siehe Beispiel 4.16. Folglich sollten
Eigenwerte nicht über das charakteristische Polynom berechnet werden.
(Das charakteristische Polynom ist nützlich für theoretische Zwecke.)
Im Gegenteil, Nullstellen von Polynomen sollten als Eigenwerte der Begleitmatrix (engl. companion matrix ) des Polynoms gesucht werden, siehe
Lemma 4.17.
24. März 2016
53
Ü 4.15 Bestimmen Sie Eigenwerte der 10 × 10-Matrix

10 10

9 10


.

8 ..

Aε = 
..

. 10


2 10
ε
1










für ε = 0, eps = 10−6 und eps = 10−5 . Wie hängen die Eigenwerte von der
Größe der Störung ε ab?
Bsp 4.16 (Kondition der Nullstellensuche)
• λn = 0 hat n-fache Nullstelle λ = 0
• λn − 10−n = 0 hat n verschiedene Nullstellen λi mit |λi | = 10−1 .
Lem 4.17 Die Nullstellen des Polynoms
p(λ) = a0 + a1 λ + · · · + an−1 λn−1 + λn
sind die Eigenwerte der Begleitmatrix





A=



0
0
1
0
0
1
0
0
..
.
···
···
..
−a0 −a1 · · ·
···
0
0
.
···
1
−an−1





.



Beweis Durch Nachrechnen von det(λI − A) = p(λ).
q.e.d.
4.3. Der QR-Algorithmus
4.3.1. Grundform
Alg 4.18 Sei A ∈ Rn×n eine gegebene Matrix.
• Initialisierung: Wähle eine orthogonale Matrix Q(0) ∈ Rn×n , zum Beispiel
Q(0) = I, und berechne T (0) = (Q(0) )T AQ(0) .
• Iteriere bis zur Konvergenz
– Berechne QR-Faktorisierung von T (k−1) = Q(k) R(k) .
– Setze T (k) = R(k) Q(k) .
54
24. März 2016
Folgerung 4.19 Es gilt also
T (k) = (Q(k) )T T (k−1) Q(k) .
Das ist eine (orthogonale) Ähnlichkeitstransformation, d. h., die Eigenwerte
bleiben erhalten.
Satz 4.20 (Konvergenz der QR-Methode) [GvL07] Sei A ∈ Rn×n mit ausschließlich reellen Eigenwerten |λ1 | > |λ2 | > · · · > |λn |. Dann gilt


λ1 t12 t13 · · ·
t1n

λ2 t23 · · ·
t2n 




..
..
..
lim T (k) = 
.
.
.
.


k→∞

λn−1 tn−1,n 
λn
(k)
Die Konvergenz der Subdiagonalelemente ti,i−1 erfüllt für k → ∞
(k) ti,i−1 = O
!
λi k
λi−1 ,
i = 2, . . . , n.
Bem 4.21
• Ist A symmetrisch, dann auch alle T (k) . Folglich konvergiert T (k) gegen
diag (λ1 , . . . , λn ).
• Hat A komplexe Eigenwerte kann der Grenzwert keine obere Dreiecksmatrix T sein, denn T ist reell.
• Auch bei reellen Eigenwerten λi = −λj entsteht im Allgemeinen keine
obere Dreiecksmatrix [QSS02, S. 227, Bsp. 5.9]. Es entsteht aber eine
Block-Dreiecksmatrix.
Bsp 4.22 [QSS02, Bsp. 5.4, S. 217] Wir wenden den QR-Algorithmus auf die
Matrix


4 3 2 1
 3 4 3 2 

A=
 2 3 4 3 
1 2 3 4
an.
c l e a r ; format s h o r t e ; format compact ;
A = [ 4
3
2
1
3
4
3
2
2
3
4
3
1
2
3
4
;
;
;
]
...
...
...
f o r i =1:20
[ Q,R] = qr (A ) ;
A = R∗Q
end
24. März 2016
55
Nach 20 Iterationen ist die

1,1099 e+01
 6,4460 e−10
A=
 1,7436 e−21
2,3185 e−25
Matrix

6,4460 e−10 −8,1781 e−16 2,7788 e−16
3,4142 e+00 1,4286 e−11 2,4560 e−15 

1,4286 e−11 9,0098 e−01 1,1603 e−04 
2,6769 e−15 1,1603 e−04 5,8579 e−01
entstanden. Die tatsächlichen Eigenwerte sind 11,0990195, 3,4142136, 0,9009805
und 0,5857864.
E 4.23 (Aufwandsbetrachtungen)
• Für A = AT benötigt die QR-Faktorisierung mit Householder rund 23 n3
und das Produkt R · Q rund 23 n3 Operationen. Das ist sehr teuer.
• Abhilfe 1: Weniger rechenintensive Iterationen durch Transformation auf
Hessenbergform
T (0) = (Q(0) )T AQ(0) = @
@
@
@
@
@
(0)
(tij = 0 für j ≤ i − 2) bzw. bei symmetrischen Matrizen auf Tridiagonalform
(0)
(0) T
(0)
@
@
T = (Q ) AQ = @
@
@
@ ,
@
@
@
@
siehe Abschnitt 4.3.2. Dann kann die QR-Faktorisierung mit n − 1 GivensRotationen ausgeführt werden. Das ergibt für symmetrische Matrizen rund
10n Operationen pro Iteration.
• Abhilfe 2: Schnellere Konvergenz der Iteration durch Spektralverschiebung,
siehe Abschnitt 4.3.3.
4.3.2. Ähnlichkeitstransformation auf Hessenberg- bzw. Tridiagonalform
E 4.24 (Householder-Reflexion) [Householder 1958] Householder-Matrizen haben die Form
vv T
H =I −2 T .
v v


a1
 .. 
Was bewirken diese? Sei a =  .  ein Vektor. Wählt man
an



v = a + sgn(a1 )kak 

1
0
..
.



,

0
56
24. März 2016
dann ist
v T v = (a1 + sgn(a1 )kak)2 +
n
X
a2i
= 2kak2 + 2sgn(a1 )kak · a1 ,
i=2
v T a = (a1 + sgn(a1 )kak) · a1 +
n
X
a2i
= kak2 + sgn(a1 )kak · a1 ,
i=2
und man erhält

Ha = a − 2

vT a
vv T

a
=
a
−
2
v
=
a
−
v
=

T
T
v v
v v

−sgn(a1 )kak
0
..
.



.

0
Desweiteren können wir zeigen, dass H für v ∈ Rn \ {0} eine orthogonale Matrix
ist: Mit H = H T gilt
vv T
vv T
T
HH = I − 2 T
I −2 T
v v
v v
T
T
T
vv vv
vv
=I −4 T +4 T T =I
v v
v vv v
E 4.25 (Transformation mittels Householder-Spiegelung) Laut Erklärung
4.24 gilt: Für jeden Vektor a gibt es eine orthogonale Matrix Q mit
 
∗
 0 
 
Qa =  .  .
 .. 
0

Folglich gibt es eine orthogonale Matrix Q1 = 
1
0
0
Q̃1

, so dass bei der
Multiplikation mit QT1 von links die erste Zeile erhalten bleibt und das 3. bis
letzte Element der ersten Spalte zu Null transformiert wird,


a11 a12 · · · a1n
 ∗
∗ ···
∗ 



∗ ...
∗ 
QT1 A =  0
.
 ..
..
.. 
 .
.
. 
∗
0
∗
...
Zur Vervollständigung der Ähnlichkeitstransformation multipliziert man nun
mit Q1 von rechts, wobei die erste Spalte erhalten bleibt, so dass


a11 ∗ ∗ · · · ∗
 ∗





QT1 AQ1 =  0
.
 ..

 .

A
2
0
24. März 2016
57
Wir transformieren nun die Restmatrix A2 analog mit



Q2 = 

1
0
0
1
0
Q̃2
0



,

so dass





T
Q2 A = 




a11 a12 · · ·
a21 a22 · · ·
0
∗
..
.
0
..
..
.
.
0
0
···
···

a1n
a2n 













T
und Q2 AQ2 = 




∗
∗
0
..
.
..
.
0
∗ ···
∗ ···
∗
0
..
.
···
···
A3

∗
∗ 








0
Die Fortsetzung dieses Verfahrens ergibt die Hessenbergform.
Bem 4.26 Ist A symmetrisch, dann auch QT AQ. Es werden also auch in den
Zeilen Nullen erzeugt, was zu einer Tridiagonalform führt.
4.3.3. Spektralverschiebung
E 4.27 (Spektralverschiebung) Die Idee ist, im Algorithmus anstelle von A die
Matrix A = A − µI zu verwenden. Es gilt: Ist (λ, v) Eigenpaar von A, dann ist
(λ − µ, v) Eigenpaar von A. Sind die Eigenwerte von A berechnet, kennt man
auch die Eigenwerte von A.
Wozu soll das gut sein? Entsprechend Satz 4.20 spielen für die Konvergenzgeschwindigkeit des QR-Verfahrens die Verhältnisse aufeinander folgender Eigenwerte eine wesentliche Rolle. Durch Spektralverschiebung versucht man diese
Verhältnisse so groß wie möglich zu bekommen.
Alg 4.28 (QR-Algorithmus mit Spektralverschiebung) Iteriere bis zur Konvergenz
• Faktorisiere T (k−1) − µk−1 I = Q(k) R(k)
• Berechne T (k) = R(k) Q(k) + µk−1 I
Folgerung 4.29 Es gilt
T (k) = (Q(k) )T T (k−1) − µk−1 I Q(k) + µk−1 I
= (Q(k) )T T (k−1) Q(k) ,
aber das Q(k) ist anders als bei der Grundform.
58
24. März 2016
E 4.30 (QR-Algorithmus mit Wilkinson-Verschiebung) Man wählt als µk−1
denjenigen Eigenwert von


(k−1)
(k−1)
tn−1,n−1 tn−1,n

,
(k−1)
(k−1)
tn,n−1
tn,n
(k−1)
der am nächsten an tn,n
(k−1)
liegt. Im Fall A = AT konvergiert tn,n−1
• beweisbar quadratisch,
• praktisch kubisch oder schneller
gegen Null und tn,n ist eine ausreichend genaue Approximation für einen Eigenwert λj von T . Die Iterierte hat dann das folgende Aussehen:



T (k) = 

T̃ (k)


0 
0 λj
Man hat einen Eigenwert berechnet und kann den Algorithmus mit der kleineren
Matrix T̃ (k) fortsetzen.
E 4.31 (Aufwand) Der Aufwand für eine symmetrische Matrix A kann wie folgt
geschätzt werden:
• Pro Eigenwert sind im Mittel 17 QR-Schritte nötig (experimentell, siehe
[SK91, S. 111].
• Die Berechnung aller Eigenwerte von T (0) benötigt etwa 10n2 Operationen.
• Der Gesamtaufwand (incl. Tridiagonalisierung) beläuft sich auf etwa 23 n3
Operationen.
• Wenn auch Eigenvektoren mit berechnet werden sollen, erhöht sich der
3
Aufwand auf etwa 11
3 n Operationen.
Bem 4.32
1. Bei Paaren konjugiert-komplexer Eigenwerte einer reellen Matrix gilt:
• Mit einem reellen Shift bleibt die Rechnung reell, konvergiert aber
nicht gegen eine Tridiagonalmatrix.
• Mit einem echt komplexen Shift konvergiert das Verfahren gegen
eine Tridiagonalmatrix, der Preis ist das Rechnen mit komplexer
Arithmetik (höherer Aufwand).
• Bei der Doppel-Verschiebungsstrategie nach Francis verschiebt man
abwechslend mit den Shifts µk und µk . Man kann zwei Iterationsschritte zusammenfassen und mit reeller Artihmetik implementieren,
siehe Literatur.
2. Die effektive Implementierung des QR-Algorithmus ist aufwendig.
24. März 2016
59
E 4.33 (Bewertung) Der QR-Algorithmus mit Tridiagonalisierung und Wilkinson-Spektralverschiebung ist ein gutes Verfahren für Eigenwertprobleme kleiner
und mittlerer Dimension, vor allem wenn alle Eigenwerte berechnet werden
sollen.
Sind nur wenige Eigenwerte von Interesse, gibt es bessere Verfahren.
Wenn A eine dünn besetzte Matrix ist, dann ist i. Allg. sowohl die entsprechende
Hessenberg-Matrix als auch der Faktor Q einer QR-Zerlegung voll besetzt. Der
QR-Algorithmus ist dann nicht empfehlenswert.
4.4. Die Potenzmethode und verwandte Verfahren
4.4.1. Vektoriteration (Potenzmethode, von-Mises-Iteration)
E 4.34 Der sehr einfache Grundalgorithmus besteht in fortgesetzter Multiplikation mit der Matrix A,
yk+1 = Ayk ,
k = 1, 2, 3, . . .
E 4.35 (Idee des Verfahrens) Seien (λj , uj ), j = 1, . . . , n, die Eigenpaare der
Matrix A, wobei
|λ1 |>|λ2 | ≥ · · · ≥ |λn |.
Wir wollen ferner annehmen, dass es keine defektiven Eigenwerte gibt, so dass
die Eigenvektoren uj eine Basis des Rn bilden. Der Startvektor y0 der Iteration
kann in dieser Basis dargestellt werden,
y0 = α1 u1 + α2 u2 + · · · + αn un
Wir benötigen eine weitere Voraussetzung, nämlich α1 6= 0. (Das ist bei einem
zufälligen Startvektor mit Wahrscheinlichkeit 1 erfüllt.) Dann ergibt sich
yk
Ak y 0
α1 λk1 u1 + α2 λk2 u2 + · · · + αn λkn un
=
=
k
k
λ1
λ1
λk1
= α1 u1 + α2
λk2
λkn
k→∞
u
+
·
·
·
+
α
un −→ α1 u1
2
n
λk1
λk1
und damit
yk −→ (λk1 α1 )u1
Die Iteration konvergiert gegen den Eigenvektor zum betragsgrößten
Eigenwert,
und die Konvergenzgeschwindigkeit wird durch den Faktor λλ12 bestimmt.
Wegen
kyk k ≈
α1 λk1
(
→0
falls |λ1 | < 1
→ ∞ falls |λ1 | > 1
muss bei dieser Grundvariante des Algorithmus mit Unter- order Überlauf
gerechnet werden. Abhilfe schafft die Normierung des Vektors in jedem Schritt.
60
24. März 2016
Bem 4.36 Die Matrix A wird nicht explizit benötigt. Es genügt ein Unterprogramm, das die Multiplikation mit einem Vektor ausführt.
E 4.37 (Eigenwerte) Die besten Eigenwertnäherungen ρk erhält man bei einer
symmetrischen (hermiteschen) Matrix A über den Rayleigh-Quotienten von yk ,
ρk =
(yk )T Ayk
,
(yk )T yk
denn dieser minimiert kAyk − ρyk k2 : Es gilt
kAy − ρyk2 = (Ay − ρy)T (Ay − ρy) = y T AT Ay − 2ρy T Ay + ρ2 y T y
∂
!
kAy − ρyk2 = −2y T Ay + 2ρy T y = 0
∂ρ
und folglich ρ = y T Ay/y T y.
Ist die Matrix A nicht symmetrisch, muss der Rayleigh-Quotient anders definiert
werden, wobei auch die Linkseigenvektor(näherung)en benötigt werden.
Alg 4.38 (Potenzmethode)
1: Wähle y0 mit ky0 k = 1
2: for k = 0(1)maxiter do
3:
wk+1 = Ayk
4:
ρk = (yk )T wk+1 {neue Eigenwertnäherung}
5:
yk+1 = wk+1 /kwk+1 k {neue Eigenvektornäherung}
6: end for
Bsp 4.39 [SK91, S. 113] Wir

2

 −1
A=


wenden die Potenzmethode auf die Matrix

−1

..
..

.
.
,
n = 10,

..
..
.
. −1 
−1 2
an, wobei wir zwei verschiedene Startvektoren testen. Die Ergebnisse sind in
den folgenden Tabellen zusammengefasst, wobei yk,1 die erste Komponente der
k-ten Eigenvektornäherung bezeichnet:
k
1
2
3
5
10
20
35
50
65
80
110
yk,1
0.316 227 766
0.248 281 767
0.213 882 337
0.177 816 820
0.142 512 794
0.124 216 323
0.120 405 213
0.120 156 902
0.120 133 207
0.120 131 328
0.120 131 167
y0 =
%k
3.0
3.863 013 70
3.896 138 48
3.907 840 65
3.916 112 28
3.918 887 69
3.918 985 33
3.918 985 94
3.918 985 94
3.918 985 95
3.918 985 95
√1 [1, −1, . . . , 1, −1]T
10
24. März 2016
k
1
2
60
75
90
150
300
400
500
600
700
yk,1
0.316 227 766
0.707 106 781
0.230 530 174
0.230 530 022
0.230 530 018
0.230 529 983
0.230 115 837
0.011 083 507
−0.119 867 700
−0.120 130 644
−0.120 131 165
y0 =
%k
0.2
2.0
3.682 507 06
3.682 507 06
3.682 507 07
3.682 507 07
3.682 509 86
3.859 980 42
3.918 985 64
3.918 985 95
3.918 985 95
√1 [1, 1, . . . , 1, 1]T
10
61
Im ersten Fall hat ρk bereits nach ca. 35 Iterationen 6 genaue Stellen erreicht,
yk erst nach etwa 80 Schritten.
Im zweiten Fall ist der Startvektor orthogonal zu u1 ist. Folglich übernimmt
λ2 die Rolle von λ1 . Das Verfahren hat nach ca. 60 Schritten konvergiert; ρk
approximiert λ2 , yk approximiert u2 .
Wird das Verfahren fortgesetzt, beginnt allmählich der Einfluss von Rundungsfehlern wirksam zu werden. Nach 300 bis 400 Iterationen kippt“ das Verfahren
”
und konvergiert schließlich gegen λ1 und u1 .
4.4.2. Teilraumiteration und Rayleigh-Ritz-Prozess
E 4.40 (wesentliche Merkmale der Teilraumiteration)
• Die Teilraumiteration ist die natürliche Verallgemeinerung der Vektoriteration wobei p Vektoren simultan iteriert werden. Die p (1 ≤ p ≤ n)
orthonormal en Startvektoren werden in eine Matrix V (0) ∈ Rn×p zusammengefasst.
• Die p Vektoren müssen aber nicht nur normiert, sondern auch orthogonalisiert werden, da sie sonst alle gegen u1 konvergieren.
Alg 4.41 (Teilraumiteration)
1: Wähle p ∈ {1, . . . , n} und V0 ∈ Rn×p
2: for k = 0(1)maxiter do
3:
W (k+1) = AV (k)
4:
W (k+1) = Q(k+1) R(k+1) {Orthogonalisierung}
5:
V (k+1) = Q(k+1)
6: end for
Dabei ist Q(k+1) ∈ Rn×p und R(k+1) ∈ Rp×p .
E 4.42 Ist |λ1 | > |λ2 | > · · · > |λp | > |λp+1 | ≥ · · · ≥ |λn |, dann konvergieren die
Spalten von V (k) gegen die Eigenvektoren zu λ1 , . . . , λp .
Gilt nur |λ1 | ≥ |λ2 | ≥ · · · ≥ |λp | > |λp+1 | ≥ · · · ≥ |λn |, so konvergiert das
Verfahren in dem Sinne, dass der durch die Spalten von V (k) aufgespannte Teilraum gegen den Teilraum geht, der durch die Eigenvektoren zu den dominanten
Eigenwerten aufgespannt wird.
E 4.43 Zur Berechnung der Eigenwerte wird der Rayleigh-Quotient verallgemeinert. Gute Eigenwertnäherungen sind die Eigenwerte der Matrix
P (k+1) = (Q(k+1) )T AQ(k+1) ∈ Rp×p
(4.26)
Da p klein ist, können diese zum Beispiel mit dem QR-Verfahren berechnet
werden. Wir erhalten die Matrizen X (k+1) , Λ̃(k+1) ∈ Rp×p , und es gilt
P (k+1) = (X (k+1) )T Λ̃(k+1) X (k+1) .
62
24. März 2016
E 4.44 Bessere Näherungen für die Eigenwerte ergeben sich durch die Modifikation
V (k+1) = Q(k+1) X (k+1) ,
so dass das folgende Verfahren entsteht. Dabei berechnet man die auch für
P (k+1) benötigte Hilfsgröße Z (k+1) := AQ(k+1) und kann statt mit V (k+1) =
Q(k+1) X (k+1) die Berechnung von W (k+2) über
W (k+2) = AV (k+1) = AQ(k+1) X (k+1) = Z (k+1) X (k+1)
mit Z (k) := AQ(k)
durchführen, wodurch pro Iteration eine Matrix-Multiplikation gespart werden
kann.
Alg 4.45 (Teilraumiteration)
1: Wähle V (0) ∈ Rn×p und initialisiere Z (0) = AV (0) , X (0) = I.
2: for k = 0(1)maxiter do
3:
W (k+1) = Z (k) X (k)
4:
W (k+1) = Q(k+1) R(k+1)
5:
Z (k+1) = AQ(k+1)
6:
P (k+1) = (Q(k+1) )T Z (k)
7:
P (k+1) = (X (k+1) )T Λ̃(k+1) X (k+1)
8: end for
Λ̃(k+1) enthält die Eigenwertnäherungen; die Eigenvektornäherungen erhält man
durch V (k+1) = Q(k+1) X (k+1) im Postprocessing.
4.4.3. Inverse Iteration und Rayleigh-Quotienten-Iteration
E 4.46 (Inverse Iteration) Seien λ1 , . . . , λn die Eigenwerte der Matrix A und
u1 , . . . , un die zugehörigen Eigenvektoren. Weil
Auj = λj uj
⇔
−1
λ−1
j uj = A uj
ist, hat die Matrix A−1 die Eigenwerte λ−1
j und dieselben Eigenvektoren uj . Der
betragskleinste Eigenwert von A kann durch Vektoriteration mit A−1 ermittelt
werden, sofern |λ1 | < |λ2 | ≤ · · · ≤ |λn |.
Das bedeutet, dass die Anweisung wk+1 = Ayk durch wk+1 = A−1 yk zu ersetzen
ist, d. h., wk+1 ist die Lösung des Gleichungssystems
Awk+1 = yk .
Das so modifizierte Verfahren heißt Inverse Iteration .
E 4.47 (Inverse Teilraum-Iteration) Analog ist bei der inversen TeilraumIteration die neue Matrix W (k+1) die Lösung des Gleichungssystems
AW (k+1) = V (k) .
24. März 2016
63
Bem 4.48 Man beachte, dass bei inverser (Teilraum-)Iteration die Matrix des
zu lösenden Gleichungssystems in jeder Iteration diesselbe ist. Es kann also aus
Varianten einer LU-Faktorisierung Nutzen gezogen werden.
E 4.49 (Spektralverschiebung) Für die Spektralverschiebung gilt: Ã = A − µI
hat die Eigenwerte λ̃j = λj − µ, j = 1, . . . , n. Ist µ eine gute Näherung für λj ,
dann ist λ̃j betragskleinster Eigenwert von Ã. Der verschobene Eigenwert λ̃j
kann dann durch inverse Iteration mit Ã bestimmt werden.
˜
˜
E 4.50 (Rayleigh-Quotienten-Iteration) Seien λ̃1 und λ̃2 die beiden betragskleinsten Eigenwerte von Ã. Die Konvergenzgeschwindigkeit der inversen Itera˜ tion hängt von λ̃1 ab. D.h., je besser µ den Eigenwert λj approximiert, desto
˜
λ̃2
˜ kleiner λ̃˜ 1 , desto besser die Konvergenz.
λ̃2
Daraus ergibt sich folgende Idee: Wir arbeiten nicht mit fester Verschiebung µ,
sondern mit dem Rayleigh-Quotienten ρk als aktuell bester Eigenwertnäherung.
Diese Variante wird als Rayleigh-Quotienten-Iteration bezeichnet:
1:
2:
3:
4:
5:
6:
Wähle y1 mit ky1 k = 1
for k=1(1)maxiter do
ρk = ykT Ayk {siehe Erklärung 4.37}
Berechne wk+1 aus (A − ρk I)wk+1 = yk
yk = wk+1 /kwk+1 k
end for
Das Verfahren hat folgende Eigenschaften:
• Die Systemmatrix ändert sich bei jeder Iteration.
• Das Verfahren konvergiert sogar kubisch.
• Es ist nicht klar, welches Eigenpaar approximiert wird.
Besser ist es, zunächst einige Schritte inverse Iteration mit fester Spektralverschiebung auszuführen.
4.5. Das Lanczos-Verfahren
E 4.51 (wesentliche Merkmale des Lanczos-Verfahrens)
• Das Verfahren setzt eine symmetrische (hermitesche) Matrix A voraus.
Die Erweiterung auf den nichtsymmetrischen Fall heißt Arnoldi-Verfahren.
• Das Lanczos-Verfahren ist zunächst ein Verfahren zur Tridiagonalisierung,
A = QT QT , wobei im k-ten Iterationsschritt die Matrix


α1 β1


 β1 α2 . . .

(k)

 ∈ Rk×k .
T =

.
.
.
.

.
. βk−1 
βk−1 αk
entsteht und T = T (n) ist.
64
24. März 2016
• Die Bedeutung des Verfahrens liegt darin, dass die extremalen (d.h. die
kleinsten und größten, die am Rand des Spektrums) Eigenwerte von A
bereits für kleine k überraschend gut durch die extremalen Eigenwerte
von T (k) approximiert werden. Da außerdem zur Berechnung der T (k) nur
Matrixmultiplikationen benötigt werden, ist das Verfahren vor allem für
große, schwach besetzte Matrizen gut geeignet.
E 4.52 (Herleitung des Algorithmus)
• Wir schreiben die Zerlegung in der Form AQ = QT mit

α1 β1

..
 β1 α2
.
Q = q1 q2 · · · qn ,
T =

..
..

.
. βn−1
βn−1 αn



.


Diese Zerlegung bedeutet spaltenweise:
Aqk = βk−1 qk−1 + αk qk + βk qk+1 ,
k = 1, . . . , n
(4.27)
wobei β0 = βn = 0 und q0 = qn+1 = 0 gesetzt wird, damit (4.27) auch für
k = 1 und k = n gilt.
• Anschließend multiplizieren wir von links mit qkT . Das liefert wegen der
vorausgesetzten Orthonormalität der Vektoren qk die Beziehung
qkT Aqk = αk .
• Wir lösen dann (4.27) nach βk qk+1 auf:
βk qk+1 = Aqk − βk−1 qk−1 − αk qk =: rk .
Mit dieser Definition von rk gilt
krk k = kβk qk+1 k = |βk |.
• Falls rk 6= 0 ist βk = ±krk k =
6 0 und damit
qk+1 =
1
rk .
βk
Beim Lanczos-Verfahren entscheidet man sich für βk = +krk k.
Damit ergibt sich folgender Algorithmus.
24. März 2016
65
Alg 4.53
1: Setze q0 = 0, β0 = 0
2: Wähle q1 mit kq1 k = 1
3: for k = 1(1)n do
4:
pk = Aqk
5:
αk = qkT pk
6:
rk = pk − αk qk − βk−1 qk−1
7:
if k = n then stop
8:
βk = krk k
9:
if βk = 0 then stop
10:
qk+1 = β1k rk
11: end for
E 4.54 Zweckmäßig ist die Kopplung mit dem Rayleigh-Ritz-Algorithmus. Man
beachte, dass T (k) = (Q(k) )T AQ(k) gerade die Matrix P aus (4.26) ist. Die Matrix
T (k) kann nach Zeile 5 des Algorithmus aus T (k−1) , αk und βk−1 aufgebaut
werden. Mit dem QR-Verfahren können die Eigenwerte von T (k) und damit
Eigenwertnäherungen von A berechnet werden.
Bem 4.55 In Computerarithmetik können die Spalten qj von Q stark von der
Orthogonalität abweichen, und das Verfahren kann zusammenbrechen. Für
Stabilisierungen siehe Spezialliteratur.
4.6. Symmetrische verallgemeinerte Eigenwertprobleme
Def 4.56 (verallgemeinertes Eigenwertproblem) Ein vom Nullvektor verschiedener Vektor x heißt Eigenvektor des geordneten Paares (A, B), wenn es eine
(reelle oder komplexe) Zahl λ gibt, so dass
Ax = λBx.
Die Zahl λ heißt Eigenwert.
Ü 4.57 Man überlege sich, dass die Eigenwerte des verallgemeinerten Eigenwertproblems gerade die Nullstellen des charakteristischen Polynoms
p(A,B) (λ) := det(λB − A)
sind.
Satz 4.58 Seien A, B ∈ Rn×n symmetrische Matrizen und sei B zusätzlich positiv definit. Dann hat das Paar (A, B) nur reelle Eigenwerte und die zugehörigen
Eigenvektoren sind B-orthogonal,
xTi Bxj = 0 für i 6= j.
Beweis Weil die Matrix B symmetrisch und positiv definit ist, existiert die
Cholesky-Faktorisierung, B = LLT . Wir können nun das verallgemeinerte
66
24. März 2016
Eigenwertproblem Axj = λj Bxj in ein spezielles Eigenwertproblem für eine
symmetrische Matrix umformen:
Axj = λj LLT xj
⇔
L−1 Axj = λj LT xj
⇔
L−1 AL−T LT xj = λj LT xj
⇔
Ãx̃j = λj x̃j
mit Ã = L−1 AL−T und x̃j = LT xj .
Wegen
ÃT = (L−1 AL−T )T = L−1 AT L−T = L−1 AL−T = Ã
ist die Matrix Ã symmetrisch. Daraus folgt, dass die Eigenwerte λj reell sind
und die zugehörigen Eigenvektoren x̃j orthogonal zueinander, das heißt, es gilt
für i 6= j
0 = x̃Ti x̃j = (LT xi )T LT xj = xTi LLT xj = xTi Bxj .
Damit ist alles gezeigt.
q.e.d.
Bem 4.59 (Diskussion der Voraussetzungen von Satz 4.58)
1. In den Anwendungen in der Dynamik ist A häufig eine Steifigkeitsmatrix
und B eine Massematrix. Diese erfüllen die Voraussetzungen des Satzes
4.58, es sind sogar beide symmetrisch und positiv definit.
2. Die positive Definitheit von B ist wesentlich für die Aussage des Satzes
4.58. Parlett hat gezeigt, dass es zu jeder beliebigen Matrix C ∈ Rn×n
symmetrische Matrizen A, B ∈ Rn×n gibt, so dass C = AB −1 gilt. Damit
gilt
Cx = λx
⇔
⇔
AB
−1
x = λx = λBB −1 x
Ax̂ = λB x̂
mit x̂ = B −1 x.
Die Matrizen A und B sind symmetrisch, aber die Eigenwerte können
komplex sein.
3. Unter den Voraussetzungen des Satzes 4.58 hat das verallgemeinerte Eigenwertproblem genau n Eigenwerte. Wenn die Matrix B nicht invertierbar
ist, gilt dies nicht mehr. So ist für
1
0
1 0
A=
, B=
0 −1
0 0
das charakteristische Polynom p(A,B) (λ) = λ − 1, und λ = 1 ist der einzige
Eigenwert.
Für
A=B=
1 0
0 0
ist p(A,B) (λ) ≡ 0, und alle λ ∈ C sind Eigenwerte.
24. März 2016
67
E 4.60 Wenn B invertierbar ist, dann ist das verallgemeinerte Eigenwertproblem
Ax = λBx äquivalent zum speziellen Eigenwertproblem B −1 Ax = λx und auch
zum speziellen Eigenwertproblem AB −1 x̃ = λx̃ mit x̃ = Bx. Es gibt aber
Gründe, diese Umformung nicht auszuführen:
• Wenn A und B symmetrisch sind, dann ist AB −1 bzw. B −1 A nicht symmetrisch.
• Wenn A und B dünn besetzt sind, dann ist AB −1 bzw. B −1 A im Allgemeinen nicht dünn besetzt.
• Wenn B schlecht konditioniert ist, aber die Eigenwerte von Ax = λBx
gut konditioniert sind, dann sind die Eigenwerte von AB −1 bzw. B −1 A
schlecht konditioniert [Wat10].
Zum Glück ist die Massematrix bei einer Finite-Elemente-Diskretisierung gut
konditioniert, so dass die dritte Eigenschaft nicht ins Gewicht fällt. Die erste
Eigenschaft kann man umgehen, indem man die Cholesky-Faktorisierung B =
LLT benutzt und wie im Beweis von Satz 4.58 zum speziellen Eigenwertproblem
Ãx̃ = λx̃
mit Ã = L−1 AL−T und x̃ = LT x.
übergeht. Der Nachteil des Auffüllens des Bandes (“Fill-In”) bei den CholeskyFaktoren bleibt aber bestehen.
Um die Matrix Ã nicht berechnen zu müssen, arbeitet man mit den einzelnen Faktoren, was besonders gut bei Verfahren wie der Potenzmethode oder
dem Lanczos-Verfahren geht, bei denen in jeder Iteration nur mit der Matrix
multipliziert wird.
4.7. Zusammenfassung
E 4.61 (Wahl des Algorithmus)
• Wenn die Matrix voll besetzt ist und kleine bis mittlere Dimension hat:
– Sind mehr als ein Viertel aller Eigenwerte zu berechnen, verwendet
man eine Variante des QR-Verfahrens.
– Sind nur wenige Eigenwerte zu berechnen, kann inverse Iteration
empfohlen werden.
• Wenn die Matrix schwach besetzt ist und hohe Dimension hat, greift man
zu Iterationsverfahren.
Bem 4.62 Eigenwertprobleme für nichtsymmetrische Matrizen unterscheiden
sich in verschiedenen Punkten vom symmetrischen Eigenwertprobem:
•
•
•
•
68
Statt einer Tridiagonalisierung entsteht die Hessenbergform.
Beim Lanczos-Verfahren entsteht keine 3-stufige Rekursion mehr.
Es können komplexe Eigenwerte und Eigenvektoren auftreten.
Es können defektive Eigenwerte (Eigenwerte, bei denen die geometrische Vielfachheit kleiner als die algebraische Vielfachheit ist, also keine
vollständige Basis aus Eigenvektoren existiert) auftreten. Mir ist kein
Algorithmus bekannt, der damit gut zurechtkommt.
24. März 2016
A. Grundlagen Eigenwertprobleme
Def A.1 Ein vom Nullvektor verschiedener Vektor x heißt Eigenvektor der
Matrix A, wenn es eine (reelle oder komplexe) Zahl λ gibt, so dass
Ax = λx.
Die Zahl λ heißt der zum Eigenvektor x gehörende Eigenwert von A, das Paar
(λ, x) heißt Eigenpaar von A.
E A.2 Nehmen wir an, ein Eigenwert λ der Matrix A wäre bekannt. Wegen
Ax = λIx
⇐⇒
(A − λI)x = 0
ist der zugehörige Eigenvektor x Lösung des homogenen Gleichungssystems mit
der Systemmatrix A − λI. Ist die Determinante det(A − λI) 6= 0, dann hat
das LGS nur die triviale Lösung x = 0. Eigenvektoren sind jedoch nichttriviale
Lösungen, die genau dann existieren, wenn die Beziehung
pA (λ) := det(A − λI) = 0
gilt. Ist A eine n × n-Matrix, dann ist pA (λ) ein Polynom n-ten Grades in λ.
Def A.3 Sei A eine n × n-Matrix. Das Polynom n-ten Grades
pA (λ) := det(A − λI)
heißt charakteristisches Polynom von A. Die Gleichung
pA (λ) = 0
heißt charakteristische Gleichung.
Satz A.4 Die Eigenwerte von A sind genau die Nullstellen des charakteristischen
Polynoms bzw. die Lösungen der charakteristischen Gleichung. Zählt man jede
Nullstelle mit ihrer Vielfachheit, so besitzt pA (λ) genau n Nullstellen (reelle und
konjugiert komplexe Paare) und damit besitzt A genau n Eigenwerte.
E A.5 (Berechnung von Eigenpaaren aus der Definition) Eine mögliche Vorgehensweise zur Lösung des speziellen Eigenwertproblems besteht in der Ausführung folgender Schritte:
Schritt 1 Bestimmen der Eigenwerte durch Lösen der charakteristischen Gleichung pA (λ) = 0,
Schritt 2 Bestimmen der zugehörigen Eigenvektoren durch Lösen des Gleichungssystems (A − λI)x = 0 der Reihe nach für alle Eigenwerte λ.
Dieser Algorithmus mag für Matrizen bis zur Ordnung 4 brauchbar sein. Für
größere Matrizen verwendet man jedoch numerische Verfahren.
24. März 2016
69
Satz A.6 (Eigenwerte und Eigenvektoren spezieller Matrizen)
1. Ist A eine Diagonalmatrix oder eine obere bzw. untere Dreiecksmatrix,




a11
0
a11
0

 a21 a22




.
..
A=
oder
A= .


.
.
.
.
.

 .
.
.
0
ann
an1 an2 · · · ann
dann sind die Eigenwerte die Elemente in der Hauptdiagonalen, λi = aii ,
i = 1, . . . , n.
2. Ist (λ, x) ein Eigenpaar einer invertierbaren Matrix A, dann ist ( λ1 , x) ein
Eigenpaar von A−1 .
3. Ist (λ, x) ein Eigenpaar einer Matrix A, dann ist (λk , x) ein Eigenpaar von
Ak , k ∈ N.
4. AT hat dieselben Eigenwerte wie A, aber im Allgemeinen andere Eigenvektoren.
Beweis
1. In den genannten Fällen gilt
pA (λ) = (a11 − λ)(a22 − λ) · · · (ann − λ),
woraus sofort die Behauptung folgt.
2. Wir multiplizieren Ax = λx mit A−1 von links und erhalten
x = A−1 (λx) = λA−1 x
und damit die Behauptung
A−1 x =
1
x.
λ
3. Es gilt
Ax = λx
A2 x = A · Ax = A · λx = λAx = λ2 x,
womit wir den Fall k = 2 bewiesen hätten. Für allgemeines k ∈ N verwendet
man vollständige Induktion.
4. Wegen I = I T gilt det(AT −λI) = det(A−λI)T . Die Behauptung folgt, weil
die Determinante einer Matrix gleich der Determinante der Transponierten
ist.
q.e.d.
Ü A.7 Man überlege sich, dass die dritte Eigenschaft sogar für alle k ∈ Z gilt.
Wie müsste A0 definiert sein, damit diese Eigenschaft für k = 0 gilt?
Ü A.8 Beweisen Sie: Ist (λ, x) ein Eigenpaar einer Matrix A, dann ist (αλ, x)
ein Eigenpaar von αA, α ∈ C.
70
24. März 2016
Entsprechend der Definition gehört zu jedem Eigenwert mindestens ein Eigenvektor. Wir wollen nun die Menge aller Eigenvektoren zu einem Eigenwert
charakterisieren.
Satz A.9 Sind x und y Eigenvektoren zum Eigenwert λ, dann auch jede Linearkombination αx + βy, wobei α und β beliebige komplexe Zahlen sind. Damit
spannen die Eigenvektoren zu einem Eigenwert λ einen linearen Raum auf, den
Eigenunterraum.
Beweis Der Beweis erfolgt durch Ausnutzen der Rechenregeln für Matrizen:
A(αx + βy) = αAx + βAy = αλx + βλy = λ(αx + βy).
q.e.d.
Def A.10 (Algebraische und geometrische Vielfachheit) Ist die Zahl λ eine kfache Nullstelle des charakteristischen Polynoms, dann sagt man, der Eigenwert
λ hat die algebraische Vielfachheit k.
Die maximale Anzahl linear unabhängiger Eigenvektoren zum Eigenwert λ (die
Dimension des Eigenunterraums) heißt geometrische Vielfachheit von λ.
Satz A.11 Die geometrische Vielfachheit ist stets kleiner oder gleich der algebraischen Vielfachheit.
Bsp A.12 (ungleiche geometrische und algebraische Vielfachheit) Die Matrix
A=
0 1
0 0
besitzt den Eigenwert λ = 0 mit der algebraischen Vielfachheit 2. Wir berechnen
nun die zugehörigen Eigenvektoren. Es gilt
(A − λI)
x1
x2
=
0 1
0 0
x1
x2
=
x2
0
!
=
0
0
.
Folglich ist x2 = 0 und x1 beliebig, d. h., es gibt nur den eindimensionalen
Eigenunterraum
Span
1
0
:=
α
0
, α ∈ C beliebig .
Def A.13 Die Spur Sp(A) einer Matrix A ist als die Summe ihrer Hauptdiagonalelemente definiert,
Sp(A) :=
n
X
aii = a11 + a22 + . . . + ann .
i=1
24. März 2016
71
Satz A.14 (Eigenschaften von Eigenwerten und Eigenvektoren)
1. Die Summe aller Eigenwerte einer Matrix A ist gleich der Spur von A,
n
X
λi = λ1 + λ2 + λ3 + . . . + λn = Sp(A).
i=1
2. Das Produkt aller Eigenwerte ist gleich der Determinante von A,
n
Y
λi = λ1 · λ2 · . . . · λn = det A.
i=1
3. Die zu verschiedenen Eigenwerten λi von A gehörenden Eigenvektoren xi
sind linear unabhängig.
Beweis Als Beweisskizze für die ersten beiden Eigenschaften sei angegeben,
dass
pA (λ) = det(A − λI)
= (−λ)n + (a11 + a22 + . . . + ann )(−λ)n−1 + . . . + det A
und andererseits aber auch
pA (λ) = (λ1 − λ)(λ2 − λ) · · · (λn − λ)
= (−λ)n + (λ1 + λ2 + λ3 + . . . + λn )(−λ)n−1 + . . . + λ1 λ2 · · · λn
gilt. Durch Koeffizientenvergleich erhält man die Behauptung. Einen Beweis der
dritten Eigenschaft findet man z. B. in [MV03, Kap. 6, Satz 6.7].
q.e.d.
Folgerung A.15 Sind alle Eigenwerte von A voneinander verschieden, so gehört
zu jedem Eigenwert (bis auf eine multiplikative Konstante) genau ein Eigenvektor.
Die n Eigenvektoren sind linear unabhängig und spannen den Rn auf.
Allgemein gilt: Stimmen algebraische und geometrische Vielfachheit aller Eigenwerte einer Matrix überein, dann gibt es n linear unabhängige Eigenvektoren,
die den Raum Rn aufspannen.
Def A.16 (Ähnlichkeitstransformation) Sei A eine gegebene Matrix. Die mit
einer beliebigen regulären (invertierbaren) Matrix Y gebildete Transformation
B = Y −1 AY
(A.28)
heißt Ähnlichkeitstransformation. Eine Matrix B heißt ähnlich zu A, wenn es
eine reguläre Matrix Y gibt, so dass die Gleichung (A.28) gilt.
Satz A.17 Ähnliche Matrizen haben dieselben Eigenwerte. (Die Eigenwerte sind
invariant unter einer Ähnlichkeitstransformation.)
72
24. März 2016
Beweis Ist λ Eigenwert von B und x der zugehörige Eigenvektor, dann ist
Bx = λx
Y
x = λx
A(Y x) = λ(Y x)
−1 AY
d. h., λ ist ein Eigenwert von A mit dem Eigenvektor Y x.
q.e.d.
Satz A.18 (Schur) Jede quadratische Matrix ist ähnlich zu einer oberen Dreiecksmatrix.
E A.19 Die Bedeutung dieses Satzes liegt darin, dass es numerische Verfahren
gibt (den QR-Algorithmus, siehe Abschnitt 4.3), die diese obere Dreiecksmatrix
(näherungsweise) bestimmt. Die Eigenwerte stehen dann (näherungsweise) auf
der Hauptdiagonale dieser Dreiecksmatrix.
Wenn algebraische und geometrische Vielfachheit aller Eigenwerte übereinstimmen,
kann man das bessere Resultat aus den nächsten Satz zeigen.
Satz A.20 (Diagonalisierung einer Matrix) Stimmen algebraische und geometrische Vielfachheit aller Eigenwerte einer Matrix überein, dann ist A ähnlich
zur Diagonalmatrix


λ1


..
Λ := diag (λ1 , . . . , λn ) = 
.
.
λn
Dabei sind die Spalten der Transformationsmatrix die Eigenvektoren von A: Sei
X := x1 . . . xn
mit
Axi = λi xi .
Dann ist
AX = XΛ,
d.h.
Λ = X −1 AX.
Beweis Es bleibt nachzurechnen, dass AX = XΛ. Dies ist eine Matrixgleichung
und man überprüft leicht die Gleichheit der i-ten Spalte, i = 1, . . . , n. Die i-te
Spalte von AX ist Axi = λi xi und die i-te Spalte von XΛ ist ebenfalls xi λi .
q.e.d.
E A.21 Die Diagonalisierung eignet sich z. B. zum Berechnen von Matrixfunktionen. Ist die Zerlegung A = XΛX −1 mit Λ = diag (λ1 , . . . , λn ) bekannt, dann
kann man Ak über
Ak = XΛX −1 · XΛX −1 · · · XΛX −1 = XΛk X −1
berechnen, wobei Λk = diag (λk1 , . . . , λkn ) wieder eine Diagonalmatrix ist.
Eine andere Anwendung der Diagonalisierung ist die Hauptachsentransformation.
24. März 2016
73

Download Report