Grundlagen der Numerischen Mathematik · InhaltsverzeichnisII Pseudo-Kode mit Skalierter Pivotsuche...

Grundlagen der Numerischen Mathematik

a.o.Univ.Prof. Mag.Dr. Stephen Keelinghttp://math.uni-graz.at/keeling/

Literatur:Numerical Analysis von R.L. Burden und J.D. Faires

Unterlagen:http://math.uni-graz.at/keeling/teaching.html

1

Inhaltsverzeichnis IEinfuhrung

Motivierendes Beispiel: BildverarbeitungApproximation von Ableitungen und IntegralenMittelwertsatzOptimalitatsbedingung furs BildverarbeitungsbeispielRandwertproblem furs BildverarbeitungsbeispielFehleranalyseSatz von TaylorO-NotationZwischenwertsatzIterative Bestimmung einer NullstelleKonvergenzgeschwindigkeitStabilitatDarstellung von Zahlen im ComputerAufrundung und AbrundungSignifikante ZiffernFehler der FließkommadarstellungAusloschungHorner Algorithmus

Lineare Gleichungssysteme: Direkte MethodenGaußsche EliminationRuckwarts SubstitutionDiskretisierung einer PDGAutomatisierung der Gaußschen EliminationPseudo-Kode der Gaußschen Eliminationflops der Gaußschen EliminationRekursive Losungen: EvolutionsgleichungLU ZerlegungRekursive Losungen mit LUPivot StrategienPseudo-Kode fur Gaußsche Elimination mit Pivot-SucheLU Zerlegung mit PermutationSkalierte Pivotsuche2

Inhaltsverzeichnis IIPseudo-Kode mit Skalierter PivotsucheTotale PivotsucheInverse und DeterminanteBesondere Eigenschaften einer MatrixGerschgorin SatzPseudo-Kode fur den Cholesky AlgorithmusBandmatrizenSparse Speicherformat

Lineare Gleichungssysteme: Iterative MethodenVektornormenAquivalenz von NormenQualitative Eigenschaften der NormenMatrixnormenCharakterisierung bekannter MatrixnormenSpektralradiusKonvergente MatrizenZerlegungen fur IterationenJacobi MethodePseudo-Kode der Jacobi MethodeGauß-Seidel MethodeSymmetrische Gauß-Seidel MethodeApproximierte InversenSOR MethodeKonvergenz von Iterativen VerfahrenKonjugierte GradientenPrakonditionierungFehler Abschatzungen

Eigenwerte und EigenvektorenVektoriterationInverse VektoriterationBerechnung aller EigenwerteDie Householder MethodeHessenberg Matrizen3

Inhaltsverzeichnis IIIPseudo-Kode zur Transformation auf Hessenberg oder Tridiagonale FormQR AlgorithmusGivens TransformationenKonvergenz des QR-AlgorithmusPseudo-Kode fur den QR-Algorithmus

AusgleichsproblemeLineare RegressionPolynome RegressionSingularwert ZerlegungEigenschaften der SWZPseudoinverseBerechnung der SWZBerechnung der SWZVereinfachung auf Bidiagonale FormPseudo-Kode zur Transformation auf Bidiagonale FormDiagonalisierende IterationFolge der Bidiagonalen MatrizenPseudo-Kode fur den QR Algorithmus zur Singularwert-Zerlegung

InterpolationLagrange Interpolierende PolynomeGenauigkeit globaler InterpolationIterierte InterpolationNevilles AlgorithmusDividierte DifferenzenHermite InterpolationStuckweise Polynome InterpolationStuckweise Kubisch Hermite InterpolationSplinesDie Kanonischen SplinesGlattheitsbedingungenInterpolation mit Kubischen SplinesTensor Produkte fur 2D Interpolation

Nicht Lineare Gleichungen4

Inhaltsverzeichnis IVPseudo-Kode fur das BisektionsverfahrenFixpunktiterationPseudo-Kode fur eine FixpunktiterationVergleich zweichen Bisektions- und FixpunktiterationNewton VerfahrenPseudo-Kode fur das Newton VerfahrenPseudo-Kode fur Sekant VerfahrenVergleich zweichen Newton und BisektionsiterationAsymptotische KonvergenzrateSysteme von Nicht Linearen GleichungenAbstiegsverfahrenImplizite FixpunktiterationNewton Verfahren fur SystemeRichtungsfelder der LosungsverfahrenEntrauschen Ergebnisse

Numerisches Differenzieren und IntegrierenApproximation der AbleitungenRichardson ExtrapolationNumerische IntegrationTrapez-RegelSimpson-RegelGenauigkeit einer geschlossenen Newton-Cotes FormelGenauigkeit einer offenen Newton-Cotes FormelMittelpunkt-RegelZusammengesetzte Mittelpunkt-RegelZusammengesetzte Simpsons-RegelnZusammengesetzte Trapez-RegelnRomberg IntegrationGauß Quadratur

Gewohnliche DifferentialgleichungenDas Eulersche VerfahrenWarmegleichungKonservative Verfahren5

Inhaltsverzeichnis VWellengleichungCrank-Nicholson VerfahrenRunge-Kutta MethodenRunge-Kutta-Fehlberg MethodenKonsequenz, Stabilitat, KonvergenzAbsolute StabilitatSteife SystemeRandwertproblemeDiskretisierung eines Sturm-Liouville RandwertproblemsRayleigh-Ritz Methode

6

Was ist Numerische Mathematik?

I Experiment zu Hause:I Im Taschenrechner: 2 eingeben, dann 10×

√· und 10×x2.

I Subtrahiere 2. Ergebnis 6= 0. Warum?I Zahlen werden nicht genau gespeichert!

I Wie soll man Ax = b, f (x) = 0, u′′(x) = f (x) losen?Antwort: Naherungsweise.

I Ein gutes numerisches Verfahren soll nicht erlauben, dassFehler sich schlecht aufbauen!

7

Motivierendes Thema: Signal- und Bildverarbeitung

I Dieses Beispiel wird mehrmals in verschiedenenKontexten analysiert.

I v(x) ist ein gegebenes verrauschtes Signal.I u?(x) ist das gewunschte Signal.I u(x) ist eine Abschatzung von u?(x).

8


I Wie bestimmt man u ≈ u??I Ansatz: Ω = (0,1), u? ≈ u = argminwJ(w),

J(u) =

∫Ω

(u − v)2dx + µ

∫Ω|u′|2dx

I Qualitativ:I µ→∞⇒ u = Mittelwert von v (ubergeglattet)I µ→ 0⇒ u = v (untergeglattet)

Das beste µ? ist irgendwo in der Mitte.I Optimalitatsbedingung zur Bestimmung von u?

Diskreter Ansatz...

9

Motivierendes Thema: Signal- und BildverarbeitungI Die Daten werden ublicherweise sowieso diskret gegeben:

vi = v(xi) auf einem Gitter

xi = ih, i = 0, . . . ,N, h = 1/N

x0

0

xN

1

xi

h

I Fragestellung: Es gelten

hN∑

i=0

[u(xi)− v(xi)]2h→0−→

∫Ω|u(x)− v(x)|2dx

µhN∑

i=1

[u(xi)− u(xi−1)

h

]2h→0−→ µ

∫Ω|u′(x)|2dx

unter welchen Bedingungen?10

Motivierendes Thema: Signal- und BildverarbeitungSatz (Mittelwertsatz): Wenn f ∈ C[a,b] differenzierbar in (a,b)ist, dann ∃c ∈ (a,b) sodass f ′(c) = [f (b)− f (a)]/[b − a].I Also ∃ξi ∈ (xi−1, xi) sodass

u(xi)− u(xi−1)

h=

u(xi)− u(xi−1)

xi − xi−1= u′(ξi)

und

µhN∑

i=1

[u(xi)− u(xi−1)

h

]2

= µ

N∑i=1

[u′(ξi)]2h ≈ µ∫

Ω|u′(x)|2dx

Satz (Mittelwertsatz fur Integrale): Wenn f ∈ C[a,b] undg(x) ≥ 0, ∀x ∈ [a,b], dann ∃c ∈ (a,b) sodass∫ b

af (x)g(x)dx = f (c)

∫ b

ag(x)dx

I Themen: Approximation von Ableitungen und Integralen.11

Motivierendes Thema: Signal- und BildverarbeitungI Fur Daten vi und Approximationen ui ≈ u(xi) seien

u = 〈u0, . . . ,uN〉T und v = 〈v0, . . . , vN〉T. Wir minimieren:

Jh(u) = hN∑

i=0

(ui − vi)2 + µh

N∑i=1

[ui − ui−1

h

]2

I Optimalitat: ∇Jh(u) = ∇uJh(u) = 0.∂Jh

∂uk︸︷︷︸0<k<N

= 2h(uk − vk ) + µh/h2[2 (uk − uk−1)︸︷︷︸i=k

−2 (uk+1 − uk )︸︷︷︸i=k+1

]

∂Jh

∂u0= 2h(u0 − v0) + µ/h[ −2 (u1 − u0)︸︷︷︸

i=k+1=1

]

∂Jh

∂uN= 2h(uN − vN) + µ/h[2 (uN − uN−1)︸︷︷︸

i=k=N

]

12


I oder ∇Jh(u) = 2h(u − v) + 2(µ/h)Du wobei

D =

1 −1 0−1 2 −1

. . .. . .

. . .

−1 2 −10 0 −1 1

I Optimalitat:

0 =1

2h∇Jh(u) = u − v +

µ

h2 Du

oder das lineare Gleichungssystem:[ µh2 D + I

]u = v

13


I Optimalitatssystem komponentenweise:

−µui+1 − 2ui + ui−1

h2 + ui = vi , i = 1, . . . ,N − 1

−µu1 − u0

h2 + u0 = v0, i = 0

−µ −uN + uN−1

h2 + uN = vN , i = N

wobei in der zweiten und dritten Zeilen,“u0 = u−1” beziehungsweise “uN = uN+1”

(virtuelle Randbedingungen u′ = 0) zu sehen sind.I Kontinuum-Ansatz fur Optimalitat: Richtungableitung sind

Null fur alle Storungen w :

0 =δJδu

(u; w) = limε→0

ddε

J(u + εw)

14


I Gilt wenn u das Randwertproblem erfullt:−µu′′ + u = v , Ω

u′ = 0, ∂Ω

I Vergleich mit dem diskreten Ergebnis zeigt:

u′′(xi) ≈ui+1 − 2ui + ui−1

h2

Approximation der zweiten Ableitung.I Hausaufgabe: Schreibe eine Approximation Jh(u) fur

J(u) =

∫Ω

(u − v)2dx + µ

∫Ω

√|u′|2 + ε2dx

und leite die Optimalitatsbedingung ∇Jh(u) = 0 her.

15

Zusammenfassung unserer Fragestellungen

I Wie leitet man Formeln zur Approximation eine Integralsher?

I Formeln fur Ableitungen?I Wann wissen wir, dass diese Formeln konvergieren wenn

h→ 0? Und zwar schnell?I Wie lost man die entstehenden Gleichungssysteme?I Wenn das Integral

∫Ω |u

′|2dx in J mit∫

Ω |u′|dx ersetzt

wird,I ist das Entrauschen-Ergebnis viel besser, aberI das Gleichungssystem ist nicht linear.I Wie lost man das nicht lineare Problem?

Dies werden wir spater mit dem Newtonschen Verfahrenlosen.

16

FehleranalyseI Fehler heissen Abbruchsfehler, wenn sie durch eine

Diskretisierung und trotz exakter Arithmetik entstehen, z.B.der Rest in

u′′(xi)−1h2 [u(xi+1)− 2u(xi) + u(xi−1)]

wird nicht berechnet sondern abgebrochen.I Fehler heissen Rundungsfehler, wenn Zahlen einer

Rechnung nicht exakt gespeichert werden konnen.I Diese konnen ublicherweise so grafisch dargestellt

werden:

17

FehleranalyseI Wie konnen wir die Abbruchsfehler abschatzen?

Satz (Taylor): Sei u ∈ Ck [a,b] und u(k+1) existiert in [a,b]. Seix0 ∈ [a,b]. ∀x ∈ [a,b], ∃ξ zwischen x und x0 sodass

u(x) =k∑

m=0

u(m)(x0)

m!(x − x0)m +

u(k+1)(ξ)

(k + 1)!(x − x0)k+1

I Anwendung fur die zweite Ableitung: xi+1 − xi = h

u(xi+1) = u(xi) + u(1)(xi)h + u(2)(xi)h2

2+ u(3)(xi)

h3

6+ u(4)(ξi+1)

h4

24

u(xi−1) = u(xi)− u(1)(xi)h + u(2)(xi)h2

2− u(3)(xi)

h3

6+ u(4)(ξi−1)

h4

24

Summe:

u(xi+1)− 2u(xi) + u(xi−1)

h2 = u′′(xi) + F (h)

18

FehleranalyseI Der Abbruchsfehler ist:

F (h) = [u(4)(ξi+1) + u(4)(ξi−1)]h2

24

I Man sagt, F (h) = O(h2).

Def: F (x) = O(G(x)) fur x → x0 (x0 ∈ R ∪ ±∞) wennlimx→x0 |F (x)/G(x)| <∞. F (x) = O(G(x)) wennlimx→x0 |F (x)/G(x)| = 0.I Wenn die exakte Losung erfullt |u(4)(x)| ≤ c, ∀x , dann gilt:

|F (h)| ≤ |u(4)(ξi+1) + u(4)(ξi−1)|h2

24

≤ (|u(4)(ξi+1)|+ |u(4)(ξi−1)|) h2

24

≤ 2ch2

24≤ ch2

12⇒ lim

h→0|F (h)/h2| ≤ c

12

Ahnlich gilt F (h) = O(h).19

KonvergenzBemerkung: Wir haben gezeigt, der Abbruchsfehler dieserApproximation der zweiten Ableitung konvergiert zu Null mitVerfeinerung. Wir haben aber noch nicht gezeigt,[ µ

h2 D + I]

u = v . . .

−µu′′ + u = v , Ω

u′ = 0, ∂Ω

dass die Losung des linearen Gleichungssystems zur exaktenLosung des Randwertproblems mit Verfeinerung konvergiert,

maxi |ui − u(xi)| → 0, h→ 0aber der Taylorsatz gehort zum Werkzeug.

Bemerkung: Die Optimalitatsbedingung ∇Jh(u) = 0 fur

J(u) =

∫Ω

(u − v)2dx +

∫Ω

√|u′|2 + ε2dx

ist nicht linear. Wir brauchen ein Verfahren, um die Nullstelle zubestimmen...20

Nullstellen

I Einfaches Beispiel: f : R→ R, finde die Nullstelle x0. Wirhaben durch Versuch und Irrtum x1 und x2 gefunden,sodass f (x1)f (x2) < 0. Gibt es eine Nullstelle zwischen x1und x2?

Satz (Zwischenwert): Wenn f ∈ C[a,b] und k zwischen f (a)und f (b) liegt, dann existiert c ∈ (a,b) sodass f (c) = k gilt.

I Gegenbeispiele:

21

Nullstellen

I Bekannte Verfahren zur Bestimmung einer Nullstelle:Bisektion

f (a)f (b) < 0 . . . (a + b)/2 = c →

a, f (c)f (b) < 0b, f (a)f (c) < 0

und Newton

xk+1 = xk − f (xk )/f ′(xk ), k = 1,2, . . .

I Unsere Fragen:I Ist ein Verfahren konvergent?I Was soll das Abbruchskriterium sein?I Was ist die Geschwindigkeit der Konvergenz?

22

Konvergenz und Stabilitat

Def: Die Folge xk konvergiert zu x∗ mitKonvergenzgeschwindigkeit O(yk ) wenn xk = x∗ +O(yk ).

I Beispiel: xk = (k + 1)/k2 = 0 +O(1/k):limk→∞ |xk/(1/k)| = limk→∞(k + 1)/k = 1.

Ahnlich yk = (k + 3)/k3 = 0 +O(1/k2).

I Falls (Rundungs)fehler in einer Iteration auftauchen,konnen sie sich schlecht aufbauen?

I D.h. ist das Verfahren stabil?I Wir stellen die Fehler als Anfangsstorung E0 dar...

23

Konvergenz und StabilitatDef: Sei Ek der Fehler eines Verfahrens nach k Operationen,wenn der Anfangsfehler E0 ist. Wenn gilt

|Ek | ≤ ck |E0| 0 < c 6= c(k)sagt man, das Fehlerwachstum ist linear. Wenn nur gezeigtwerden kann,

|Ek | ≤ cak |E0| 0 < c 6= c(k),a > 1sagt man, das Fehlerwachstum ist exponentiell. DaFehlerwachstum in der Regel unvermeidlich ist, sagen wir, einVerfahren mit hochstens linearem Fehlerwachstum ist stabil.I Beispiel: Wir berechnen die Folge 3−k .

I Methode 1: p0 = 1, pk = 13 pk−1. Probiere sie! Sie ist stabil!

I Methode 2: p0 = 1, p1 = 13 , pk = 10

3 pk−1 − pk−2, z.B.p2 = 10

313 − 1 = 1

32 (mit exakter Arithmetik). Probiere sie!Sie ist nicht stabil!

I Die Fehler in diesem Beispiel sind Rundungsfehler. Wieentstehen sie genau?

24

Darstellung von Zahlen im Computer

I Beispiel: IEEE Fließkommazahl einfacher Genauigkeit(32 bits, 4 bytes).

e0· · · →e7 m1 · · · · · · m23

8bits︷︸︸︷ 8bits︷︸︸︷ 8bits︷︸︸︷ 8bits=1byte︷︸︸︷σ︸︷︷︸

e︸︷︷︸

m

σ →Vorzeichen, e→Exponent, m→Mantissa

e = 〈e0, . . . ,e7〉, m = 〈m1, . . . ,m23〉, ek ,mk ∈ 0,1

V = (−1)σ, E = −127 +7∑

k=0

ek2k , M = 1 +23∑

k=1

mk2−k

Fließkommazahl: x = V · 2E ·M

25


I Beispiel: x = 1

x = (1)10 = (1.0)2 × 20 (M ∈ [1,2))M = (1)2 ⇒ mk = 0, k = 1, . . . ,230 = E = −(127)10 + e⇒ e = (127)10

Wie bringt man (127)10 auf die binare Darstellung?

Algorithmus fur die dezimalen Ziffer von 127:12710 = 12 + 7

10 ,1210 = 1 + 2

10 ,1

10 = 0 + 110

Algorithmus fur die binaren Ziffer von 127:127

2 = 63 + 12 ,

632 = 31 + 1

2 ,312 = 15 + 1

2 ,152 = 7 + 1

2 ,72 = 3 + 1

2 ,32 = 1 + 1

2 ,12 = 0 + 1

2

Also e = (127)10 = (1111111)2

x (IEEE): 00111111 10000000 00000000 0000000026


I Beispiel: x = 0.375

Algorithmus fur die dezimalen Ziffer von 0.375:0.375× 10 = 3 + 0.75,0.750× 10 = 7 + 0.50, 0.50× 10 = 5 + 0Algorithmus fur die binaren Ziffer von 0.375:0.375× 2 = 0 + 0.75,0.750× 2 = 1 + 0.50, 0.50× 2 = 1 + 0Alsox = (0.375)10 = (0.011)2 = (1.1)2 × 2−2 (M ∈ [1,2))M = (1.1)2 ⇒ m1 = 1,mk = 0, k = 2, . . . ,23e − (127)10 = E = −(2)10 ⇒ e = (125)10 = (01111101)2

x (IEEE): 00111110 11000000 00000000 00000000

27


I Beispiel: x = 12.375.

x = (12)10 + (0.375)10x = (1100)2 + (0.011)2 = (1100.011)2

x = (1.100011)2 × 23, (M ∈ [1,2))M = (1.100011)2 ⇒ m1 = 1,m5 = 1,m6 = 1, sonst mk = 0e − (127)10 = E = (3)10 ⇒ e = (130)10 = (10000010)2

x (IEEE): 01000001 01000110 00000000 00000000

28

Darstellung von Zahlen im ComputerI Beispiel: x = −0.1

Algorithmus fur die binaren Ziffer von 0.1:0.1× 2 = 0 + 0.2, 0.2× 2 = 0 + 0.4,0.4× 2 = 0 + 0.8, 0.8× 2 = 1 + 0.6,0.6× 2 = 1 + 0.2→ zuruck zum Anfang

Also x = −(0.1)10 = −(0.00011)2 = −(1.10011)2 × 2−4

Probe: x = −∑∞

k=1 2−4k − 12∑∞

k=1 2−4k = − 110

M = (1.10011)2 ⇒ m1 = 1, m4k = m4k+1 = 15k=1,sonst mk = 0?

e − (127)10 = E = −(4)10 ⇒ e = (123)10 = (1111011)2

Mit Aufrundung (m23 = 1)x (IEEE): 10111101 11001100 11001100 11001101Mit Abrundung (m23 = 0)x (IEEE): 10111101 11001100 11001100 11001100

29

Darstellung von Zahlen im ComputerI Fortsetzung des Beispiels: x = −0.1

Mit Aufrundung (m23 = 1)x (IEEE): 10111101 11001100 11001100 11001101

→ (0.100000001490116119384765625)10Mit Abrundung (m23 = 0)x (IEEE): 10111101 11001100 11001100 11001100

→ (0.0999999940395355224609375)10Mit einfacher Genauigkeit gibt es keine speicherbarenZahlen zwischen diesen!

I IEEE Fließkommazahl doppelter Genauigkeit (64 bits, 8bytes)

x = V · 2E ·M, E = −1023 +10∑

k=0

ek2k

V = (−1)σ, M = 1 +52∑

k=1

mk2−k

30

Fehler in der Darstellung?

Def: Wenn p eine Approximation zu p ist, ist der absoluteFehler |p − p|. Der relative Fehler ist |p − p|/|p| wenn p 6= 0.

Def: p approximiert p zu t signifikanten Ziffern wenn t ∈ N diegroßte Zahl ist, bei der der relative Fehler kleiner als 5× 10−t

ist.

I Beispiel: x = −0.1, x = −0.10000000149 · · · (Aufrundung)

5× 10−9 <|x − x ||x |

=1.49 · · · × 10−9

10−1 < 5× 10−8 ⇒ t = 8

31

Darstellung von Zahlen im ComputerI Relativer Fehler der Fließkommazahl-Darstellung fl(x), mit

Abrundung und p-Ziffer Genauigkeit M =∑p−1

k=0 mk2−k ,

|x − fl(x)||x |

=|V · 2E (

∑∞k=0)mk2−k − V · 2E (

∑p−1k=0)mk2−k |

|V · 2E (∑∞

k=0)mk2−k |

=

∣∣∣∣∣∣∞∑

k=p

mk2−k

∣∣∣∣∣∣/ ∣∣∣∣∣

∞∑k=0

mk2−k

∣∣∣∣∣≥1

≤∞∑

k=p

mk2−k

≤∞∑

k=p

2−k = 2−p∞∑

k=0

2−k =2−p

1− 2−1 = 2−(p−1)(< tol)

I IEEE einfache Genauigkeit, p = 24,2−p+1 ≈ 1.19× 10−7, t = 7

I IEEE doppelte Genauigkeit, p = 53,2−p+1 ≈ 2.22× 10−16, t = 16

32


I Verlust von signifikanten Ziffern: Ausloschung.Bildlich gesehen:

x1 = (−1)σ2EM1

x2 = (−1)σ2EM2

gleich anders

x2 − x1 E #######000 · · · 000

(−1)σ2E (M2 −M1) E #######000 · · · 000

Nur die Ziffer ####### uberleben Subtraktion. In exakterArithmetik ware der letzte Bereich 000 · · · 000 mit richtigenZiffern besetzt.

33


I Beispiel: Berechne die Nullstellen vonp(x) = x2 + 62.10x + 1

mit 4-Dezimal-zifferiger Genauigkeit.

−62.10 + [(62.10)2 − 4.000 · 1.000 · 1.000]12

2.000

→ −62.10 + [3856− 4.000]12

2.000→ −0.04000

2.000= −0.02000

6= x1 = −0.01610723

I Besserer Weg: A− B = (A2 − B2)/(A + B)⇒

−b +√

b2 − 4ac2a

=−b2 + (b2 − 4ac)

2a[b +√

b2 − 4ac]→ −0.01610

I Auch anfallig fur Ausloschung: Skalarprodukt

34


I Auch anfallig fur Ausloschung: PolynomauswertungI Beispiel: An der Stelle x = 4.71 das Polynom

p(x) = x3 − 6x2 + 3x − 0.149mit 3-Dezimal-zifferiger Genauigkeit auswerten.exakt:104.487111− 133.1046 + 14.13− 0.149 = −14.6364893-zifferig:105− 133 + 14.1− 0.149→ −14.0, 4% relativer Fehler

I Horner Algorithmus:p(x) = x [x2 + 6x + 3]− 0.149 = x [x(x − 6) + 3]− 0.149

p ← x − 6, p ← x · p, p ← p + 3, p ← x · p, p ← p − 0.149

Ergebnis: -14.5, 0.25% relativer Fehler

35

Lineare Gleichungssysteme: Direkte MethodenGaußsche Elimination zur Losung von Ax = b. Beispiel:I Das System:

E1E2

E3

3x1 + 6x2 9x3 = −122x1 + 5x2 − 8x3 = −11

x1 − 4x2 − 7x3 = −10

I E2 ← E2 − 23E1, E3 ← E3 − 1

3E1,

3x1 + 6x2 + 9x3 = −12x2 − 14x3 = −3

− 6x2 − 10x3 = −6

I E3 ← E3 + 61E2,

3x1 + 6x2 + 9x3 = −12x2 − 14x3 = −3− 94x3 = −24

36

Gaußsche Elimination und Ruckwarts SubstitutionI Ruckwarts Substitution:

3x1 + 6x2 + 9x3 = −12 ⇒ x1 = 13(−12− 6x2 − 9x3) = −278

47x2 − 14x3 = −3 ⇒ x2 = −3 + 14x3 = 27

47 ↑− 94x3 = −24 ⇒ x3 = 12

47 ↑

