Numerik 1 Eine Zusammenfassung Steffi, Balthasar, Julian und Julian WS 2009/2010 - M. Bause

Contents 1 Fehleranalyse 1.1 Kondition eines Problems . . . . . . . 1.1.1 Definition einer Norm . . . . . 1.1.2 Relative & absolute Kondition 1.1.3 Fehler der Ausgabe . . . . . . . 1.1.4 Lineare Abbildungen Satz 2.19 1.2 Rundungs- Maschinenegenauigkeit . . 1.2.1 relative Rundungsfehler . . . . 1.3 Stabilitaet eines Algorithmus . . . . . 1.3.1 Leitlinien . . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

3 3 3 3 3 3 4 4 4 4

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

2 Direkte L¨ osungsverfahren f¨ ur lineare Gleichungssysteme 2.1 Zeilenskalierung . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Dreiecksmatrizen . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Vorw¨ arts- und R¨ uckw¨artseinsetzen . . . . . . . . . . 2.3 Gauß-Elimination und LR-Zerlegung . . . . . . . . . . . . . 2.3.1 Gauß-Elimination ohne Pivotisierung . . . . . . . . . 2.3.2 LR-Zerlegung . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Gauß-Elimination mit Spaltenpivotisierung . . . . . 2.3.4 Anwendungen der LR-Zerlegung . . . . . . . . . . . 2.3.5 Cholesky-Zerlegung . . . . . . . . . . . . . . . . . . 2.4 QR-Zerlegung . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Givens-Rotationen . . . . . . . . . . . . . . . . . . . 2.4.2 Householder-Transformationen . . . . . . . . . . . .

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

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

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

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

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

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

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

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

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

5 . 5 . 5 . 6 . 6 . 6 . 7 . 7 . 7 . 8 . 9 . 9 . 10

3 Lineare Ausgleichsrechnung 3.1 Kondition des Problems . . . . . . . . . . 3.2 Numerische L¨ osung . . . . . . . . . . . . . 3.2.1 L¨ osung u ¨ ber Normalengleichungen 3.2.2 L¨ osung u ¨ ber QR-Zerlegung . . . . 3.3 Singul¨ arwertzerlegung . . . . . . . . . . . 3.3.1 Berechnung der Singul¨ arwerte . . . 3.3.2 Pseudoinverse . . . . . . . . . . . . 3.3.3 Bestimmung des Rangs . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

11 11 12 12 12 13 13 14 14

4 L¨ osungsverfahren f¨ ur nichtlineare Gleichungssysteme 4.1 Kondition . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Fixpunktiteration . . . . . . . . . . . . . . . . . . . . . . 4.3 Konvergenzordnung und Fehlersch¨ atzung . . . . . . . . 4.4 Nullstellenberechnung skalarer Gleichung . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

14 14 14 15 16

1

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

2

CONTENTS

4.5

4.4.1 Bisektion . . . . . . . . . . . . . . . . . . . 4.4.2 Newton-Verfahren . . . . . . . . . . . . . . 4.4.3 Sekanten-Verfahren . . . . . . . . . . . . . . Newton-Verfahren f¨ ur Systeme . . . . . . . . . . . 4.5.1 Vereinfachtes Newton-Verfahren . . . . . . 4.5.2 Numerische Berechnung der Jacobi-Matrix 4.5.3 Ged¨ampftes Newton-Verfahren . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

16 16 16 16 17 17 18

5 Nichtlineare Ausgleichsrechnung 18 5.1 Gauß-Newton-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5.2 Levenberg-Marquardt-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . 19 6 Interpolation 6.1 Lagrange-Polynominterpolation . . . . . . . . . . . . . . . . . . . . 6.1.1 Auswertung an bestimmten Stellen: Neville-Aitken-Schema 6.2 Explizite Darstellung . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Potenzform . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Newton’sche Interpolationsformel . . . . . . . . . . . . . . . 6.2.3 Fehler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Hermite-Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Numerische Differentiation . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

20 20 20 21 21 21 22 22 22

3

1 FEHLERANALYSE

1

Fehleranalyse

1.1

Kondition eines Problems

• unabhaengig vom Algorithmus • gibt an welche Genauigkeit (bei exater Rechnung) mit gestoerten Eingabedaten erwartet werden kann Eingabedaten x ˜eX ⇒ behaftet mit Fehler: ∆x = x˜ − x 1.1.1

→ Problem, Prozess → f

Ausgabedaten y = f (˜ x)eY mit Fehlern: ∆y = f (˜ x) − f (x)

Definition einer Norm

(N1) ||v|| ≥ 0veV und ||v|| = 0 impliziert v = 0 (N2) Fuer alle aeK , veV gilt: ||aV || = |a|||V || (N3) Fuer alle v, w e V gilt die Dreiecksungleichung ||v + w|| ≤ ||v|| + ||w|| Bei einem durch eine Lipschitz-stetige Abbildung beschriebenen Problem laesst sich also die absolute Kondition des Problems gleichmaesig durch die Lipschitz-Konstante beschraenken 1.1.2

Relative & absolute Kondition

absolute Fehler: ||∆x||X , ||∆y||Y ||∆y||Y ||∆x||X ,δy = relative Fehler: δx = ||x||X ||y||Y δy ⇒ relative Kondition: δx ||∆y||Y ⇒ absolute Kondition: ||∆x|X je kleiner die relative Kondition ( mit δx → 0) desto besser ist das Problem konditioniert 1.1.3

Fehler der Ausgabe f (˜ x) − f (x) . Pn = j=1 φj (x) · f (x)

x ˜j −xj xj

∂f (x) xj mit Fehlerverstaerkung: φj (x) = · ∂xj f (x) f (˜ Pn x˜ − x) − f (x) ≤ κrel (x) j=1 jxj j mit κrel (x) = κ∞ ⇒ rel (x) = max |φj (x)| f (x) vereinfachter Fall: f ist nur von einer Variablen abhaengig: f (˜x)−f (x) . f (x) = κrel (x) x˜−x x

1.1.4

x mit κrel (x) := f ′ (x) f (x)

Lineare Abbildungen Satz 2.19

Eine Abbildung L : X → Y ist linear, wenn fuer alle x, z ∈ X und α, β ∈ R gilt: L(αx+βz) = αL(x) + βL(z) ||L(x)||y → ||L||x→y := sup||x||X=1 ||L(x)||y ≡ ||L||x→y := ||x||X ⇒ ||L||y ≤ ||L||X→Y ||x||X , ∀x ∈ X