I Zur Losung von Ax = b kann man die Schritte auf dieerweiterte Matrix [A|b] durchfuhren:

[A|b] =

1 1 0 | 42 1 −1 | 13 −1 −1 | −13

∈ R3×4 (Rn×(n+1))

I Wenn es mehrere rechte Seiten gibt,

Ax i = bi , i = 1,2, . . . ,m oderAX = B, X = [x1, . . . ,xm], B = [b1, . . . ,bm]

kann man die Schritte gleichzeitig auf die erweiterte Matrix[A|B] durchfuhren. Aber nicht wenn bi von x i−1 abhangt!

37

Diskretisierung einer PDGI Die Diskretisierung [ µh2 D + I]u = v des Randwertproblems,

−µu′′ + u = v , Ωu′ = 0, ∂Ω

wird sehr groß mit Verfeinerung, aber sie ist trotzdemimmer tridiagonal; deswegen ist es ziemlich billig,Gaußsche Elimination fur dieses System durchzufuhren.

I Bessere Motivation: Ω = (0,1)2, ∆u = uxx + uyy−µ∆u + u = v , Ω

∂u/∂n = 0, ∂Ω

1 · · · N+1

N+2 · · · 2N+2

· · ·· · ·· · ·

(N+1)2

uxx + uyy ≈ui+1,j − 2uij + ui−1,j

h2 +ui,j+1 − 2uij + ui,j−1

h2

38

Diskretisierung einer PDG

I Gleichungssystem:

. . .

. . .. . .

. . .. . . 0

− µ

h2 − µ

h24µh2 + 1 − µ

h2 − µ

h2

0. . .

. . .. . .

. . .. . .

u0,0

...uN,0

...u0,N

...uN,N

=

v0,0

...vN,0

...v0,N

...vN,N

-

O(N) -

O(N)

I Obwohl die Matrix dunn besetzt ist (sparse, d.h. ganzwenig nicht triviale Elemente), fuhrt Gaußsche Eliminationzur Auffullung (fill-in) der Zwischendiagonalen.

39

Automatisierung der Gaußschen Elimination

I Um die Losungsschritte der Gaußschen Elimination zuautomatisieren, mussen wir diese symbolisch schreiben:

A ← [A|b]

2. Zeile:

A21 = A21/A11 (Multiplikator: m21)A2j = A2j − A21 · A1j , j = 2, . . . ,n(+1)

3. Zeile:


...

i . Zeile:

Ai1 = Ai1/A11 (Multiplikator: mi1)Aij = Aij − Ai1 · A1j , j = 2, . . . ,n(+1)

i = 2, . . . ,n

I Ergebnis:

A11

m21

mn1

...

Pivot

neues b

neue Zielspalte

-

sonst Null40

Automatisierung der Gaußschen Elimination

I Mit der nachsten Zielspalte:

3. Zeile:


...

i . Zeile:

Ai2 = Ai2/A22 (Multiplikator: mi2)Aij = Aij − Ai2 · A2j , j = 3, . . . ,n(+1)

i = 3, . . . ,n

I Ergebnis:

A11

A22m21

mn1

...

m32

mn2

...

...

Pivot neues b

neue Zielspalte

-

sonst Null

I Wird mit Zielspalten fortgesetzt, bis man ein Dreieck fur dieMultiplikatoren hat.

41

Automatisierung der Gaußschen EliminationI Pseudo-Kode zur Gaußschen Elimination:A <- [A,b] % Erweiterungfor k=1,...,n-1 % Pivot-Index

for i=k+1,...,n % Zeilen-IndexA ik = A ik / A kk % Multiplikatorfor j=k+1,...,n+1 % Spalten-Index

A ij = A ij - A ik*A kjend

endend

I Pseudo-Kode zur Ruckwarts Substitution:x n = A n,n+1 / A nnfor i=n-1,...,1

s=0for j=i+1,...,n

s = s + A ij*x jendx i = (A i,n+1 - s) / A ii

end42

flops der Gaußschen EliminationI Wie viele Operationen kostet Gaußsche Elimination?

I Divisionenk = 1, ...,n− 1

i = k + 1, ...,n → (n − k)

n−1∑k=1

(n − k)

I Multiplikationen

k = 1, ...,n− 1i = k + 1, ...,n → ×(n − k)

j = k + 1, ...,n(+1) → (n + 1− k)

n−1∑k=1

(n−k)(n+1−k)

I Subtraktionen

...ebenso :n−1∑k=1

(n − k)(n + 1− k)

I Nutzliche Summen: (vgl.∫ m

0 kp−1dk = mp/p)

m∑k=1

1 = m,m∑

k=1

k =m(m + 1)

2,

m∑k=1

k2 =m(m + 1)(2m + 1)

643

flops der Gaußschen EliminationI Kosten fur Gaußsche Elimination:

I Divisionen

n−1∑k=1

(n−k) = nn−1∑k=1

1−n−1∑k=1

k = n(n−1)− (n − 1)n2

=n(n − 1)

2

I Multiplikationen

n−1∑k=1

(n − k)(n + 1− k)︸︷︷︸n(n+1)−(2n+1)k+k2

= (n2 + n)n−1∑k=1

1− (2n + 1)n−1∑k=1

k +n−1∑k=1

k2 =

(n2 + n)(n − 1)− (2n + 1)(n − 1)n

2+

(n − 1)n[2(n − 1) + 1]

6=

n(n2 − 1)

3

I Subtraktionen: auch n(n2 − 1)/3I Gesamte Operationen:

n(n − 1)/2 + 2n(n2 − 1)/3 = n(n − 1)(4n + 7)/6 = O(n3)

44

flops der Gaußschen Elimination

I Kosten fur Ruckwarts Substitutionen:I Divisionen: nI Subtraktionen: n − 1I Additionen:

i = n− 1, ...,1j = i + 1, ...,n → (n − i)

n−1∑i=1

(n − i) =n(n − 1)

2

I Multiplikationen: auch n(n − 1)/2I Gesamte Operationen:

2n − 1 + 2n(n − 1)/2 = n2 + n − 1 = O(n2)

I Botschaft: Ruckwarts Substitution ist viel billiger!I Hausaufgabe: Einen Pseudo-Kode fur Vorwarts

Substitution schreiben und die gesamten Operationenausrechnen.

45

Rekursive Losungen: EvolutionsgleichungI Wenn man sehr viele Systeme losen muss,

Ax i = bi , i = 1,2, . . . ,mist es vorteilhaft wenn man die folgenden Schrittedurchfuhren kann:I Gaußsche Elimination nur einmal mit Kosten O(n3)

machen,I das Ergebnis irgendwie speichern undI dann m-Mal Ruckwarts und Vorwarts Substitutionen mit

Kosten O(mn2) machen.Sonst sind die Kosten O(mn3)!

I Beispiel: Evolutionsgleichung,u′ = Bu

u(0) = u0

zeitlich

−→diskretisieren

u(tk+1)− u(tk )

tk+1 − tk = Bu(tk+1)

oder mit τ = tk+1 − tk ,

Auk+1 = [I − τB]uk+1 = uk , k = 0,1, . . .46

LU Zerlegung

I Wie kann man die O(n3) Arbeit fur Gaußsche Eliminationgunstig speichern?

I Definiere A(1) = A. Ergebnis nach der ersten Zielspalte istA(2) = M(1)A(1) wobei:

A(2) =

0

0

...M(1) =

0· · ·01−m21

−mn1

... I(n−1)

. . .

. . .

Probe: Sei A =

z1...

zn

, M(1)A =

z1z2 −m21z1

...

47

LU Zerlegung

I Ergebnis nach der zweiten Zielspalte ist A(3) = M(2)A(2)

wobei:

A(3) =

0

0

...0

0

...M(2) =

0· · ·0

0· · ·0

1

10

0

...

−m32

−mn2

... I(n−2)

. . .

. . .

I Am Ende von Gaußscher Elimination:A(n) = M(n−1)A(n−1) = M(n−1) · · ·M(1)A

wobei A(n) eine obere Dreiecksmatrix ist.I Definiere U = A(n). Bemerke:

A = M(1)−1 · · ·M(n−1)−1UI Was sind die Inversen von M(i) explizit?

48

LU ZerlegungI Bemerke:

1 0 · · · 0. . .

. . .

I(n−1)

m21

mn1

...∗

0· · ·01-m21

-mn1

... I(n−1)

. . .

. . .

= M(1)−1M(1) = I(n)

I Es gelten

M(i)−1 =

001. . .0 1

mi+1,i

mn,i

0 ... I(n−i)

. . .

. . .

i = 1, . . . ,n − 1

I Bemerke:1 0 · · · 0. . .

. . .

I(n−1)

m21

mn1

...∗

0· · ·01

1 0 · · · 0m32

mn2

... I(n−2)

. . .

. . .

=

0· · ·01

1 0 · · · 0m21

mn1

m32

mn2

... ... I(n−2)

. . .

. . .49

LU ZerlegungI Es gilt:

0· · ·01

1

1

01mn,n−1

0 · · · 0m21

mn1

m32

mn2 · · ·

. . .. . .

. . ....... ...

M(1)−1 · · ·M(n−1)−1 = = L

I Wir haben gezeigt:

Satz: Wenn alle Pivot-Elemente a(k)kk nicht Null sind und alle

Multiplikatoren mjk = a(k)jk /a(k)

kk wohl definiert sind, kann dieMatrix mit A = LU faktorisiert werden, wobei

U = uij, uij =

a(n)

ij , 1 ≤ i ≤ n, i ≤ j ≤ n0, 2 ≤ i ≤ n, 1 ≤ j < i

L = `ij, `ij =

mij , 2 ≤ i ≤ n, 1 ≤ j < i

1, 1 ≤ i ≤ n, j = i0, 1 ≤ i ≤ n, i < j ≤ n

50

Rekursive Losungen mit LU

I Mit dieser Zerlegung kann man mehrere sequentielleSysteme gunstig losen: Ax i = bi ⇔ Ly i = bi ,Ux i = y i .

I Beispiel: Evolutionsgleichung,

u′ = Bu

u(0) = u0

zeitlich

−→diskretisieren

u(tk+1)− u(tk )

tk+1 − tk = Bu(tk+1)

mit τ = tk+1 − tk ,

Auk+1 = [I − τB]uk+1 = uk , k = 0,1, . . .

I Durch die LU Zerlegung A = LU, O(n3),

Vorwarts Substitution: Luk+ 12 = uk , O(n2)

Ruckwarts Substitution: Uuk+1 = uk+ 12 , O(n2)

k = 0,1, . . .

51

Pivot StrategienPivot Strategien:I Beispiel:[

0.003000 59.145.291 −6.130

] [x1x2

]=

[59.1746.78

]I Exakte Losung:

x∗1 = 10.00, x∗2 = 1.000

I Mit 4-zifferiger Arithmetik,

m21 =5.291

0.003000= 1763.66→ 1764

I Erster Schritt von Gaußscher Elimination,[0.003000 59.14 59.17

(1764) −104300 −104400

]I In exakter Arithmetik ware es gewesen,[

0.003000 59.14 59.17(1763.66) −104309.3766 −104309.3766

]52

Pivot Strategien

I Nun Ruckwarts Substitution mit 4-zifferiger Arithmetik

x2 =−104400−104300

→ 1.001 ≈ x∗2

x1 =59.17− (59.14)(1.001)

0.003000→ −10.00 6= 10.00 = x∗1

I Hausaufgabe: Zeige, wenn die Zeilen getauscht werden,[5.291 −6.130 46.78

0.003000 59.14 59.17

]bekommt man

x2 → 1.000, x1 → 10.00

53

Pivot Strategien

Tauschmit I(·)

k →p →

I Pseudo-Kode zur Gaußschen Elimination mit Pivot-Suche:A <- [A,b] % ErweiterungI(i)=i, i=1,...,n % Zeiger Initialisierungfor k=1,...,n-1 % Pivot-Index

p=argmax k<=i<=n |A I(i),k| % Großtes Pivot-Elementif A I(p),k = 0 then stop % Test fur RegularitatIt = I(k), I(k) = I(p), I(p) = It % Virtueller Zeilentauschfor i=k+1,...,n % Zeilen-Index

A I(i),k = A I(i),k / A I(k),k % Multiplikatorfor j=k+1,...,n+1 % Spalten-Index

A I(i),j = A I(i),j - A I(i),k*A I(k),jend

endend

I Pseudo-Kode zur Ruckwarts Substitution:x n = A I(n),n+1 / A I(n),nfor i=n-1,...,1

s=0for j=i+1,...,n

s = s + A I(i)j*x jendx i = (A I(i),n+1 - s) / A I(i),i

end54

LU Zerlegung mit PermutationI Wie kann man mit Pivot-Suche die Matrix schon zerlegen?

Def: Eine Permutationsmatrix P hat genau einen nicht trivialenEintrag (mit dem Wert 1) in jeder Spalte und in jeder Zeile, undes gilt PPT = I.I Beispiel:

P =

[0 11 0

]A =

[0 13 2

]PA =

[3 20 1

]Bemerkung: I(i) : i = 1, . . . ,n ist eine gegebenePermutation von den Zahlen 1, . . . ,n. Die Permutationsmatrix

Pi,j = δI(i),j erfullt A = PA wobei Aij = AI(i),j .

I Wegen dem obigen Satz haben wir A = LU von demAlgorithmus mit Pivot-Suche, und PA = LU wobei

Pi,j = δI(i),j Lij =

A(n)

I(i),j , i > j1, i = j0, i < j

Uij =

A(n)

I(i),j , i ≤ j0, i > j

55

Skalierte Pivotsuche

Weitere Pivotsuche-StrategienI Das 2× 2 Beispiel hatten wir so umschreiben konnen:[

30.00 5914005.291 −6.130

] [x1x2

]=

[59170046.78

](← ×104)

I Erster Schritt von Gaußscher Elimination ohne Tauschen,

m21 =5.29130.00

→ 0.1764[

30.00 591400 591700(0.1764) −104300 −104400

]fuhrt zur gleichen Losung wie vorher ohne Pivotsuche:

x2 → 1.001, x1 → −10.00 6= 10.00 = x∗1

I Neue Strategie: Man sollte die Zeilen skalieren.

56

Skalierte PivotsucheI Pseudo-Kode zur G.E. mit skalierter Pivotsuche:

A <- [A,b] % Erweiterungs i=max j |A i,j|, i=1,...,n % ZeilenskalenI(i)=i, i=1,...,n % Zeiger Initialisierungfor k=1,...,n-1 % Pivot-Index

p=argmax k<=i<=n |A I(i),k|/s I(i) % Großtes Pivot-Elementif A I(p),k = 0 then stop % Test fur RegularitatIt = I(k), I(k) = I(p), I(p) = It % Virtueller Zeilentauschfor i=k+1,...,n % Zeilen-Index

A I(i),k = A I(i),k / A I(k),k % Multiplikatorfor j=k+1,...,n+1 % Spalten-Index

A I(i),j = A I(i),j - A I(i),k*A I(k),jend

endend

I Pseudo-Kode zur Ruckwarts Substitution:x n = A I(n),n+1 / A I(n),nfor i=n-1,...,1

s=0for j=i+1,...,n

s = s + A I(i)j*x jendx i = (A I(i),n+1 - s) / A I(i),i

end57

Totale PivotsucheBemerkung: Kosten fur Tauschen und Skalieren sind O(n2)

Weitere Pivotsuche-Strategien: Totale PivotsucheI Virtueller Zeilentausch mit einem Zeiger I(·).I Virtueller Spaltentausch mit einem Zeiger J(·).I Bestimmung des Pivot-Elements:

(p,q)=argmax (k≤i,j≤n) |A I(i),J(j)|I Zusatzliche Arbeit ist O(n3).I Geeignet wenn das Problem sehr verschiedene Skalen

hat, d.h. wenn die Matrix A sehr steif ist.I Verwendung:

Zeilentausch: Ax = b → P1Ax = LUx = P1bSpaltentausch: Ax = b → AP2P2x = LUP2x = b

Beides: (P1AP2)P2x = LUP2x = P1b

Ly = P1b, Uz = y , x = PT2 z

58

Inverse und DeterminanteBemerkung: Falls die Inverse notwendig ist, kann man siegunstig durch die Losung des folgenden Systems berechnen:

AX = I ⇒ X = A−1

Die direkte Berechnung ist aber viel teurer.

Bemerkung: Die Determinante,

δk = det

a11 · · · a1k...

...ak1 · · · akk

kann man mit Gaußscher Elimination so berechnen,

δk = a(1)11 · a

(2)22 · · · a

(k)11

falls keine Pivotsuche notwendig ist. Die direkte Berechnung istaber viel teurer.59

Diagonal Dominante MatrizenWann wird keine Pivotsuche notwendig?

Def: Eine Matrix aij = A ∈ Rn×n ist streng diagonal dominantwenn

|aii | >n∑

i 6=j=1

|aij |, i = 1, . . . ,n

Bei Gleichheit ist sie schwach diagonal dominant.

I Beispiel: Wir haben die folgende Diskretisierung,−µu′′ + u = v , Ω

u′ = 0, ∂Ω−→ Au =

[ µh2 D + I

]u = v

wobei

A =

µh2 + 1 − µ

h2

− µh2

2µh2 + 1 − µ

h2

. . .. . .

. . .

− µh2

2µh2 + 1 − µ

h2

− µh2

µh2 + 1

60

SPD MatrizenI A ist streng diagonal dominant:

(x0, xN) i = 0,N : µh2 + 1 > | − µ

h2 |√

(xi) 0 | − µ

h2 |+ | − µh2 |

√

Satz: Fur eine streng diagonal dominante Matrix, kannGaußsche Elimination ohne Pivotsuche stabil durchgefuhrtwerden.

Def: Eine symmetrische Matrix A ∈ Rn×n ist positiv definit(SPD) wenn xTAx > 0, ∀x ∈ Rn, x 6= 0.

Bemerkung: Wenn A ∈ Rn×n symmetrisch ist, sind dieEigenwerte reel und es gilt

minλ ∈ σ(A) ≤ xTAxxTx

≤ maxλ ∈ σ(A)

Wenn minλ ∈ σ(A) > 0 gezeigt werden kann, ist A SPD.61

Gerschgorin SatzI Dies kann man fur das obige Beispiel mit dem folgenden

Satz zeigen.

Satz (Gerschgorin): Fur aij = A ∈ Rn×n gilt

σ(A) ⊂n⋃

i=1

Ri wobei Ri =

z ∈ C : |z − aii | ≤n∑

i 6=j=1

|aij |

Weiters liegen genau k Eigenwerte in ∪k

j=1Rij wenn dieseMenge einen leeren Schnitt mit den anderen Scheiben hat.I Beispiel: Wo liegen die Eigenwerte fur das obige A?

R0 = B( µh2 + 1, µh2 )

Ri = B(2µh2 + 1, 2µ

h2 )

RN = B( µh2 + 1, µh2 )

C

&%'$

1

µ

h2 + 12µh2 + 1

4µh2 + 1

@

@62

Notwendigkeit einer PivotsucheI Das obige A ist symmetrisch und es gilt

1 ≤ minλ ∈ σ(A) ≤ xTAxxTx

also ist A SPD.Satz: Fur eine SPD Matrix kann Gaußsche Elimination ohnePivotsuche stabil durchgefuhrt werden.

Botschaft: Fur viele Matrizen in Anwendungen ist einePivotsuche nicht notwendig, aber es ist sicher nicht immer so.

Folgesatz: Eine Matrix A ∈ Rn×n ist SPD genau dann wenn∃L ∈ Rn×n eine untere Dreiecksmatrix sodass A = LLT.I Mit der Gleichung fur LLT = A = aij, L = `ij,

aij =

min(i,j)∑k=1

`ik`jk

a11 = `211a21 = `21`11a12 = `11`21

. . .

landet man auf dem folgenden Algorithmus:63

Pseudo-Kode fur den Cholesky AlgorithmusL 11 =

√A 11

for i=2,...,nL i1 = A i1 / L 11

endfor j=2,...,n-1

s = 0for k=1,...,j-1

s = s + L jk * L jkendL jj =

√A jj− s

for i=j+1,...,ns = 0for k=1,...,j-1

s = s + L ik * L jkendL ij = [A ij - s]/L jj

endends = 0for k=1,...,n-1

s = s + L nk * L nkendL nn =

√A nn− s

I Hausaufgabe: Mit demobigen Gleichungssystemleite diese Methode her.

I Bemerke: Der CholeskyAlgorithmus ist nur fur SPDMatrizen.

I Hausaufgabe: Schatze dieflops des Kodes ab.

64

BandmatrizenI Die Bandbreite einer solchen Matrix,

A =

a1,1 · · · a1,q+1 0 · · · 0...

. . .. . .

. . .. . .

...

ap+1,1. . .

. . .. . .

. . . 0

0. . .

. . .. . .

. . . an−q,n...

. . .. . .

. . .. . .

...0 · · · 0 an,n−p · · · an,n

∈ Rn×n

is p + q + 1, die kleinste Anzahl der benachbartenDiagonalen, zu denen die nicht trivialen Elementeeingeschrankt sind.

I Bemerke:

p ≤ n − 1, q ≤ n − 1, Bandbreite ≤ 2n − 1.

Wenn p + q + 1 << 2n − 1 ist A eine Bandmatrix.65

Bandmatrizen

I Um Speicherplatz zu reduzieren, wird A als Liste derDiagonalen gespeichert:

B =

0 · · · 0 a1,1 · · · · · · a1,q+1... . .

.. .. ...

......

...

0 . .. ...

......

... an−q,n

ap+1,1...

......

... . ..

0...

......

... . ..

. .. ...

an,n−p · · · · · · an,n 0 · · · 0

∈ Rn×(p+q+1)

I Bemerke: Die nicht trivialen Elemente von A = ai,jlassen sich bezuglich der Elemente von B = bα,β sodarstellen:

ai,j = bi,j−i+p+1

66

Sparse SpeicherformatI Hausaufgabe: Schreibe einen Pseudo-Kode zur

Implementierung des Gauß Algorithmus fur A bezublich Bmit moglichst wenig arithmetischen Operationen.

I Es gibt auch sparse format:

A = sparse(i,j,w,m,n) ⇒ A ∈ Rm×n, Ai(k),j(k) = w(k)

I Beispiel:

A =

µh2 + 1 − µ

h2

− µh2

2µh2 + 1 − µ

h2

. . .. . .

. . .

− µh2

2µh2 + 1 − µ

h2

− µh2

µh2 + 1

A = spdiags([ µh2 + 1; (2µ

h2 + 1)ones(N − 1,1); µh2 + 1],0,N + 1,N + 1)

− spdiags( µh2 ones(N + 1,1),+1,N + 1,N + 1)

− spdiags( µh2 ones(N + 1,1),−1,N + 1,N + 1)67

Lineare Gleichungssysteme: Iterative MethodenIterative Verfahren zur Losung von Ax = b.I Um Ax = b zu losen, werden wir solche Iterationen

formulieren,

xk+1 = T xk + c, k = 0,1, . . .

I Um Konvergenz zu untersuchen, mussen wir eine gewisseInfrastruktur aufbauen.

Def: Eine Vektornorm ist eine Funktion ‖ · ‖ : Rn → R mit denEigenschaften:

1. ‖x‖ ≥ 0, ∀x ∈ Rn (∑n

i=1 xi keine Norm)2. ‖x‖ = 0⇔ x = 0, (|x1| keine Norm fur n > 1)3. ‖αx‖ = |α|‖x‖, ∀x ∈ Rn, ∀α ∈ R (

∑ni=1 x2

i keine Norm)4. ‖x + y‖ ≤ ‖x‖+ ‖y‖, ∀x ,y ∈ Rn

(f (√∑n

i=1 x2i ) keine Norm wenn f nicht konvex)

68

VektornormenSatz: Die folgenden sind Normen:

‖x‖`p =

[n∑

i=1

|xi |p] 1

p

,1 ≤ p <∞, ‖x‖`p = max1≤i≤n

|xi |, p =∞

Satz: ∀x ,y ∈ Rn gilt

|x · y | ≤n∑

i=1

|xiyi | ≤ ‖x‖`p‖y‖`q

wobei 1/p + 1/q = 1.&%'$ ‖x‖2=1-‖x‖∞=1 @

@

@@

-‖x‖1=1

Beweis (fur p = q = 2): Mit x ,y ∈ Rn beliebig, seienx = 〈|x1|, . . . , |xn|〉 und y = 〈|y1|, . . . , |yn|〉. Dann gilt

‖x‖22 + 2|∑n

i=1 xiyi |+ ‖y‖22 ≤‖x‖22 + 2

∑ni=1 |xiyi |+ ‖y‖22 = ‖x‖22 + 2x · y + ‖y‖22

= ‖x + y‖22 ≤ [‖x‖2 + ‖y‖2]2

= ‖x‖22 + 2‖x‖2‖y‖2 + ‖y‖2269

Aquivalenz von Normen

Satz: In Rn sind alle Normen aquivalent, d.h. fur 2 beliebigeNormen ‖x‖A und ‖x‖B auf Rn ∃c1, c2 3

c1‖x‖A ≤ ‖x‖B ≤ c2‖x‖A ∀x ∈ Rn

Insbesondere:

Satz: ∀x ∈ Rn,

‖x‖∞ ≤ ‖x‖2 ≤√

n‖x‖∞

Beweis: Fur x ∈ Rn sei |xk | = ‖x‖∞. Dann gilt

‖x‖2∞ = |xk |2 ≤n∑

i=1

|xi |2 = ‖x‖22 ≤n∑

i=1

|xk |2 = n|xk |2 = n‖x‖2∞

Hausaufgabe: Zeige Aquivalenz der Normen ‖ · ‖1 und ‖ · ‖2.

70

Qualitative Eigenschaften der NormenI Soll man einen Fehler ‖x − x∗‖ mit der `1-, `2- oder`∞-Norm messen? Mit einer anderen Norm?

I Antwort: Es hangt vom Kontext ab!I Beispiel: Die gemessenen pH-Werte einer chemischen

Losung sind f = 〈a,b, . . . ,b〉 ∈ Rn+1.

Dieses Beispiel ubertreibt die Situation, in der ein einziger Wert einstatistischer Ausreißer ist, wahrend die anderen Werte naher beieinander liegen.

I In einem Blick ist b die beste Abschatzung vom pH.I Hausaufgabe: Fur e = 〈1, . . . ,1〉 zeige

a + nb1 + n

= arg minc∈R‖f−ce‖22 wahrend b = arg min

c∈R‖f−ce‖1

I Also laßt die ‖ · ‖2-Norm sich von Ausreißern beeinflussen.I Botschaft: Die ‖ · ‖1-Norm ist statistisch robuster.

71

MatrixnormenAuch mussen wir den Abstand zwischen Matrizen messen.

Def: Eine Matrixnorm ist eine Funktion ‖ · ‖ : Rn×n → R mit denEigenschaften:

1. ‖A‖ ≥ 0, ∀A ∈ Rn×n

2. ‖A‖ = 0⇔ A = 03. ‖αA‖ = |α|‖A‖, ∀A ∈ Rn×n, ∀α ∈ R4. ‖A + B‖ ≤ ‖A‖+ ‖B‖, ∀A,B ∈ Rn×n

5. ‖AB‖ ≤ ‖A‖‖B‖, ∀A,B ∈ Rn×n

Def: Fur eine Vektornorm ‖ · ‖V ist die entsprechendeinduzierte (oder naturliche) Matrixnorm ‖ · ‖M definiert durch:

‖A‖M = max‖x‖V=1

‖Ax‖V

Hausaufgabe: Zeige dass diese Funktion ‖ · ‖M eineMatrixnorm ist.

72

Charakterisierung bekannter MatrixnormenSatz: Die Matrixnormen ‖ · ‖1, ‖ · ‖2 und ‖ · ‖∞ sind fur eineMatrix Aij = A ∈ Rn×n explizit so gegeben:

‖A‖1 = max‖x‖1=1

‖Ax‖1 = max1≤j≤n

n∑i=1

|aij |

‖A‖2 = max‖x‖2=1

‖Ax‖2 =[maxλ ∈ σ(ATA)

] 12

‖A‖∞ = max‖x‖∞=1

‖Ax‖∞ = max1≤i≤n

n∑j=1

|aij |

Def: Eine Matrixnorm ‖ · ‖M ist mit einer Vektornorm ‖ · ‖Vkompatibel (oder vertraglich) wenn gilt

‖Ax‖V ≤ ‖A‖M‖x‖V, ∀x ∈ Rn, ∀A ∈ Rn×n

Hausaufgabe: Zeige fur eine Matrix Aij = A ∈ Rn×n, dieFrobenius Matrixnorm,

‖A‖F =

n∑i,j=1

|Aij |2 1

2

ist mit der Vektornorm ‖x‖2 kompatibel.73

SpektralradiusDef: Sei Rn mit einer gewissen Vektornorm ‖ · ‖V versehen. Furxk ,x ∈ Rn gilt limk→∞ xk = x genau dann wenn ∀ε > 0, ∃K , 3∀k ≥ K , ‖x − xk‖V < ε.

Satz: Fur xk = 〈(x1)k , . . . , (xn)k 〉,x = 〈x1, . . . , xn〉 ∈ Rn giltlimk→∞ xk = x genau dann wenn limk→∞(xi)k = xi .

Beweis: Hausaufgabe.

Def: Die Eigenwerte einer Matrix A ∈ Rn×n sind die Nullstellendes charakteristischen Polynoms

p(λ) = det(A− λI), σ(A) = λ : p(λ) = 0

und die entsprechenden Eigenvektoren sind die nicht trivialenLosungen (x 6= 0) zur Gleichung Ax = λx .

Def: Der Spektralradius einer Matrix A ∈ Rn×n ist

ρ(A) = max|λ| : λ ∈ σ(A)74

SpektralradiusI ρ(A) ist der großte Vergroßerungsfaktor.

Beispiel: A = SΛS−1, S = [x1,x2],

Λ =

[1 00 2

],

Λe1 = e1, Ax1 = x1Λe2 = 2e2, Ax2 = 2x2 ←

Satz: Sei ‖ · ‖? eine beliebige Matrixnorm auf Rn×n, die von derVektornorm ‖ · ‖? auf Rn induziert ist. Dann gilt

ρ(A) ≤ ‖A‖?, ∀A ∈ Rn×n

Beweis: Wahle λ, x so aus, dass Ax = λx , |λ| = ρ(A) und‖x‖? = 1 gelten. Dann gilt

ρ(A) = |λ| = |λ|‖x‖? = ‖λx‖? = ‖Ax‖? ≤ ‖A‖?‖x‖? = ‖A‖?.

I Wie kann man eine Losungsiteration

xk+1 = T xk + c

von Ax = b herleiten? Zerlegung: A = M + N, wobei75

Konvergente MatrizenSysteme mit M leicht zu losen sind:

Mxk+1 = b − Nxk

so c = M−1b und T = −M−1N. Die exakte Losung x? = A−1berfullt

Mx? = b − Nx?,

Subtraktion,(xk+1 − x?) = T (xk − x?)

Wenn gilt ‖T‖ < 1 fur eine induzierte (oder vertragliche)Matrixnorm, dann gilt

‖xk+1−x?‖ ≤ ‖T‖‖xk−x?‖, ‖xk+1−x?‖ ≤ ‖T‖k‖x0−x?‖ k→∞−→ 0

Def: T ∈ Rn×n ist konvergent wenn gilt

limk→∞

(T k )ij = 0, ∀1 ≤ i , j ≤ n

76

Konvergente MatrizenI Beispiele:

A =

[2 00 1

], Ak =

[2k 00 1

]−→×

[0 00 0

]B =

[ 12 00 1

3

], Bk =

[ 12k 00 1

3k

]−→

[0 00 0

]C =

[ 12 014

12

], Ck =

[ 12k 0k

2k+112k

]−→

[0 00 0

]D =

[−3

2 −165

65

52

], Dk =

[1 23 4

]−1 [( 9

10)k 00 ( 1

10)k

] [1 23 4

]−→

[0 00 0

]I ‖D‖1 = 5.7, ‖D‖2 = 4.4922, ‖D‖∞ = 4.7 und ρ(D) = 0.9.I Gibt es eine Matrixnorm ‖ · ‖M wobei ‖D‖M < 1?

Satz: Sei A ∈ Rn×n gegeben. ∀ε > 0 gibt es eine induzierteMatrixnorm ‖ · ‖Mε wobei gilt ρ(A) + ε > ‖A‖Mε .I Bemerkung: Diese Matrixnorm ‖ · ‖Mε hangt von A ab.

77

Zerlegungen fur IterationenSatz: Die folgenden sind aquivalent:

1. T ∈ Rn×n ist konvergent2. limk→∞ ‖T k‖ = 0 fur eine induzierte Matrixnorm3. limk→∞ ‖T k‖ = 0 fur alle induzierten Matrixnormen4. ρ(T ) < 15. limk→∞ T kx = 0, ∀x ∈ Rn

I In der Zerlegung A = M + N,

Mxk+1 = b − Nxk

welches M ist geeignet? Besonders wenn A diagonaldominant ist, ist die folgende Zerlegung naturlich:

A = D + L + U, M = D,N = L + U

wobei D eine diagonale Matrix, L ist eine streng untereDreiecksmatrix und U ist eine streng obere Dreiecksmatrix.

78

Jacobi Methode

I Die Jacobi Methode: TJ = −D−1(L + U)

xk+1 = D−1[b − (L + U)xk ]

= D−1[b − (L + U + D)= Axk + Dxk ]

= xk + D−1[b − Axk ]

= xk + D−1rk

wobei r = b − Ax das Residuum fur gegebenes x ist.I Die Jacobi Methode, komponentenweise,

∑ji

aijxkj =

n∑j=1

aijxj = bi , 1 ≤ i ≤ n

I Im folgenden Pseudo-Kode wird xi sequentiell in einerSchleife i = 1, . . . ,n berechnet:

79

Pseudo-Kode der Jacobi Methoded i = A ii, i=1,...,nx i = b i, i=1,...,nfor it=1,...,itmax

for i=1,...,nr i = b ifor j=1,...,n

r i = r i - A ij * x jenddx i = r i / d i

endfor i=1,...,n

x i = x i + dx iendif (||dx|| < tol * ||x||)

breakend

end

% Abbruchsparameter: itmax O(n)% tol = relative Fehlertoleranz

80

Jacobi MethodeI Wenn ein neuer Wert x i+dx i im Pseudo-Kode verfugbar

ist, warum nicht in x i sofort speichern? Somit waren dieZeilen im Pseudo-Kode

dx i = r i / d ix i = x i + dx i

neben einander innerhalb einer einzigen i-Schleife.I Da im Pseudo-Kode die Schleife i=1,...,n so lauft,

wurde die neue Iteration so aussehen:∑ji

aijx(k)j =

n∑j=1

aijxj = bi , 1 ≤ i ≤ n

oder

x (k+1)i =

1aii

bi −∑ji

aijx(k)j

I In Matrixnotation,

81

Gauß-Seidel MethodeI In Matrixnotation,

xk+1 = D−1[b − Lxk+1 − Uxk ]

oder die Gauß-Seidel Methode,

(D + L)xk+1 = b − Uxk

I Die Zerlegung A = M + N ist M = D + L, N = U undTGS = −M−1N = −(D + L)−1U.

I Umgeschrieben,

xk+1 = (D + L)−1︸︷︷︸Vorwarts Substitution

[b − Uxk ]

= (D + L)−1[b − (D + L + U)︸︷︷︸=A

xk + (D + L)xk ]

= xk + (D + L)−1[b − Axk ] = xk + (D + L)−1rk

82

Symmetrische Gauß-Seidel Methode

I Es gibt eine gewisse Asymmetrie bei der Gauß-SeidelMethode.

I Wir hatten die Schleife ebensogut i=n,...,1 laufenlassen konnen.

I Mit der Schleife in dieser entgegengesetzten Richtungbekommt man xk+1 = (D + U)−1[b − Lxk ].

I Mit beiden Schleifen bekommt man die SymmetrischeGauß-Seidel Methode,

xk+ 12 = (D + L)−1[b − Uxk ]

xk+1 = (D + U)−1[b − Lxk+ 12 ]

I Hausaufgabe: Finde die Iterationsmatrix TSGS fur dieSymmetrische Gauß-Seidel Methode.

83

Approximierte Inversen

Bemerkung: Die Zerlegung A = M + N kann soumgeschrieben werden:

xk+1 = −M−1Nxk + M−1b= [−M−1(M + N) + I]xk + M−1b= [I −M−1A]xk + M−1b

und so wird M−1 eine Approximierte Inverse benannt, da‖I −M−1A‖ klein sein sollte.

Hausaufgabe: Oben sind genannt worden: MJ = D undMGS = (D + L).I Leite die Approximierte Inverse MSGS fur die Symmetrische

Gauß-Seidel Methode.I Zeige, wenn A symmetrisch ist, ist MSGS symmetrisch,

wahrend MGS im allgemeinen nicht symmetrisch ist.

84

SOR MethodeI Gauß-Seidel noch einmal komponentenweise,

x (k+1)i =

1aii

bi −∑j

aijx(k)j

ω+x (k)i (1−ω)

Man konnte das alte und das neue so gewichten.I In Matrix-Notation,

xk+1 = D−1[b − Lxk+1 − Uxk ]ω + xk (1− ω)

Dxk+1 = [b − Lxk+1 − Uxk ]ω + Dxk (1− ω)

(D + ωL)xk+1 = [b − Uxk ]ω + Dxk (1− ω)

xk+1 = (D + ωL)−1[(1− ω)D − ωU]xk + ω(D + ωL)−1b

I Diese Methode heisst Successive Over-Relaxation (SOR)I 0 < ω < 1⇒ under relaxed, d.h. sie konvergiert wegen

Dampfung, wenn sie bei ω = 1 moglicherweise nichtkonvergiert.

I ω > 1⇒ over relaxed, d.h. sie konvergiert schneller wegenExtrapolation, wenn sie bei ω = 1 schon zuverlassigkonvergiert.

85

Konvergenz von Iterativen VerfahrenI Die Symmetrische Successive Over-Relaxation (SSOR)

Methode:xk+ 1

2 = (D + ωL)−1[(1− ω)D − ωU]xk + ω(D + ωL)−1bxk+1 = (D + ωU)−1[(1− ω)D − ωL]xk+ 1

2 + ω(D + ωU)−1b

Konvergenz

Satz: Wenn fur das System Ax? = b gilt x? = T x? + c, dannerfullt das iterative Verfahren xk+1 = T xk + c die Konvergenzlimk→∞ xk = x?, genau dann wenn ρ(T ) < 1.

Bemerkung: Wenn fur eine induzierte Norm ‖T‖ < 1 gezeigtwerden kann, folgt ρ(T ) < 1.

Satz: Wenn A streng diagonal dominant ist, giltlimk→∞ xk = x?, ∀x0 ∈ Rn, fur die Jacobi und Gauß-SeidelMethoden.86

Konvergenz von Iterativen VerfahrenI Wie wahlt man ω fur SOR aus? Sodass ρ(TSOR) < 1,

wobei TSOR = (D + ωL)−1[(1− ω)D − ωU].

Satz: Wenn fur aij = A ∈ Rn gilt aii 6= 0, 1 ≤ i ≤ n, folgtρ(TSOR) ≥ |ω − 1|.

Deswegen gilt ρ(TSOR) < 1 nur fur ω ∈ (0,2) gilt. Umgekehrt:

Satz (Ostrowski-Reich): Wenn A SPD ist und 0 < ω < 2 gilt,folgt fur SOR die Konvergenz limk→∞ xk = x?, ∀x0 ∈ Rn.

I Beispiel: Fur die Diskretisierung Au = [ µh2 D + I]u = v desRandwertproblems,

−µu′′ + u = v , x ∈ (0,1)u′ = 0, x ∈ 0,1

haben wir schon gezeigt, A ist streng diagonal dominantsowohl SPD. Deswegen konvergieren Jacobi undGauß-Seidel sowohl SOR, ∀ω ∈ (0,2).

87

Konjugierte GradientenDie Methode der Konjugierten Gradienten

Satz: Sei f ∈ C2(Ω). Wenn x? ∈ Ω erfullt ∇f (x?) = 0 und∇2f (x?) SPD ist, ist x? ein lokales Minimum fur f .

Explizite Bespiele:I Fur f (x) = bTx gelten ∇f (x) = b und ∇2f (x) = 0.

∂

∂xk

n∑i=1

bixi =n∑

i=1

bi∂xi

∂xk=

n∑i=1

biδik = bk

I Fur f (x) = xTMx , M ∈ Rn×n, gelten ∇f (x) = (M + MT)xund ∇2f (x) = M + MT.

∂

∂xk

n∑i=1

n∑j=1

xiMijxj =n∑

i=1

n∑j=1

(∂xi

∂xkMijxj + xiMij

∂xj

∂xk

)=

n∑i=1

n∑j=1

(δikMijxj + xiMijδjk

)=

n∑j=1

Mkjxj +n∑

i=1

Mikxi

88

Konjugierte GradientenSatz: Sei A ∈ Rn×n SPD. Dann lost x? das System Ax = bgenau dann wenn x? die Funktion g(x) = 1

2xTAx − bTxminimiert.

Beweis: (⇒) Aus Ax? = b folgt ∇g(x?) = Ax? − b = 0. Weil ASPD ist, hat g genau einen kritischen Punkt. Da ∇2g(x) = ASPD ∀x ∈ Rn ist, ist g uberall konvex, und x? ist ein lokalesMinimum. Aus der konvexen Analysis ist x? ein globalesMinimum. (⇐) Mit einem Minimum g(x?) von g ∈ C∞(Rn) giltnotwendigerweise ∇g(x?) = 0. Da A SPD ist, ∃!x? 3∇g(x?) = Ax? − b = 0.

Vorausgesetzt ist A SPD, und Ax = b soll gelost werden.

Entwicklung der Methode der Konjugierten Gradienten:

I Sei x0 eine Approximation der Losung x?.I Sei v1 eine gegebene Suchrichtung zur Minimierung von g.I Wo ist das Minimum von g entlang v1 von x0 weg?89

Konjugierte Gradienten

I Antwort durch die Ketten-Regel,

ddt

g(x0+tv1) = ∇g(x0+tv1)·v1 = [A(x0+tv1)−b]·v1 = 0

und das Minimum ist in x1 = x0 + t1v1 wobei

t1 =(−Ax0 + b) · v1

v1 · Av1 =r0 · v1

v1 · Av1

I Eigenschaft des Gradienten:∇g(x1) ⊥ x : g(x) = c, c = g(x1)

I Orthogonalitat −r1 · v1 = ∇g(x1) · v1

= [A(x0 + t1v1)− b] · v1 = 0bedeutet v1 ‖ x : g(x) = g(x1).

I Wie soll man die Suchrichtung v2 auswahlen?I Antwort: Minimiere g uber spanx1 + r1,x1 + v1.

90

Konjugierte Gradienten

I Die Funktion (t , σ) 7→ g(x1 + tr1 + σv1) wird minimiert:

ddσ

g(x1 + tr1 + σv1) = ∇g(x1 + tr1 + σv1) · v1 =σ=ts· · · =

[A(x1 + tr1 + tsv1)− b] · v1 = (Ax1 − b) · v1︸︷︷︸−r1·v1=0

+t [A (r1 + sv1)︸︷︷︸=v2

] · v1 = 0

I Nimm v2 = r1 + s1v1 wobei

s1 = − r1 · Av1

v1 · Av1

I Bemerke: v2 · Av1 = (r1 + s1v1) · Av1 = 0,d.h. v1 und v2 sind A-orthogonal oder A-konjugiert.

Iddt

g(x1 + tr1 + ts1v1) =ddt

g(x1 + tv2) = 0 fuhrt zu t = t2.

I Die obigen Formeln sind allgemein fur jede Iteration...91

Konjugierte GradientenZusammenfassung der Formeln:

x0 sei gegeben Fur k = 1,2, . . .v1 = r0 = b − Ax0 tk = rk−1 · vk/vk · Avk

xk = xk−1 + tkvk

rk · vk = 0 rk = b − Axk

vk+1 · Avk = 0 sk = −rk · Avk/vk · Avk

vk+1 = rk + skvk

I Vereinfachungen:

rk · vk = 0⇒ tk =rk−1 · [rk−1 + sk−1vk−1]

vk · Avk =rk−1 · rk−1

vk · Avk

Alsork = b − Axk = b − A(xk−1 + tkvk ) = rk−1 − tkAvk

⇒ rk · rk−1 = rk−1 · rk−1 − tk (vk − sk−1vk−1) · Avk = 0⇒ rk · rk = rk · rk−1 − tk rk · Avk = −tk rk · Avk

undsk = − rk · Avk

vk · Avk =rk · rk/tk

rk−1 · rk−1/tk =rk · rk

rk−1 · rk−192

Konjugierte GradientenDie Methode der Konjugierten Gradienten:

x0sei gegeben Fur k = 1,2, . . .v1 = r0 = b − Ax0 tk = rk−1 · rk−1/vk · Avk

xk = xk−1 + tkvk

rk = rk−1 − tkAvk

sk = rk · rk/rk−1 · rk−1

vk+1 = rk + skvk

I Hausaufgabe: Schreibe einen effizienten Pseude-Kode zurImplementierung dieser Methode.

Vorausgesetzt dass A SPD ist:

Satz: Es gelten rk · r l = 0, 0 ≤ l < k und vk · Av l = 0, 1 ≤ l < k .

Satz: Mit exakter Arithmetik gilt Axk = b fur ein k ≤ n.93

PrakonditionierungI Angenommen gilt C2 ≈ A aber Systeme Cy = c mit SPD

C sind leichter zu losen als Systeme Ax = b.I Das Problem Ax = b kann so umgeschrieben werden:

(C−1AC−1)(Cx) = C−1b

I Die Methode der Konjugierten Gradienten fuhrt schnellerzur Konvergenz, wenn sie auf das hergeleitete Systemangewendet wird, falls die SPD Matrix (C−1AC−1) besserkonditioniert ist als die Matrix A.

Def: Sei ‖ · ‖? eine induzierte Matrixnorm. Die Konditionszahleiner Matrix A bezuglich dieser Norm ist κ?(A) = ‖A‖?‖A−1‖?.

Bemerkung: 1 = ‖I‖? = ‖AA−1‖? ≤ ‖A‖?‖A−1‖? = κ?(A).Beispiel:

A =

[1 1− ε

1− ε 1

]ε∈(0,1)

σ(A) = ε,2− εκ2(A) = (2− ε)/ε

94

Fehler Abschatzungen

Satz: Eine gegebene Vektornorm und die entsprechendeinduzierte Matrixnorm sei mit ‖ · ‖ bezeichnet, und κ ist dieentsprechende Konditionszahl. Dann fur A ∈ Rn×n undb,x? = A−1b, x ,∈ Rn gilt:

‖x? − x‖‖x?‖

≤ κ(A)‖b − Ax‖‖b‖

Beweis: Aus Ax? = b folgt

‖b‖ = ‖Ax?‖ ≤ ‖A‖‖x?‖ ⇒ 1‖x?‖

≤ ‖A‖‖b‖

Mit e = x? − x und r = b − Ax = Ax? − Ax = Ae gilt

‖e‖ = ‖A−1Ae‖ ≤ ‖A−1‖‖Ae‖ = ‖A−1‖‖r‖

und daher‖e‖‖x?‖

≤ ‖A−1‖‖A‖ ‖r‖‖b‖

95

Fehler AbschatzungenI Ist (‖b−Axk‖ < τ · ‖b‖) ein gutes Abbruchskriterium einer

Iteration xk k→∞−→ x??I Wegen der obigen Abschatzung bleibt es offen, dass‖x? − xk‖/‖x?‖ sehr groß sein kann. Pessimistisch?

I Beispiel:

A =

[1 1− ε

1− ε 1

]ε∈(0,1)

σ(A) = ε,2− εκ2(A) = (2− ε)/ε

x? =

[cos(θ)sin(θ)

]θ=π

4

b = Ax? =

[cos(θ) + (1− ε) sin(θ)sin(θ) + (1− ε) cos(θ)

]‖x?‖2 = 1, ‖b‖22 = 1+(1−ε)2+4 cos(θ) sin(θ)︸︷︷︸

2 sin(2θ)=2

(1−ε) = (2−ε)2

x = x? +zδ, z =

1√2

[1−1

], Az = εz, ‖z‖2 = 1

‖x? − x‖2 =‖z‖2δ

=1δ, ‖b − Ax‖2 =

‖Az‖2δ

=ε

δ96

Fehler Abschatzungen

I Zusammenfassung: mit δ =√ε,

‖x? − x‖2‖x?‖2︸︷︷︸

1/δδ→0−→∞

= κ2(A)︸︷︷︸(2−δ2)/δ2δ→0−→∞

‖b − Ax‖2‖b‖2︸︷︷︸

δ/(2−δ2)δ→0−→0

d.h. es kann doch sein, dass das Residuum sehr kleinwird, wahrend der Fehler groß bleibt.

I Abbruchskriterium? Zu empfehlen ist:

(‖xk+1 − xk‖ < τ · ‖xk+1‖)

wobei die Toleranz τ rein von der Genauigkeit derFließkommazahl-Darstellung gewahlt werden kann.

97

Eigenwerte und EigenvektorenZuerst eine gewisse Infrastruktur:

Satz: Hat A ∈ Rn×n die verschiedenen Eigenwerte λiki=1,(k ≤ n), mit den entsprechenden Eigenvektoren x iki=1, dannsind diese Eigenvektoren linear unabhangig. Wenn k = n gilt,dann bilden die Eigenvektoren eine Basis fur Rn.

Def: P ∈ Rn×n ist orthogonal wenn P−1 = PT gilt.

Def: A,B ∈ Rn×n sind ahnlich wenn ∃S ∈ Rn×n 3 A = S−1BS .

Satz: A = S−1BS ⇒ σ(A) = σ(B) undAx = λx ⇒ BSx = λSx .

Satz: Ist A ∈ Rn×n symmetrisch und D = diag(λini=1),λini=1 = σ(A), dann ∃P ∈ Rn×n orthogonal wobei D = PTAP.

Satz: Ist A ∈ Rn×n symmetrisch, gilt σ(A) ⊂ R.

Satz: Ist A ∈ Rn×n symmetrisch, dann bilden die Eigenvektoreneine orthonormale Basis fur Rn.

Satz: A ist SPD⇔ A ist symmetrisch und σ(A) ⊂ R+.98

VektoriterationZur Bestimmung eines dominanten EigenwertsI Der Plan: Fur A ∈ Rn×n und einen Startvektor x0, soll Akx0

einigermaßen zu einem Eigenvektor konvergieren, derzum dominanten Eigenwert gehort.

I Angenommen hat A ∈ Rn×n die Eigenwerte λini=1 mitden entsprechenden Eigenvektoren v ini=1, die eine Basisfur Rn bilden.

I Zusatzlich gelten

|λ1| > |λ2| ≥ · · · ≥ |λn| ≥ 0.

d.h. λ1 ist dominant.I Ein beliebiger Startvektor fur die Vektoriteration kann so

dargestellt werden:x0 =

n∑i=1

αiv i

I Angenommen gilt α1 6= 0 und x0 ist normalisiert sodass‖x0‖∞ = 1 gilt.

99

VektoriterationI Die Vektoriteration: Fur k = 1,2, . . .

yk = Axk−1

xk = yk/ykpk

µk = ykpk−1

wobei|yk

pk| = ‖yk‖∞

d.h. ‖xk‖∞ = 1

I Behauptungen:

µkk→∞−→ λ1 und skxk k→∞−→ v1

‖v1‖∞wobei s2

k = 1

I Zuerst eine Darstellung von xk ,

xk =yk

ykpk

=Axk−1

ykpk

=Ayk−1

ykpk

yk−1pk−1

= · · · =Ak−1y1

ykpk· · · y1

p1

=Akx0∏ki=1 y i

pi

I Dann eine Darstellung fur µk ,

µk =yk

pk−1

xk−1pk−1

(= 1)=

(Axk−1)pk−1

xk−1pk−1

=(Akx0)pk−1/

∏k−1i=1 y i

pi

(Ak−1x0)pk−1/∏k−1

i=1 y ipi

=(Akx0)pk−1

(Ak−1x0)pk−1

100

VektoriterationI Weil λ1 dominant ist, gilt

µk =

λk1

[α1v1 +

n∑i=2

(λi

λ1

)k

↓0

v i

]pk−1

λk−11

[α1v1 +

n∑i=2

(λi

λ1

)k−1

↓0

v i

]pk−1

k→∞−→ λ1

I Wegen der obigen Darstellung von xk gilt

yk = ykpk

xk = ykpk

Akx0∏ki=1 y i

pi

=λk

1∏k−1i=1 y i

pi︸︷︷︸Q=1/(α1νk )

[α1v1 +

n∑i=2

αi

(λi

λ1

)k

↓0

v i

]

Mit |v1p | = ‖v1‖∞ gilt

|ykpk|

|ykp |

=‖yk‖∞|yk

p |=‖νkyk‖∞|νkyk

p |k→∞−→ ‖v1‖∞

|v1p |

= 1

101

Inverse VektoriterationI Also konvergieren die Vektoren xk wie behauptet:

ykpk

ykp

∣∣∣∣∣→±1

xk =yk

ykp

=Axk−1

(Axk−1)p=

Akx0/∏k−1

i=1 y ipi

(Akx0)p/∏k−1

i=1 y ipi

=Akx0

(Akx0)pk

=

λk1

[α1v1 +

n∑i=2

αi

(λi

λ1

)k

↓0

v i

]

λk1

[α1v1 +

n∑i=2

αi

(λi

λ1

)k

↓0

v i

]p

k→∞−→ v1

v1p

= sign(v1p )

v1

‖v1‖∞

Inverse VektoriterationI Um λk zu berechnen, wird die Vektoriteration auf

(A− qI)−1 angewendet, wobei |λk − q| < |λi − q|,∀i 6= k .I Hausaufgabe: Es gilt µk

k→∞−→ (λk − q)−1 fur die Iteration,(A− qI)yk = xk−1

xk = yk/ykpk

µk = ykpk−1

wobei|yk

pk| = ‖yk‖∞

d.h. ‖xk‖∞ = 1

102

Berechnung aller EigenwerteI Die Vektoriteration und die Inverse Vektoriteration sind nur

zur Berechnung bestimmter Eigenwerte geeignet.I Nun suchen wir Verfahren, die alle Eigenwerte und ihre

entsprechenden Eigenvektoren liefern.I Nach dem Abel-Ruffini Satz gibt es fur n ≥ 5 keine

allgemeine endliche Formel mit Arithmetik und Wurzeln zurBestimmung der Nullstellen des charakteristischenPolynoms.

I Daher mussen Eigenwerte iterativ berechnet werden.I Der Plan vom QR-Algorithmus: A(1) = A (SPD)

Fur k = 1,2 . . . A(k) = Q(k)R(k), A(k+1) = R(k)Q(k) wobei

Q(k) orthogonal ist und R(k) eine obere Dreiecksmatrix ist.I Aber diese Iteration ist viel billiger, wenn zuerst so

transformiert wird: A = PAPT, wobei P orthogonal und ASPD tridiagonal sind.

103

Die Householder MethodeI Zur Bestimmung der Transformation A = PAPT, A SPD,

wobei P orthogonal und A SPD tridiagonal sind:I Mit P = P(n−2) · · ·P(1), A(1) = A, soll gelten

A(2) = P(1)A(1)P(1) =

0 · · · 0 + − − +0 | |... | |0 + − − +

I und im nachsten Schritt

A(3) = P(2)A(2)P(2) =

0 · · · 0 0 · · ·0 + − +... 0 | |

0... + − +

I usw bis A = A(n−1).

104

Householder TransformationenDef: Sei w ∈ Rn mit w · w = 1. Dann heisst P = I − 2wwT

eine Householder Transformation.

Satz: Eine Householder Transformation ist symmetrisch undorthogonal.

PTP = P2 = (I−2wwT)(I−2wwT) = I−4wwT+4wwTwwT = I

Satz: Seien 0 6= u ∈ Rn und θ = 12‖u‖

22. Dann ist

P = I − uTu/θ eine Householder Transformation.

w = u/‖u‖2 ⇒ P = I − 2wTw

Satz: Seien x ∈ Rn und σ = ±‖x‖2. Angenommen x 6= −σe(1).Seien u = x + σe(1) und θ = 1

2‖u‖22. Dann ist P = I − uTu/θ

eine Householder Transformation und es gilt Px = −σe(1).

Beweis: x 6= −σe(1) ⇒ u = x + σe(1) 6= 0. Nach dem obigenSatz ist P eine Householder Transformation. Dann gilt105

Householder Transformationenθ = 1

2‖u‖22 = 1

2(x + σe(1))T(x + σe(1)

)

= 12(xTx + 2σxTe(1)

+ σ2) (σ2 = ‖x‖22)

= 12(σ2 + 2σxTe(1)

+ σ2) = σ2 + σx1

Daher gilt(I − uuT/θ)x = x − uuTx/θ

= x − (x + σe(1))

(x + σe(1))Tx

σ2 + σx1

= x − (x + σe(1))

(xTx + σx1)

σ2 + σx1= −σe(1)

Bemerkung: Um Ausloschung zu vermeiden, soll σ mit demgleichen Vorzeichen wie x1 ausgewahlt werden. P wird auchnicht verandert, wenn x zur Sicherheit skaliert wird.

Algorithmus: Hausaufgabe: θ = 12‖u‖

22

v = x/‖x‖∞ θ = ρu1ρ = s(v1)‖v‖2 σ = ρ‖x‖∞u = v + ρe(1) P = I − uTu/θ

s(t) =

1, t ≥ 0−1, sonst

Px = −σe(1)106

Transformation auf Tridiagonale FormI Zur Bestimmung von P(1) sei

a(1) =

a(1)

11a(1)

21...

a(1)n1

=

a(1)

11

a(1)

P = I − uTu/θ, u ∈ Rn−1

wobei Pa(1) = −σe(1) ∈ Rn−1

I Dann

P = I − uTu/θ, u =

[0u

]∈ Rn

erfullt

Pa(1) =

1 0 · · · 00 + − +... | P |0 + − +

a(1)11

a(1)

=

a(1)

11−σ0...

107

Transformation auf Tridiagonale FormI Zur Bestimmung von P(2) sei

a =

a(2)

12a(2)

22...

a(2)n2

=

a(2)

12a(2)

22

a(2)

P = I − uTu/θ, u ∈ Rn−2

wobei Pa(2) = −σe(1) ∈ Rn−2

I Dann

P = I − uTu/θ, u =

00u

∈ Rn

erfullt

Pa(2) =

1 0 0 · · · 00 1 0 · · · 00 0 + − +...

... | P |0 0 + − +

a(2)

12a(2)

22

a(2)

=

a(2)

12a(2)

22−σ0...

108

Hessenberg MatrizenI Dann sehen die Matrizen P(i) · · ·P(1)A so aus:

P(1)A =

+ − − + | |0 | |... | |0 + − − +

P(2)P(1)A =

· · · · · · + − +0 | |... 0 | |

0... + − +

I Und P(n−2) · · ·P(1)A ist eine obere Hessenberg Matrix:

P(n−2) · · ·P(1)A =

· · · · · · · · ·

. . .

...

0. . .

. . ....

.... . .

. . .. . .

0 · · · 0

109

Transformation auf Tridiagonale FormI Wenn A(1) = A SPD ist, gilt A(2) = P(1)A(1)P(1) =

1 0 · · · 00 + − +... | P |0 + − +

a(1)11 a(1)T

+ − +

a(1) | A |+ − +

1 0 · · · 00 + − +... | P |0 + − +

=

a(1)

11 −σ 0 · · ·−σ + − +

0 | PAP |... + − +

I P(2) wird von A(2) = P(1)A(1)P(1) wie vorher bestimmt, und

ahnlich ist A(3) = P(2)A(2)P(2) tridiagonal.I Dann mit

P = P(n−2) · · ·P(1)

ist PAPT tridiagonal.110

Pseudo-Kode zur Transformation auf Hessenbergoder Tridiagonale Form

H = A, Q = I, ASym = (A = A’)for k=1,...,n-2

h i = H k+i,k, i=1,...,n-k(theta,sigma,u) = Householder(h) % (I-uu’/theta)h=-sigma e1for l=k,...,n % H <- (I-uu’/theta)H, l=k Verbesserung?

sm = sum(H k+i,l * u i: i=1,...,n-k)H k+i,l = H k+i,l - sm*u i/theta, i=1,...,n-k

endfor l=1,...,n % Q <- Q(I-uu’/theta)

sm = sum(Q l,k+i * u i: i=1,...,n-k)Q l,k+i = Q l,k+i - sm*u i/theta, i=1,...,n-k

endIf ASym

for l=k,...,n % H <- H(I-uu’/theta)sm = sum(H l,k+i * u i: i=1,...,n-k)H l,k+i = H l,k+i - sm*u i/theta, i=1,...,n-k

endend

end% ASym = false => A=Q*H, ASym = true => A=Q*H*Q’111

QR AlgorithmusI Fur die Iteration, A(k) = Q(k)R(k), A(k+1) = R(k)Q(k), gilt

A(k+1) = R(k)Q(k) = Q(k)TA(k)Q(k) = · · · =

Q(k)T · · ·Q(1)T︸︷︷︸QT

A Q(1) · · ·Q(k)︸︷︷︸Q

also sind die Matrizen A(k) alle ahnlich zu A.I Unter gewissen Bedingungen gilt

A(k) k−→∞−→ Λ = diag(σ(A)).I Der erste Bedarf ist, ein Algorithmus zur QR-Zerlegung.I Householder Transformationen funktionieren: Fur Q(1)

a(1) =

a(1)

11a(1)

21...

a(1)n1

Q(1) = I − uTu/θ, u ∈ Rn

wobei Q(1)a(1) = −σe(1) ∈ Rn

112

QR AlgorithmusI Fur Q(2)

a(2) =

a(2)

12a(2)

22...

a(2)n2

=

a(2)

12

a(2)

Q(2) = I − uTu/θ, u ∈ Rn−1

wobei Q(2)a(2) = −σe(1) ∈ Rn−1

Q(2) = I − uTu/θ, u =

[0u

]I Es gelten

Q(1)A=

−σ1 + − − +

0 | |... | |0 + − − +

Q(2)Q(1)A=

−σ1 + − +

0 −σ2 | |... 0 | |

0... + − +

I Nach n − 1 solchen Transformationen ergibt sich eine

obere Dreiecksmatrix,

QTA = Q(n−1) · · ·Q(1)A = R113

Givens TransformationenI Aber wenn A schon tridiagonal ist, kann Q billiger

konstruiert werden:

A =

a1 b1 0 · · · 0

b1 a2 b2. . .

...

0. . .

. . .. . . 0

.... . . bn−2 an−1 bn−1

0 · · · 0 bn−1 an

Def: Eine Givens Transformation G = Gij(θ) hat die Form

GGT = I, G =

1. . . 1qii qij

1. . . 1qji qjj

1. . . 1

qii = cos(θ) = qjjqij = sin(θ) = −qjifur irgendwelche

i , j und θ

114

Givens TransformationenI Fur x1 = a1, y1 = b1, soll θ so ausgewahlt werden, dass

cos(θ) sin(θ)− sin(θ) cos(θ)

11. . . 1

x1 y1b1 a2 b2

. . .. . .

. . .

bn−1 an

=

g1 d1 e10 x2 y2

b2 a3 b3. . .

. . .. . .

bn−1 an

I Mit s1 = sin(θ1), c1 = cos(θ1),

−s1x1 + c1b1 = 0s2

1 + c21 = 1

s1 = b1/√

b21 + x2

1

c1 = x1/√

b21 + x2

1

I Hausaufgabe: Fur g1 > 0 ist Q eindeutig bestimmt.115

Die QR-ZerlegungI Initialisierung: R = A, Q = I.I Fur k = 1, . . . ,n − 1,

Q = Gk ,k+1(θk )Q, R = Gk ,k+1(θk )R

wobei sin(θk ) = bk/√

b2k + x2

k , cos(θk ) = xk/√

b2k + x2

k .I Es gilt

1. . . 1

cos(θk ) sin(θk )− sin(θk ) cos(θk )

1. . . 1

. . .. . .

. . .

gk−1 dk−1 ek−1

0 xk yk 0bk ak+1 bk+1

. . .. . .

. . .

. . .. . .

. . .

gk dk ek

0 xk+1 yk+1 0bk+1 ak+2 bk+2

. . .. . .

. . .

I Am Ende setze Q = QT und es gilt A = QR. Eindeutig?

116

QR AlgorithmusSatz: Sei A ∈ Rn×n mit Rang(A) = n. Dann gibt es eineZerlegung A = QR, wobei Q orthogonal ist und R eine obereDreiecksmatrix ist. Weiters ist die Zerlegung eindeutig, wenndie diagonalen Eintrage von R positiv sind.

Bemerkung: Fur die Konvergenz des QR-Algorithmus zurBestimmung der Eigenraume von A ist diese eindeutigeZerlegung nicht notwendig.

I Nun haben wir eine QR-Zerlegung A(k) = Q(k)R(k).I Zu zeigen ist, dass

A(k+1) = R(k)Q(k) = Q(k)TA(k)Q(k)

wieder tridiagonal ist.I Zuerst wird gezeigt, Q(k) ist obere Hessenberg.

117

Tridiagonale Form Bleibt ErhaltenI Bemerke:

1c2 s2−s2 c2

1. . . 1

c1 s1−s1 c1

11. . . 1

=

c1 s1

11. . . 1

I Induktiv:

1. . . 1

ck sk−sk ck

1. . . 1

. . .

. . ....

. . .. . .

. . ....

. . .. . .

· · · · · · 1

1. . . 1

=

. . .

. . ....

. . .. . .

. . ....

. . .. . .

.... . .

· · · · · · · · · 1. . . 1

I Also ist Gk ,k+1(θk ) · · ·G1,2(θ1) untere Hessenberg.118

Tridiagonale Form Bleibt ErhaltenI Daher haben wir folgende Struktur fur A(k+1) = R(k)Q(k) =

. . .

. . .. . .

. . .. . . . . .

· · · · · ·

. . .

. . ....

. . .. . .

. . ....

. . .. . .

und A(k+1) ist

eine obereHessenberg

Matrix. =

· · · · · ·

. . .

. . ....

. . .. . .

. . ....

. . .. . .

I Aber A(k+1) = Q(k)TA(k)Q(k) ist symmetrisch und

deswegen tridiagonal!119

Konvergenz des QR-AlgorithmusSatz: Sei A ∈ Rn×n symmetrisch mit Eigenwerten λini=1, dieerfullen |λ1| > |λ2| > · · · > |λn| > 0. Seien v ini=1 dieentsprechenden Eigenvektoren, die eine orthonormale Basisbilden. Fur die Matrix Q = [v1, . . . ,vn] soll es eineLU-Zerlegung QT = LU geben, ohne dass QT mit einerPermutationsmatrix multipliziert werden muss. Dann mitA(1) = A und einer (beliebigen) QR-Zerlegung Ak = QkRk ,konvergiert die Folge A(k+1) = R(k)Q(k) zur diagonalen MatrixΛ = diag(λini=1).

Bemerkung: Mit A(k) = a(k)ij gilt

|a(k)i+1,i | = O

(∣∣∣∣λi+1

λi

∣∣∣∣k)

oder |a(k)i+1,i | = O

(∣∣∣∣λi+1 − qλi − q

∣∣∣∣k)

wenn man eine Verschiebung macht,

A(k) − qI = Q(k)R(k), A(k+1) = R(k)Q(k) + qI120

Pseudo-Kode fur den QR-Algorithmusa i = H ii, i=1,...,n, b i = H i,i+1 = H i+1,i, i=1,...,n-1for k=1,...,kmax

x 1 = a 1, y 1 = b 1for j=1,...,n-1 % compute R, Q (storing c’s & s’s)

g j =√x jˆ2 + b jˆ2

c j = x j / g j, s j = b j / g jd j = c j * y j + s j * a j+1x j+1 = c j * a j+1 - s j * y jif j 6= (n-1)

e j = s j * b j+1, y j+1 = c j * b j+1end

endg n = x na 1 = s 1 * d 1 + c 1 * g 1, b 1 = s 1 * g 2for j=2,...,n-1 % compute R * Q, overwriting H

a j = s j * d j + c j * c j-1 * g jb j = s j * g j+1

enda n = c n-1 * g nif (||b|| < tol * ||a||) then break

end% a i, i=1,...,n ≈ Eigenwerte von H121

AusgleichsproblemeBeispiel (Lineare Regression): Gegeben sind (xi , yi)m−1

i=0 .Was ist die beste Gerade durch diese Daten?

I Es soll minimiert werden:

E(a,b) =m−1∑i=0

[yi − (a + bxi)]2 V = V (x) = [e,x ]

= ‖y − Vc‖22 c = 〈a,b〉T (ei = 1)

= cTV TVc − 2cTV Ty + yTy

I E(c?) = min⇒ 0 = ∇cE(c?) = 2V TVc? − 2V TyI Hausaufgabe: xim−1

i=0 verschieden⇒ V TV ist SPD.I Wenn c? die normalen Gleichungen lost,

V TVc = V Ty

folgt E(c?) = min.122

Lineare und Polynome RegressionI Hausaufgabe: Zeige, fur lineare Regression ist die Losung

explizit so gegeben:

a? = y−b?x, b? =xy − x yx2 − x2 wobei z.B. xy =

1m

m−1∑i=0

xiyi

I Polynome Regression: (V=Vandermonde Matrix)

E(c) =m−1∑i=0

[yi − P(xi ; c)]2, P(x ; c) =n−1∑k=0

ckxk

= ‖y − Vc‖22 V = [e,x , . . . ,xn−1]

= cTV TVc − 2cTV Ty + yTy c = 〈c0, . . . , cn−1〉T

I Losung gegeben durch die normalen Gleichungen,

V TVc = V Ty

aber V TV kann fur hohes n schlecht konditioniert sein.123

Singularwert ZerlegungI Ein alternativer Ansatz ist durch die Singularwert

Zerlegung gegeben.

Satz (SWZ): Fur A ∈ Rm×n ∃ orthogonale Matrizen U ∈ Rm×m,V ∈ Rn×n sodass UTAV = Σ = diagσ1, . . . , σp ∈ Rm×n wobeip = minm,n und σ1 ≥ σ2 ≥ · · · ≥ σp ≥ 0.

I Also UTAV = Σ⇒ A = UUTAVV T = UΣV T.

Def: A = UΣV T ist die Singularwert Zerlegung (SWZ) von A,und σipi=1 sind die Singularwerte von A.

Satz: Sei A = UΣV T die SWZ fur A ∈ Rm×n, wobei m > n undσ1 ≥ · · · ≥ σr > σr+1 = · · · = σn = 0 gelten. Dann wird‖Ax − b‖2 durch x? = V Σ†UTb ∈ Rn minimiert, wobeiΣ† = diagσ−1

1 , . . . , σ−1r ,0, . . . ,0 ∈ Rn×m. Wenn x ∈ Rn auch

minimierend ist und x 6= x? gilt, dann gilt ‖x?‖2 < ‖x‖2.124

Eigenschaften der SWZ

Satz (Eigenschaften der SWZ): Sei A = UΣV T die SWZ furA ∈ Rm×n mit p = minm,n. Es gelten:I U = u1, . . . ,um, V = v1, . . . ,vn ⇒

Av i = σiui , ATui = σiv i , i = 1, . . . ,p.I σ1 ≥ · · · ≥ σr > σr+1 = · · · = σn = 0⇒N (A) = spanv ini=r+1, R(A) = spanuiri=1.

I ‖A‖2 = σ1.I ‖A‖2F = σ2

1 + · · ·+ σ2p.

I κ?(A) = max‖x‖?=1

‖Ax‖?/ min‖x‖?=1

‖Ax‖?m≥n⇒ κ2(A) = σ1/σn.

I ATA = V ΣTΣV T und AAT = UΣΣTUT.I σ2

i pi=1 sind die eindeutigen Nullstellen von det(λI − ATA)

fur p = n oder die eindeutigen Nullstellen vondet(λI − AAT) fur p = m.

125

PseudoinverseDef: Sei A = UΣV T die SWZ fur A ∈ Rm×n. Dann ist istA† = V Σ†UT die Pseudo-Inverse von A.

Satz: Fur A ∈ Rm×n, N (A) = 0⇒ A† = (ATA)−1AT.

Beweis: Es gilt xTATAx = ‖Ax‖22 ≥ 0 und zusatzlich folgt‖Ax‖22 = 0⇔ x = 0 sowohl m ≥ n aus N (A) = 0. Daher istATA SPD. Aus ATA = (UΣV T)T(UΣV T) = V ΣTUTUΣV T =V ΣTΣV T folgt, dass die positiven Eigenwerte durch 0 < ΣTΣ =diagσ2

1, . . . , σ2n ∈ Rn×n gegeben sind. Es folgt (ATA)−1AT =

V [ΣTΣ]−1V TV ΣTUT = V [ΣTΣ]−1ΣTUT = V Σ†UT = A†.

Satz: Fur A ∈ Rm×n und orthogonale Matrizen P ∈ Rm×m undQ ∈ Rn×n haben A und PAQ die gleichen Singularwerte.

Beweis: Fur m > n sind die Eigenwerte σ2i

ni=1 von ATA =

V ΣTΣV T die Eigenwerte von (PAQ)T(PAQ) = (QTV )ΣTΣ(V TQ).Fur n > m sind die Eigenwerte σ2

i mi=1 von AAT = UΣΣTUT die

Eigenwerte von (PAQ)T(PAQ) = (QTU)ΣΣT(UTQ).126

Konditionierung und die SWZI Fur A ∈ Rm×n, m > n, mit SWZ, A = UTΣV ,

Σ =

(Σ0

)∈ Rm×n, Σ = diagσini=1, σi > 0

sei die Projektion P = ΣΣ† von Rm in Rn definiert.I Fur ein gegebenes b sei die Zerlegung definiert:

Ub = 〈c,d〉T, c ∈ Rr , d ∈ Rm−r

PUb = 〈c,0〉T, (I − P)Ub = 〈0,d〉TI Das Ausgleichsproblem kann so umformuliert werden:

‖Ax − b‖22 = ‖UTΣVx − b‖22 = ‖ΣVx − Ub‖22= ‖P(ΣVx − Ub)‖22 + ‖(I − P)(ΣVx − Ub)‖22

+2〈P(· · · ), (I − P)(· · · )〉=0 = ‖ΣVx − c‖22 + ‖d‖22I Da κ2(Σ) = κ2(ΣV ) gilt, erfullen die Losungen der

Systeme ΣVx = c und ΣV x = c die Fehlerabschatzung‖x − x‖2‖x‖2

≤ κ2(Σ)‖c − c‖2‖c‖2

I Vergleiche mit κ2(ATA) = κ2(Σ)2 fur die normalenGleichungen.

127

Berechnung der SWZBerechnung der SWZ (zur Losung eines Ausgleichsproblems)I Da die Singularwerte durch die Eigenwerte von ATA

(m > n) gegeben sind, konnten wir den QR-Algorithmusauf ATA anwenden, um ΣTΣ zu bestimmen.

I Der folgende Algorithmus mit bidiagonalen Matrizen isteffizienter.

I Der Plan:I Vereinfachung PAQ = 〈BT,0〉T, A ∈ Rm×n, P ∈ Rm×m und

Q ∈ Rn×n orthogonal und B ∈ Rn×n obere bidiagonal.I Mit B(1) = B ist B(l)TB(l) tridiagonal. (warum?)I Setze B(l)TB(l) = Q(l)R(l). (QR-Zerlegung)I Dann ist B(l)Q(l) untere bidiagonal. (warum?)I Berechne P(l) sodass B(l+1) = P(l)B(l)Q(l)

obere bidiagonal ist. (wie?)I Zeige, es gilt B(l+1)TB(l+1) = R(l)Q(l)!I Im Endeffekt wird ein Schritt eines QR-Algorithmus effizient

durchgefuhrt, und es folgt |B(l)| → diagσini=1, l →∞.

128

Vereinfachung auf Bidiagonale FormI Fur die Vereinfachung PAQ = (B,0)T setze

A = A(1) = (d (1)|C(1)), d (1) ∈ Rm×1, C(1) ∈ Rm×(n−1)

I Fur k = 2, . . . ,n − 1,A(k) = P(k−1)A(k−1)Q(k−1) =

(B(k) 0

a(k) 0

0 d (k) C(k)

)wobei

k > 2

B(k) ∈ R(k−1)×(k−1) a(k) ∈ R1

C(k) ∈ R(m−k+1)×(n−k) d (k) ∈ R(m−k+1)×1

I Im nachsten Schritt,

P(k) =

(I(k−1) 0

0 P(k)

)Q(k) =

(I(k) 00 Q(k)

)wobei durch Householder Transformationen

P(k)d (k) = −σe(1) ∈ Rm−k+1 c(k) = 1.Zeile von P(k)C(k)

Q(k)Tc(k)T = −τ e(1) ∈ Rn−k a(k+1) = −τ129

Pseudo-Kode zur Transformation auf BidiagonaleForm

d i = A i,1, C i,j = A i,j+1, i=1,...,m, j=1,...,n-1P = I (m×m), Q = I (n×n)for k=1,...,n

% d is (m-k+1)×1(u,sigma,theta) = Householder(d) % (I-uu’/theta)d=-sigma e1B k,k = -sigma % B growsif (k>1) then B k-1,k = afor l=1,...,n-k % C <- (I-uu’/theta)C

sm = sum(C i,l * u i: i=1,...,m-k+1)C i,l = C i,l - sm*u i/theta, i=1,...,m-k+1

endfor l=1,...,m % P <- (I-uu’/theta)P

sm = sum(P k-1+i,l * u i: i=1,...,m-k+1)P k-1+i,l = P k-1+i,l - sm*u i/theta, i=1,...,m-k+1

endif (k < n-1)

c j = C 1,j, j=1,...,n-k% c is (n-k)×1

130

Pseudo-Kode zur Transformation auf BidiagonaleForm

(u,tau,theta) = Householder(c) % (I-uu’/theta)c=-tau e1a = -taufor l=2,...,m-k+1 % C <- C(I-uu’/theta)

sm = sum(C l,i * u i: i=1,...,n-k)C l,i = C l,i - sm*u i/theta, i=1,...,n-k

endfor l=1,...,n % Q <- Q(I-uu’/theta)

sm = sum(Q l,k+i * u i: i=1,...,n-k)Q l,k+i = Q l,k+i - sm*u i/theta, i=1,...,n-k

endendif (k = n-1) then a = C 1,1if (k < n) % d & C shrink

d i = C i+1,1, i=1,...,m-kC i,j = C i+1,j+1, i=1,...,m-k, j=1,...,n-k-1

endend

% P*A*Q = (BˆT,0)ˆT, B n×n

131

Diagonalisierende Iteration

I Bei k = n − 1,A(n) = P(n−1)A(n−1)

und Q(n−1) ist nicht notwendig, weil C(n) nur eine einzigeSpalte hat, und c(n−1) ∈ R1 ist schon a(n).

I Bei k = n,A(n+1) = P(n)A(n)

und Q(n) ist nicht notwendig, weil C(n) nicht existiert. Dad (n) ∈ R(m−n+1)×1, ist P(n) aber notwendig.

I Dann ist B = B(n+1) bidiagonal.I Fur die Diagonalisierung von B setze B(1) = B.I Sei die obere Hessenberg Matrix

Q(l) = G(1,2)(θ1) · · ·G(n−1,n)(θn−1)

durch Givens Transformationen definiert, damit

132

Diagonalisierende Iteration

B(l)Q(l) =

. . .. . .

. . .

Q(l) =

0

. . .

. . .

. . .. . . 0

eine untere bidiagonale Matrix ist.

I Dann durch Multiplikation links mit B(l)T

R(l)T = B(l)T(B(l)Q(l)) =

. . .

. . .. . .

. . .

. . .. . .

ist R(l) eine obere tridiagonale Matrix.

133

Folge der Bidiagonalen MatrizenI Daher ist

B(l)TB(l) = Q(l)R(l)

eine QR Zerlegung von der tridiagonale Matrix B(l)TB(l),

. . .

. . .

. . .. . .

= B(l)TB(l) = Q(l)R(l) =

· · · · · ·

. . .

. . ....

. . .. . .

. . ....

. . .. . .

. . .

. . .. . .

. . .. . . . . .

134

Folge der Bidiagonalen MatrizenI Sei die untere Hessenberg Matrix

P(l) = G(n−1,n)(θn−1) · · ·G(1,2)(θ1)

durch Givens Transformationen definiert, damit

B(l+1) =P(l)(B(l)Q(l))=P(l)

. . .

. . .. . .

=

0. . .

. . .

. . .. . . 0

eine obere bidiagonale Matrix ist.

I Schliesslich,

B(l+1)TB(l+1) = (P(l)B(l)Q(l))T (P(l)B(l)Q(l))

= (Q(l)TB(l)TP(l)T) (P(l)B(l)Q(l))

= Q(l)T(B(l)TB(l))Q(l) = Q(l)T(Q(l)R(l))Q(l)

= R(l)Q(l)

und ein Schritt des QR Verfahrens fur BTB ist vollstandig.135

Pseudo-Kode fur den QR Algorithmus zurSingularwert-Zerlegung

Q = I, P = I, Bk = Bfor k=1,...,kmax

a i = Bk i,i, i=1,...,nb i = Bk i,i+1, i=1,...,n-1x 1 = a 1for j=1...,(n-1) % compute Qk by storing c’s & s’s


c j = x j / g js j = -b j / g jd j = -s j * a j+1x j+1 = c j * a j+1

endg n = x n % Bk*Qk has (lower) diags d,g

136


Qk = Ifor j=1,...,(n-1)

q1 i = Qk i,j * c j - Qk i,j+1 * s j, i=1,...,nq2 i = Qk i,j+1 * c j + Qk i,j * s j, i=1,...,nQk i,j = q1 i, Qk i,j+1 = q2 i, i=1,...,n

end % (Bk*Qk)’*Bk = RkQ = Q*Qka i = (Bk*Qk) i,i = g i, i=1,...,nb i = (Bk*Qk) i+1,i = d i, i=1,...,n-1x 1 = a 1for j=1...,(n-1)


c j = x j / g js j = b j / g jd j = s j * a j+1x j+1 = c j * a j+1

endg n = x n % Tk has (upper) diags g,d

137

Pk = Ifor i=1,...,(n-1)

p1 j = Pk i,j * c i + Pk i+1,j * s i, j=1,...,np2 j = Pk i+1,j * c i - Pk i,j * s i, j=1,...,nPk i,j = p1 j, Pk i+1,j = p2 j, j=1,...,n

end % Pk*(Bk*Qk) = TkP = Pk*PBk i,i = g i, i=1,...,nBk i,i+1 = d i, i=1,...,(n-1) % Bk’*Bk = Rk*Qkif (||d|| < tol * ||g||) then break

end% |g i|, i=1,...,n ≈ Singularwerte von BSn i,i = sign(g i), i=1,...,nS i,i = |g i|, i=1,...,nP=Sn*P% S = P*B*Q

138

InterpolationI Fur Polynome Regression haben wir das

Ausgleichsproblem gehabt:

E(c) =m−1∑i=0

[yi − P(xi ; c)]2, P(x ; c) =n−1∑k=0

ckxk

mit m > n. Nun n = m.I Die explizite Losung ist gegeben durch die Lagrange

interpolierenden Polynome.I Beispiel: Gegeben sind (x0, y0) und (x1, y1). Die Lagrange

interpolierenden Polynome sind:

L0(x) =x − x1

x0 − x1L1(x) =

x − x0

x1 − x0Bemerke:

L0(x0) = 1 L1(x0) = 0L0(x1) = 0 L1(x1) = 1

oder Li(xj) = δij

und

P(x) =1∑

i=0

yiLi(x) erfullt P(xj) =1∑

i=0

yiLi(xj) =1∑

i=0

yiδij = yj

139

Lagrange Interpolierende PolynomeI Beispiel: Gegeben sind (x0, y0), (x1, y1) und (x2, y2). Die

Lagrange interpolierenden Polynome sind:

L0(x) =(x−x1)(x−x2)

(x0−x1)(x0−x2)L1(x) =

(x−x0)(x−x2)

(x1−x0)(x1−x2)L2(x) =

(x−x0)(x−x1)

(x2−x0)(x2−x1)

und

P(x) =2∑

i=0

yiLi(x) erfullt P(xj) =2∑

i=0

yiLi(xj) =2∑

i=0

yiδij = yj

I Im allgemeinen werden diese Polynome so definiert:

LN,i(x) =N∏

i 6=j=0

x − xj

xi − xj

PN(x) =N∑

i=0

yiLN,i(x)

PN(xj) =N∑

i=0

yiLN,i(xj) =N∑

i=0

yiδij = yj

I Was ist die Qualitat dieser Interpolation?140

Genauigkeit globaler InterpolationI Beispiel: Die Funktion f (x) = x2/(1 + 20x2) wird mit

n = 10 oder m = 100 Punkten folgendermasseninterpoliert:I P(xi ) = f (xi ), xi = −1 + 2i/n, i = 0, . . . ,nI Q(tj ) = f (tj ), tj = cos(π(j + 1/2)/(N + 1)), j = 0, . . . ,nI∑m

k=0 |R(yk )− f (yk )|2 = min, yk = −1 + 2k/m, k = 0, . . . ,m

(Gleichmaßig)

(Tchebychev)

(Ausgleich)

I Die Ergebnisse sehen so aus:

wobei die theoretischen Fehler rechts gezeigt werden.141

Genauigkeit globaler InterpolationSatz: Wenn xiNi=0 ⊂ [a,b] und f ∈ CN+1([a,b]) dann∀x ∈ [a,b], ∃ξ(x) ∈ [a,b] 3

f (x) = PN(x) +f (N+1)(ξ(x))

(n + 1)!ψN+1(x ; x0, . . . , xN)

wobei

PN(x) =N∑

i=0

f (xi)LN,i(x) ψN+1(x ; x0, . . . , xN) =N∏

i=0

(x − xi)

Insbesondere gilt

|f (x)− PN(x)| ≤ M|ψN+1(x ; x0, . . . , xN)|wenn |f (N+1)(ξ)| ≤ M(N + 1)!, ∀ξ ∈ [a,b].

I Wie kann man |ψN+1(x ; x0, . . . , xN)| minimieren?I Diese Funktion ist oben grafisch dargestellt fur die Falle:

I Gleichmaßig verteilte Stutzstellen undI Tchebychev Stutzstellen.

I Das zweite Ergebnis ist klar besser, aber man kann nichtimmer die Stutzstellen auswahlen.

142

Iterierte InterpolationDef: Gegeben sind die Daten (xi , fi)Ni=0. Mit verschiedenenmjkj=1, 0 ≤ mj ≤ N, wird das Polynom hier mit Pm1,...,mk

bezeichnet, das die Daten (xmj , fmj )kj=1 interpoliert.

I Beispiel: Fur xi = i , i = 0,1,2,

P0,2(x) =(x − x2)

(x0 − x2)f0 +

(x − x0)

(x2 − x0)f2

Satz: Gegeben sind die Daten (xi , fi)Ni=0, P0,...,j−1,j+1,...,N undP0,...,i−1,i+1,...,N . Dann gilt

P0,...,N(x) =(x − xj)P0,...,j−1,j+1,...,N(x)− (x − xi)P0,...,i−1,i+1,...,N(x)

xi − xj

I Dieser Satz fuhrt zu einem Algorithmus zur Berechnungdes Polynoms P0,...,N .

143

Iterierte Interpolation

I Beispiel: Gegeben seien (xi , fi)2i=0 und eineAuswertungsstelle x fur P012. Man bekommt P012(x) ausder Tabelle:

x0 P0(x) = f0x1 P1(x) = f1 P01(x) = (x−x1)P0−(x−x0)P1

x0−x1

x2 P2(x) = f2 P12(x) = (x−x2)P1−(x−x1)P2x1−x2

P012(x) = (x−x2)P01−(x−x0)P12x0−x2

I Wenn ein neuer Datenpunkt verfugbar wird, dann kann dieTabelle einfach erganzt werden:

x0 P0(x) = f0

x1 P1(x) = f1 P01(x) =(x−x1)P0−(x−x0)P1

x0−x1

x2 P2(x) = f2 P12(x) =(x−x2)P1−(x−x1)P2

x1−x2P012(x) =

(x−x2)P01−(x−x0)P12x0−x2

x3 P3(x) = f3 P23(x) =(x−x3)P2−(x−x2)P3

x2−x3P123(x) =

(x−x3)P12−(x−x1)P23x1−x3

P0123(x) =(x−x3)P012−(x−x0)P123

x0−x3

144

Nevilles AlgorithmusAlgorithmus (Neville): Eingaben sind Daten (xi , fi)Ni=0 unddie Auswertungsstelle x , Ausgabe ist der Wert P0,...,N(x).

Qi,0 = fi , i = 0, . . . ,Nfor i = 1, . . . ,N

for j = 1, . . . , i

Qij =(x − xi−j)Qi,j−1 − (x − xi)Qi−1,j−1

xi − xi−jend

endreturn QNN

Bildlich gesehen:

x0 Q00 Qij = Pi−j,...,i...

... Q11 QNN = PN−N,...,N...

.... . .

xN QN0 QN1 · · · QNN145

Dividierte DifferenzenDef: PN(x) bezeichnet hier das Polynom, das die Daten(xi , fi)Ni=0 interpoliert.I Nun wollen wir nicht nur den Wert PN(x) an einer

Auswertungsstelle x berechnen, sondern auchKoeffizienten eines Polynoms bestimmen.

I PN(x) kann so dargestellt werden:PN(x) = a0 + a1(x−x0) + · · ·+ aN(x−x0)×· · ·× (x−xN−1)

I Die Koeffizienten aiNi=0 finden wir so:

f0 = PN(x0) = a0 =: f [x0]f1 = PN(x1) = a0 + a1(x1 − x0) = f0 + a1(x1 − x0)

⇒ a1 =f1 − f0x1 − x0

=: f [x0, x1]

I Die restlichen Koeffizienten werden induktiv durch ahnlicheDifferenzen gegeben:

ak = f [x0, . . . , xk ], f [xi , . . . , xi+k ] =f [xi+1, . . . , xi+k ]− f [xi , . . . , xi+k−1]

xi+k − xi146

Dividierte DifferenzenI Das interpolierende Polynom ist so gegeben:

PN(x) = f [x0] +N∑

k=1

f [x0, . . . , xk ]k−1∏l=0

(x − xl)

Satz: Wenn f ∈ CN([a,b]) und xiNi=0 ⊂ [a,b] verschiedensind, dann ∃ξ ∈ (a,b) 3

f [x0, . . . , xN ] =f (N)(ξ)

N!

Def: Pk (Ω) bezeichnet die Menge der Polynome eines Gradeshochstens k auf einem Gebiet Ω.

Bemerkung: Seien PL,PD ∈ PN die durch Lagrangebeziehungsweise Dividierte Differenzen gegebenen Polynome,die die Daten (xi , fi)Ni=0 interpolieren. Dann giltPL(x) ≡ PD(x), weil Q = PL − PD ∈ PN , und Q hat die N + 1Nullstellen xiNi=0.147

Hermite InterpolationHermite InterpolationI Nun seien Daten (xi , fi)Ni=0 und (xi , fi)Ni=0 gegeben, und

eine interpolierende Funktion H(x) soll konstruiert werden,die erfullt H(xi) = fi , H ′(xi) = fi , i = 0, . . . ,N.

I Seien die Basis Funktionen definiert

HN,i(x) = [1− 2(x − xi)L′N,i(xi)]L2N,i(x)

HN,i(x) = (x − xi)L2N,i(x)

148

Hermite Interpolation

I Bemerke:HN,i(xj) = δij HN,i(xj) = 0H ′N,i(xj) = 0 H ′N,i(xj) = δij

und

H2N+1(x) =N∑

i=0

[fiHn,i(x) + fi Hn,i(x)

]erfullt

(?) H2N+1(xj) = fj , H ′2N+1(xj) = fj , i = 0, . . . ,N

Satz: Wenn f ∈ C1([a,b]) und xiNi=0 ⊂ [a,b] verschieden sind,dann ist H2N+1 ∈ P2N+1([a,b]) das eindeutige Polynomkleinsten Grades, das (?) erfullt mit f (xi) = fi , f ′(xi) = fi ,i = 0, . . . ,N.

149

Stuckweise Polynome InterpolationStuckweise Polynome InterpolationI Seien die Basis Funktionen definiert,

`i(x) =