4

1 FEHLERANALYSE

Identitaet I:

X → X ⇒ ||I||X := ||I||X→X = sup||x||X=1 ||Ix||X = 1

relative Kondition einer linearen Abbildung: ||˜ x − x||X ||L(˜ x) − L(x)||Y ≤ κ(L) ||L(x)||Y ||x||X sup||x||X=1 ||L(x)||Y ||L(x)||X→Y = mit κ(L) ist das Verhaeltnis von inf ||x||X=1 ||L(x)||Y inf ||x||X=1 ||L(x)||Y maximaler Dehnung zur staekstmoeglichen Stauchung der Einheitskugel K||·||X (0, 1) unter der Abbildung L : x → L(x) (jeweils gemessen in der Bildnorm || · ||Y mit κ(L) =

1.2 1.2.1

Rundungs- Maschinenegenauigkeit relative Rundungsfehler b−m e b f l(x) − x b1−m ≤ 2 = =: eps x b−1 be 2

b: Basis des Zahlensystems e: Exponent f: Mantisse m: Mantissenlaenge fl(x): gerundete Zahl eps: relative Maschinengenauigkeit

1.3

Stabilitaet eines Algorithmus

Stabilitaet: Ein Algorithmus heisst gutartig oder stabil, wenn die durch ihn im Laufe der Rechnung erzeugten Fehler in der Groessenordnung des durch die Kondition des Problems bedingten unvermeidbaren Fehlers bleiben. exakte Daten x

Rueckwaertsfehler

f (exakt)

f˜(n u

f (x)

me ris ch

Fehler im Resultat )

f˜(x) = f (˜ x)

exakte Daten x f (exakt)

1.3.1

Leitlinien

• Kenntisse ueber die Kondition eines Problems sind oft fuer die Interpretation oder Bewertung der Ergebnisse von entscheidender Bedeutung • Multiplikation/Division sind fuer alle Eingangsdaten gut konditioniert Subtraktion zweier annaehernd gleicher Zahlen sind schlecht konditioniert (Ausloeschung). Rundungsfehler koennen enorm verstaerkt werden. • In einem Loesungsverfahren sollen (wegen der Stabilitaet) Ausloeschungsefekte vermieden werden • Bei einem stabilen Loesungsverfahren bleiben die im Laufe der Rechung erzeugten Rundungsfehler in der Groessenordnung der durch die Kondition des Problems bedingten vermeidbaren Fehler.

¨ ¨ LINEARE GLEICHUNGSSYSTEME 2 DIREKTE LOSUNGSVERFAHREN FUR

5

• Im Rahmen der Gleitpunktarithmetik sollen Abfragen, ob eine Groesse gleich Null ist oder ob zwei Zahlen gleich sind, vermieden werden. Auf Grund der Aufloesbarkeit bis auf Maschinengenauigkeit ist diese Frage nicht entscheidbar im Sinne einer essentiellen Voraussetzung fuer weitere Entscheidungen und Schritte • Vorteilhafte Reihenfolgen bei der Summenbildung sollen beruecksichtigt werden

2

Direkte L¨ osungsverfahren f¨ ur lineare Gleichungssysteme

H¨ aufig gilt es Probleme der Form Ax = b mit A ∈ Rn×n , x, b ∈ Rn zu l¨osen. Die Verfahren dazu werden hier zusammengefasst. Verfahren zu einer nichtquadratischen Matrix A werden in der linearen Ausgleichsrechnung behandelt.

2.1

Zeilenskalierung

Die Kondition eines linearen Gleichungssystems wird durch die Konditionszahl κ(A) = kAkkA−1 k beschrieben. D.h. beim L¨ osen mit einer Maschinengenauigkeit eps tritt ein unvermeidbarer Fehler in der Gr¨oßenordnung κ(A) · eps auf. Es ist also von Vorteil, eine m¨oglichst geringe Konditionszahl zu erreichen. Dazu verwendet man eine Diagonalmatrix D, die jede Zeile der Matrix A skaliert:   d11 0 . . . 0   a11 . . . a1n  ..   0 d22 . . .  .   . ..  .. DA =  . .   .  ·  .. . . .. ..  .. 0  an1 . . . ann 0 ... 0 dnn

Dabei ist die Dz diejenige Diagonalmatrix, f¨ ur die gilt:

κ∞ (Dz A) ≤ κ∞ (DA) Damit ist die Kondition des Problems optimal und man nennt die Matrix Dz A (zeilenweise) ¨ aquilibriert.

2.2

Dreiecksmatrizen

Besonders einfach lassen sich Gleichungssysteme der Form Ax = b l¨osen, wenn die Matrix A eine obere- oder untere Dreiecksmatrix ist:

¨ ¨ LINEARE GLEICHUNGSSYSTEME 2 DIREKTE LOSUNGSVERFAHREN FUR

obere Dreiecksmatrix

untere Dreiecksmatrix



  R=  

r11 0 .. .

 0 l11   l21 L=  .  ..

ln1

6

 r1n ..  r22 .   ..  .. .. . . .  ... 0 rnn  0 ... 0 ..  .. . l22 .    .. . 0  . . . ln,n−1 lnn r21

...

Diagonalmatrix obere- und untere Dreiecksmatrix Stehen auf der Diagonalen der Dreiecksmatrix nur Einsen, spricht man von einer normierten Dreiecksmatrix.

2.2.1

Vorw¨ arts- und R¨ uckw¨ artseinsetzen

Mit Dreiecksmatrizen lassen sich die Eintr¨age xj des L¨ osungsvektors x durch Vorw¨arts- bzw. R¨ uckw¨artseinsetzen l¨osen.

R¨ uckw¨ artseinsetzen bei oberen Dreiecksmatrizen xj = (bj −

n X

k=j+1

rj,k · xk )/rj,j

j = n, . . . , 1

Vorw¨ artseinsetzen bei unteren Dreicksmatrizen xj = (bj − Aufwand

2.3

1 2 2n

j−1 X

k=1

lj,k · xk )/lj,j

j = 1, . . . , n

Operationen

Gauß-Elimination und LR-Zerlegung

Die L¨ osung von Ax = b ¨ andert sich nicht, wenn man Vielfache einer Gleichung (Zeile) von einer anderen Gleichung (Zeile) abzieht. F¨ ur k = 1, . . . , n−1 zieht man ein Vielfaches der k-ten Zeile der Matrix A von allen folgenden Zeilen ab, so dass in der k-ten Spalte alle Eintr¨age unterhalb der k-ten Zeile Null werden. Auf diese Weise u ¨ berf¨ uhrt man die Matrix A in eine obere Dreiecksmatrix R. Die rechte Seite b muss dabei ebenfalls entsprechend angepasst werden! Wird ein Eintrag aj j (so genanntes Pivotelement ) auf der Hauptdiagonalen w¨ahrend der Umformungen Null, so versagt der Algorithmus 2.3.1

Gauß-Elimination ohne Pivotisierung

¨ 1. Uberf¨ uhre Ax = b in Rx = c 2. L¨ ose Rx = c

¨ ¨ LINEARE GLEICHUNGSSYSTEME 2 DIREKTE LOSUNGSVERFAHREN FUR

2.3.2

7

LR-Zerlegung

¨ Notiert man sich w¨ ahrend der Uberf¨ uhrung die Koeffizienten jeder Zeilenumformung in einer Matrix, ergibt sich eine untere Dreiecksmatrix   1 0 ... 0  ..  ..  l21 . 1 .   L= .  . . .. ..  .. 0 ln1 . . . ln,n−1 1

wobei der Eintrag lkj den Koeffizient zur Umformung der k-ten Zeile darstellt, so dass die j-te Spalte Null wird. Auf diese Weise erh¨ alt man eine Faktorisierung der Matrix A: A = LR 2.3.3

Gauß-Elimination mit Spaltenpivotisierung

Verschwindet w¨ ahrend des Eliminationsverfahrens ein Pivotelement, so ist eine Zeilenvertauschung notwendig. Dazu wird einfach die j-te Zeile durch eine k-te Zeile k > j ersetzt. Ebenfalls ist es auf diese Weise m¨oglich Rundungsfehler zu reduzieren, wenn immer das Spaltenmaximum als Pivotelement verwendet wird. Die Spaltenvertauschung l¨asst sich durch eine Pivotmatrix P ausdr¨ ucken: P A = LR Beispiel:



1 0 P A = 0 0 0 1

  0 1 1 1 · 2 2 0 3 3

  1 1 2 = 3 3 2

 1 1 3 3 2 2

Vorgehen: 1. Bestimme Diagonalmatrix D, so dass DA zeilenweise ¨aquilibriert ist. 2. Wende Gauß-Elimination mit Spaltenpivotisierung an, d.h. suche f¨ ur jeden Schritt das betragsm¨aßig gr¨ oßte Element in der ersten Spalte der Restmatrix und tausche gegebenenfalls. Faktorisierung: P DA = LR Aufwand:

1 3 3n

Totalpivotisierung Um noch g¨ unstigere Pivotelemente zu bekommen, l¨asst sich auch noch die Reihenfolge der Unbekannten (d.h. Spalten) ¨andern. Allerdings ist der Aufwand daf¨ ur sehr groß. 2.3.4

Anwendungen der LR-Zerlegung

L¨ osen eines Gleichungssystems Die LR-Zerlegung erlaubt, ein lineares Gleichungssystem Ax = b ⇔ P Ax = LRx = P b in zwei Schritten zu l¨osen: 1. Ly = P b durch Vorw¨ artseinsetzen 2. Rx = y durch R¨ uckw¨artseinsetzen

¨ ¨ LINEARE GLEICHUNGSSYSTEME 2 DIREKTE LOSUNGSVERFAHREN FUR

Aufwand

1 3 n + 3 |{z}

Zerlegung

8

2n2 |{z}

Vorw¨ arts- und R¨ uckw¨ artseinsetzen

Berechnung der Inversen Sei A−1 = (a1 a2 . . . an ), d.h. ai die i-te Spalte der Inversen. Mit AA−1 = I gilt Aai = ei . 1. Bestimme LR-Zerlegung P A = LR 2. L¨ ose die Gleichungssysteme LRai = P ei Aufwand

4 3 3n

Berechnung von Determinanten Mit P A = LR gilt: n Y det L det R #Zeilenvertauschungen det A = R = (−1) · det P j=1

Nachiteration ˜ R ˜ Da die begrenzte Genauigkeit der Rechnerarithmetik nur fehlerbehaftete Ergebnisse x ˜, L, liefert, gilt f¨ ur das Residuum r = b − A˜ x 6= 0

Betrachtet man x0 = x ˜ und r0 = b − Ax0 erh¨ alt man den Fehler δ 0 = x − x0 durch L¨ osung des Defektsystems Aδ 0 = Ax − Ax0 = b − Ax0 = r0 Dies ist z.B. durch eine bereits bekannte Faktorisierung A = LR einfach m¨oglich. Dadurch erh¨ alt man ein genaueres x1 = x0 + δ 0 . Wiederholt durchgef¨ uhrt kann man so die Genauigkeit der L¨ osung erh¨ ohen. 2.3.5

Cholesky-Zerlegung

Ist die Matrix A symmetrisch positiv definit, so existiert die Choleskyzerlegung A = LDLT mit einer normierten unteren Dreiecksmatrix L und einer Diagonalmatrix D mit di,i > 0, i = 1, . . . , n. Diese l¨asst sich einfach bestimmen durch dk,k = ak,k −

li,k = (ai,k − Aufwand

1 3 6n

k−1 X j=1

k−1 X

2 lk,j dj,j

j=1

li,j dj,j lk,j )/dk,k

¨ ¨ LINEARE GLEICHUNGSSYSTEME 2 DIREKTE LOSUNGSVERFAHREN FUR

9

Symmetrisch positiv definit Eine Matrix A ∈ Rn×n heißt symmetrisch positiv definit (s.p.d.), falls AT = A

(symmetrisch)

und xT Ax > 0

(positiv definit)

f¨ ur alle x ∈ Rn , x 6= 0 gilt. F¨ ur A gilt: • A ist invertierbar und A−1 ist s.p.d. • Alle Eigenwerte λi > 0 • Jede Hauptuntermatrix von A ist s.p.d. • det(A) > 0 • Der betragsgr¨oßte Eintrag der Matrix befindet sich auf der Hauptdiagonalen

2.4

QR-Zerlegung

Orthogonalit¨ at Eine Matrix Q ∈ Rn×n ist orthogonal, falls QT Q = I d.h. die Spalten von Q bilden eine Orthonormalbasis. F¨ ur Q gilt: 1. Q−1 = QT ist ebenfalls orthogonal 2. kQxk2 = kxk2 3. kAk2 = kQAk2 = kAQk2 4. κ2 (Q) = 1 5. κ2 (A) = κ2 (QA) = κ2 (AQ) ˜ sei orthogonal, dann ist QQ ˜ orthogonal 6. Q Idee 1. Zerlege A ∈ Rn×n in A = QR Q: orthogonal, R: obere Dreiecksmatrix 2. L¨ ose Ax = b ⇔ QRx = b ⇔ Rx = QT b Daf¨ ur ist lediglich R¨ uckw¨artseinsetzen n¨ otig! Die QR-Zerlegung kann dabei f¨ ur beliebige rechteckige Matrizen durchgef¨ uhrt werden! Zur Durchf¨ uhrung stehen zwei Methoden zur Verf¨ ugung: Givens-Rotationen und HouseholderTransformationen. 2.4.1

Givens-Rotationen

Givens-Rotationen nehmen eine Drehung der Spalten von A vor. Mit einer Folge von Rotationen kann man eine Matrix auf obere Dreiecksgestalt transformieren: GTiN ,kN · . . . GTi1 ,k1 · A = R

¨ ¨ LINEARE GLEICHUNGSSYSTEME 2 DIREKTE LOSUNGSVERFAHREN FUR

10

Vorgehen Nullsetzen des Eintrags (i, j) der Matrix A q 1. Berechne rij = sign(aij ) a2jj + a2ij wobei sign(0) = 1 2. Berechne cij =

ajj rij

und sij =

aij rij

3. Stelle die Givens-Matrix auf:   1    0      j→    Gij =  ..  .     i→      

↓j

0 .. . .. .

..

...

↓i



.

1 0

0 c 0 .. .

.. .

0 −s

0 1 .. .

... ... 0 .. .

0 ... .. .

s 0 .. .

.. .

... 0

0 ...

1 0

0 c

0

0

0

...

0

0

...

1 .. .

..

.

..

.

0

0             ..  .           0  1

4. Wiederhole f¨ ur n¨ achsten Eintrag bis A obere Dreiecksgestalt besitzt. Vorsicht: F¨ ur die Wiederholung verwende die neu erzeugte Matrix A′ = Gij · A! Damit erh¨ alt man eine Faktorisierung A = G2,1 · · · · · Gn,n−1 · R = QR Zusammenfassung • sehr stabil • Ber¨ ucksichtigung vorhandener Nulleintr¨age ist m¨oglich • Aufwand 43 n3 f¨ ur m ≈ n, sonst 2mn2 2.4.2

Householder-Transformationen

Hosueholder-Transformationen nehmen eine Spiegelung der Matrix A an (Hyper-)Ebenen durch den Ursprung vor, um diese auf obere Dreicksform zu bringen. Dies erfolgt durch orthogonale Householder-Matrizen Qv mit v als Normalenvektor der (Hyper-)Ebene: Qv = I −

2 vv T vT v

Den zugeh¨origen Normalenvektor zur i-ten Spalte erh¨ alt man mit v i = ai + αi · ei wobei αi = sign(aii )kai k2 mit sign(0) = 1.

11

3 LINEARE AUSGLEICHSRECHNUNG

Vorgehen Nullsetzen der Eintr¨age unterhalb der Hauptdiagonalen in der i-ten Spalte der Matrix A 1. Bestimme αi und v i 2. Berechne die erste Spalte der neuen Matrix A(i+1) mit Qi ai = −αe1 3. Berechne die u ¨ brigen Spalten A′ : A = ai Qi A′ = A′ − 4. Fahre mit der Untermatrix A˜(i+1) :  ∗ 0  Qi A(i) =  .  ..

∗


A′ mit

2

T vi vi

... A˜(i+1)

0

fort, bis A obere Dreicksgestalt besitzt.

T

· v i · v i A′

∗



   = A(i+1) 

Damit erh¨ alt man eine Faktorisierung A = QT1 · . . . · QTn−1 · R = QR Zusammenfassung • sehr stabil • Aufwand 23 n3 f¨ ur m ≈ n, sonst mn2

3

Lineare Ausgleichsrechnung

ur das gilt: Problemstellung Bestimme ein x∗ ∈ Rn zu A ∈ Rm×n und b ∈ Rm f¨ kAx∗ − bk2 = minx∈Rn kAx − bk2 mit m: Anzahl der Messungen und n: Anzahl der unbekannten, wobei i. A. m ≫ n.

3.1

Kondition des Problems

Bez¨ uglich St¨ orungen in b ergibt sich: κ2 (A) k˜b − bk2 k˜ x − x∗ k2 ≤ · kx∗ k2 cos θ kbk2 ∗

∗

k2 k2 kAxk2 mit cos θ = kAx und κ2 (A) = maxx6=0 kAx kbk2 kxk2 /minx6=0 kxk2 . Bez¨ uglich St¨ orungen in A ergibt sich:

k˜ x − x∗ k2 kA˜ − Ak2 ≤ (κ2 (A) + κ2 (A)2 tan θ) ∗ kx k2 kAk2 ∗

k2 mit tan θ = kb−Ax kAx∗ k2 . Falls die Norm des Residuums klein ist gegen¨ uber der Norm der Eingabe, ist das Problem ¨ahnlich einem linearen Gleichungssystem konditioniert.

12

3 LINEARE AUSGLEICHSRECHNUNG

3.2 3.2.1

Numerische L¨ osung L¨ osung u ¨ ber Normalengleichungen

Das Problem der linearen Ausgleichsrechnung l¨asst sich umformen zu AT Ax∗ = AT b AT A ist s.p.d., d.h. u ¨ ber Cholesky-Zerlegung l¨osbar: 1. Berechne AT A, AT b 2. F¨ uhre Cholesky-Zerlegung durch: LDLT = AT A 3. L¨ ose Ly = AT b LT x∗ = D−1 y Aufwand:

1 1 n2 mn2 + n3 + |{z} 2 6 | {z } |{z} 3. 1.

Nachteile

2.

• F¨ ur große m aufw¨andig • Gefahr von Genauigkeitsverlusten durch Ausl¨oschung bei Skalarprodukten • Fehlerverst¨ arkung durch κ2 (AT A) = (κ2 (A))2 , d.h. m¨oglicherweise nicht stabil 3.2.2

L¨ osung u ¨ ber QR-Zerlegung

Suche die Orthogonalmatrix Q ∈ Rm×m so dass gilt:     ˜ b R und Qb = 1 QA = R = b2 0 ˜ ∈ Rn×n . mit oberer Dreicksmatrix R Vorgehen 1. Bestimme QR-Zerlegung von A:   ˜ R QA = 0 2. Berechne Qb =

˜ ∈ Rn×n ,R

  b1 b2

˜ = b1 3. L¨ ose Rx Aufwand

1 2 2mn + n2 mn |{z} + |{z} 2 |{z} 1. 2. 3.

Zusammenfassung

• Aufwand etwa um Faktor 2 h¨ oher als bei L¨ osung u ¨ ber Normalengleichungen ˜ • Kondition und damit Stabilit¨ at wesentlich besser, da κ2 (A) = κ2 (R)

13

3 LINEARE AUSGLEICHSRECHNUNG

3.3

Singul¨ arwertzerlegung

Zu jeder Matrix A ∈ Rm×n existieren orthogonale Matrizen u ∈ Rm×m , V ∈ Rn×n und eine Diagonalmatrix Σ = diag(σ1 , . . . , σp ) ∈ Rm×n , p = min{m, n} mit den Singul¨ arwerten σ1 ≥ σ2 ≥ . . . ≥ σp ≥ 0, so dass U T AV = Σ Damit ergeben sich folgende Eigenschaften ui , vi seien die i-ten Spalten der Matrizen U, V • Avi = σi ui

und AT ui = σi vi

mit i = 1, . . . , p

• Rang(A) = r • Bild(A) = span{u1 , . . . , ur }

und Kern(A) = span{vr+1 , . . . , vn }

• kAk2 = σ1 • Konditionszahl, falls Rang(A) = n ≤ m: κ∗2 (A) = κ2 (A) =

maxkxk2 =1 kAxk2 minkxk2 =1 kAxk2

p • Singul¨ arwerte: {σi |i = 1, . . . , r} = { λi (AT A)|i = 1, . . . , r} wobei λi die Eigenwerte von AT A darstellen 3.3.1

Berechnung der Singul¨ arwerte

Die Multiplikation mit Orthogonalmatrizen ¨andert die Singul¨ arwerte einer Matrix nicht. Vorgehen 1. Umformung von A mit Householdertransformationen auf Bidiagonalmatrix: 

ˆ1 · . . . · Q ˆ n−2 Qm−1 · Qm−2 · . . . · Q1 · A · Q

∗

 0       =B= .  ..     

0

∗

0

∗ .. .

∗ .. .

...

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

 0 ..  .   0   ∗   ∗  0  ..  . 0

ˆi Qi setzt dabei die Eintr¨ age unterhalb der Bidiagonalen in der i-ten Spalte auf 0, Q setzt die Eintr¨ age rechts der Bidiagonalen in der i-ten Zeile auf 0. 2. Berechne die Eigenwerte der Tridiagonalmatrix B T B p 3. Bestimme die Singul¨ arwerte σi = λi (B T B)

Die Berechnung der Eigenwerte einer Tridiagonalmatrix f¨allt wesentlich leichter als die direkte Berechnung der Eigenwerte von AT A, deshalb der Umweg u ¨ ber die Transformationen.

¨ ¨ NICHTLINEARE GLEICHUNGSSYSTEME 4 LOSUNGSVERFAHREN FUR

3.3.2

14

Pseudoinverse

Die Gleichung AT Ax∗ = AT b l¨asst sich umschreiben zu x∗ = (AT A)−1 AT b. Die Matrix A+ = (AT A)−1 AT heißt Pseudoinverse von A. F¨ ur m = n gilt: A+ = A−1 . Die Pseudoinverse ist explizit bestimmbar u ¨ ber A+ = V Σ+ U T wobei Σ+ = diag(σ1−1 , . . . , σr−1 , 0, . . . , 0) ∈ Rn×m und σ1 ≥ . . . σr > σr+1 = . . . = σp = 0, p = min{m, n}. 3.3.3

Bestimmung des Rangs

Da die Anzahl der Singul¨ arwerte ungleich Null dem Rang der Matrix entspricht, ist prinzipiell die Bestimmung des Rangs u ¨ ber die Singul¨ arwerte m¨oglich. Jedoch ist im Allgemeinen aufgrund von Rundungsfehlern der Rang so nicht bestimmbar, da z.B. σk = 0 nicht entscheidbar ist. Jedoch ist es m¨oglich den numerischen Rang zu bestimmen: √ ˜ = min{1 ≤ k ≤ p|˜ Rangnum (A) σk+1 ≤ σ ˜1 mn · eps} wobei A˜ ∈ Rm×n , σ ˜p+1 = 0 und eps die Maschinengenauigkeit darstellt

4

L¨ osungsverfahren f¨ ur nichtlineare Gleichungssysteme

Besteht kein linearer Zusammenhang f¨ ur die Eingangsvariablen, handelt es sich um folgende Problemstellung Suche x∗ = (x∗1 , . . . , x∗n )T ∈ Rn zu gegebenem f = (f1 . . . fn )T : Rn → Rn mit f1 (x∗1 , . . . , x∗n ) = 0 .. .. .. . . .

fn (x∗1 , . . . , x∗n ) = 0 oder kurz f (x∗ ) = 0

4.1

Kondition

Brauchen wir das?

4.2

Fixpunktiteration f (x∗ ) = 0

l¨asst sich umschreiben zu x∗ = x∗ − Mx∗ f (x∗ )

alt man das Fixmit Mx∗ ∈ Rn×n als von x∗ abh¨angige, invertierbare Matrix. Daraus erh¨ punktproblem x∗ = Φ(x∗ ) mit Φ(x) = x − Mx f (x)

1. W¨ ahle Startwert x0 (in Umgebung von x∗ ) 2. Bilde xk+1 = Φ(xk ) Die Fixpunktiteration beschreibt der

k = 0, 1, 2, . . .

¨ ¨ NICHTLINEARE GLEICHUNGSSYSTEME 4 LOSUNGSVERFAHREN FUR

15

Banach’sche Fixpunktsatz X sei ein linearer normierter Raum mit der Norm k · k und Φ eine Selbstabbildung auf E ⊆ X: Φ:E→E sowie eine Kontraktion auf E: kΦ(x) − Φ(y)k ≤ Lkx − yk

∀x, y ∈ E mit L < 1

Dann gilt: 1. Es existiert genau ein Fixpunkt x∗ von Φ auf E 2. F¨ ur beliebiges x0 ∈ E konvergiert xk+1 = Φ(xk ),

k = 0, 1, 2, . . .

gegen x∗ 3. A-priori-Fehlerabsch¨ atzung: kxk − x∗ k ≤

Lk kx1 − x0 k 1−L

4. A-posteriori-Fehlerabsch¨ atzung (genauer als a-priori-Absch¨ atzung): kxk − x∗ k ≤

L kxk − xk−1 k 1−L

F¨ ur die notwendige Anzahl k an Iterationen, um eine Genauigkeit ε zu erreichen, gilt: Lk kx1 − x0 k ≤ ε 1−L k≥

4.3

ε(1−L) ) log( kx 1 −x0 k

log L

Konvergenzordnung und Fehlersch¨ atzung

¨ Ziel ist es, eine Ann¨aherung an x∗ m¨oglichst schnell zu erreichen. Uber die Geschwindigkeit der Ann¨aherung einer Folge {xk }k∈N gibt die Konvergenzordnung p Auskunft: kxk+1 − x∗ k ≤ ckxk − x∗ kp f¨ ur alle k ≥ k0 ∈ N gilt, wobei 0 < c < 1 falls p = 1. Abh¨angig von der Konvergenzordnung l¨asst sich eine vereinfachte a-posteriori-Fehlersch¨ atzung vornehmen: skalare Folgen Ak (xk − xk−1 ) p = 1: x∗ − xk ≈ 1−A k xk −xk−1 wobei Ak = xk−1 −xk−2 etwa konstant sein sollte.

p > 1: x∗ − xk ≈ xk+1 − xk Vektorfolgen p = 1: keine allgemeine Formel bekannt p > 1: kxk − x∗ k ≈ kxk+1 − xk k

¨ ¨ NICHTLINEARE GLEICHUNGSSYSTEME 4 LOSUNGSVERFAHREN FUR

4.4

16

Nullstellenberechnung skalarer Gleichung

Ein h¨ aufig zu l¨osendes Problem, ist die Nullstellenfindung einer stetigen Funktion f , also eine Problemstellung der Form f (x∗ ) = 0 4.4.1

Bisektion

• Zwei Werte ak < bk seien bekannt, f¨ ur die f (ak ) · f (bk ) < 0 gilt, d.h. eine Nullstelle findet sich in [ak , bk ]. • Berechne

1 (ak + bk ) 2 und bestimme so die neuen Ausgangswerte ak+1 , bk+1 mit [ak+1 , bk+1 ] = [ak , xk ] oder [ak+1 , bk+1 ] = [xk , bk ], so dass gilt: f (ak+1 ) · f (bk+1 ) ≤ 0

Konvergenzordnung: 4.4.2

xk =

p=1

Newton-Verfahren

Zur Nullstellensuche kann man eine Fixpunktiteration xk+1 = Φ(xk ) aufstellen: xk+1 = xk −

f (xk ) , f ′ (xk )

k = 0, 1, 2, . . .

Geometrisch stellt das Newton-Verfahren die Konstruktion von Tangenten im Punkt xk dar, deren Nullstelle xk+1 den Ausgangspunkt f¨ ur den n¨ achsten Iterationsschritt darstellt. Konvergenzordnung: p = 2, jedoch nur in einer Umgebung von x∗ , d.h. ein geeigneter Startwert muss vorliegen. 4.4.3

Sekanten-Verfahren

¨ Ahnlich dem Newton-Verfahren n¨ ahert das Sekanten-Verfahren die Nullstelle einer Funktion u ¨ ber die Nullstellen von Geraden an. Jedoch wird dort jeweils eine Sekante durch die Punkte xk , xk−1 gelegt und deren Nullstelle xk+1 gew¨ahlt: xk+1 = xk − f (xk )

xk − xk−1 f (xk ) − f (xk−1 )

Konvergenzordnung: p ≈ 1.6, jedoch entf¨ allt die h¨ aufig teure Berechnung der Ableitung f ′ und das Verfahren ist somit meist sehr effizient. Nachteilig ist die Notwendigkeit von zwei Startwerten.

4.5

Newton-Verfahren f¨ ur Systeme

Das Newton-Verfahren l¨asst sich f¨ ur mehrdimensionale Probleme verallgemeinern auf die Form xk+1 = xk − [f ′ (xk )]−1 f (xk )

¨ ¨ NICHTLINEARE GLEICHUNGSSYSTEME 4 LOSUNGSVERFAHREN FUR

17

Vorgehen Zu einem Startwert x0 f¨ ur k = 0, 1, 2, . . . 1. Berechne f (xk ), f ′ (xk ) 2. L¨ ose das Lineare Gleichungssystem f¨ ur sk : f ′ (xk )sk = −f (xk ) 3. Setze xk+1 = xk + sk wobei

 ∂f

 f ′ (x) =  

1 (x) ∂x1

.. .

∂fn (x) ∂x1

... .. . ...

∂f1 (x) ∂xn

.. .

∂fn (x) ∂xn

   

Jacobi-Matrix

Konvergenz Sei Ω ∈ Rn , f : Ω → Rn , f ′ (x) eine invertierbare Jacobi-Matrix und k[f ′ (x)]−1 k ≤ β f¨ ur alle x ∈ Ω sowie f ′ Lipschitz-stetig auf Ω: kf ′ (x) − f ′ (y)k ≤ γkx − yk,

x, y ∈ Ω

Existiert eine L¨ osung x∗ von f (x) = 0 und ein Startwert x0 ∈ Kω (x∗ ) = {x ∈ Rn |kx∗ − xk < ω} mit 2 ω≤ βγ Dann konvergiert {xk }∞ k=0 quadratisch: kxk+1 − x∗ k ≤ 4.5.1

βγ k kx − x∗ k2 , 2

k = 0, 1, 2, . . .

Vereinfachtes Newton-Verfahren

Die Berechnung der Jacobi-Matrix ist sehr aufwendig, deshalb erlaubt das vereinfachte Newton-Verfahren, diese nur im ersten Schritt zu berechnen und auch in den folgenden zu verwenden: f ′ (xk ) = f ′ (x0 ) Nachteil: keine quadratische Konvergenz mehr, deshalb wird in der Praxis h¨ aufig die Jacobi-Matrix alle 3 bis 5 Schritte neu berechnet. 4.5.2

Numerische Berechnung der Jacobi-Matrix

Bei der Berechnung der Jacobi-Matrix kann man die partiellen Ableitungen ann¨ahern u ¨ ber einen Differenzenquotienten: fi (x + h · ej ) − fi (x) ∂fi (x) ≈ ∂xj h ej : j-ter Basisvektor h: muss abh¨angig von der Aufgabe gew¨ahlt werden; ein kleinerer Wert verbessert die Ann¨aherung, provoziert jedoch auch Ausl¨oschung.

5 NICHTLINEARE AUSGLEICHSRECHNUNG

4.5.3

18

Ged¨ ampftes Newton-Verfahren

Mit der Schrittweite sk =

−f (xk ) f ′ (xk )

kann es zugunsten einer bessern Ann¨aherung sinnvoll sein, nur einen Teil des Schrittes zu machen, also xk+1 = xk + λ · sk

5

Nichtlineare Ausgleichsrechnung

¨ Ahnlich wie bei der Linearen Ausgleichsrechnung, kann es auch bei nichtlinearen Problemen n¨ otig sein, zu einer Reihe von Messwerten eine L¨ osung anzun¨ ahern. Problemstellung Bestimme x∗ ∈ Rn so, dass kF (x∗ )k2 = minx∈Rn kF (x)k2 mit f : Rn → Rm , oder ¨ aquivalent

Fi (x) = y(ti , x) − bi , i = 1, . . . , m Φ(x∗ ) = minx∈Rn Φ(x)

wobei Φ : Rn → R, Φ(x) = 12 kF (x)k22 = 12 F (x)T F (x) Dabei gilt: ∇Φ(x) = F ′ (x)T F (x)

19

5 NICHTLINEARE AUSGLEICHSRECHNUNG

Φ′′ (x) = F ′ (x)T F ′ (x) +

m X

Fi (x)Fi′′ (x)

i=1

5.1

Gauß-Newton-Verfahren

Ann¨ahern der L¨ osung des nichtlinearen Problems u ¨ ber eine Reihe geeigneter linearer Probleme: W¨ahle einen Startwert x0 . F¨ ur k = 0, 1, 2, . . . 1. Berechne F (xk ), F ′ (xk ) 2. L¨ ose f¨ ur sk das lineare Ausgleichsproblem kF ′ (xk )sk + F (xk )k2 = mins∈Rn kF ′ (xk )s + F (xk )k2 3. Setze xk+1 = xk + sk Als Abbruchkriterium formuliert man h¨ aufig kF ′ (xk )T F (xk )k2 ≤ ε da die Ableitung von Φ in diesem Punkt (fast) Null wird. Die Gauß-Newton-Methode ist eine Fixpunktiteration. Eigenschaften • Wenn Gauß-Newton-Methode konvergiert, tut sie das i. A. nicht schneller als linear • Kritische Punkte x∗ , die ein Maximum oder Sattelpunkt darstellen (d.h. deren Ableitung ebenfalls Null ist), sind abstoßend f¨ ur das Newton-Verfahren, k¨onnen also nicht f¨alschlich gefunden werden • lokale Konvergenz, wenn ̺(k)kF (x∗ )k2 < 1 (̺(k): Spektralradius ̺(k) = max|λk |)

5.2

Levenberg-Marquardt-Verfahren

Eine dem ged¨ ampften Gauß-Verfahren ¨ahnliche Methode ist das Levenberg-MarquardtVerfahren. Problemstellung Finde sk ∈ Rn , so dass kF ′ (xk )sk + F (xk )k22 + µ2 ksk k22 = min wobei µ > 0 oder ¨ aquivalent

 ′ k   

F (x ) k F (xk )

= min s +

µI 0 2

Somit besitzt das Problem immer vollen Rang und eine eindeutige L¨ osung. Die Korrektur in jedem Schritt betr¨ agt somit ksk k2 ≤

kF (xk )k2 µ

kann also durch µ ged¨ ampft werden. Auch das Levenberg-Marquardt-Verfahren ist eine Fixpunktiteration.

20

6 INTERPOLATION

Vorgehen W¨ ahle Startwert x0 und einen Anfangswert f¨ ur µ F¨ ur k = 0, 1, 2, . . . 1. Berechne F (xk ), F ′ (xk ) 2. L¨ ose das lineare Ausgleichsproblem

 ′ k   

F (x ) k F (xk )

= min s +

µI 0 2

3. Teste, ob die Korrektur xk akzeptabel ist. Wenn nicht: passe µ an und wiederhole Schritt 2. 4. Setze xk+1 = xk + sk

6

Interpolation

Problemstellung Die Berechnung von beliebigen Funktionswerten einer Funktion, die nur an diskreten Stellen bekannt ist, durch das Finden einer Funktion, die an den bekannten Stellen mit den gegebenen Werten u ¨ bereinstimmt.

6.1

Lagrange-Polynominterpolation

Aufgabe

Finde zu Daten f (x0 ), f (x1 ), . . . , f (xn ) ein Polynom Pn ∈ Πn mit Pn (xj ) = f (xj ),

j = 0, 1, . . . , n

Das Lagrange-Interpolationspolynom ist stets eindeutig l¨osbar zu gegebenen, beliebigen Daten f (x0 ), f (x1 ), . . . , f (xn ) als Pn (x) =

n X

f (xj )ljn (x)

j=0

wobei ljn (x) =

n Y x − xk xj − xk

k=0 k6=j

die so genannten Lagrange-Fundamentalpolynome sind. 6.1.1

Auswertung an bestimmten Stellen: Neville-Aitken-Schema

H¨ aufig sind nur wenige interpolierte Werte n¨ otig, d.h. keine explizite Bestimmung des Polynoms. Diese zu bestimmen erlaubt das Neville-Aitken-Schema:

mit Pi,k = Pi,k−1 + Aufwand

n2

Pi,0 f (x0 ) f (x1 ) f (x2 ) .. .

Pi,1

Pi,2

x0 x1 x2 .. .

P1,1 P2,1

P2,2

xn

f (xn )

Pn,1

x−xi xi −xi−k (Pi,k−1

...

Pi,n

..

− Pi−1,k−1 )

Pn,2

. ...

Pn,n

21

6 INTERPOLATION

6.2

Explizite Darstellung

Will man das Polynom explizit bestimmen, gibt es zwei M¨ oglichkeiten: die simple Potenzform oder die Newtonform. 6.2.1

Potenzform

Die Darstellung in der klassischen Potenzform P (f |x0 , . . . , xn )(x) = a0 + a1 x + a2 x2 + . . . + an xn ist einfach u ¨ ber ein lineares Gleichungssystem zu bestimmen:       f (x0 ) 1 x0 x20 . . . xn0 a0 1 x1 x21 . . . xn1     f (x1 )  a1       .. ..  .. .. ..  ·  . . . =   .  .  . . . 2 n a n f (xn ) 1 xn x ... x n

n

Dieses ist mit den bekannten Verfahren l¨osbar. Allerdings ist die Kondition des Systems schlecht und somit Rundungsfehler sehr groß. F¨ ur numerische Berechnungen ist diese Darstellung nicht geeignet. Liegt dennoch ein Polynom in dieser Form vor, so l¨asst es sich mit dem Horner-Schema auswerten: • Setze bn = an • Berechne f¨ ur k = n − 1, n − 2, . . . , 0 bk = ak + xbk+1 • Dann ist p(x) = b0 .

Aufwand 6.2.2

Gegen¨ uber naivem Vorgehen spart das Horner-Schema ca. n Operationen.

Newton’sche Interpolationsformel

Eine numerisch stabilere Darstellung ist die Newton-Form: P (f |x0 , . . . , xn ) = [x0 ]f + (x − x0 )[x0 , x1 ]f + (x − x0 )(x − x1 )[x0 , x1 , x2 ]f + . . . + (x − x0 ) · · · (x − xn−1 )[x0 , . . . , xn ]f wobei [x0 , . . . , xi ]f =

f (xi ) − Pi−1 (xi ) [x1 , . . . xi ]f − [x0 , . . . , xi−1 ]f = (x − x0 ) · · · (x − xi−1 ) xi − x0

und [x0 ]f = f (x0 ). Die praktische Berechnung ist mit dem Schema der dividierten Differenzen m¨oglich:

x0

[xi ]f [x0 ]f

x1

[x1 ]f

x2

[xi , xi+1 ]f >

[x0 , x1 ]f

>

[x1 , x2 ]f

>

[x2 ]f

[x3 ]f .. .

[x2 , x3 ]f .. .

[xi , xi+1 , xi+2 , xi+3 ]f

[x0 , x1 , x2 ]f >

> >

x3 .. .

[xi , xi+1 , xi+2 ]f

[x1 , x2 , x3 ]f .. .

[x0 , x1 , x2 , x3 ]f .. .

Die gesuchten Koeffizienten treten am oberen Rand der Tafel auf.

...

22

6 INTERPOLATION

Aufwand

ca.

1 2 2n

Divisionen und n(n + 1) Additionen

Eigenschaften dividierter Differenzen 1. [x0 , . . . , xn ]f ist eine symmetrische Funktion, die Reihenfolge der St¨ utzstellen ist deshalb egal (z.B. [x0 , x1 , x2 ]f = [x1 , x0 , x2 ]f ) 2. F¨ ur Q ∈ Πk−1 gilt: [x0 , . . . , xk ]Q = 0 6.2.3

Fehler

Sei a = min{x0 , . . . , xn }, b = max{x0 , . . . , xn }, x ∈ R und I = [min{a, x}, max{b, x}]. Dann existiert f¨ ur f ∈ C n+1 (I) ein ξ ∈ I so, dass f (x) − P (f |x0 , . . . , xn )(x) = (x − x0 ) · · · (x − xn ) ·

f (n+1) (ξ) (n + 1)!

und maxx∈[a,b] |f (x) − P (f |x0 , . . . , xn )(x)| ≤ maxx∈[a,b] |

n Y

j=0

(x − xj )|maxx∈[a,b]

|f (n+1) (x)| (n + 1)!

Zusammenfassung • Eine h¨ ohere Anzahl an St¨ utzstellen reduziert den Fehler im Intervall, erh¨ oht jedoch die Oszillation und den Fehler am Rand des Intervalls • Eine geringere Schrittweite verbessert die Interpolation • Die Qualit¨ at der Interpolation ist also begrenzt. Eine bessere Ann¨aherung als eine polynomielle Funktion, ist eine Approximation, die st¨ uckweise polynomial ist. Dies erfolgt bei Splinefunktionen.

6.3

Hermite-Interpolation

Die Lagrange-Interpolation ist nur dann wohldefiniert, wenn die St¨ utzstellen paarweise verschieden sind und f stetig ist. Bei mehrfachen St¨ utzstellen l¨asst sich das Interpolationspolynom u ¨ ber die Hermite-Interpolation bestimmen: Zu x0 ≤ x1 ≤ . . . ≤ xn definiere f¨ ur j = 0, 1, . . . , n µj (Pn ) = µj (f ) = f (lj ) (xj ), j = 0, 1, . . . , n wobei lj = max{r|xj = xj−r } Die Bestimmung des Polynoms erfolgt ebenfalls u ¨ ber dividierte Differenzen, wobei ( ([x0 , . . . , xi−1 , xi+1 , . . . , xk ]f − [x0 , . . . , xj−1 , xj+1 , . . . , xk ]f )/(xj − xi ) [x0 , . . . , xk ]f = f (k) (x0 ) k!

falls xi 6= xj

falls x0 = . . . = xk

Fehler F¨ ur den auftretenden Fehler gilt das gleiche wie bei der Lagrange-Interpolation.

6.4

Numerische Differentiation

F¨ ur Funktionen, deren Werte nur an diskreten Stellen bekannt sind, ist eine numerische Ann¨aherung der Ableitungen u ¨ ber Differenzenquotienten m¨oglich: f (x + h) − f (x) h 1 1 f (x + 2 h) − f (x − 2 h) f ′ (x) = h f (x + h) − 2f (x) + f (x − h) f ′′ (x) = h2 f ′ (x) =

−

h ′′ · f (ξ) 2

h2 ′′′ · f (ξ) 24 h2 − · f (4) (ξ) 12 −

(Vorw¨artsdifferenzen) (zentrale Differenzen)

23

6 INTERPOLATION

Fehler Bei der numerischen Differentiation treten sowohl Rundungsfehler als auch Diskretisierungsfehler auf, die sich addieren: ˜ h − f ′′ (x)| ≤ |∆ ˜ h − ∆h | + |∆ | {z } Rundungsfehler

|∆h − f ′′ (x)| | {z }

≤ 4εh−2 + ch2

Diskretisierungsfehler

wobei im Allgemeinen der Rundungsfehler kleiner als der Diskretisierungsfehler sein sollte!