x − xi−1

xi − xi−1, xi−1 ≤ x ≤ xi

xi+1 − xxi+1 − xi

, xi ≤ x ≤ xi+1

0, sonst

I Es gilt `i(xj) = δij .150

Stuckweise Kubisch Hermite InterpolationStuckweise Kubisch Hermite InterpolationI Nun seien Daten (xi , fi)Ni=0 und (xi , fi)Ni=0 gegeben, und

eine interpolierende Funktion H(x) soll konstruiert werden,die erfullt H(xi) = fi , H ′(xi) = fi , i = 0, . . . ,N.

I Seien die Basis Funktionen φiNi=0, ψiNi=0 definiert,

die stuckweise kubisch sind und erfullen:φi(xj) = δij ψi(xj) = 0φ′i(xj) = 0 ψ′i (xj) = δij152

Stuckweise Kubisch Hermite Interpolation

I Hausaufgabe: Konstruiere φiNi=0, ψiNi=0 fur einregelmaßiges Gitter.

I Die Funktion

H(x) =N∑

i=0

[fiφi(x) + fiψi(x)

]erfullt

H(xj) = fj , H ′(xj) = fj , i = 0, . . . ,N

Satz: Sei f ∈ C4([a,b]) und xiNi=0 mit h = maxi(xi+1 − xi).Dann gelten

maxa≤x≤b

|Dk [f (x)− H(x)]| ≤ ckh4−k maxa≤x≤b

|D4f (x)|

mit c0 = 1/384, c1 =√

3/216, c2 = 1/12 und c3 = 1/2.153

SplinesI Spline Funktionen sind stuckweise Polynome.I Obwohl die gewunschte Glattheit eingebaut werden kann,

werden hier nur die B-Splines mit maximaler Glattheituntersucht.

I Gezeigt wird eine Basis siN−1i=−k fur die Menge

Sk (x) = s ∈ Ck−1([x0, xN ]) : s ∈ Pk ([xi−1, xi ]),1 ≤ i ≤ N

fur k = 0, . . . ,3.I Basis fur stuckweise konstante Funktionen:

S0(x) 3 f (x) =N−1∑i=0

αisi(x)

154

SplinesI Basis fur stuckweise lineare Funktionen in C0([a,b]):

S1(x) 3 f (x) =N−1∑i=−1

αisi(x)

I Basis fur stuckweise quadratische Funktionen in C1([a,b]):

S2(x) 3 f (x) =N−1∑i=−2

αisi(x)

155

SplinesI Basis fur stuckweise kubische Funktionen in C2([a,b]):

S3(x) 3 f (x) =N−1∑i=−3

αisi(x)

I Um jede Basis zu erzeugen,

πk (x) =∫ +∞

−∞πk−1(x − y)π0(y)dy

k∑i=1

πk (i) = 1156

Die Kanonischen Splines

I Die Kanonischen Splines explizit:

π0(x) =

1 x ∈ [0,1]0 sonst

π1(x) =

x x ∈ [0,1]2− x x ∈ [1,2]0 sonst

=

∫ 1

0π0(x − y)dy

π2(x) =

12x2 x ∈ [0,1]34 − (x − 3

2)2 x ∈ [1,2]12(x − 3)2 x ∈ [2,3]0 sonst

=

∫ 1

0π1(x − y)dy usw

π3(x) =

16x3 x ∈ [0,1]16 [1 + 3(x − 1) + 3(x − 1)2 − 3(x − 1)3] x ∈ [1,2]16 [1 + 3(3− x) + 3(3− x)2 − 3(3− x)3] x ∈ [2,3]16(4− x)3 x ∈ [3,4]0 sonst

157

GlattheitsbedingungenI Wenn das Gitter gleichmaßig ist, x = xiNi=0, h = xi − xi−1,

ist eine Basis fur Sk (x) explizit gegeben durch:

sk ,i(x) = πk ((x − xi)/h), i = −k , . . . ,N − 1

I Im allgemeinen werden die Basis Funktionen durchGlattheitsbedingungen bestimmt.

I Fur die linearen Splines: 4 Unbekannte, 4 Bedingungen,

I Hausaufgabe: Schreibe das Gleichungssystem fur dieKoeffizienten der Splines Grades 1, 2 und 3.

158

GlattheitsbedingungenI Fur die quadratischen Splines: 9+9

I Fur die kubischen Splines: 16+16

159

Interpolation mit Kubischen SplinesInterpolation mit Kubischen SplinesI Gegeben seien die Daten (xi , fi)Ni=0 mit xiNi=0 =

x ∈ RN+1.I Zur Vereinfachung wird angenommen, dass das Gitter

regelmaßig ist, obwohl es nicht notwendig ist.I Die Dimension von S3(x) ist N + 3.I Im Datensatz fehlen noch 2 Informationsstucke, um eine

Interpolantes(x) =

N−1∑i=−3

αisi(x) ∈ S3

zu bestimmen.I Eine mogliche Erganzung sind die Randbedingungen:

0 = s′(x0) = s′−3(x0)︸︷︷︸=−1/(2h)

α−3+s′−2(x0)︸︷︷︸=0

α−2+s′−1(x0)︸︷︷︸=1/(2h)

α−3 =α−1 − α−3

2h

und s′(xN) = 0 oderα−3 = α−1 und αN−1 = αN−3.

160

Interpolation mit Kubischen SplinesI Mit den Randbedingungen:

0 = s′′(x0) = s′′−3(x0)︸︷︷︸=1/h2

α−3+s′′−2(x0)︸︷︷︸=−2/h2

α−2+s′′−1(x0)︸︷︷︸=1/h2

α−3 =α−1 − 2α−2 + α−3

h2

und s′′(xN) = 0 bekommt manα−3 = 2α−2 − α−1 und αN−1 = 2αN−2 − αN−3.

Durch

f0 = s(x0) = s−3(x0)α−3 + s−2(x0)α−2 + s−1(x0)α−3= 1

6α−3 + 23α−2 + 1

6α−3= 1

6(2α−2 − α−1) + 23α−2 + 1

6α−3 = α−2

sieht man, fur diese Randbedingungen gelten

s′(x0) =α−1 − α−3

2h=α−1 − α−2

h=α−1 − f0

h≈ f ′(x0)

und ahnlich s′(xN) ≈ f ′(xN).161

Interpolation mit Kubischen SplinesI Um die Interpolante s ∈ S3(x) zu bestimmen, werden die

Koeffizienten der Basis Funktionen durch die Losung deslinearen Gleichungssystems berechnet:

fj = s(xj) = α−3s−3(xj)+N−2∑i=−2

αisi(xj)+αN−1sN−1(xj), j = 0, . . . ,N

mit Randbedingungen, z.B.

α−3 = α−1 oder α−3 = 2α−2 − α−1αN−1 = αN−3 αN−1 = 2αN−2 − αN−3

I Bemerke:

si(xj) = 0, j < i+1; si(xi+1), si(xi+2), si(xi+3) 6= 0; sk (xj) = 0, j > i+3

Die Matrixsi(xj) : −3 ≤ i ≤ N − 1,−1 ≤ j ≤ N + 1

ist tridiagonal!I Die entsprechende Matrix fur lineare Splines ist diagonal!

162

Interpolation mit Kubischen SplinesI Das System sieht so aus:

f0......

fN

=

1− θ θ

16

23

16

. . .. . .

. . .16

23

16

θ 1− θ

α−2......

αN−2

wobei θ = 1

3 fur die Randbedingungen s′(x0) = 0 = s′(xN)und θ = 0 fur die Randbedingungen s′′(x0) = 0 = s′′(xN).

Satz: Sei f ∈ C1([a,b]). Dann ∃!s ∈ S3 3 s(xi) = f (xi),0 ≤ i ≤ N, und s′(x0) = f ′(x0), s′(xN) = f ′(xN). Furf ∈ C4([a,b]) und h = max

0≤i≤N+1(xi+1 − xi) gilt

maxa≤x≤b

|Dr [f (x)− s(x)]| ≤ cr h4−r maxa≤x≤b

|D4f (x)|, 0 ≤ r ≤ 3

wobei c0 = 5/384, c1 = 1/24, c2 = 3/8 und c3 = 1/2(β + 1/β),β = h/ min

0≤i≤N−1(xi+1 − xi).

163

Tensor Produkte fur 2D InterpolationTensor Produkte fur 2D InterpolationI Gegeben seien Daten (xi , yj , fij)Ni,j=0 auf dem 2D Gitter

xi = ih, i = 0, . . . ,N, yi = jh, j = 0, . . . ,N, h = 1/N.I Die Interpolante ist

s(x , y) =N−1∑i=−3

N−1∑j=−3

αijsi(x)sj(y)

wobei die Koeffizienten αijN−1i,j=−3 durch Losung des

Gleichungssystems bestimmt werden:

fkl = s(xk , yl) =N−1∑i=−3

N−1∑j=−3

αijsi(xk )sj(yl), 0 ≤ k , l ≤ N

mit z.B. Randbedingungen, θ = 1 oder 2, φ = 1− θ,

α−3,j = θα−2,j + φα−1,j , αN−1,j = θαN−2,j + φαN−3,j , −3 ≤ j ≤ N − 1αi,−3 = θαi,−2 + φαi,−1, αi,N−1 = θαi,N−2 + φαi,N−3, −3 ≤ i ≤ N − 1

I Solche Tensor Produkte konnen fur andere BasisFunktionen ahnlich definiert werden, wie z.B. StuckweiseKubisch Hermite Interpolanten.164

Nicht Lineare Gleichungen

I Das Ziel ist, eine Nullstelle zu finden,f (x) = 0

wobei f nicht linear ist.I Wenn f eine Ableitung ist, F ′(x) = f (x), kann das Ziel sein,

F zu minimieren oder maximieren.I Gegeben seien a und b mit f (a)f (b) < 0. Nach dem

Zwischenwertsatz ∃ξ ∈ (a,b) mit f (ξ) = 0 wennf ∈ C([a,b]).

I Die einfachste Methode ist das Bisektionsverfahren:

f (a)f (b) < 0 . . . (a + b)/2 = c →

a, f (c)f (b) < 0b, f (a)f (c) < 0

I In einem Kode soll c entweder a oder b uberschreiben.I Solang f (a)f (b) < 0 gilt, ist es egal welcher Wert positiv

oder negativ ist.

165

BisektionsverfahrenI Pseudo-Kode fur das Bisektionsverfahrenfa = f(a), fb = f(b)if (fa * fb > 0)

error: sign(fa) = sign(fb)endfor k=1,...,kmax

c = (a + b)/2fc = f(c)if (|b-a| < tol*|c|)

return cendif (fa * fc > 0)

a = c % sign(fa) = sign(fc)else

b = c % sign(fb) = sign(fc)end

enderror: failed to converge to relativedifference tol after kmax iterations

166

BisektionsverfahrenSatz: Sei f ∈ C([a,b]) mit f (a)f (b) < 0. Die Iterierten ck desBisektionsverfahrens konvergieren zu einer Nullstellex? ∈ (a,b), f (x?) = 0, und zwar mit folgender Geschwindigkeit:

|x? − ck | ≤ 2−k |b − a|, k ≥ 1

Beweis: Sei [ak ,bk ] das k te Bisektionsintervall mitck = (ak + bk )/2 und [a1,b1] = [a,b]. Es gilt

(b − a) = (b1 − a1) = 2(b2 − a2) = · · · = 2k−1(bk − ak ).Aus ak+1 ∈ ak , ck und bk+1 ∈ ck ,bk folgen

(ak+1 − ak ) ≤ 12(bk − ak ) und (bk − bk+1) ≤ 1

2(bk − ak )und daher

|ck+1 − ck | = 12 |(ak+1 − ak )− (bk − bk+1)| ≤

12 [|ak+1 − ak |+ |bk − bk+1|] ≤ 1

2(bk − ak ) = 2−k (b − a).Also fur k > l gilt|ck − cl | ≤ |ck − ck−1|+ · · · |cl+1 − cl | =

∑k−1m=l |cm+1 − cm| ≤∑k−1

m=l 2−m(b − a) = (b − a)(2−l − 2−k )/(1− 2−1)k ,l→∞−→ 0.

167

BisektionsverfahrenDeswegen ist ck eine Cauchy Folge. Sei x? der Limes.Aus ck ∈ (ak ,bk ) folgen|x? − ak | ≤ |x? − ck |+ |ck − ak | ≤ |x? − ck |+ |bk − ak |

k→∞−→ 0.und|x? − bk | ≤ |x? − ck |+ |ck − bk | ≤ |x? − ck |+ |ak − bk |

k→∞−→ 0.Daher aus f ∈ C([a,b]) folgen

limk→∞

f (ak ) = f (x?) = limk→∞

f (bk ).

Aus f (ak )f (bk ) ≤ 0, ∀k , folgtlim

k→∞f (ak )f (bk ) ≤ 0.

Diese Limiten implizieren0 ≤ f (x?)2 = lim

k→∞f (ak )f (bk ) ≤ 0

und es folgt f (x?) = 0.

Bemerkung: Obwohl Konvergenz fur das Bisektionsverfahrengarantiert ist, ist die Geschwindigkeit niedriger als die fur dasNewton Verfahren, das aber nur lokal konvergiert.

168

FixpunktiterationI Die Nullstelle in der Gleichung f (x) = 0 kann bezuglich

eines Fixpunkts x = g(x) einer Funktion g(x) inverschiedenen Weisen umgeschrieben werden:

g(x) = x − f (x), g(x) = x − f (x)/f ′(x)I Dann gibt es eine naturliche Fixpunktiteration:

xk+1 = g(xk ), k = 0,1,2, . . .I Wann konvergiert diese Iteration? Grafische Darstellung:

I Wichtige Eigenschaft: |g′(x?)| < 1.169

FixpunktiterationSatz: Sei g ∈ C([a,b]) mit g([a,b]) ⊂ [a,b]. Dann ∃x? ∈ [a,b] 3g(x?) = x?. Falls g ∈ C1([a,b]) und ∃γ 3 |g′(x)| ≤ γ < 1,∀x ∈ [a,b], dann ist x? ein eindeutiger Fixpunkt.

Beweis: g([a,b]) ⊂ [a,b]⇒ f (x) := g(x)− x erfulltf (a) = g(a)− a ≥ a− a = 0 und f (b) = g(b)− b ≤ b − b = 0.Wenn f (a) = 0 gelten sollte, ist x? = a ein Fixpunkt. Ahnlichwenn f (b) = 0 gelten sollte. Nimm an, f (b) < 0 < f (a). Nachdem Zwischenwertsatz ∃x? ∈ [a,b] 3 0 = f (x?) = g(x?)− x?,da f stetig ist. Wenn es einen zweiten Fixpunkt x ∈ [a,b] gabe,dann wurde ein Widerspruch durch den Mittelwertsatzentstehen:

|x? − x | = |g(x?)− g(x)| = |g′(ξ)(x? − x)|≤ γ|x? − x |< |x? − x |

wobei ξ zwischen x? und x in [a,b] liegt.170

FixpunktiterationSatz: Sei g ∈ C1([a,b]) mit g([a,b]) ⊂ [a,b]. Wenn ∃γ 3|g′(x)| ≤ γ < 1, ∀x ∈ [a,b], dann ∀x0 ∈ [a,b] konvergiert dieFixpunktiteration xk+1 = g(xk ) zum eindeutigen Fixpunkt x?.

Beweis: g([a,b]) ⊂ [a,b]⇒ xk+1(= g(xk )) ∈ [a,b], ∀k . Nachdem Mittelwertsatz gilt

|xk+1 − x?| = |g(xk )− g(x?)| = |g′(ξk )(xk − x?)|≤ γ|xk − x?| ≤ γ2|xk−1 − x?| ≤ · · ·≤ γk+1|x0 − x?| k→∞−→ 0

wobei ξk zwischen x? und xk in [a,b] liegt.

Satz: Unter den Bedingungen des letzten Satzes gilt

|x? − xk | ≤γk

1− γ|x1 − x0|, k ≥ 1

Beweis: Hausaufgabe.171

FixpunktiterationI Pseudo-Kode fur eine Fixpunktiterationx0 gegebenfor k=1,...,kmax

x = g(x0)if (|x-x0| < tol*|x|)

return xendx0 = x


I Welche Iteration ist schneller, Bisektion oder Fixpunkt?Hangt von der Funktion ab: g′(x?) = 0 favorisiert Fixpunkt.

I Beispiel: g(x) = xer(1−x), r ≥ 1, f (x) = x − g(x).Entspricht einem diskreten Populations-Modell.

I Hausaufgabe: Zeige, fur r ∈ [1,2),I ∃ε > 0 3 g(B(1, ε)) ⊂ B(1, ε) undI ∃γ ∈ (0,1) 3 |g′(x)| ≤ γ, ∀x ∈ B(1, ε).

172

Vergleich zweichen Bisektions- und Fixpunktiteration

I Fur die erste Grafik (r = 1,g′(1) = 0) verlangt Bisektion 12Iterationen und Fixpunkt 6 Iterationen, um x = 0.99995bzw. x = 0.999999995 zu berechnen.

I Fur die zweite Grafik (r = 1.99,g′(x?) = −0.99) verlangtBisektion 12 Iterationen und Fixpunkt 542 Iterationen, umx = 0.99995 bzw. x = 1.0005 zu berechnen.

173

Newton Verfahren

I Die Strategie: Eine Nullstelle wird durch eine Folge vonlinearen Approximationen der Funktion bestimmt:

I Fur ein gegebenes x0 ist

y − f (x0) = f ′(x0)(x − x0)

die zu f tangente Geradean der Stelle (x0, f (x0)).

I Die Nullstelle dieser Gerade erfullt−f (x0) = f ′(x0)(x − x0)oderx = x0 − f (x0)/f ′(x0)

I Im allgemeinen werden die Newton Iterierten so gegeben:xk+1 = xk − f (xk )/f ′(xk ), k = 0,1,2, . . .

I Wichtig ist, dass f ′(x?) 6= 0.

174

Newton Verfahren

I Pseudo-Kode fur das Newton Verfahrenf(x), fp(x) = f’(x), x0 gegebenfor k=1,...,kmax

f0 = f(x0)fp0 = fp(x0)if (fp0 =/= 0)

x1 = x0 - f0/fp0else

x1 = x0endif (|x1-x0| < tol*|x1|)

return x1endx0 = x1


175

Sekant VerfahrenI Pseudo-Kode fur das Sekant Verfahrenf(x), x0 =/= x1 gegebenf1 = f(x0)for k=1,...,kmax

f0 = f1f1 = f(x1)if (f1 =/= f0)

x2 = x1 - f1 * (x1 - x0)/(f1 - f0)else

x2 = x1endif (|x2-x1| < tol*|x2|)

return x2endx0 = x1x1 = x2


176

Vergleich zweichen Newton und Bisektionsiteration

I Fur die erste Grafik (r = 1) verlangt Bisektion 12Iterationen und Newton 5 Iterationen, um x = 0.99995bzw. x = 1.0000000000006 zu berechnen.

I Fur die zweite Grafik (r = 1.99) verlangt Bisektion 12Iterationen und Newton 4 Iterationen, um x = 0.99995bzw. x = 1.0! zu berechnen.

177

Konvergenz des Newton VerfahrensSatz: Sei f ∈ C2([a,b]) mit f (x?) = 0 und f ′(x?) 6= 0 furx? ∈ (a,b). Dann ∃ε > 0 3 ∀x0 ∈ B(x?, ε) die Newton Iteriertenxn zu x? konvergieren.

Beweis: Die Funktion g(x) = x − f (x)/f ′(x) erfullt g(x?) = x?.Da f ′(x?) 6= 0 und f ′ ∈ C([a,b]), ∃δ > 0 3 f ′(x) 6= 0,∀x ∈ B(x?, δ). Deswegen ist g(x) in B(x?, δ) wohl definiert undstetig. Aus der Rechnung

g′(x) = 1− (f ′(x)2 − f (x)f ′′(x))/f ′(x)2 = f (x)f ′′(x)/f ′(x)2

folgt g′ ∈ C(B(x?, δ)) und g′(x?) = 0. Deswegen ∃ε > 0 undγ ∈ (0,1) 3 |g′(x)| ≤ γ,∀x ∈ B(x?, ε). Nun sei x ∈ B(x?, ε).Nach dem Mittelwertsatz ∃ξ zwischen x und x? 3

|g(x)− g(x?)| = |g′(ξ)(x − x?)|≤ γ|x − x?| < ε

Es folgt, g(B(x?, ε)) ⊂ B(x?, ε), und die Bedingungen desobigen Fixpunkt-Konvergenz-Satzes sind erfullt.

178

Asymptotische KonvergenzrateBemerkung: Mit g′(x?) = 0 ergibt sich die bestmogliche lokaleKonvergenz der Fixpunktiteration.

Bemerkung: Das Newton Verfahren konvergiert nur lokal, d.h.nur wenn man nah genug beim Fixpunkt x? startet.

Bemerkung: Um die Vorteile des Bisektions- und des NewtonVerfahrens zu kombinieren, kann man mit Bisektion beginnenund dann auf Newton umschalten. Wenn die Abstandezwischen Newton Iterierten nicht kleiner werden, kann manwieder kurzfristig auf Bisektion welchseln und spater aufNewton zuruckschalten.

Def: Sei xn gegeben mit xmk→∞−→ x?. Wenn gilt

limk→∞

|x? − xk+1||x? − xk |α

= λ ∈ (0,1]

hat die Folge Konvergenzordnung α mit asymptotischerFehlerkonstante λ.

179

Asymptotische KonvergenzrateI Beispiele:α = 1⇒ aymptotisch lineare Konvergenzrate,

z.B. xk+1 = g(xk ), 0 < |g′(x?)| < 1α = 2⇒ aymptotisch quadratische Konvergenzrate,

z.B. xk+1 = g(xk ), |g′(x?)| = 0

Satz: Sei g ∈ C1([a,b]) mit g([a,b]) ⊂ [a,b], g(x?) = x? ∈ (a,b)und |g′(x)| ≤ γ < 1, ∀x ∈ [a,b]. Wenn g′(x?) 6= 0, konvergiertxk+1 = g(xk ) asymptotisch nur linear zu x?.

Beweis: Nach dem Mittelwertsatz, ∃ξk zwischen xk und x? in[a,b] 3

|xk+1 − x?| = |g(xk )− g(x?)| = |g′(ξk )||xk − x?|Nach einem obigen Konvergenzsatz 171 gilt xk

k→∞−→ x?, unddaher gilt auch ξk

k→∞−→ x?. Dann g ∈ C1([a,b])⇒ g′(ξk )k→∞−→

g′(x?). Daher gilt|x? − xk+1||x? − xk |

= |g′(ξk )| k→∞−→ |g′(x?)| 6= 0.180

Asymptotische KonvergenzrateSatz: Sei g ∈ C2([a,b]) mit g(x?) = x? ∈ (a,b) und g′(x?) = 0.Dann ∃ε > 0 3 ∀x0 ∈ B(x?, ε) konvergiert xk+1 = g(xk )mindestens asymptotisch quadratisch zu x?.

Beweis: Wahle ε > 0 so aus, dass |g′(x)| ≤ γ < 1,∀x ∈ B(x?, ε) ⊂ [a,b]. Seien xlkl=0 ⊂ B(x?, ε). Nach demMittelwertsatz ∃ξk zwischen xk und x? in B(x?, ε) 3|x? − xk+1| = |g(x?)− g(xk )| = |g′(ξk )||x? − xk | < |x? − xk |

Es folgt, xk+1 ∈ B(x?, ε) und daher xk∞k=0 ⊂ B(x?, ε). Nachdem Taylorsatz ∃ηk zwischen xk und x? in B(x?, ε) 3

xk+1 = g(xk ) = g(x?)︸︷︷︸=x?

+ g′(x?)︸︷︷︸=0

(xk − x?) + 12g′′(ηk )(xk − x?)2

Nach einem obigen Konvergenzsatz gilt xkk→∞−→ x?, und daher

gilt auch ηkk→∞−→ x?. Daher folgt

|x? − xk+1||x? − xk |2

=12|g′′(ηk )| k→∞−→ 1

2|g′′(x?)| = λ <∞.

Bemerkung: Mit g′′(x?) = 0 ist α noch hoher.181

Nullstellen hoherer MultiplizitatI Beispiel: f ′(x?) 6= 0⇒ g(x) = x − f (x)/f ′(x) erfullt

g′(x) = 1− f ′(x)2 − f (x)f ′′(x)

f ′(x)2 =f (x)f ′′(x)

f ′(x)2x→x?−→ 0

Was ist wenn f ′(x?) = 0?Def: Die Nullstelle, f (x?) = 0, hat Multiplizitat m, wenn m diehochste Zahl ist, die erfullt limx→x? |f (x)|/|x − x?|m ∈ (0,∞).Bei m = 1 ist die Nullstelle einfach.I Kann Newton asymptotisch quadratisch konvergieren,

wenn x? keine einfache Nullstelle fur f ist? Im allgemeinennicht.

I Beispiel: f (x) = ex − x − 1, x? = 0, m = 2,g(x) = x − f (x)/f ′(x) aber g′(0) = 1

2 6= 0!I Wenn die Nullstelle f (x?) = 0 Multiplizitat m > 1 hat,

verwende das modifizierte Newton Verfahren:gm(x) = x −mf (x)/f ′(x)

I Hausaufgabe: Zeige asymptotisch quadratischeKonvergenz fur gm.

182

Systeme von Nicht Linearen GleichungenI Beispiel: Minimiere

Jh(u) =h2

N∑i=0

(ui − vi)2 + µh

N∑i=1

[(ui − ui−1

h

)2

+ ε2

] 12

≈ 12

∫Ω

(u − v)2dx + µ

∫Ω

√|u′|2 + ε2dx

I Fur ∇Jh(u) = 0 sei Di(u)=

[(ui−ui−1

h

)2

+ε2

]− 12

=Di

∂Jh

∂uk︸︷︷︸0<k<N

= h(uk − vk ) + µ/h[(uk − uk−1)Dk︸︷︷︸i=k

− (uk+1 − uk )Dk+1︸︷︷︸i=k+1

]

∂Jh

∂u0= h(u0 − v0) + µ/h[ − (u1 − u0)D1︸︷︷︸

i=k+1=1

]

∂Jh

∂uN= h(uN − vN) + µ/h[(uN − uN−1)DN︸︷︷︸

i=k=N

]

183

OptimalitatssystemI Die Optimalitatsbedingung ist

0 =1h∇Jh(u) = u − v +

µ

h2 D(u)u =: F (u)

wobei

D(u) =

D1 −D1−D1 D1 + D2 −D2

. . .. . .

. . .

−DN−1 DN−1 + DN −DN−DN DN

I Gesucht wird eine Nullstelle F (u) = 0.I Ein Ansatz ist, eine Fixpunktiteration uk+1 = G(uk ) zu

verwenden, wobei

G(u) = v − µ

h2 D(u)u

184

Konvergenz der FixpunktiterationSatz: Fur B ⊂⊂ Rn sei G ∈ C(B,Rn) mit G(B) ⊂ B. Dann∃u? ∈ B 3 G(u?) = u?. Wenn G ∈ C1(B,Rn) und es giltρ(G′(u)) ≤ γ < 1, ∀u ∈ B, ist u? ein eindeutiger Fixpunkt in B.Wenn fur eine induzierte Matrixnorm ‖ · ‖∗ gilt‖G′(u)‖∗ ≤ γ < 1, ∀u ∈ B, erfullt die Fixpunktiterationuk+1 = G(uk ) die Abschatzung,

‖u? − uk‖∗ ≤γk

1− γ‖u1 − u0‖∗ k ≥ 1

I Beispiel: Es wird fur das obige Beispiel gezeigt, dieBedingungen dieses Satzes werden mit B = B∞(v , 2µ

h )

erfullt, und es gilt ‖G′(u)‖∞ ≤ 4µεh2 , ∀u ∈ B.

I Um G(B) ⊂ B zu zeigen,

|[G(u)−v ]i |=∣∣∣− µ

h2 [D(u)u]i

∣∣∣= µ

h

∣∣∣∣Diui − ui−1

hδi>0−Di+1

ui+1 − ui

hδi<N

∣∣∣∣≤ 2µh

Also ‖G(u)− v‖∞ ≤ 2µh , ∀u ∈ RN+1.

185

Konvergenz der FixpunktiterationI Um ‖G′(u)‖∞ ≤ 4µ

εh2 zu zeigen, wird ∇Gi(u) berechnet,

Gi = vi −µ

h

[Di

ui − ui−1

hδi>0 − Di+1

ui+1 − ui

hδi<N

]Man sieht beispielsweise mit Di(u)=

[(ui−ui−1

h

)2

+ε2

]− 12

∂Gi

∂ui−1= −µ

h

[D3

i

(ui − ui−1

h

)2 1h− Di

1h

]δi>0

=µ

h2 D3i

[D−2

i −(

ui − ui−1

h

)2]δi>0 =

µε2

h2 D3i δi>0

und ahnlicherweise,

∂Gi

∂ui= −µε

2

h2 (D3i δi>0 + D3

i+1δi<N),∂Gi

∂ui+1=µε2

h2 D3i+1δi<N

I Da εDi ≤ 1 gilt, folgt ‖G′(u)‖∞ ≤ 4µεh2 , ∀u ∈ RN+1.

186

Abstiegsverfahren

I Nach dem Konvergenzsatz ∃!u? ∈ B 3 G(u?) = u? und esgilt

‖u? − uk‖∞ ≤γk

1− γ‖u1 − u0‖∞

solang 4µ/(εh2) ≤ γ < 1.

Abstiegsverfahren

I Man sucht in die Richtung −∇Jh(u)

uk+1 = uk − αh∇Jh(uk )

= uk − αuk + αv − αµh2 D(uk )uk

= (1− α)uk + αG(uk )=: GA(uk )

I Hausaufgabe: Zeige fur α ∈ (0,1), es gilt GA(B) ⊂ B furB = B∞(v , 2µ

h ) und ‖G′A(u)‖∞ ≤ 1− α + 4αµεh2 , ∀u ∈ B.

187

Implizite FixpunktiterationImplizite FixpunktiterationI Die Optimalitatsbedingung,

0 =1h∇Jh(u) = u − v +

µ

h2 D(u)u

kann auch so umgeschrieben werden:[I +

µ

h2 D(u)]

u = v

I Dann ist die Iteration[I +

µ

h2 D(uk )]

uk+1 = v

eine implizite Fixpunktiteration mit der Fixpunktfunktion:

GI(u) =[I +

µ

h2 D(u)]−1

v

I Hausaufgabe: Zeige es gilt GI(B) ⊂ B fur B = B2(0, ‖v‖2)und ‖G′I(u)‖2 ≤ 16µ

(ε3h4)‖v‖22, ∀u ∈ B.

188

Newton Verfahren fur SystemeI Was ist mit Newton?

F ′(uk )(uk+1 − uk ) = −F (uk ), k = 0,1,2, . . .Aus

F (u) = u−v− µ

h2 D(u)u = u−G(u), F ′(u) = I−G′(u)

folgt

F ′(u)= I+µε2

h2 E(u)= I+µε2

h2

D3

1 −D31

−D31 D3

1 + D32 −D3

2. . .

. . .. . .

−D3N−1 D3

N−1 + D3N −D3

N−D3

N D3N

I F ′(u) = I + µε2

h2 E(u) ist SPD ∀u ∈ RN+1:I Die Gerschgorin Scheiben in C sind B(zi , ri) wobei

z0 = 1 + µε2

h2 D31 r0 = µε2

h2 D31

zi = 1 + µε2

h2 (D3i + D3

i+1) ri = µε2

h2 (D3i + D3

i+1)

zN = 1 + µε2

h2 D3N rN = µε2

h2 D3N

I Da εDi ≤ 1 gilt, folgt 1 ≤ zi ± ri ≤ 1 + 4µ/(εh2), ∀i .189

Newton Verfahren fur SystemeSatz: Fur B ⊂⊂ Rn sei G ∈ C2(B,Rn) mit G(B) ⊂ B. Dann∃u? ∈ B 3 G(u?) = u?. Wenn gilt G′(u?) = 0, dann ∃ε > 0 3∀u0 ∈ B(u?, ε) die Fixpunktiteration uk+1 = G(uk )quadratisch zu u? konvergiert.

I Fur die Newton Iteration sei

GN(u) = u − F ′(u)−1F (u)

fur F ausreichend glatt, und es gilt

G′N(u) = I − [∂uF ′(u)−1] F (u)︸︷︷︸→0,u→u?

−F ′(u)−1F ′(u)u→u?−→ 0

Somit werden die Bedingungen des letzten Satzes erfullt.I Fur das obige Beispiel ist F ′(u) SPD ∀u ∈ RN+1, und

daher erfullt die Newton Iteration die Bedingungen desobigen Satzes.

190

Vergleich der Iterativen VerfahrenI Wir vereinfachen das obige Beispiel, um die verschiedenen

Ansatze zu vergleichen und zu visualisieren.I Nimm N = 1, h = 1, u = (u0,u1) und v = (v0, v1),

J(u) = 12(u0 − v0)2 + 1

2(u1 − v1)2 + µ√

(u1 − u0)2 + ε2

I Das Abstiegsverfahren ist uk+1 = GA(uk ) oder

uk+1 = (1− α)uk + α(

v − µ

h2 D(uk )uk

)I Die Implizite Fixpunktiteration ist uk+1 = GI(uk ) oder[

I +µ

h2 D(uk )]

uk+1 = v

I Die Newton Iteration ist uk+1 = GN(uk ) oder[I +

µε2

h2 E(uk )

](uk+1 − uk ) = −

[uk − v +

µ

h2 D(uk )uk

]I Die Implizite Fixpunktiteration funktioniert hier am besten.

191

Richtungsfelder der LosungsverfahrenI Fur N = 1, h = 1, v0 = 0.1, v1 = 0.5, ε = 10−2 und µ = 1

die Richtungsfelder uk+1 − uk , uk ∈ (0,1)2, der jeweiligenMethoden sind hier grafisch dargestellt:

I Der Punkt markiert v .I Fur ε klein und µ groß, soll das Minimum u sich naher bei

u0 = u1 d.h. dem Punkt befinden.I Nur das Richtungsfeld fur das Implizite Fixpunkt Verfahren

zeigt uberall zuverlassig in die Richtung des Minimums.192

Entrauschen ErgebnisseI Nun N = 128, ε = 10−3, µ = 7 · 10−3 und v ist rot:

I Zusammenfassung der (notwendigen) Parameter:Methode Dampfung Iterationen Fehler ZeitAbstieg: α = 10−3 104 (max) 2 · 10−3 2.3Implizit: 0 22 9 · 10−5 0.02Newton: ω = 10−2 104 (max) 2 · 10−3 8.5

wobei:Fehler = ‖uk+1 − uk‖2/‖uk+1‖2, Zeit = SekundenNewton Dampfung: uk+1 ← ωuk+1 + (1− ω)uk193

Numerisches Differenzieren und IntegrierenWir beginnen mit Approximationen der Ableitungen.

I Fur die erste Ableitung, Vorwarts Differenzen:

f ′(x) =f (x + h)− f (x)

h+O(h)

wobei mit dem Taylorsatz ∃ξ ∈ [x , x + h] 3

f (x + h) = f (x) + f ′(x)h + 12 f ′′(ξ)h2︸︷︷︸

hO(h)

I Ruckwarts Differenzen:

f ′(x) =f (x)− f (x − h)

h+O(h)

wobei mit dem Taylorsatz ∃η ∈ [x − h, x ] 3

f (x − h) = f (x)− f ′(x)h + 12 f ′′(η)h2︸︷︷︸

hO(h)

194

Approximation der AbleitungenI Zentrale Differenzen:

f ′(x) =f (x + h)− f (x − h)

2h+O(h2)

wobei mit dem Taylorsatz ∃η ∈ [x − h, x ], ξ ∈ [x , x + h] 3

f (x + h) = f (x) + f (1)(x)h + 12 f (2)(x)h2 + 1

6 f (3)(ξ)h3

f (x − h) = f (x)− f (1)(x)h + 12 f (2)(x)h2 − 1

6 f (3)(η)h3

undf (x + h)− f (x − h)

2h−f ′(x) =

112

f (3)(ξ)h2 +1

12f (3)(η)h2︸︷︷︸

O(h2)I Ahnlich fur die zweite Ableitung, wie oben 18 gezeigt,

f ′′(x) =f ′(x + h

2 )− f ′(x − h2 )

h+ O(h2)

=1h

[f (x + h)− f (x)

h− f (x)− f (x − h)

h+O(h3)

]+ O(h2)

=f (x + h)− 2f (x) + f (x − h)

h2 + O(h2)

195

Approximation der AbleitungenI Fur die vierte Ableitung,

f (4)(x) =f (2)(x + h)− 2f (2)(x) + f (2)(x − h)

h2 +O(h2)

=1h2

[f (x + 2h)− 2f (x + h) + f (x)

h2

−2f (x + h)− 2f (x) + f (x − h)

h2

+f (x)− 2f (x − h) + f (x − 2h)

h2 +O(h4)

]+O(h2)

=f (x − 2h)− 4f (x − h) + 6f (x)− 4f (x + h) + f (x + 2h)

h4 +O(h2)

I Fur Ableitungen hoherer Ordnung,

δnh f (x)

hn =1hn

n∑k=0

(−1)k(

nk

)f (x+(n/2−k)h) = f (n)(x)+O(h2)

(Bei n ungerade δnh f (x ± h

2 ) mitteln, um f auf hZ auszuwerten.)I Wie bestimmt man die Koeffizienten im allgemeinen?

196

Approximation der AbleitungenI Beispiel: Wir bestimmen die Koeffizienten a−1,a0,a1

einer Approximation der ersten Ableitung,

f ′(x) = a−1f (x − h) + a0f (x) + a1f (x + h) +O(hp)

um die Genauigkeit O(hp) zu maximieren.I Mit dem Taylorsatz,

f (x − h) = f (x)− f (1)(x)h + 12 f (2)(x)h2 − 1

6 f (3)(ξ(x))h3

f (x + h) = f (x) + f (1)(x)h + 12 f (2)(x)h2 + 1

6 f (3)(η(x))h3

bilden wir die gewichtete Summe,

a−1f (x − h) + a0f (x) + a1f (x + h) = (a−1 + a0 + a1)f (x)

+(−a−1 + a1)f (1)(x)h + (a−1 + a1)12 f (2)(x)h2

−a−116 f (3)(ξ(x))h3 + a1

16 f (3)(η(x))h3

mit der ein System fur die Koeffizienten hergeleitet wird: 1 1 1−1 0 1

1 0 1

a−1a0a1

=

0h−1

0

⇒ a−1

a0a1

=1

2h

−101

197

Approximationen der AbleitungenI Aus

−a−1

6f (3)(ξ(x))h3+

a1

6f (3)(η(x))h3 =

112

[f (3)(ξ(x))+f (3)(η(x))]h2

folgt die Genauigkeit O(h2) oder p = 2.I Ahnlich konnen wir die Koeffizienten a0,a1,a2 dieser

Approximation bestimmen,

f ′(x) = a0f (x) + a1f (x + h) + a2f (x + 2h) +O(hq)

fur ein moglichst hohes q.I Mit dem Taylorsatz,

f (x + h) = f (x) + f (1)(x)h + 12 f (2)(x)h2 + 1

6 f (3)(ξ(x))h3

f (x + 2h) = f (x) + f (1)(x)2h + 12 f (2)(x)4h2 + 1

6 f (3)(η(x))8h3

bilden wir die gewichtete Summe,

a0f (x) + a1f (x + h) + a2f (x + 2h) = (a0 + a1 + a2)f (x)

+(a1 + 2a2)f (1)(x)h + (a1 + 4a2)12 f (2)(x)h2

+a116 f (3)(ξ(x))h3 + a2

16 f (3)(η(x))8h3

198

Approximationen der AbleitungenI mit der ein System fur die Koeffizienten hergeleitet wird: 1 1 1

0 1 20 1 4

a0a1a2

=

0h−1

0

⇒ a0

a1a2

=1

2h

−34−1

I Aus

a1

6f (3)(ξ(x))h3+

a2

6f (3)(η(x))8h3 =

13

[f (3)(ξ(x))−2f (3)(η(x))]h2

folgt die Genauigkeit O(h2) oder q = 2.I Hausaufgabe: Leite ein 5× 5 System fur die Koeffizientena−2, . . . ,a2 her,

f ′′(x) =2∑

k=−2

ak f (x + kh) +O(hr )

um die Genauigkeit r = 4 zu erreichen.199

Richardson ExtrapolationI Angenommen es ist eine Formel N0(h) zur Approximation

einer Große M verfugbar. Man mochte die Genauigkeitbillig erhohen.

I Wenn giltM = N0(h) + a0h2 +O(h4)

dann folgt

M = N0(h/2) + a0(h/2)2 +O(h4)

Die gewichtete Summe ist

M =4N0(h/2)− N0(h)

4− 1︸︷︷︸=:N1(h)

+O(h4)

und N1(h) ist eine Approximation hoherer Genauigkeit.I Wenn allgemeiner gilt

M = N0(h) + a0hk0 + a1hk1 + · · ·200

Richardson Extrapolation

dann gilt

M = Nj(h) +O(hkj ), j = 1,2, . . .

wobei fur t > 1 (z.B. t = 2)

Ni+1(h) =tki Ni(h/t)− Ni(h)

tki − 1

I Beispiel: Mit

N0(h) =f (x0 + h)− f (x0 − h)

2h

gilt

f ′(x0) = N0(h)− h2

2!f (3)(x0)− h4

4!f (5)(x0) · · ·

Nimm f (x) = xex , x0 = 2, h = 0.2.201

Richardson ExtrapolationSchon bekannt: f ′(x0) = 22.167168 · · ·

I Mit Extrapolation,

N0(0.2) = [f (2.2)− f (1.8)]/[2(0.2)] = 22.414160N0(0.1) = [f (2.1)− f (1.9)]/[2(0.1)] = 22.228786

N0(0.05) = [f (2.05)− f (1.95)]/[2(0.05)] = 22.182564

konvergiert mit O(h2).I Erweiterte Tabelle:

N0(0.2) = 22.414160N0(0.1) = 22.228786 N1(0.2) = [4N0(0.1)− N0(0.2)]/3 = 22.166995

N0(0.05) = 22.182564 N1(0.1) = [4N0(0.05)− N0(0.1)]/3 = 22.167157

I Noch weiter:N0(0.2) = 22.414160N0(0.1) = 22.228786 N1(0.2) = 22.166995

N0(0.05) = 22.182564 N1(0.1) = 22.167157 N2(0.2) = [16N1(0.1)− N1(0.2)]/15= 22.167168

I Die Methode wird auch fur numerische Integrationverwendet.202

Numerische IntegrationI Die meisten Integrale konnen nicht analytisch

ausgerechnet werden, z.B.

erf(z) =2√π

∫ z

0e−t2

dt =?

I Integrale werden mit gewichteten Summen vonFunktionswerten approximiert:∫ b

af (x)dx ≈

n∑i=0

ai f (xi)

I Wie leitet man die Gewichte ai und die Quadraturpunktexi her?

I Von Interpolation bekommen wir eine Newton-CotesFormel:

f (x) ≈n∑

i=0

f (xi)Li(x) ⇒∫ b

af (x)dx ≈

n∑i=0

f (xi)

∫ b

aLi(x)dx︸︷︷︸=:aiaber xi =?

203

Trapez-RegelI Der Approximationsfehler

E(f ) =

∫ b

af (x)dx −

n∑i=0

f (xi)

∫ b

aLi(x)dx

kann mit dem Satz 142 so dargestellt werden:

E(f ) =1

(n + 1)!

∫ b

af (n+1)(ξ(x))

n∏i=0

(x − xi)dx

und wir mochten den Fehler abschatzen.I Beispiel: Trapez-Regel. Gegeben sind x0, x1, n = 1,

h = x1 − x0 und

f (x) ≈ f (x0)x − x1

x0 − x1+ f (x1)

x − x0

x1 − x0

Dann mit fi = f (xi) gilt∫ x1

x0

f (x)dx = f0∫ x1

x0

x − x1

x0 − x1dx + f1

∫ x1

x0

x − x0

x1 − x0dx + E(f )

204

Trapez-Regeloder∫ x1

x0

f (x)dx =f0

x0 − x1

(x − x1)2

2

∣∣∣∣x1

x0

+f1

x1 − x0

(x − x0)2

2

∣∣∣∣x1

x0

+ E(f )

=f0 + f1

2(x1 − x0) + E(f )

I Der Approximationsfehler ist:

E(f ) =12

∫ x1

x0

f (2)(ξ(x))(x − x0)(x − x1)dx

Nach dem Mittelwertsatz 11 fur Integrale ∃x ∈ [x0, x1]wobei

E(f ) =f (2)(ξ(x))

2

∫ x1

x0

(x − x0)(x − x1)dx = −h3

6f (2)(ξ(x))

I Fur die Trapez-Regel gilt∫ x1

x0

f (x)dx =f0 + f1

2(x1 − x0) +O(h3)

205

Simpson’s RegelI Beispiel: Simpson-Regel. Gegeben sind x0, x1, x2,

n = 2, x1 = (x0 + x2)/2, h = (x2 − x0)/2 und

f (x) ≈ f (x0)(x − x1)(x − x2)

(x0 − x1)(x0 − x2)+f (x1)

(x − x0)(x − x2)

(x1 − x0)(x1 − x2)+f (x2)

(x − x0)(x − x1)

(x2 − x0)(x2 − x1)

Dann mit fi = f (xi) gilt∫ x2

x0

f (x)dx =x2 − x0

6[f0 + 4f1 + f2] +O(h5)

I Bemerkung: Es gibt einen Sprung in Genauigkeit vonTrapez mit n = 1 (ungerade) und O(h3) zu Simpson mitn = 2 (gerade) und O(h5).

Def: Der Genauigkeitsgrad einer Integrationsformel ist n, wobeiE(P) = 0, ∀P ∈ Pn aber E(P) 6= 0 fur ein P ∈ Pn+1.

Def: Eine Newton-Cotes Formel heisst geschlossen wenna = x0 und b = xn gelten. Sonst heisst sie offen.206

Genauigkeit einer Newton-Cotes FormelSatz: Sei eine geschlossene Newton-Cotes Formel gegeben,∫ b

af (x)dx =

n∑i=0

ai f (xi) + E(f ), ai =

∫ b

aLi(x)dx

mit xi = a + ih, i = 0, . . . ,n, h = (b − a)/n. DannI Wenn n ungerade ist und f ∈ Cn+1([a,b]), ist der

Genauigkeitsgrad n und ∃ξ ∈ [a,b] 3

E(f ) = hn+2 f (n+1)(ξ)

(n + 1)!

∫ n

0t(t − 1) · · · (t − n)dt

I Wenn n gerade ist und f ∈ Cn+2([a,b]), ist derGenauigkeitsgrad n + 1 und ∃η ∈ [a,b] 3

E(f ) = hn+3 f (n+2)(η)

(n + 2)!

∫ n

0t2(t − 1) · · · (t − n)dt

207

Geschlossene Newton-Cotes FormelnI Beispiel: Es gibt einen Sprung in Genauigkeit von Trapez

mit n = 1 (ungerade) und O(h3) zu Simpson mit n = 2(gerade) und O(h5).

I Botschaft: Es ist vorteilhaft wenn n gerade ist.I Fur die folgenden Beispiele,

xi = a + ih, fi = f (xi), i = 0, . . . ,n, h = (b − a)/n.

und es gibt einen Sprung in Genauigkeit von n = 3 mitO(h5) zu n = 4 mit O(h5).

I Beispiel: n = 3,∫ x3

x0

f (x)dx =3h8

[f0 + 3f1 + 3f2 + f3] +O(h5)

I Beispiel: n = 4,∫ x4

x0

f (x)dx =2h45

[7f0 + 32f1 + 12f2 + 32f3 + 7f4] +O(h7)

208

Genauigkeit einer offenen Newton-Cotes FormelSatz: Sei eine offene Newton-Cotes Formel gegeben,∫ b

af (x)dx =

n∑i=0

ai f (xi) + E(f ), ai =

∫ b

aLi(x)dx

mit xi = a + (i + 1)h, i = 0, . . . ,n, h = (b − a)/(n + 2). DannI Wenn n ungerade ist und f ∈ Cn+1([a,b]), ist der

Genauigkeitsgrad n und ∃ξ ∈ [a,b] 3

E(f ) = hn+2 f (n+1)(ξ)

(n + 1)!

∫ n+1

−1t(t − 1) · · · (t − n)dt

I Wenn n gerade ist und f ∈ Cn+2([a,b]), ist derGenauigkeitsgrad n + 1 und ∃η ∈ [a,b] 3

E(f ) = hn+3 f (n+2)(η)

(n + 2)!

∫ n+1

−1t2(t − 1) · · · (t − n)dt

209

Mittelpunkt-RegelI Beispiel: Mittelpunkt-Regel, n = 0, h = (x1 − x−1)/2,∫ x1

x−1

f (x)dx = 2hf (x0) +O(h3)

I Direkt zeigen:

f (x) = f (x0) + (x − x0)f ′(x0) +12

(x − x0)2f ′′(ξ(x))

∫ x1

x−1

f (x)dx =

∫ x1

x−1

f (x0)dx︸︷︷︸=(x1−x−1)f (x0)

+f ′(x0)

∫ x1

x−1

(x − x0)dx︸︷︷︸=0

+f ′′(η)

2

∫ x1

x−1

(x − x0)2dx︸︷︷︸=2h3/3

210

Offene Newton-Cotes FormelnI Weitere offene Newton-Cotes Formeln:

n = 1 ∫ x2

x−1

f (x)dx =3h2

[f0 + f1] +O(h3)

n = 2 ∫ x3

x−1

f (x)dx =4h3

[2f0 − f1 + 2f2] +O(h5)

n = 3 ∫ x4

x−1

f (x)dx =5h24

[11f0 + f1 + f2 + 11f3] +O(h5)

I Hier gibt es einen Sprung in Genauigkeit von n = 1 mitO(h3) zu n = 2 mit O(h5).

I Auch keinen Sprung von n = 2 zu n = 3, die beiden mitGenauigkeit O(h5).

211

Zusammengesetzte Mittelpunkt-RegelFur Zusammengesetzte Integrationsregeln,I wird ein Intervall zerlegt:

[a,b] = ∪mj=1[zj−1, zj ]

I Das Integral wird dementsprechend zerlegt:∫ b

af (x)dx =

m∑j=1

∫ zj

zj−1

f (x)dx

I Eine Grundregel wird in jedem Teilintervall angewendet:∫ b

af (x)dx =

m∑j=1

∫ zj

zj−1

f (x)dx =m∑

j=1

[n∑

i=0

aij f (zj−1 + ih) +O(hp)

]Fur die Mittelpunkt-Regel,I Seien z0 = a, zm = b, zj = x2j , h = (b − a)/(2m),∫ b

af (x)dx =

m∑j=1

∫ x2j

x2j−2

f (x)dx =m∑

j=1

[2hf (x2j−1) +O(h3)

]212

Zusammengesetzte Simpsons-Regelnoder mit mh3 = h2(b − a)/2 = O(h2),∫ b

af (x)dx = 2h

m∑j=1

f (x2j−1) +O(h2)

Satz: E(f ) = O(h2) gilt fur die ZusammengesetzteMittelpunkt-Regel wenn f ∈ C2([a,b]).

Fur die Simpsons-Regel,I Seien z0 = a, zm = b, zj = x2j , h = (b − a)/(2m),∫ b

af (x)dx =

m∑j=1

∫ x2j

x2j−2

f (x)dx =m∑

j=1

[h3

(f2j−2 + 4f2j−1 + f2j) +O(h5)

]=

h3

m∑j=1

f2j−2 +4h3

m∑j=1

f2j−1 +h3

m∑j=1

f2j +O(mh5)

=h3

f0 +m−1∑j=1

f2j

+4h3

m∑j=1

f2j−1 +h3

m−1∑j=1

f2j + f2m

+O(h4)

213

Zusammengesetzte Trapez-Regelnoder∫ b

af (x)dx =

h3

f0 + 4m∑

j=1

f2j−1 + 2m−1∑j=1

f2j + f2m

+O(h4)

wobei O(mh5) = O(h4) folgt aus h = (b − a)/(2m) odermh = (b − a)/2.

Satz: E(f ) = O(h4) gilt fur die ZusammengesetzteSimpsons-Regel wenn f ∈ C4([a,b]).I Hausaufgabe: Leite die Zusammengesetzte Trapez-Regel

her:∫ b

af (x)dx =

h2

f (a) + 2m−1∑j=1

f (xj) + f (b)

+O(h2)

wobei xj = a + jh, j = 0, . . . ,m, h = (b − a)/m.Satz: E(f ) = O(h2) gilt fur die ZusammengesetzteTrapez-Regel wenn f ∈ C2([a,b]).214

Romberg IntegrationI Mit der zusammengesetzten Trapez-Regel,∫ b

af (x)dx = N0(h) +O(h2)

wobei h = (b − a)/m und

N0(h) =h2

f (a) + 2m−1∑j=1

f (a + jh) + f (b)

wird Richardson Extrapolation durchgefuhrt,

Nj(h) =4jNj−1(h/2)− Nj−1(h)

4j − 1, j = 1,2, . . .

I In der ersten Spalte der Richardson Tabelle sind WerteN0(h/2k ) : k = 0,1, . . . , wobei die Halfte derSummanden in N0(h/2k ) schon in N0(h/2k−1) dabei sind.

215

Romberg IntegrationI Diese erste Spalte kann durch die folgende Rekusion

gunstig berechnet werden:

N0(h/2) =h4

[f (a) + 2

2m−1∑j=1

f (a + jh2

) + f (b)

]

=h4

[f (a) + 2

m−1∑k=1

f (a + kh)︸︷︷︸j=2k

+2m∑

l=1

f (a + (l − 12)h)︸︷︷︸

j=2l−1

+f (b)

]

=12

[N0(h) + h

m∑l=1

f (a + (l − 12)h)

]I Im allgemeinen,

N0(h/2k ) =12

[N0(h/2k−1)+h/2k−1

m2k−1∑l=1

f (a+(l−12)h/2k−1)

]I Die Tabelle,N0(h)

N0(h/2) N1(h) = [4N0(h/2)− N0(h)]/3N0(h/4) N1(h/2) = [4N0(h/4)− N0(/2)]/3 N2(h) = [16N1(h/2)− N1(h)]/15· · ·216

Romberg IntegrationI Beispiel:

∫ π0 sin(x)dx = 2, f (x) = sin(x).

Nimm a = 0, b = π, m = 1, h = (b − a)/m = π.

k = 0N0(h) = h

2 [f (a) + f (b)] = π2 [sin(0) + sin(π)] = 0

k = 1N0(h/2) = 1

2 [N0(h) + h∑1

l=1 f (a + (l − 12)h︸︷︷︸

π/2

)] ≈ 1.571

k = 2N0(h/4) = 1

2 [N0(h/2) + h/2∑2

l=1 f (a + (l − 12)h/2︸︷︷︸

π/4,3π/4

)] ≈ 1.896

k = 3N0(h/8) = 1

2 [N0(h/4) + h/4∑4

l=1 f ( a + (l − 12)h/4︸︷︷︸

π/8,3π/8,5π/8,7π/8

)] ≈ 1.974

I Die Tabelle: 01.571 2.0941.896 2.004 1.9991.974 2.000 2.000 2.000· · ·217

Gauß QuadraturI Nun wird untersucht, ob man die Quadraturpunkte xi

vorteilhaft auswahlen kann.I Das Resultat: Der Genauigkeitsgrad einer Quadraturformel∫ b

af (x)dx ≈

n∑i=0

ai f (xi)

ist maximal bei 2n + 1, wenn die Quadraturpunkte xierfullen Pn+1(xi) = 0, wobei Pn+1 ∈ Pn+1 zu einer Gruppevon orthogonalen Polynomen gehort.

I Sei φknk=0 eine Menge von Polynomen mit den folgendenEigenschaften:• φk ∈ Pk

• ∃a,b ∈ R, w = w(x) wobei w(x) ≥ 0, ∀x ∈ [a,b] und∫ b

aw(x)φk (x)φl(x)dx = δkl , 0 ≤ k , l ≤ n

Diese Polynome sind orthogonal auf [a,b] bezuglich desGewichts w .218

Gauß Quadratur

Satz: Sei φknk=0 orthogonale Polynome auf [a,b] bezuglichdes Gewichts w . Dann hat φk genau k Nullstellen in (a,b).

I Beispiel: [a,b] = [−1,1], w(x) = 1, die LegendrePolynome:

P0(x) = 1, P1(x) = x , P2(x) = 12(3x2 − 1)

und im allgemeinen,

Pk (x) =1

2kn!

dk

dxk [(x2 − 1)n]

I Beispiel: [a,b] = [−1,1], w(x) = 1/√

1− x2, dieTchebyshev Polynome:

T0(x) = 1, T1(x) = x , T2(x) = 2x2 − 1und im allgemeinen,

Tk (x) = cos[k cos−1(x)]

219

Gauß Quadratur

Def: Die Nullstellen xini=0 des Legendre Polynoms Pn+1 sinddie Gauß-Legendre Punkte.

Satz: Wenn xini=0 die Gauß-Legendre Punkte sind, ist derGenauigkeitsgrad der Quadraturformel∫ 1

−1f (x)dx ≈

n∑i=0

ai f (xi), ai =

∫ b

aLi(x)dx

maximal bei 2n + 1.

Bemerkung: Dieser Genauigkeitsgrad ist maximal weil 2n + 2Freiheitsgrade, aini=0 und xini=0, in der Quadraturformel zuverfugung stehen, und es gibt 2n + 2 unabhangigeBedingungen in∫ 1

−1P(x)dx =

n∑i=0

aiP(xi), ∀P ∈ P2n+1

weil alle P ∈ P2n+1 2n + 2 Koeffizienten haben.220

Gauß QuadraturBemerkung: Um diese Quadraturformel fur ein allgemeinesIntegral anzuwenden,∫ b

af (t)dt =

∫ 1

−1g(x)dx ≈

n∑i=0

aig(xi)

muss man eine Koordinatentransformation durchfuhren:t = a + (b − a)(x + 1)/2

dt = dx(b − a)/2g(x) = f (a + (b − a)(x + 1)/2)(b − a)/2

I Beispiel:∫ 1.5

1e−t2

dt = 0.1093643

Newton-Cotes:

n: 0 1 2 3 4geschl: − 0.1183197 0.1093104 0.1093404 0.1093643offen: 0.1048057 0.1063473 0.1094116 0.1093971 −

221

Gauß Quadratur

Romberg-Trapez:0.11831970.1115627 0.10931040.1099114 0.1093610 0.10936430.1095009 0.1093641 0.1093643 0.1093643

Gauß-Legendre:∫ 1.5

1e−t2

dt =

∫ 1

−1g(x)dx , g(x) = 1

4 exp[− (x+5)2

16 ]

n = 1 x0 = −0.5773502692 a0 = 1x1 = +0.5773502692 a1 = 1

1∑i=0

aig(xi) = 0.1094003

n = 2 x0 = −0.7745966692 a0 = 5/9x1 = 0 a1 = 8/9x2 = +0.7745966692 a2 = 5/9

2∑i=0

aig(xi) = 0.1093642

222

Gewohnliche Differentialgleichungen

Zu losen ist das Anfangswertproblem: Finde y : R→ Rn 3

(?)

y ′(t) = f (t ,y(t)), t ∈ [t0,T ]y(t0) = y0

Satz: Sei f ∈ C(D,Rn), D = [t0,T ]×Ω, Ω ⊂ Rn. Wenn ∃L > 0 3

‖f (t ,y2)− f (t ,y1)‖ ≤ L‖y2 − y1‖, ∀(t ,y1), (t ,y2) ∈ D

dann ∃! Losung y zum Anfangswertproblem (?).

I Analytische Losung? Selten.I Diskretisierung: Zur Vereinfachung nimm ein regelmaßiges

Gitter in Zeit:

tkMk=0 ⊂ [t0,T ], tk = t0+kτ, τ = (T−t0)/M, k = 0,1, . . . ,M

223

Das Eulersche VerfahrenI Wegen dem Taylorsatz ∃ξk ∈ (tk , tk+1) 3

y(tk+1) = y(tk ) + (tk+1 − tk )y ′(tk ) + 12(tk+1 − tk )2y ′′(ξk )

oder mit τ = tk+1 − tk ,

f (tk ,y(tk )) = y ′(tk ) =y(tk+1)− y(tk )

τ+O(τ)

I Das Vorwarts Euler Verfahren ist:

yk+1 = yk + τ f (tk ,yk ), k = 0,1, . . . ,M − 1

wobei yk ≈ y(tk ).I Beispiel: f (t , y) = −y ,

y0 = 1, t0 = 0, T > 0,

y ′ = −y , t ∈ [0,T ]

y(0) = 1

⇒ y(t) = e−t .

224

Das Eulersche VerfahrenSatz: Nimm an, f erfullt die Bedingungen des Satzes 223 , unddie eindeutige Losung y zum (?) erfullt ‖y ′′(t)‖ ≤ B,∀t ∈ [t0,T ]. Dann erfullt das Vorwarts Euler Verfahren:

‖y(tk )− yk‖ ≤τB2L

[exp(L(tk − t0))− 1]

und dahermax

0≤k≤M‖y(tk )− yk‖ ≤

τB2L

[exp(L(T − t0))− 1] = O(τ)

I DasVorwarts Euler Verfahren konvergiert,yk+1 = yk + τ f (tk ,yk ), k = 0,1, . . . ,M − 1

aber man sieht durch die nachsten Beispiele, dass dieZeitschrittweite τ oft sehr klein sein muss, aber

I das Ruckwarts Euler Verfahren (auch O(τ)) ist vorteilhaft:yk+1 = yk + τ f (tk+1,yk+1), k = 0,1, . . . ,M − 1

I Die Vorwarts Methode ist explizit, weil yk+1 bezuglichschon berechneter Großen explizit gegeben ist.

I Die Ruckwarts Methode ist implizit, weil yk+1 nur implizitdefiniert ist. Sie kann eine iterative Losung kosten!225

WarmegleichungI Beispiel: Eine raumliche Diskretisierung der

Warmegleichung,ut = κuxx , Ω× (0,∞)u′ = 0, ∂Ω× (0,∞)u = u0, Ω× t = 0

Ω = (0,1)

ist mit ui(t) ≈ u(xi , t), xi = ih, i = 0, . . . ,N, h = 1/N:

u′ = Au

u(0) = u0A = − κ

h2 D, D =

1 −1 0−1 2 −1

. . .. . .

. . .

−1 2 −10 −1 1

I Die Losung u erfullt:

limt→∞

u(t) = Mittelwert(u0)e, e = 〈1,1 . . . ,1〉T

Laut Gerschgorin gilt σ(A) ≤ 0, und daher gilt12

ddt‖u‖2 = uTu′ = uTAu ≤ 0 ⇒ ‖u(t)‖2 ≤ ‖u(0)‖2

226

WarmegleichungI Vorwarts Euler:

uk = uk−1+τAuk−1 =(

I − κτ

h2 D)

uk−1 = · · · =(

I − κτ

h2 D)k

u0

I Um Instabilitat zu vermeiden, wollen wirσ(I − κτ

h2 D) ⊂ B(0,1).Die Gerschgorin Scheiben in C sind B(zi , ri) wobei

z0 = 1−κτh2 , r0 =

κτ

h2 ; zi = 1−2κτh2 , ri =

2κτh2 ; zN = 1−κτ

h2 , rN =κτ

h2

I Verlangt wird:

−1 < 1− 2κτh2 < 1 und − 1 < 1− 4κτ

h2 < 1

oder die umstandliche Bedingung,

1N2 = h2 > 2κτ = 2κ

T − t0M

⇒ M = O(N2)

227

Konservative VerfahrenI Ruckwarts Euler:

uk = uk−1 + τAuk ⇒(

I +κτ

h2 D)

uk = uk−1

Laut Gerschgorin gilt σ(D) ≥ 0, und daher

‖uk‖2 = uTk uk ≤ uT

k

(I +

κτ

h2 D)

uk = uTk uk−1 ≤ ‖uk‖‖uk−1‖

I Es folgt ‖uk‖ ≤ ‖uk−1‖ ohne Bedingungen auf τ !I Wenn z.B. BT = −B gilt, folgtu′ = Bu

u(0) = u0⇒ 1

2ddt‖u‖2 = uTBu =

12(uTBu + uTBTu

)= 0

und ‖u(t)‖ bleibt erhalten bei ‖u(0)‖.

Fur ein solches Problem soll das numerische Verfahrenkonservativ sein:

‖uk‖ = ‖uk−1‖.228

WellengleichungI Beispiel: Eine raumliche Diskretisierung der

Wellengleichung,utt = ν2uxx , Ω× (0,∞)u = 0, ∂Ω× (0,∞)u = u0, Ω× t = 0ut = u1, Ω× t = 0

(νuxut

)t

=

(0 ν∂xν∂x 0

)(νuxut

)

ist mit xi = ih, i = 0, . . . ,N, h = 1/N,

Ui(t) ≈ ν∂xu(xi+ 12, t) i = 0, . . . ,N − 1

UN+i(t) ≈ ∂tu(xi , t) i = 0, . . . ,N

(und ahnlich U0 ≈ (ν∂xu0,u1)) gegeben durch:

U ′ = BU

U(0) = U0B =

ν

h

(0 Q−QT 0

)Q =

−1 +1 0. . .

. . .

0 −1 +1

I Die Losung U erfullt: ‖U(t)‖ = ‖U(0)‖.

229

Crank-Nicholson VerfahrenI Das Crank-Nicholson Verfahren fur y ′ = f (t ,y) ist:

yk+1 − yk

τ=

f (tk ,yk ) + f (tk+1,yk+1)

2(‖y(tk )−yk‖ = O(τ2))

I Fur die raumlich diskretisierte Wellengleichung,

Uk+1 − Uk

τ=

BUk + BUk+1

2⇒

[I − τ

2B]

Uk+1 =[I +

τ

2B]

Uk

I Wegen BT = −B gelten

UTBU = 12(UTBU + UTBTU) = 0

UTk+1Uk+1 = UT

k+1[I − τ2 B]Uk+1 = UT

k+1[I + τ2 B]Uk

= UTk+1Uk + τ

2 UTk+1BUk

und

UTk Uk = UT

k [I + τ2 B]Uk = UT

k [I − τ2 B]Uk+1

= UTk Uk+1 − τ

2 UTk BUk+1 = UT

k+1Uk − τ2 UT

k+1BTUk

und daher gilt ‖Uk+1‖ = ‖Uk‖.230

Runge-Kutta MethodenI Zur numerischen Losung des Anfangswertproblems,

(?)

y ′(t) = f (t ,y(t)), t ∈ [t0,T ]y(t0) = y0

ist eine Runge-Kutte Methode durch das System definiert:yk ,i = yk + τ

q∑j=1

aij f (tk + θjτ,yk ,j)

yk+1 = yk + τ

q∑i=1

bi f (tk + θiτ,yk ,i)

I Die Methode wird mit einer Tabelle dargestellt:

A θ

bT wobei A ∈ Rq×q,b,θ ∈ Rq

I Die Methode heisst explizit wenn gilt aij = 0, i ≤ j , und sieheisst implizit wenn ∃i 3 aii 6= 0.

231

Runge-Kutta MethodenI Eine wohl bekannte RK-Methode ist:

A θ

bT →

0 0 0 0 012 0 0 0 1

20 1

2 0 0 12

0 0 1 0 116

13

13

16

oder

yk ,1 = ykyk ,2 = yk + τ

2 f (tk ,yk ,1)

yk ,3 = yk + τ2 f (tk + τ

2 ,yk ,2)

yk ,4 = yk + τ f (tk + τ2 ,yk ,3)

yk+1 = yk + τ [16 f (tk ,yk ,1) + 1

3 f (tk + τ2 ,yk ,2)

+13 f (tk + τ

2 ,yk ,3) + 16 f (tk + τ,yk ,4)]

Die Methode erfullt ‖y(tk )− yk‖ = O(τ4) wenn y ∈ C(5)([t0,T ]).232

Runge-Kutta-Fehlberg MethodenI Nun verwenden wir kein regelmaßiges Gitter mehr.I Die Zeitschrittweite wird automatisch ausgewahlt.I Eine grobere RK-Methode (A) wird in eine genauere

RK-Methode (A) eingebettet:

A =

0 0...

. . .. . .

· · · 0

A =

[A 0

aT+ 0

]θ =

[θθ+

]b = bi 63 bi = b

I Die grobere MethodeA θ

bT kann so dargestellt werden:

yk+1 = yk + τφ(tk ,yk , τ)

wobei φ die gesamte Abbildung (tk ,yk , τ)→ [yk+1 − yk ]/τdarstellt.

233

Runge-Kutta-Fehlberg MethodenI Ahnlich kann die genauere Methode

A θ

bT so dargestellt

werden:

yk+1 = yk + τ φ(tk ,yk , τ)

aber mit nur wenig zusatzlicher Arbeit bekommt manφ(tk ,yk , τ) wenn φ(tk ,yk , τ) verfugbar ist.

I yk+1 und yk+1 sollen verwendet werden, um dasgeeignete τ zu bestimmen.

I Seien die jeweiligen Abbruchsfehler definiert:

Tk+1(τ) = [y(tk+1)− y(tk )]/τ − φ(tk ,y(tk ), τ) = O(τn)

Tk+1(τ) = [y(tk+1)− y(tk )]/τ − φ(tk ,y(tk ), τ) = O(τn+1)

I yk ≈ y(tk )⇒y(tk+1)− yk+1 = y(tk+1)− yk − τφ(tk ,yk , τ)

≈ y(tk+1)− y(tk )− τφ(tk ,y(tk ), τ)= τTk+1(τ)

I Ahnlich gilt y(tk+1)− yk+1 ≈ τ Tk+1(τ).234

Runge-Kutta-Fehlberg MethodenI Die Abbruchsfehler konnen so kombiniert werden:

O(τn) = Tk+1(τ) ≈y(tk+1)− yk+1

τ

=y(tk+1)− yk+1

τ+

yk+1 − yk+1

τ

≈ Tk+1(τ)︸︷︷︸O(τn+1)

+yk+1 − yk+1

τ︸︷︷︸signifikanter

I MitCτn ≈ Tk+1(τ) ≈

yk+1 − yk+1

τund

Tk+1(qτ) ≈ C(qτ)n = qnCτn ≈ qn

τ[yk+1 − yk+1]

wahlen wir q so aus,

qn|yk+1 − yk+1| ≤ ετ

um die Genauigkeit |Tk+1(qτ)| ≤ ε zu erreichen.235

Runge-Kutta-Fehlberg MethodenDie RKF-Strategie:I Die Fehlertoleranz ε wird eingegeben.I Ein Max τmax und ein Min τmin fur die Zeitschrittweite

werden eingegeben, und τ0 = τmax.I Fur k ≥ 0 wird ein τk vom letzten Zeitschritt gegeben und

verwendet, um yk+1 und yk+1 zu berechnen.I Wenn |yk+1 − yk+1| ≤ ετ gilt, akzeptiere dieses τk .

(Falls |yk+1 − yk+1| << ετ , τk vergroßern?)I Wenn |yk+1 − yk+1| ≤ ετ nicht gilt, nimm das neue

τk ← min(τmax,max(τmin,qτk ))

wobei

q ≤(

ετk

|yk+1 − yk+1|

)1/n

und berechne yk+1 und yk+1 frisch.I Setzte τk+1 = τk fur den nachsten Zeitschritt.

236

Konsequenz, Stabilitat, KonvergenzI Fur das Anfangswertproblem,

(?)

y ′(t) = f (t ,y(t)), t ∈ [t0,T ]y(t0) = y0

sei die Einzelschritt-Methode gegeben:


wobei φ die gesamte Abbildung (tk ,yk , τ)→ [yk+1 − yk ]/τdarstellt.

Def: Die Methode heisst konsequent mit (?) wenn

φ(t ,y ,0) = f (t ,y)

Def: Die Genauigkeitsordnung ist n wenn

y(t + τ)− y(t)τ

− φ(t ,y(t), τ) = O(τn)

237

Konsequenz, Stabilitat, KonvergenzI Beispiel: Das Vorwarts Euler Verfahren hat

Genauigkeitsordnung n = 1:φ(t ,y , τ) = f (t ,y)

undy(t + τ)− y(t)

τ−φ(t ,y(t , τ))︸︷︷︸

f (t ,y(t))=y ′(t)

=y(t + τ)− y(t)

τ−y ′(t) = O(τ)

I Hausaufgabe: Zeige, die Genauigkeitsordnung vomCrank-Nicholson Verfahren ist n = 2.

Def: Mit der Einzelschritt-Methode, yk+1 = yk + τφ(tk ,yk , τ),τ = (T − t0)/M, seien ykMk=0 und zkMk=0 gegeben mitAnfangswerten y0 beziehungsweise z0. Die Methode heisststabil wenn ∃C, τ0 > 0 wobei ∀τ ∈ (0, τ0),

max0≤k≤M

‖yk − zk‖ ≤ C‖y0 − z0‖.

238

Konsequenz, Stabilitat, KonvergenzSatz: Fur das Anfangswertproblem (?) sei eineEinzelschritt-Methode gegeben:


Wenn ∃τ0 > 0 sodass φ auf D = [t0,T ]× Rn × [0, τ0] erfullt

‖φ(t ,y1, τ)−φ(t ,y2, τ)‖ ≤ L‖y1−y2‖, ∀(t ,y1, τ), (t ,y2, τ) ∈ D

dann:1. Die Methode ist stabil.2. Die Methode ist konvergent genau dann wenn sie

konsequent ist.3. Wenn ∃T (τ) sodass

Tk (τ) =y(tk )− y(tk−1)

τ−f (tk−1,y(tk−1)) erfullt max

1≤k≤M|Tk (τ)| ≤ T (τ)

dann folgt

‖y(tk )− yk‖ ≤T (τ)

Lexp[L(tk − t0)], k = 1,2, . . . ,M.

239

Bedeutung der StabilitatFur das Anfangswertproblem (?) seien ykMk=0 und zkMk=0gegeben durch das Vorwarts Euler Verfahren mitAnfangswerten y0 beziehungsweise z0. Es gelten:

[yk − zk ] = [yk−1 − zk−1] + τ [f (tk−1,yk−1)− f (tk−1, zk−1)]

und

‖yk − zk‖ ≤ ‖yk−1 − zk−1‖+ τL‖yk−1 − zk−1‖≤ (1 + Lτ)k‖y0 − z0‖ ≤ (1 + Lτ)M‖y0 − z0‖= (1 + L(T − t0)/M)M‖y0 − z0‖≤ exp(L(T − t0))‖y0 − z0‖

und das Eulersche Verfahren ist stabil weil ∃C > 0 3

max0≤k≤M

‖yk − zk‖ ≤ C‖y0 − z0‖

aber die Konstante C = exp((T − t0)L) kann riesig sein.240

Bedeutung der Stabilitat

I Beispiel: [t0,T ] = [0,1], y0 = 0 undf (t , y) = −1000y + 1000t2 + 2t

Die genaue Losung ist y(t) = t2, t ∈ [0,1]. Sei zkMk=0 dienumerische Losung mit der Anfangsstorung z0 ≈ y0.Es gelten:

[yk − zk ] = [yk−1 − zk−1]− 1000τ [yk−1 − zk−1]

|yk − zk | ≤ (1− 1000τ)k |y0 − z0| ≤ (1− 1000τ)M |y0 − z0|

I Mit τ < 0.002 gilt yM ≈ 1 = y(tM).I Mit τ = 0.01 werden Fehler bei jedem Schritt mit einem

Faktor |1− 1000 · 0.01| = 9 vergroßert.I Die obige Fehler-Abschatzung wird nicht verletzt, aber

C = e1000.

Solche Beispiele motivieren die folgende Definition.241

Absolute StabilitatDef: Ein numerisches Verfahren zur Losung von (?) heisstabsolut stabil, wenn die numerische Losung yk des Problemsy ′ = λy , <λ < 0, erfullt limk→∞ yk = 0. Die Menge der Werteτλ, in der eine Methode absolut stabil ist, heisst die Mengeder absoluten Stabilitat.

I Die Menge der absoluten Stabilitat fur das VorwartsEulersche Verfahren ist z ∈ C : |z + 1| < 1:

yk = yk−1+τλyk−1 = · · · = (1+τλ)ky0

und

ykk→∞−→ 0 ⇔ |1 + τλ| < 1

oder mit x = <[τλ], y = =[τλ],

[(x + 1)2 + (y − 0)2] < 1

C

·&%'$

-1

242

A-StabilitatI Die Menge der absoluten Stabilitat fur das Ruckwarts

Eulersche Verfahren ist z ∈ C : <z < 0:

yk+1 = yk +τλyk+1 ⇒ (1−τλ)yk+1 = yk ⇒ yk+1 = yk/(1−τλ)

und

yk = y0/(1−tλ)k k→∞−→ 0 ⇔ |1−tλ|−1 < 1 ⇔ 1 < |1−tλ|

oder mit x = <[τλ], y = =[τλ],

x < 0 und1 < [(x − 1)2 + (y − 0)2]

⇒ x < 0

C

·&%'$

+1

I Wenn die Menge der absolutenStabilitat so maximal ist,heisst ein Verfahren A-Stabil.

I Solche Stabilitat ist geeignet fur folgende Probleme.243

Steife SystemeI Ein Prozess heisst steif wenn mindestens zwei sehr

unterschiedliche Zeitskalen dabei sind.I Physikalisches Beispiel: Stromung einer Flussigkeit, in der

chemische Reaktionen stattfinden. Es ist ublich, dass dieReaktionen viel schneller als die Stromung laufen.

I Mathematisches Beispiel:

x ′(t) = Ax(t), AS = SΛ, Λ = diagλ1, λ2

λ1, λ2 < 0 aber |λ1| >> |λ2|y(t) = S−1x(t) erfullt

y ′(t) = S−1x ′(t) = S−1Ax(t) = S−1ASy(t) = Λy(t)

y(t) = exp(Λt)y(0)⇒ x(t) = S exp(Λt)S−1x(0) oder

x1(t) = Aeλ1t + Beλ2t , x2(t) = Ceλ1t + Deλ2t

Jedes xi(t) enthalt eine Komponente eλ1t , die viel schnellerzerfallt als die Komponente eλ2t .

244

Steife SystemeMit S =

1√2

[1 1−1 1

], Λ =

[−10 0

0 −1

], x(0) =

[01

],

gilt x1(t) = (e−t − e−10t )/2:

x1 hat eine schnellereKomponente e−10t undeine langsamere Kom-ponents e−t .

Die numerischen Ergebnisse mit [t0,T ] = [0,1], M = 7,27:

I A-stabile Methoden verlangen auch ein kleines τ furGenauigkeit in schnellen Bereichen eines steifen Problems.245

RandwertproblemeI Zu losen ist das Sturm-Liouville Randwertproblem:

(‡)−(pu′)′ + qu = f , x ∈ Ωαu + ωDnu = g, x ∈ ∂Ω

Ω = (a,b)

wobei Dnu = −u′|x=a,u′|x=b und

p = p(x) ≥ p0 > 0, q = q(x) ≥ q0 > 0

I Fur verschiedene Werte von α und ω hat man folgendeArte von Randbedingungen:

ω = 0 ⇒ Dirichlet RBα = 0 ⇒ Neumann RB

α, ω 6= 0 ⇒ Robin RBI Ergibt sich durch Minimierung des Funktionals,

J(u) =

∫Ω

[(mu − r)2 + p|u′|2]dx

wobei q = m2 und f = mr , und die Randbedingung ist eineNebenbedingung.

246

Randwertprobleme

I Das Problem (‡) kann durch Finite Differenzenfolgendermaßen diskretisiert werden.

I Sei das Gitter: xi = a + ih, i = 0, . . . ,N, h = (b − a)/N,

x0

a

xN

bh

xi− 12

xi+ 12

I Mit der Notation vi ≈ v(xi), vi± 12≈ v(xi± 1

2), wird die

Differentialgleichung in (‡) so diskretisiert:(pu′)′i ≈

(pu′)i+ 12− (pu′)i− 1

2

h≈ 1

h

[pi+ 1

2

ui+1 − ui

h− pi− 1

2

ui − ui−1

h

](qu)i ≈ qiui

247

RandwertproblemeI Mit der Neumann Randbedingung u′ = 0, ∂Ω gelten

(pu′)′0 ≈ p0+ 12

u1 − u0

h2 − p0− 12

[u0 − u−1]=0

h2 = p 12

u1 − u0

h2

(pu′)′N ≈ pN+ 12

[uN+1 − uN ]=0

h2 − pN− 12

uN − uN−1

h2 = −pN− 12

uN − uN−1

h2

I Das algebraische System zur numerischen Losung von (‡)ist:

h−2p 12

+q0 −h−2p 12

−h−2p 12

h−2[p 12

+p 32

]+q1 −h−2p 32

. . .. . .

. . .

−h−2pN− 3

2h−2[p

N− 32

+pN− 1

2]+qN−1 −h−2p

N− 12

−h−2pN− 1

2h−2p

N− 12

+qN

u0.........

uN

=

f0.........

fN

I Hausaufgabe: Schreibe das System fur die DirichletRandbedingung u = 0, ∂Ω, und fur die RobinRandbedingung σu + Dnu = 0, ∂Ω.

248

Rayleigh-Ritz Methode

I Das Randwertproblem (‡) ergibt sich durch Minimierungdas Funktionals,

J(u) =

∫Ω

[(mu − r)2 + p|u′|2]dx

I Die Richtungsableitungen des Funktionals an der Stelle uin die Richtung einer Storung v sind gegeben durch:

δJδu

(u; v) = limε→0

J(u + εv)− J(u)

ε

= 2∫

Ω[(mu − r)mv + pu′v ′]dx

= 2∫

Ω[(mu − r)m − (pu′)′]vdx + 2pu′v |∂Ω

I Wenn J(u) = min gilt, muss δJ/δu(u; v) = 0 fur alleStorungen v gelten.

249

Rayleigh-Ritz MethodeI Fur ein beliebiges B(x0, δ) ⊂ Ω, sei

vδ0 (x) =

1/(2δ), x ∈ B(x0, δ)

0, sonstNach dem Mittelwertsatz fur Integrale 11 ∃ξ ∈ B(x0, δ) sodass:

0 =12δJδu

(u; vδ0 ) = [(mu−r)m−(pu′)′](ξ)

∫Ω

vδ0dx︸︷︷︸=1

δ→0−→ [(mu−r)m−(pu′)′](x0)

oder m2u − rm − (pu′)′ = 0, Ω, wobei q = m2 und f = mr ,d.h. die Differentialgleichung in (‡) ergibt sich.

I Fur ein beliebiges x0 ∈ ∂Ω, sei

vδ0 (x) =

1, x ∈ B(x0, δ) ∩ Ω0, sonst

Nach dem Mittelwertsatz fur Integrale ∃ξ ∈ B(x0, δ) ∩ Ω sodass

0 =12δJδu

(u; vδ0 ) = [(mu−r)m−(pu′)′](ξ)

∫Ω

vδ0dx︸︷︷︸=δ/2

+pu′v |∂Ωδ→0−→ [pu′](x0)

oder u′ = 0, ∂Ω, d.h.die naturliche Neumann Randbedingung in (‡) ergibt sich.250

Rayleigh-Ritz MethodeI Wenn andere Randbedingungen, z.B. Dirichlet,

vorgegeben werden, darf es keine Storung am Randgeben, und so muss v = 0, ∂Ω gelten. Dann ist dieRechnung des letzten obigen Integrals hinfallig.

I Nachdem Ableitungen 1. Ordnung der Test-Funktionen vauf Ableitungen 2. Ordnung der Losung u verschobenwerden, ergibt sich die starke Formulierung (‡) derOptimalitatsbedingung fur J, fur die hohere Regularitat vonu angenommen wird.

I Wenn die gleiche Glattheit fur Test-Funktionen und fur dasminimierende u verwendet wird, ergibt sich die schwacheFormulierung der Optimalitatsbedingung fur J:∫

Ω[quv + pu′v ′]dx =

∫Ω

fvdx , ∀v ∈ H1(Ω)

wobei q = m2 und f = mr undH1(Ω) = v ∈ L2(Ω) : v ′ ∈ L2(Ω), L2(Ω) = v :

∫Ω |v |

2dx <∞.251

Rayleigh-Ritz MethodeI Diskretisierung der schwachen Formulierung: Mit dem

Gitter x = xi, xi = a + hi , i = 0, . . . ,N und h = (b − a)/Nwerden die Splines S1(x) als Basis Funktionen verwendet:

uh(x) =N∑

i=0

aisi(x)

I Die diskrete schwache Formulierung von (‡) ist∫Ω

[quhs + pu′hs′]dx =

∫Ω

fsdx , ∀s ∈ S1(x)

252

Rayleigh-Ritz MethodeI Die diskrete schwache Formulierung kann explizit

umgeschrieben werden, um die Koeffizienten ai derSpline-Darstellung von u zu bestimmen:

N∑i=0

ai

∫Ω

[qsisj + ps′i s′j ]dx =

∫Ω

fsjdx , j = 0, . . . ,N

I Definiere die Matrix und die Vektoren,

A =

∫Ω

[qsisj + ps′i s′j ]dx

N

i,j=0, a = aiNi=0, b =

∫Ω

fsjdx

und die Koeffizienten werden durch Losung des linearenGleichungssystems Aa = b bestimmt.

I Bespiel: Das Gleichungssystem wird explizit fur denfolgenden Fall bestimmt: Ω = (0,1), p = 1 = q,f = fi−1 in [xi−1, xi), i = 1, . . . ,N.

253

Rayleigh-Ritz MethodeI Zuerst merkt man,∫ 1

0sisjdx =

∫ 1

0s′i s′j dx = 0, |i − j | > 1

I Fur |i − j | = 1,∫ 1

0sisjdx =

∫ h

0s0s1dx =

∫ h

0

(1− x

h

)(xh

)dx =

h6

und∫ 1

0s′i s′j dx =

∫ h

0s′0s′1dx =

∫ h

0

(−1

h

)(+

1h

)dx = −1

hI Fur |i − j | = 0,∫ 1

0s2

i dx =

∫ 2h

0s2

1dx =

∫ h

0

(xh

)2dx+

∫ 2h

h

(2− x

h

)2dx =

2h3

und∫ 1

0|s′i |

2dx =

∫ 2h

0|s′1|dx =

∫ h

0

(1h

)2

dx+

∫ 2h

h

(−1

h

)2

dx =2h254

Rayleigh-Ritz MethodeI Das Integral von f ist:∫ 1

0fsjdx =

N∑i=1

∫ xi

xi−1

fisjdx = fj−1

∫ xj

xj−1

sjdx + fj∫ xj+1

xj

sjdx

= fj−1

∫ h

0s1dx + fj

∫ h

0s0dx

= fj−1

∫ h

0

xh

dx + fj∫ h

0

(1− x

h

)dx = h

fj−1 + fj2

fur 1 ≤ j ≤ N − 1.I Fur j = 0, ∫ 1

0fs0dx = f0

∫ h

0

(1− x

h

)= h

f02

I Ahnlich fur j = N,∫ 1

0fsNdx = fN−1

∫ h

0s1 = fN−1

∫ h

0

(xh

)= h

fN−1

2255

Rayleigh-Ritz Methode

I Das Gleichungssystem ist Aa = b mit:

A =

1h + h

3 −1h + h

6−1

h + h6

2h + 2h

3 −1h + h

6. . .

. . .. . .

−1h + h

62h + 2h

3 −1h + h

6−1

h + h6

1h + h

3

=1h

D+hK

wobei

D =

1 −1 0−1 2 −1

. . .. . .

. . .

−1 2 −10 0 −1 1

K =16

2 1 01 4 1

. . .. . .

. . .

1 4 10 0 1 2

256

Rayleigh-Ritz Methodeund

a =

a0...

aN

b = h

f02

f0+f12...

fN−2+fN−12

fN−12

I Hausaufgabe: Leite das Gleichungssystem fur den

folgenden Fall her: Ω = (0,1) undp = pi−1, q = qi−1, f = fi−1 in [xi−1, xi), i = 1, . . . ,N.

I Bemerke: Fur allgemeine Funktionen p(x), q(x) und f (x),kann man Gauß Quadratur fur die Integrale verwenden.

257

Grundlagen der Numerischen Mathematik · InhaltsverzeichnisII Pseudo-Kode mit Skalierter Pivotsuche...

Documents

Transcript of Grundlagen der Numerischen Mathematik · InhaltsverzeichnisII Pseudo-Kode mit Skalierter Pivotsuche...