Dr. Alexander Veit Institut für Mathematik Universität Zürich

Numerik und wissenschaftliches Rechnen Alexander Veit Frühlingssemester 2013 Version: 11. April 2013

1

Inhaltsverzeichnis 1 Computerarithmetik 1.1 Zahlendarstellung . . . . . . . . . . . . . . . 1.2 Zahlendarstellung auf Computern . . . . . . 1.3 Rundung . . . . . . . . . . . . . . . . . . . 1.4 Rundungsfehler bei Gleitkommarechnungen 2 Approximation und Interpolation 2.1 Polynominterpolation . . . . . . . . . . . . 2.1.1 Lagrange Interpolation . . . . . . . 2.1.2 Interpolationsfehlerabschätzungen 2.1.3 Newtonsche Dividierte Differenzen 2.2 Interpolation mit linearen Splines . . . . .

. . . . .

. . . .

. . . . .

3 Numerische Differentiation und Integration 3.1 Numerische Differentiation . . . . . . . . . . . 3.2 Numerische Integration . . . . . . . . . . . . . 3.2.1 Zusammengesetzte Trapezregel . . . . 3.2.2 Zusammengesetzte Simpsonregel . . . 3.2.3 (Gewichtete) Newton-Cotes und Gauss 4 Nichtlineare Gleichungen 4.1 Bisektion . . . . . . . . . . 4.2 Lineares Eingabeln . . . . . 4.3 Das Newton-Verfahren . . . 4.4 Fixpunkt-Iterationen . . . . 4.5 Die ∆2 -Methode von Aitken

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

2

. . . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . .

. . . . .

. . . . . . . . . . . . . . . . . . . . Formeln

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

. . . . .

. . . . .

. . . . .

. . . .

4 4 5 6 7

. . . . .

10 10 11 12 14 18

. . . . .

21 21 23 24 25 26

. . . . .

36 37 38 39 44 45

Literatur • W. Gautschi. Numerical Analysis: An Introduction. • J. Stoer. Numerische Mathematik 1 • S. Sauter. Vorlesungsskript Numerik 1

3

1 Computerarithmetik Die meisten praktischen Probleme in der Mathematik lassen sich nicht explizit lösen. Die numerische Mathematik beschäftigt sich deshalb mit der Frage wie die Lösung solcher Probleme möglichst gut angenähert werden kann. Wir werden in dieser Vorlesung grundlegende Lösungsverfahren für verschiedene Problemstellungen kennenlernen. Um solche Lösungsverfahren zu beurteilen ist es wichtig die Grösse der Problemklasse zu kennen auf die sich das Verfahren anwenden lässt. Ausserdem spielt die Konvergenzgeschwindigkeit der Methode und der Implementierungsaufwand eine grosse Rolle. Oft stehen für ein gegebenes Problem mehrere numerische Verfahren zur Verfügung. Bei der Auswahl einer konkreten Methode müssen die oben stehenden Kriterien genau betrachtet und typischerweise gegeneinander abgewägt werden. Eine der wichtigsten Aufgaben der numerische Mathematik ist es, die Genauigkeit eines Rechenresultats (z.B. Ergebnis eines numerischen Verfahrens) zu beurteilen. Es gibt drei wesentliche Einflussfaktoren die diese Genauigkeit begrenzen: • Fehler in den Eingabedaten der Rechnung, • Rundungsfehler, • Approximationsfehler. Fehler in den Eingabedaten lassen sich typischerweise nicht vermeiden (z.B. Messwerte). Da auf einem Computer nur endlich viele Zahlen darstellbar sind und somit von reellen Zahlen auf Zahlen auf dem Computer übergegangen werden muss, entstehen zusätzlich Fehler, die sich nicht vermeiden lassen. Diese nennt man Rundungsfehler. Approximationsfehler sind abhängig von der verwendeten Methode die gewählt wird um ein gegebenes Problem näherungsweise zu lösen. Da hier typischerweise das gegebene Problem durch ein Problem ersetzt wird, welches sich einfacher lösen lässt, liefert selbst exakte Arithmetik nur zu einer fehlerbehafteten Annäherung. Diesen Fehler nennt man Approximationsfehler. In diesem Kapitel werden wir die Auswirkung von Eingangs- und Rundungsfehlern einer Rechnung auf das Endresultat untersuchen. Approximationsfehler lassen sich natürlich nicht in dieser Allgemeinheit untersuchen und müssen für jede Methode einzeln betrachtet werden.

1.1 Zahlendarstellung In diesem Unterkapitel werden wir zunächst die reellen Zahlen einführen. Dies kann auf verschiedenste Weise gemacht werden, im Zusammenhang mit Computern hat sich allerdings die Darstellung im Dualsystem (Binärsystem) als vorteilhaft herausgestellt. Jede Zahl x ∈ R kann demnach in der Form x = ±(bn 2n + bn−1 2n−1 + · · · + b0 + b−1 2−1 + b−2 2−2 + · · · )

(1.1)

geschrieben werden, wobei n > 0 eine natürliche Zahl ist und die bi die sogenannten Binärkoeffizienten sind, für die bi = 0 oder bi = 1 für alle i, gilt. Anstatt die Darstellung (1.1) werden wir im Folgenden die Kurzschreibweise x = ±(bn bn−1 · · · b0 . b−1 b−2 · · · )2

4

(1.2)

verwenden. Der Index 2 soll hier andeuten, dass die Zahl im Dualsystem vorliegt. Der Punkt in (1.2) wird Binärpunkt genannt und separiert den ganzzahligen Anteil links vom bruchzahligen Anteil rechts. Aufgabe 1.1 Zeigen Sie, dass die Darstellung (1.1) einer Zahl nicht eindeutig ist. Wie kann man Eindeutigkeit erreichen? Beispiel 1.1 (a) (10010.001)2 = 1·24 +0·23 +0·22 +1·21 +0·20 +0·2−1 +0·2−2 +1·2−3 = (18.125)10 P P∞ −k 1 −k = − 3 + (b) (0.01¯ 1)2 = ∞ = − 32 + 1−1/2 = (0.5)10 , wobei hier verwendet k=2 2 k=0 2 2 wurde, dass ∞ X 1 rk = 1−r k=0

falls |r| < 1 (geometrische Reihe). (c)

1 5

= (0.2)10 = (.00110011)2 .

Das letzte Beispiel zeigt, dass eine reelle Zahl mit endlicher Dezimaldarstellung eine unendliche Binärdarstellung haben kann.

1.2 Zahlendarstellung auf Computern Mit Hilfe der Binärdarstellung einer Zahl, lassen sich Gleitkommazahlen, so wie sie auf modernen Rechnern verwendet werden, einführen. Da auf Computern nur endlich viele Stellen einer Zahl gespeichert werden können, muss die maximale Anzahl an realisierbaren Stellen angegeben werden. Sei t die vom System erlaubte Anzahl von Stellen im bruchzahligen Teil der Zahl und s die zulässige Anzahl von Stellen im Exponenten. Wir bezeichnen die Menge der Gleitkommazahlen auf diesem Computer mit R(t, s), wobei x ∈ R(t, s)

x = f · 2e

⇐⇒

(1.3)

mit f = ±(. b−1 b−2 · · · b−t )2

und e = ±(bs−1 bs−2 · · · b0 .)2 .

Der bruchzahlige Teil f wird für gewöhnlich als Mantisse von x, die ganze Zahl e als Exponent von x bezeichnet. Die Zahl x in der Darstellung (1.3) heisst normalisiert, falls b−1 = 1. Es ist offensichtlich, dass die Menge der Gleitkommazahlen R(t, s) ⊂ R endlich ist. Abbildung 1.1 zeigt die Verteilung der Punkte auf der reellen Achse für die Spezialfälle t = 1, 2, 3, 4, s = 1 und Basis 4. Es wird ersichtlich, dass die Zahlen nicht gleichmässig auf der reellen Achse verteilt sind und dass sich 0 in diesem Format nicht darstellen lässt. Aus der Definition der (normalisierten) Gleitkommazahlen wird klar, dass die betragsmässig grösste und kleinste Zahl durch s −1

max |x| = (1 − 2−t )22 x∈R(t,s)

,

min |x| = 2−2

s

(1.4)

x∈R(t,s)

gegeben sind. Typische Werte für t und s sind t = 23 und s = 7. Damit ergibt sich für die betragsmässig grösste bzw. kleinste darstellbare Zahl 1.7 · 1038 bzw. 2.94 · 10−39 . Jede Zahl, die betragsmässig nicht in diesem Bereich liegt lässt sich auf einem solche Rechner nicht darstellen.

5

Abbildung 1.1: Exakt darstellbare Gleitkommazahlen der Form x = f · 4e für verschiedene Mantissenlängen t = 1, . . . , 4 und s = 1. Quelle: http: // de. wikipedia. org/ wiki/ Gleitkommazahl

Falls während einer Rechnung eine Zahl berechnet wird, die grösser als das Maximum in (1.4) ist, kommt es zu einem sogenannten Überlauf. Überlaufe sind fatal und führen typischerweise zu einem sofortigen Abbruch der Berechnung. Falls während einer Rechnung eine Zahl produziert wird, die kleiner als das Minimum in (1.4) ist, kommt es zu einem sogenannten Unterlauf. Dies ist für gewöhnlich weniger folgenreich als ein Überlauf, da die entsprechende Zahl systemintern zu 0 gerundet wird. In bestimmten Fällen kann jedoch auch das zu falschen Resultaten führen.

1.3 Rundung Wie oben beschrieben ist R(t, s) eine endliche nicht gleichmässig verteilte Menge auf der reellen Achse. Da während einer Rechnung zwangsläufig auch Zahlen produziert werden, die nicht in dieser Menge liegen (sich also nicht in der Form (1.3) darstellen lassen), müssen solche Zahlen auf vom System darstellbare Zahlen abgebildet werden. Dies erfolgt durch Rundung. Um dies genauer zu erläutern, betrachten wir eine “exakte” reelle Zahl ! ∞ X −k x ∈ R, x = ± b−k 2 2e . k=1

Diese soll auf die Zahl x? ∈ R,

t X

x? = ±

! b?−k 2−k

?

2e .

k=1

durch Rundung abgebildet werden. Eine Möglichkeit dies zu tun ist durch abschneiden: x? = chop(x),

wobei e? = e,

6

b?−k = b−k

(1.5)

für k = 1, 2, . . . , t. Der Rundungsfehler, der durch Abschneiden entsteht kann wie folgt abgeschätzt werden: ∞ ∞ X X ? −k e |x − x | = |x − chop(x)| = ± b−k 2 2 ≤ 2−k · 2e = 2−t · 2e . (1.6) k=t+1

k=t+1

Die Grösse |x − x? | ist der sogenannte absolute Fehler. Dieser Fehler hängt von der absoluten Grösse von x ab. Man geht deshalb oft zum sogenannten relativen Fehler |x − x? |/|x|, x 6= 0 über um diese Abhängigkeit zu eliminieren. Im Fall der Abschneideoperation kann dieser durch x − chop(x) 2−t · 2e 2−t · 2e ≤ P ≤ = 2 · 2−t (1.7) |± ( ∞ b−k 2−k ) 2e | −1 · 2e | x |2 k=1 abgeschätzt werden. Eine weitere Art des Rundens ist das symmetrische Runden. Dies entspricht dem gewöhnlich Auf- und Abrunden. Dies ist im Binärsystem aber einfacher, da es nur zwei Möglichkeiten gibt: Falls die (t+1)-ste Stelle eine 1 ist, wird aufgerundet, ansonsten abgerundet. Das symmetrische Runden lässt sich mit Hilfe der chop-Operation ausdrücken und lässt sich schreiben als:   x? = rd(x), wobei rd(x) = chop x + 2−(t+1) · 2e . (1.8) Der relative Fehler für das symmetrische Runden erfüllt die etwas bessere Abschätzung x − rd(x) ≤ 2−t (1.9) x Die rechte Seite in (1.9) hängt vom verwendeten Computer ab und wird als relative Maschinengenauigkeit eps := 2−t bezeichnet. Sie bestimmt die Genauigkeit aller Gleitkommarechnungen auf dem zugehörigen Rechner. So ergibt sich beispielsweise für t = 23 ein Wert von eps = 1.19 · 10−7 , was 6 bis 7 signigfikanten Dezimalstellen entspricht. Im folgenden werden wir die Auswirkung von Rundungsfehlern auf elementare mathematische Operationen betrachten. Dazu ist es vorteilhaft mit Gleichheiten zu arbeiten und nicht mit Ungleichheiten. Wir schreiben deshalb rd(x) = x(1 + ε),

|ε| ≤ eps .

1.4 Rundungsfehler bei Gleitkommarechnungen In diesem Abschnitt werden wir untersuchen wie sich Rundungsfehler auf die elementaren arithmetischen Operationen +, −, ×, / auswirken. Jede dieser Operationen kann, wenn sie auf zwei Zahlen in R(t, s) angewendet wird, ein Resultat produzieren, das nicht mehr auf dem Computer darstellbar ist (also nicht in R(t, s) enthalten ist). Da dieses Resultat dann gerundet werden muss, ist jede elementare arithmetische Operation typischerweise fehlerbehaftet. Der Fehler, den eine einzelne arithmetische Operation verursacht, ist häufig klein und kann vernachlässigt werden. Da in den meisten numerische Berechnungen jedoch eine grosse Anzahl von Operationen aneinandergereiht werden, stellt sich die Frage ob sich solche Fehler kritisch

7

anhäufen können bzw. wie sich solche Fehler fortpflanzen (Fehlerfortpflanzung). Im Folgenden betrachten wir wie Fehler in den Eingabedaten (hier durch Rundung) sich durch die einzelnen arithmetischen Operationen fortpflanzen. Wie schliessen Über- und Unterlauf aus und nehmen an, dass die arithmetischen Operationen exakt ausgeführt werden. Seien x, y ∈ R und deren interne Darstellung sei gegeben durch x ˜ = x(1 + εx ) und y˜ = y(1 + εy ) mit den durch Rundung verursachten Fehlern εx und εy . Weiterhin nehmen wir an, dass εx und εy hinreichen klein sind, sodass quadratische Terme ε2x , εx εy , ε2y (sowie Terme höherer Ordnung) vernachlässigbar klein werden. • Multiplikation Es gilt x ˜ · y˜ = x(1 + εx ) · y(1 + εy ) = x · y(1 + εx + εy + εx εy ) ≈ x · y(1 + εx + εy ). Somit ist der relative Fehler des Produkts εx·y ungefähr gegeben durch εx·y = εx + εy . Der relative Fehler des Produkts ist also die Summe der Eingabefehler. Diese Fehlerfortpflanzung ist akzeptabel und wir nennen die Multiplikation deshalb eine gutartige Operation. • Division Sei y 6= 0, dann gilt mittels Taylorentwicklung x ˜ x(1 + εx ) x x = = (1 + εx )(1 + εy + ε2y − · · · ) ≈ (1 + εx − εy ). y˜ y(1 + εy ) y y Somit ist εx/y = εx − εy und damit ist die Division ebenfalls eine gutartige Operation. • Addition und Subtraktion Die kann als ein Fall betrachtet werden, da wir die Vorzeichen von x und y nicht festgelegt haben. Es gilt   x y x ˜ + y˜ = x(1 + εx ) + y(1 + εy ) = x + y + xεx + yεy = (x + y) 1 + εx + εy , x+y x+y falls x + y 6= 0. Somit εx+y =

y x εx + εy x+y x+y

(1.10)

Wie zuvor ist der Fehler eine Linearkombination aus den einzelnen Fehlern. Im Gegensatz zu oben können hier die Koeffizienten aber beliebig gross werden. Haben x und y gleiches x Vorzeichen (es sich also um eine Addition handelt) gilt | x+y | ≤ 1 und wir können den relativen Fehler betragsmässig durch |εx+y | ≤ |εx | + |εy | abschätzen. Im Fall der Addition handelt es sich deshalb wieder um eine gutartige Operation.

8

Falls x und y jedoch unterschiedliches Vorzeichen haben und betragsmässig ungefähr gleich gross sind (d.h. |x + y| ist klein im Verhältnis zu |x| und |y|), können die Koeffizienten in (1.10) und somit der relative Fehler beliebig gross werden. Dies nennt man Auslöschung. Die einzige rundungsfehlerkritische arithmetische Operation ist demnach die Subtraktion. Der Auslöschungseffekt, der dabei auftritt kann verheerend sein und numerische Berechnungen unbrauchbar machen. Beispiel 1.2 Wir betrachten einen Computer, der 2 Dezimalstellen einer Zahl abspeichern kann und symmetrisch rundet. Die Identität (a − b)2 = a2 − 2ab + b2 , die zweifellos gilt wenn man in der Menge der reellen Zahlen rechnet, stimmt auf unserem Computer nicht mehr. Sei a = 1.8 und b = 1.7 dann gilt (a − b)2 = 0.01 auf unserem Computer und mit exakter Arithmetik, jedoch gilt a2 − 2ab + b2 = 3.2 − 6.2 + 2.9 = −0.1 auf dem Computer, was mit dem echten Ergebnis nicht einmal das Vorzeichen gemeinsam hat. Hier tritt Auslöschung auf, da die Zahlen a2 + b2 und 2ab intern auf eine Nachkommastelle gerundet werden und dann voneinander abgezogen werden.

9

2 Approximation und Interpolation In diesem Kapitel werden wir uns mit der Approximation von Funktionen beschäftigen. Die Funktionen, die wir dabei im Sinn haben sind auf einem endlichen Intervall, oder sogar nur auf einer endlichen Menge von Punkten definiert. Der erste Fall tritt in der Praxis auf, wenn spezielle Funktionen als Teil eines Algorithmus ausgewertet werden müssen. Da solche Funktionen typischerweise einen unendlichen Definitions-und Wertebereich haben, geht man dazu über sie durch eine endliche Anzahl von einfacheren Funktionen wie Polynome oder rationale Funktionen zu approximieren. Der zweite Fall tritt auf, wenn Messungen z.B. einer physikalischen Grössen als Funktion einer anderen physikalischen Grösse (wie der Zeit) vorgenommen werden. Auch hier will man häufig die unbekannte kontinuierliche Funktion mit einfachen Funktionen approximieren. Im Folgenden wollen wir eine Funktion f ∈ C ∞ ([a, b]), wobei a < b, mit Funktionen aus einem Raum Φ ⊂ C ∞ ([a, b]), der aus einfacheren Funktionen besteht, approximieren. Um zu beurteilen wie gut eine solche Approximation ist, führen wir die Unendlich-Norm kf k∞ = max |f (t)| a≤t≤b

(2.1)

auf C ∞ ([a, b]) ein. Falls f nur auf einer endlichen Menge ∆ := {xi , i = 1, . . . , N } definiert ist, ersetzen wir die kontinuierliche Norm (2.1) durch die diskrete Norm kf k∞ = max |f (ti )|. xi ∈∆

Das Ziel ist es im Allgemeinen eine approximierende Funktion φ ∈ Φ zu finden, sodass der Fehler kf − φk∞ möglichst klein ist. Dabei ist die Wahl des Approximationssystems Φ und die Wahl der Approximationsmethode massgeblich. Ein häufige Wahl für Φ sind Polynome (Kapitel 2.1), Splines (Kapitel 2.2), trigonometrische Polynome und rationale Funktionen. Bei der Wahl der Approximationsmethode werden wir uns auf Interpolation beschränken. Polynome gehören zu den am häufigsten verwendeten Funktionensystemen zur Approximationen von ‘allgemeinen’ Funktionen auf beschränkten Gebieten. Neben der einfachen Handhabung ist der Weierstrasssche Approximationssatz ein weiterer Grund dafür. Satz 2.1 Sei I ein endliches, abgeschlossenes Intervall und C 0 (I) der Raum der stetigen Funktionen auf I. Für alle f ∈ C 0 (I) und alle ε > 0 existiert ein n ∈ N0 und ein p ∈ Pn , sodass kf − pk∞ < ε, wobei Pn den Raum aller Polynome vom Maximalgrad m bezeichnet.

2.1 Polynominterpolation In diesem Abschnitt wollen wir eine gegebene Funktion f ∈ C 0 ([a, b]) mit Hilfe von Polynomen approximieren. Sei dazu ( ) n X i Pn := ϕ : ϕ(t) = ci t , ci ∈ R i=0

der Raum der Polynome vom Maximalgrad n. Wir wollen f nun so durch ein Polynom p approximieren, dass die Funktionswerte von f und p in gewissen Punkten übereinstimmen. Sei Θn := {xi , i = 0, . . . , n} ⊂ [a, b]

10

eine Menge von n + 1 verschiedenen Punkte in [a, b]. Die Interpolationsaufgabe besteht nun darin ein Polynom p ∈ Pn zu finden, sodass p(xi ) = f (xi ) für 0 ≤ i ≤ n

(2.2)

gilt. Ein Polynom, welches Bedingung (2.2) erfüllt interpoliert die Funktion f in den Punkten xi . Die Fragen, die sich nun stellen sind, ob ein solches Polynom immer existiert und wenn ja, ob dieses eindeutig ist. Ausserdem muss beurteilt werden, ob ein Polynom, das (2.2) erfüllt tatsächlich eine gute Approximation von f im Intervall [a, b] ist. Dies ist gleichbedeutend mit der Frage wie sich der Fehler der Approximation kf − pk∞ abschätzen lässt. Beispiel 2.1 Wir betrachten den einfachsten Fall n = 1. Hier wird eine gegebene Funktion f ∈ C 0 ([a, b]) mittels einer linearen Funktion approximiert. Seien x0 , x1 ∈ [a, b] zwei von einander verschiedene Punkte. Es lässt sich leicht berechnen, dass das Polynom p ∈ P1 , gegeben durch x − x1 x − x0 p(x) = f (x0 ) + f (x1 ). (2.3) x0 − x1 x1 − x0 die Funtion f in x0 und x1 interpoliert. Eine andere Darstellung von p ist gegeben durch p(x) = f (x0 ) +

f (x1 ) − f (x0 ) (x − x0 ). x1 − x0

(2.4)

Die Darstellung (2.3) entspricht der Lagrange-Darstellung (Kapitel 2.1.1), die Darstellung (2.4) entspricht der Newtonschen Darstellung (Kapitel 2.1.3) von p. Die wollen wir im Folgenden präzisieren. 2.1.1 Lagrange Interpolation Wir zeigen die Existenz des Interpolationspolynoms (2.2) durch direkte Konstruktion. Dazu seien n Y x − xj li (x) := ∈ Pn , für i = 0, 1, . . . , n. (2.5) xi − xj j=0,j6=i

Diese Polynome vom Grad n erfüllen li (xk ) = δik

( 1 falls i = k, = 0 falls i 6= k.

(2.6)

Durch Multiplizieren dieser Polynome mit den Funktionswerten fi := f (xi ) und aufsummieren erhält man das sogenannte Lagrange Interpolationspolynom pn (f, x) := pn (f, Θn , x) :=

n X

fi li (x) ∈ Pn .

(2.7)

i=0

Dieses erfüllt Bedingung (2.2), da pn (f, xk ) =

n X

fi li (xk ) =

i=0

n X

fi δik = fk .

(2.8)

i=0

Die Existenz des Interpolationspolynoms ist somit gezeigt. Weiterhin kann es in der Form (2.7) dargestellt werden.

11

Satz 2.2 Das Interpolationspolynom ist durch (2.2) eindeutig bestimmt. Beweis. Angenommen es gäbe zwei Polynome vom Grad n, die (2.2) erfüllen, d.h. ∃p, p˜ ∈ Pn mit p(xi ) = p˜(xi ) = fi für i = 0, . . . , n. Dann hätte das Polynom P = p − p˜ ∈ Pn höchstens den Grad n und mindestens n + 1 verschiedene Nullstellen xi , i = 0, . . . , n. Das ist allerdings nicht möglich ausser wenn P auf dem ganzen Intervall [a, b] verschwindet, d.h. P ≡ 0. Die ist aber gleichbedeutend mit p = p˜.  2.1.2 Interpolationsfehlerabschätzungen In diesem Abschnitt wollen wir den Interpolationsfehler f (x) − pn (f, x) für ein x ∈ [a, b] genauer untersuchen. Im Allgemeinen können Fehlerdarstellungen und Fehlerabschaätzungen numerischer Methoden nur dann bewiesen werden, wenn bestimmte Annahmen an die zu approximierende Funktion f gestellt werden. Hier werden wir annehmen, dass die Funktion f (n + 1)-mal stetig differenzierbar ist. Damit lässt sich der folgende Satz beweisen: Satz 2.3 (Interpolationsfehlerdarstellung) Sei f ∈ C n+1 ([a, b]) und x ∈ [a, b]. Dann gilt f (n+1) (ξx ) ωn (x), (n + 1)!

f (x) − pn (f, Θn , x) = mit dem Stützstellenpolynom ωn (x) :=

n Y

(x − xj )

(2.9)

(2.10)

j=0

und einer (unbekannten) Zwischenstelle ξx ∈ (a, b), die von x abhängt. Die Fehlerdarstellung (2.9) ist nicht überraschend. Da der Fehler der Interpolation an den Stützstellen xi ∈ Θn verschwindet (siehe (2.8)), ist der Faktor ωn (x) sinnvoll. Da ausserdem pn (f, x) = f (x) für alle f ∈ Pn aufgrund der Eindeutigkeit des Interpolationspolynoms gilt, rechfertigt dies weiterhin das Auftauchen der (n+1)-sten Ableitung von f . Um Satz 2.3 formal beweisen zu können, benötigen wir den folgenden Hilfssatz. Satz 2.4 (Satz von Rolle) Sei f ∈ C 0 ([a, b]) und f ∈ C 1 ((a, b)). Falls f (a) = f (b) gilt, dann existiert ein ξ ∈ (a, b) s.d. f 0 (ξ) = 0. Beweis von Satz 2.3. Falls x ∈ Θn ist die Behauptung trivialerweise wahr. Sei also x ∈ [a, b]\Θn fest gewählt. Wir definieren die Funktion F (t) := f (t) − pn (f, Θn , t) −

f (x) − pn (f, Θn , x) ωn (t) ωn (x)

Da f ∈ C n+1 ([a, b]), gilt auch F ∈ C n+1 ([a, b]). Weiterhin gilt F (xi ) = 0,

für i = 0, . . . , n und

F (x) = 0.

12

Somit hat F n + 2 Nullstellen in [a, b]. Wiederholtes Anwenden des Satz von Rolle ergibt: F0 F .. .

hat mindestens n + 1 verschiedene Nullstellen,

00

hat mindestens .. .

n .. .

verschiedene Nullstellen, .. .

F (n+1) hat mindestens

1

Nullstelle.

Wir bezeichnen mit ξx die Nullstelle von F (n+1) und differenzieren F (n + 1)-mal bezüglich t. Auswertung an ξx ergibt 0 = f (n+1) (ξx ) −

f (x) − pn (f, Θn , x) (n + 1)!, ωn (x) 

was die Aussage beweist. Korollar 2.1 Sei f ∈ C n+1 ([a, b]) und x ∈ [a, b]. Dann gilt (grobe) Fehlerschranke kf − pn (f, Θn , ·)k∞ =

kf (n+1) k∞ (b − a)n+1 . (n + 1)!

Beweis. Die Aussage folgt direkt aus Satz 2.3.



Beispiel 2.2 Im Fall der linearen Interpolation (n = 1), gilt für a < x < b und h := b − a die Fehlerdarstellung f (2) (ξx ) f (x) − p1 (f, x) = (x − a)(x − b). 2 Eine einfache Rechnung ergibt die Abschätzung kf − p1 (f, ·)k∞ ≤

kf (2) k∞ kf (2) k∞ 2 · k(x − a)(x − b)k∞ = h . 2 8

Aufgabe 2.1 Leiten Sie für die quadratische Interpolation mit äquidistanten Stützstellen x0 , x1 = x0 + h, x2 = x0 + 2h eine obere Schranke für kf − p2 (f, ·)k∞ her, welche kf (3) k∞ und h beinhaltet. (Hier sei kuk∞ = maxx0 0,

a

ist dies unmöglich.

28

• Die Wahl k = n + 1 ist also in dem Sinn optimal, dass sie eine Quadraturformel (basierend auf n + 1 Stützstellen) mit maximalem Exaktheitsgrad d = 2n + 1 liefert. Solche Formeln werden als Gausssche Quadraturformeln bezeichnet. Beweis von Satz 3.3. • Wir zeigen zunächst die Hinrichtung: Falls die Quadraturformel Exaktheitsgrad d = n+k hat, dann gilt (a) und (b). Da die Formel Exktheitsgrad n + k hat, ist sie trivialerweise interpolatorisch. Somit gilt (a). Sei nun ein p ∈ Pk−1 gegeben. Dann gilt ωn (t)p(t) ∈ Pn+k . Da die Quadraturformel nach Voraussetzung Polynome vom Grad n + k exakt integriert gilt Z b n X ωn (t)p(t)w(t)dt = wj ωn (tj )p(tj ) = 0, a

j=0

da ωn (tj ) = 0 für 0 ≤ j ≤ n. Somit gilt (b). • Wir zeigen nun die Rückrichtung: Falls (a) und (b) gilt, hat die Quadraturformel Exaktheitsgrad d = n + k. Es muss gezeigt werden, dass ein Polynom p ∈ Pn+k exakt integriert wird. Wir verwenden Polynomdivision mit Rest und schreiben p(t) = q(t)ωn (t) + r(t) mit q ∈ Pk−1 und r ∈ Pn . Damit gilt Z b Z b Z b p(t)w(t)dt = q(t)ωn (t)w(t)dt + r(t)w(t)dt a

a

a

Wegen Voraussetzung (b) verschwindet das erste Integral auf der rechten Seite. Für das zweite Integral ergibt sich mit Voraussetzung (a) und r ∈ Pn , dass Z

b

r(t)w(t)dt = a

n X

wj r(tj ) =

j=0

n X

wj [p(tj ) − q(tj )ωn (tj )] =

j=0

n X

wj p(tj ),

j=0

wobei wieder verwendet wurde, dass die tj Nullstellen des Stützstellenpolynoms sind. Somit gilt also Z b n X p(t)w(t)dt = wj p(tj ), a

j=0

d.h. das Polynom wird exakt integriert.  Beispiel 3.6 Wir möchten Newton-Cotes und Gauss Formeln für n = 2 miteinander vergleichen. Dazu betrachen wir Integrale des Typs Z 1 t−1/2 f (t)dt. 0

29

Bei den Newton-Cotes Formeln sind die Stützstellen vorgeschrieben und werden hier als Endpunkte des Intervalls gewählt. Die beiden Quadraturformeln sind also von der Form C NC NC QN 1 (f ) = w0 f (0) + w1 f (1) G G QG 1 (f ) = w0 f (t0 ) + w1 f (t1 )

Um die Gewichte der Newton-Cotes Formel zu berechnen wird wieder Formel det. Dies liefert mit w(t) = t−1/2 : Z 1 Z 1 Z 1 −1/2 t − 1 −1/2 NC t t l0 (t)dt = w0 = t−1/2 − t1/2 dt = dt = 0−1 0 0 0 Z 1 Z 1 Z 1 2 −1/2 t − 0 −1/2 G t t1/2 dt = . t l1 (t)dt = dt = w1 = 1−0 3 0 0 0

(3.6) angewen4 , 3

Somit lautet die zugehörige 2-Punkt Newton-Cotes Formel 2 4 C QN 1 (f ) = f (0) + f (1) 3 3 Aufgrund der Singularität in der Gewichtsfunktion wird der Wert f (0) doppelt gewichtet. Um die zwei Stützstellen für die Gauss Quadratur zu bestimmen, muss das zugehörige Stützstellenpolynom orthogonal auf dem Raum P1 bezüglich w sein. Sei das Stützstellenpolynom in der Form π1 (t) = t2 + p1 t + p2 gegeben. Dieses muss orthogonal zur konstanten Funktion 1, sowie zu t sein. Somit müssen die Gleichungen Z 1 Z 1 2 2 −1/2 0= t π1 (t)dt = (t3/2 + p1 t1/2 + p2 t−1/2 )dt = − p1 + 2p2 5 3 0 0 Z 1 Z 1 2 2 2 0= t−1/2 tπ1 (t)dt = (t5/2 + p1 t3/2 + p2 t1/2 )dt = − p1 + p2 7 5 3 0 0 erfüllt sein. Die führt auf das lineare Gleichungssystem 1 1 p1 − p2 = 3 5 1 1 1 p1 − p2 = , 5 3 7 dessen Lösung durch p1 = 67 und p2 = stellenpolynom ist also durch

3 35

gegeben ist. Das Orthogonalpolynom bzw. Stütz-

6 3 π1 (t) = t2 + t + 7 35 gegeben. Die Stützstellen der Gauss-Quadratur sind die Nullstellen von π1 . Somit r ! 1 6 t0 = 3−2 = 0.1155871100 7 5 r ! 1 6 t1 = 3+2 = 0.7415557471, 7 5

30

Um die Gewichte von QG 1 (f ) zu bestimmen, kann wieder Formel (3.6) verwendet werden. Es ist allerdings einfacher erneut ein lineares Gleichungssystem aufzustellen und auszunutzen, dass die Formel für f (t) = 1 und f (t) = t exakt ist. QG 1 (1)

=

w0G

! w1G =

Z

! t1 w1G =

Z

+

1

t−1/2 dt = 2

0 G QG 1 (t) = t0 w0 +

1

0

2 t−1/2 · t dt = , 3

Dies ergibt −2t1 + 32 =1+ = t0 − t1 −2t0 − 32 =1− w1G = t0 − t1

w0G

r 1 5 = 1.304290310 3 6 r 1 5 = 0.6957096903, 3 6

Für die 2-Punkt Gauss Formel ergibt sich somit r ! r !! 1 5 1 6 G Q1 (f ) = 1 + f 3−2 + 3 6 7 5

1 1− 3

r ! 5 f 6

1 7

r !! 6 3+2 5

Auch hier ist das Gewicht w0G grösser als w1G um der Singularität der Gewichtsfunktion Rechnung zu tragen. Wir illustrieren die Formeln für die Funktion f (t) = cos(0.5πt2 ). Es gilt Z I=

1

t−1/2 cos(0.5πt2 )dt = 1.5597865

0

Die Newton-Cotes und Gauss Formel liefern folgende Approximationen an I: 4 = 1.333 . . . 3 = 1.55758956,

INC = IG Für die Fehler ergibt sich

E2N C = 0.226 bzw. E2G = 0.00220. Dies zeigt die Überlegenheit von Gauss Formeln sogar für n = 2. Im Folgenden werden wir einige nützliche Eigenschaften von Gauss Quadraturformeln kennenlernen. Bemerkung 3.2 (a) Die Stützstellen tk von Gaussschen Quadraturformeln sind reel, verschieden und im offenen Intervall (a, b) enthalten. (b) Die Gewichte wk sind positiv. Beweis. (a) Seien tk , 0 ≤ k ≤ n die Stützstellen der Gauss Formel und ωn das zugehörige Stützstellenpolynom. Seien weiterhin a < tk0 < tk1 < · · · < tkl < b diejenigen Stützstelln in

31

denen ωn das Vorzeichen wechselt. Wir zeigen l = n durch Widerspruch. Angenommen l < n und sei l Y q(t) = (t − tkl ) ∈ Pl+1 . j=0

Rb Aufgrund der Orthogonalität von ωn zu q bezüglich der Gewichtsfunktion w gilt a ωn (t)q(t)w(t)dt = 0. Andererseits haben wir q so konstruiert, dass die Funktion ωn (t)q(t)w(t) ihr Vorzeichen in Rb (a, b) nicht wechselt, weshalb a ωn (t)q(t)w(t)dt 6= 0 gilt. Dies ist ein Widerspruch, sodass die Annahme l < n falsch gewesen sein muss. Somit gilt l = n und die Stützstellen sind verschieden und im Intervall (a, b) enthalten. (b) Wir bezeichnen mit li ∈ Pn wieder die Lagrange Basisfunktionen. Da die Gauss Formeln mit (n + 1) Stützstellen den Exaktheitsgrad 2n + 1 besitzen gilt b

Z

2

0
0 für i = j, dann bilden die πk ein Orthogonalsystem von Polynomen über [a, b] bezüglich der Gewichtsfunktion w. Bemerkung 3.4 Die orthogonalen Polynome π0 , . . . , πk bilden eine Basis des Pk . Deshalb gilt k−1 X a(pk , q) = a(pk , αi p i ) = 0 j=0

für alle q ∈ Pk−1 . Ein Orthogonalsystem kann mit folgender Rekursion erzeugt werden. Satz 3.5 Sei durch (3.9) ein Skalarprodukt gegeben, wobei w eine postive, beschränkte Gewichtsfunktion ist. Ein Orthogonalsystem von Polynomen πk ∈ Pk kann durch die 3-Term Rekursion π−1 (t) = 0,

π0 (t) = 1

πk+1 (t) = (t − αk )πk (t) − βk πk−1 (t),

k = 0, 1, . . .

(3.10)

erzeugt werden, wobei a(tπk , πk ) , k = 0, 1, . . . a(πk , πk ) a(πk , πk ) βk = , k = 1, 2, . . . a(πk−1 , πk−1 ) Z b β0 = w(t)dt.

αk =

a

Beweis. Wir zeigen durch Induktion, dass π0 , π1 , · · · , πk+1 ein Orthogonalsystem bilden. Zunächst gilt: a(π1 , π0 ) = a((t − α0 )π0 (t) − β0 π−1 (t), π0 ) = a(t − α0 , 1) a(t, 1) a(1, 1) = 0. = a(t, 1) − α0 a(1, 1) = a(t, 1) − a(1, 1)

34

Somit ist π1 orthogonal zu π0 . Seien nach Induktionsvoraussetzung nun also π0 , π1 , · · · , πk paarweise orthogonal. Wir müssen wir zeigen, dass dann auch πk+1 zu diesen Polynomen orthogonal ist. Sei j ≤ k − 2. Dann gilt a(πk+1 , πj ) = a((t − αk )πk − βk πk−1 , πj ) = a(tπk , πj ) − αk a(πk , πj ) − βk a(πk−1 , πj ) = a(πk , tπj ) = 0, wobei die letzte Gleichheit aus tπj ∈ Pk−1 und Bemerkung 3.4 folgt. Ausserdem a(πk+1 , πk−1 ) = a((t − αk )πk − βk πk−1 , πk−1 ) = a((t − αk )πk , πk−1 ) − a(πk , πk ) = a(πk , tπk−1 ) − αk a(πk , πk−1 ) − a(πk , πk ) = a(πk , πk + αk−1 πk−1 + βk−1 πk−1 ) − a(πk , πk ) = a(πk , αk−1 πk−1 + βk−1 πk−1 ) =0 und a(πk+1 , πk ) = a((t − αk )πk − βk πk−1 , πk ) = a(tπk , πk ) − αk a(πk , πk ) = 0. Somit sind auch π0 , π1 , · · · , πk+1 paarweise orthogonal. Dies zeigt die Behauptung.



Beispiel 3.7 Für die Gewichtsfunktion w(t) = 1 und a = −1, b = 1 ergibt die 3-Term Rekursion die Polynome P0 (t) = 1,

P1 (t) = t,

1 3 P2 (t) = t2 − , 2 2

5 3 P3 (t) = t3 − t, . . . . 2 2

Diese werden Legendre Polynome genannt. Die Nullstellen des Polynoms Pk entsprechen den Stützstellen der k-Punkt Gauss Formel im Intervall [−1, 1]. Seien t0 , . . . , tk die Stützstellen und w0 , . . . , wk die Gewichte dieser Formel. Ein Integral auf einem beliebigen Intervall [a, b] kann durch Variablentransformation auf das Interval [−1, 1] abgebildet werden, wo sich dann die Quadraturformel anwenden lässt: Z a

b

b−a f (t)dt = 2

Z

1

 f

−1

b−a b+a t+ 2 2



k

b−aX dt ≈ wk f 2 j=0



b−a b+a tk + 2 2

 .

Bemerkung 3.5 Die Nullstellen der orthogonalen Polynome aus Satz 3.5 sind gerade die Stützstellen der Gauss Quadraturformel. Es lässt sich zeigen, dass diese Nullstellen den Eigenwerten bestimmter tridiagonaler Matrizen (die von αk , βk abhängen) entsprechen. Diese Eigenwerte lassen sich sehr effizient und stabil berechnen.

35

4 Nichtlineare Gleichungen In diesem Kapitel wollen wir Gleichungen der abstrakten Form f (x) = 0

(4.1)

lösen. Im einfachsten Fall ist (4.1) eine einzelne Gleichung einer Variablen x, wobei f dann eine gegebene Funktion darstellt. Die Aufgabe besteht nun darin Werte von x zu finden, für die die Funktion f verschwindet. Solche Werte werden Nullstellen der Funktion f genannt. Im Fall f : Rd → Rd , stellt (4.1) ein System von Gleichungen dar, wobei jede Komponente von f eine Funktion von d Variablen x1 , x2 , . . . , xd ist. Ein solches System wird nichtlinear genannt, falls mindestens eine Komponenten von f nichtlinear von mindestens einer der Variablen abhängt. Da selbst die einfachsten nichtlinearen Gleichungen im Allgemeinen nicht exakt gelöst werden können, müssen iterative Methoden entwickelt werden, welche die Nullstellen der Gleichung approximieren. Dazu werden Folgen von Näherungswerten x0 , x1 , . . . , xn , . . . nach einer Iterationsvorschrift xn+1 = ϕ(xn ) berechnet, sodass die Folge gegen eine Lösung α der Gleichung (4.1) konvergiert, d.h. lim xn = α. n→∞

Fall man eine Iterationsvorschrift ϕ gefunden hat, bricht man die Iteration für ein hinreichend grosses n ab und erhält eine Approximation an α. Um eine solche iterative Methode zu beurteilen ist der Begriff der Konvergenzordnung wesentlich. Definition 4.1 (Lineare Konvergenz) Die Folge {xn }∞ n=1 konvergiert linear gegen α, falls |xn − α| ≤ n , wobei für die Folge {n }∞ n=1 n+1 = c, n→∞ n lim

mit

0 0 (4.3) n→∞ n gilt. Die Definitionen 4.1 und 4.2 decken nicht alle Arten von Konvergenz ab. Es kann beispielsweise vokommen, dass die Folge {n }∞ n=1 langsamer als linear gegen Null konvergiert, sodass (4.2) nur mit c = 1 erfüllt ist. In diesem Fall heisst die Konvergenz sublinear. Weiterhin kann (4.2) auch mit c = 0 erfüllt sein, Bedingung (4.3) jedoch für kein p > 1 gelten. In diesem Fall heisst die Konvergenz superlinear. In praktischen Berechnungen muss ein geeignetes Stopp-Kriterium die Iteration abbrechen sobald eine gewisse Genauigkeit der Approximation erreicht wurde. Im Idealfall würde man

36

die Iteration abbrechen, falls der Fehler |xn − α| kleiner als eine vorgegebene Toleranz ist. Da α jedoch nicht bekannt ist, ersetzt man den Fehler |xn − α| überlicherweise durch |xn − xn−1 | und fordert |xn − xn−1 | ≤ tol, wobei tol = |xn |r + a und r , a vorgegebene Toleranzen sind. r kontrolliert den relativen Fehler, a kontroliert den absoluten Fehler der Approximation.

4.1 Bisektion Die Bisektion ist die einfachste Möglichkeit Nullstellen einer nichtlinearen Funktion zu berechnen. Die Methode generiert eine Folge geschachtelter Intervalle [an , bn ], n = 0, 1, . . . wobei jedes der erzeugten Intervalle mindestens eine Nullstelle der Funktion enthält. Für n → ∞ strebt die Länge dieser Intervalle gegen Null, sodass im Grenzwert eine (isolierte) Nullstelle berechnet wird. Wir nehmen an, dass zwei Zahlen a, b mit a < b gegeben sind, sodass f ∈ C([a, b]),

f (a) < 0,

f (b) > 0.

(4.4)

Für die Bisektion is wesentlich, dass f and den Intervallenden unterschiedliche Vorzeichen hat. Da f stetig sein muss, garantiert Bedingung (4.4) die Existenz (mindestens) einer Nullstelle im Intervall [a, b]. Bemerkung 4.1 Der Fall f (a) > 0, f (b) < 0. kann durch Multiplikation der zu lösenden Gleichung mit -1 auf den obigen Fall zurückgeführt werden. Um eine Nullstelle der Funktion im Intervall [a, b] nun zu approximieren, geht mal wie folgt vor: Man halbiert sukzessive das Intervall und sucht jeweils in denjenigen Teilintervallen weiter, welche (4.4) erfüllen. Auf diese Weise wird eine Folge von geschachtelten Intervallen erzeugt, deren Länge jeweils halbiert wird und welche eine Nullstelle von f enthalten. Algorithmus 4.1 Es gelte (4.4), Setze: a1 = a, b1 = b. for n = 1, 2, 3, . . . do xn = 12 (an + bn ) if f (xn ) < 0 then an+1 = xn , bn+1 = bn else an+1 = an , bn+1 = xn endif endfor Nach n Schritten liefert dieser Algortihmus ein Teilintervall [an , bn ] ⊂ [a, b] zurück, in dem die unbekannte Nullstelle α enthaltens sein muss. Sei xn der Mittelpunkt dieses Intervalls, dann kann die tatsächliche Nullstelle nur um 12 (bn − an ) von xn abweichen. Somit gilt b−a 1 |xn − α| ≤ (bn − an ) = n 2 2

37

Somit konvergiert die durch Bisektion konstruierte Folge xn gegen die Nullstelle α mit n = b−a 2n . Da b−a 1 n+1 n+1 = 2b−a = für alle n, n 2 2n konvergiert die Bisektion linear gegen α. Ist nun eine absolute Fehlertoleranz tol > 0 vorgegeben, ist der Fehler der Approximation |xn − α| kleiner oder gleich der Toleranz, falls b−a ≤ tol. 2n Auflösen nach n liefert



 b−a n > log2 , tol

wobei d·e Aufrunden auf die nächste Ganzzahl bedeutet. Man weiss somit von vornherein wieviele Schritte nötig sind um eine Gewisse Genauigkeit zu erreichen.

4.2 Lineares Eingabeln Wie bei der Bisektion nehmen wir an, dass es zwei Zahlen [a, b] gibt, sodass f ∈ C([a, b]),

f (a)f (b) < 0

(4.5)

und konstruieren geschachtelte Teilintervalle [an , bn ], n = 1, 2, . . ., mit a1 = a, b1 = b, sodass die Eigenschaft f (an )f (bn ) < 0 erhalten bleibt. Im Gegensatz zum Bisektionsverfahren wird nicht der Mittelpunkt von [an , bn ] gewählt um das neue Intervall zu bestimmen, sondern die Nullstelle des linearen Interpolanten von f an den Stützstellen an und bn . D.h. xn wird so gewählt, dass p1 (f, {an , bn }, xn ) = 0. Algorithmisch geht man wie folgt vor: Algorithmus 4.2 Es gelte (4.5), Setze: a1 = a, b1 = b. for n = 1, 2, 3, . . . do −bn xn = an − f (aann)−f (bn ) f (an ) if f (xn )f (an ) > 0 then an+1 = xn , bn+1 = bn else an+1 = an , bn+1 = xn endif endfor Im Gegensatz zur Bisektionsmethode ist die Anzahl an Iterationsschritten, die benötigt wird um eine gegebene Toleranz zu erreichen, a priori nicht bekannt. In einem Computerprogramm muss deshalb die Anzahl an Iterationen beschränkt werden. Um das Konvergenzverhalten des linearen Eingabelns zu analysieren nehmen wir der Einfachheit halber an, dass die Funktion im Intervall [a, b] konvex oder konkav ist.

38

Definition 4.3 Eine Funktion f : [a, b] → R heisst konvex, wenn für alle x, y ∈ [a, b] und für alle t ∈ [0, 1] f (tx + (1 − t)y) ≤ tf (x) + (1 − t)f (y) gilt. Die Funktion heisst konkav, falls f (tx + (1 − t)y) ≥ tf (x) + (1 − t)f (y). Bemerkung 4.2 Sei f ∈ C 2 ([a, b]). Die Funktion f ist konvex (konkav) genau dann wenn f 00 (x) > 0 (f 00 (x) < 0) im Intervall [a, b]. Wir betrachten im Folgenden den Fall einer konvexen Funktion f ∈ C 2 ([a, b]) mit f 00 (x) > 0 für a ≤ x ≤ b und f (a) < 0, f (b) > 0. Die (einzige) Nullstelle von f bezeichnen wir wieder mit α. Der lineare Interpolant, der f (a) und f (b) verbindet, liegt komplett überhalb des Graphen von f und schneidet die x-Achse deshalb links von α. Da dies auch für die folgenden Sekanten der Fall sein wird, bleibt der Punkt b für alle Interationsschritte fest, während der andere Endpunkt nach jedem Schritt aktualisiert wird und eine monoton wachsende Folge von Approximationen xn+1 = xn −

xn − b f (xn ) f (xn ) − f (b)

(4.6)

bildet. Da die Folge nach oben durch α beschränkt ist, zeigt dies, dass xn → α bzw. f (xn ) → 0. Um die Konvergenzrate zu bestimmen, subrahieren wir α von beiden Seiten in (4.6) und nutzen f (α) = 0 aus: xn − b xn+1 − α = xn − α − [f (xn ) − f (α)] f (xn ) − f (b) Damit gilt xn+1 − α xn − b f (xn ) − f (α) =1− . xn − α f (xn ) − f (b) xn − α Für den Grenzwert n → ∞ erhalten wir xn+1 − α f 0 (α) = 1 − (b − α) =: c. n→∞ xn − α f (b) lim

Da 0 ≤ c ≤ 1 haben wir lineare Konvergenz des linearen Eingabelns gezeigt. Für konkave und konvexe Funktionen liefert dieses Verfahren also diesselbe Konvergenzrate wie die Bisektion. Allerdings kann für Funktionen, die im Intervall [an , bn ] nicht strikt konvex oder konkav sind superlineare Konvergenz gezeigt werden. Dies macht die Methode in einigen Fällen also effizienter als Bisektion.

4.3 Das Newton-Verfahren Ein weiteres Verfahren zur Bestimmung von Nullstellen einer Funktion ist das Newton-Verfahren (auch Newton-Raphson-Verfahren). Die grundlegende Idee ist hierbei in einem Ausgangspunkt die Tangente der Funktion zu bestimmen und die Nullstelle der Tangente als verbesserte Näherung der Nullstelle der Funktion zu verwenden. Die erhaltene Approximation kann dann als Ausgangspunkt für einen weiteren Iterationsschritt verwendet werden. Wir werden sehen, dass

39

diese Methode (lokal) schneller konvergiert als Bisektion oder lineares Eingabeln, allerdings ist auch die Glattheitsanforderung an f grösser. Wir werden in diesem Abschnitt deshalb immer annehmen, dass f ∈ C 1 ([a, b]). Formal lässt sich die Methode durch Linearisierung der Funktion f in einem Ausgangspunkt x0 herleiten. Dazu wird f in x0 in eine Taylorreihe entwickelt und nur Terme erster Ordnung berücksichtigt. Es gilt f (x) ≈ f (x0 ) + (x − x0 )f 0 (x0 ), wobei die Approximation nur in einer hinreichend kleinen Umgebung von x0 gut ist. Ersetzt man in der ursprünglich zu lösenden Gleichgung f (x) = 0, die Funktion f nun durch die obige Approximation führt dies auf die Gleichung f (x0 ) + (x − x0 )f 0 (x0 ) = 0. Wir bezeichnen die Lösung dieser Gleichung mit x1 und es gilt x1 = x0 −

f (x0 ) . f 0 (x0 )

x1 kann nun als Startwert für einen weiteren Iterationsschritt verwendet werden, sodass die allgemeine Iterationsvorschift des Newton-Verfahrens durch xn+1 = xn −

f (xn ) f 0 (xn )

(4.7)

gegeben ist. Ob das Newton-Verfahren gegen eine Nullstelle von f konvergiert, hängt stark von der Wahl des Startwertes x0 der Iteration ab. Wir werden im Folgenden einige Beispiele betrachten, die die möglichen Verhaltensweisen des Newton-Verfahrens veranschaulichen. Beispiel 4.1 Approximation von muss lautet

√

a, a > 0. Die zugehörige Gleichung, die gelöst werden f (x) = x2 − a.

Die Newton Iteration für diese Gleichung ist xn+1

x2 − a 1 = xn − n = 2xn 2



a xn + xn

 ,

n = 0, 1, . . ..

Die Iteration konvergiert für jeden positiven Startwert x0 . (Warum?) Beispiel 4.2 Wir betrachten die Gleichung sin(x) = 0 für |x|
1. Die Newton Iteration konvergiert in diesem Fall nur gegen die gesuchte Nullstelle, falls der Startwert x0 in einer hinreichend kleinen Umgebung um die Nullstelle α = 0 gewählt wird.

Abbildung 4.1: Graph der Funktion f (x) =

x x2 +1

Beispiel 4.4 Wir betrachten die Funktion f (x) = x20 −1, x > 0. Diese hat genau eine einfache Nullstelle α = 1. Die zugehörige Newton Iteration ist xn+1 = xn −

x20 19 1 n −1 = xn + , 19 20x19 20 20x n n

41

n = 0, 1, . . .

Die ist ein Beispiel dafür, dass die Newton Methode bei einer schlechten Wahl des Startwerts 19 (vorasymptotisch) sehr langsam konvergieren kann. Falls wir x0 = 12 wählen, gilt x1 = 220 ≈ 2.6 · 104 . Dies ist weit von der Nullstelle entfernt und es dauert lange bis die Iteration sich 1 wieder annähert. Um das zu sehen, bemerken wir, dass für sehr grosses xn xn+1 ≈

19 xn 20

gilt. Somit wird die Approximation in jedem Schritt nur um einen Faktor 0.95 reduziert. Falls sich die Approximation allerdings α = 1 annähert konvergiert das Newton-Verfahren sehr schnell (quadratisch), wie wir im Folgenden sehen werden. Um die (lokal) quadratische Konvergenz des Newton-Verfahrens zu zeigen, benötigen wir folgenden Satz. Satz 4.1 (Mittelwertsatz der Differentialrechnung) Sei [a, b] ein Intervall, f ∈ C([a, b]) sowie f ∈ C 1 ((a, b)). Dann existiert ein ξ ∈ (a, b), sodass f 0 (ξ) =

f (b) − f (a) b−a

gilt. Beweis. Wir definieren eine Hilfsfunktion h : [a, b] → R mit h(x) = f (x) −

f (b) − f (a) (x − a). b−a

h ist stetig in [a, b] und differenzierbar in (a, b). Ausserdem gilt h(a) = f (a) = h(b). Nach dem Satz von Rolle existiert nun eine Zwischenstelle ξ ∈ (a, b), sodass h0 (ξ) = 0. Daraus folgt die Behauptung.  Satz 4.2 Sei α eine einfache Nullstelle von f . Sei I = {x ∈ R : |x − α| ≤ } ein Intervall um α und sei weiterhin f ∈ C 2 [I ]. Für hinreichend kleines  definieren wir 00 f (s) M () = max 0 . (4.8) s,t∈I 2f (t) Falls  so klein ist, dass 2M () < 1

(4.9)

erfüllt ist, dann ist das Newton-Verfahren für jeden Startwert x0 ∈ I wohldefiniert und konvergiert quadratisch gegen die einzige Nullstelle α ∈ I . Beweis. Wir zeigen zunächst, dass α die einzige Nullstelle von f in I ist. Dazu verwenden wir eine Taylorentwicklung um α und erhalten f (x) = f (α) + (x − α)f 0 (α) +

(x − α)2 00 f (ξ), 2

mit einer Stelle ξ zwischen α und x. Da α eine Nullstelle von f ist, gilt f (α) = 0. Ausklammern ergibt deshalb   f 00 (ξ) 0 f (x) = (x − α)f (α) 1 + (x − α) 0 , (4.10) 2f (α)

42

wobei ξ ∈ I , falls x ∈ I . Falls x 6= α, ist der erste Faktor in (4.10) klar ungleich 0. Der zweite Faktor ist ebenfalls ungleich 0, wegen der Annahme, dass α eine einfache Nullstelle ist. Der dritte Faktor ist auch ungleich 0, da 00 00 (ξ) f (ξ) f (x − α) ≤  2f 0 (α) ≤ M () ≤ 2M () < 1 2f 0 (α) nach Annahme. Somit ist α die einzige Nullstelle von f im Intervall I . Als nächstes zeigen wir, dass alle Iterierten xn im Intervall I liegen und zwei aufeinanderfolgende Werte sich voneinander unterscheiden, falls xn 6= α. Dies wird duch Induktion bewiesen. Zunächst bemerken wir dabei, dass die Aussage für x1 richtig ist, da x1 = x0 −

f (x0 ) 6= x0 f 0 (x0 )

für x0 6= α. Sei xn ∈ I und xn+1 6= xn . Aufgrund von Satz 4.1 gibt es ein ξ1 ∈ [α, xn ] und ein ξ2 ∈ [ξ1 , xn ], sodass   f (xn ) 1 f (xn ) − f (α) xn+1 − α = xn − α − 0 = (xn − α) 1 − 0 f (xn ) f (xn ) (xn − α)   0   f (xn ) − f 0 (ξ1 ) f 0 (ξ1 ) = (xn − α) = (xn − α) 1 − 0 f (xn ) f 0 (xn )  0  f (xn ) − f 0 (ξ1 ) = (xn − α)(xn − ξ1 ) f 0 (xn )(xn − ξ1 ) f 00 (ξ2 ) = (xn − α)(xn − ξ1 )2 0 . 2f (xn ) Mit der Annahme aus dem Satz erhalten wir |xn+1 − α| ≤ 2 · 2 · M () ≤  und somit xn+1 ∈ I . Zuletzt muss die quadratische Konvergenzordnung gezeigt werden. Die folgt auch aus der obigen Rechnung. Da |xn − ξ1 | ≤ |xn − α| gilt: |xn+1 − α| ≤ |xn − α|2 · 2M () und somit

|xn+1 − α| ≤ 2M (). |xn − α|2

Daraus ergibt sich die quadratische Konvergenzordnung.



Beispiel 4.5 Wir betrachten die Funktion f (x) = cos(x) im Intervall [ π8 , 7π 8 ]. Um herauszufinden wann das Newton-Verfahren quadratisch gegen die einzige Nullstelle α = π2 konvergiert, wenden wir den obigen Satz an. Die Voraussetzungen von Satz 4.2 sind erfüllt. Somit müssen wir ein  finden, sodass (4.9) erfüllt ist. Es gilt 00 f (s) − cos(s) cos(s) = = 2f 0 (t) −2 sin(t) 2 sin(t)

43

Sei

und  ≤

π 7π π I = {x ∈ [ , ] : x − ≤ }. 8 8 2 3π 8 .

Dann gilt | cos(s)| ≤ cos

und |2 sin(t)| ≥ 2 sin Damit ergibt sich

π 2 π 2

 −  = sin() für s ∈ I  −  = 2 cos() für s ∈ I .

00 f (s) 1 sin() 1 max 0 ≤ = tan() s,t∈I 2f (t) 2 cos() 2

Für die quadratische Konvergenz des Newton-Verfahrens muss  also 1 2 tan() < 1 2

⇒

 < 0.8603

erfüllen.

4.4 Fixpunkt-Iterationen Nichtlineare Gleichungen sind oft als sogenannte Fixpunktgeichungen gegeben. Dabei besteht die Aufgabe darin ein x zu finden, sodass x = ϕ(x)

(4.11)

gilt. Eine Zahl x, welche die obige Gleichung erfüllt wird Fixpunkt von ϕ genannt. Sei beispielsweise eine nichtlineare Gleichung f (x) = 0 mit f 0 (x) 6= 0 gegeben. Dann kann diese Gleichung (unter anderem) auch in der Form (4.11) geschrieben werden mit ϕ(x) = x −

f (x) f 0 (x)

(4.12)

Sei α ein Fixpunkt von (4.11) und x0 eine anfängliche Approximation an α. Dann ist durch xn+1 = ϕ(xn ),

n = 0, 1, . . .

(4.13)

eine Folge von Approximationen gegeben. Eine solche Folge wird Fixpunktiteration genannt. Falls diese Iteration konvergiert (und ϕ stetig ist), konvergiert sie gegen α, den Fixpunkt von (4.11). Die Iterationsvorschrift (4.13) entspricht gerade Newton’s Methode zur Lösung der Gleichung f (x) = 0, falls ϕ wie in (4.12) definiert ist. Newton’s Methode kann deshalb als Fixpunktiteration aufgefasst werden. Für Fixpunktiterationen der Form 4.13 kann die Konvergenzordnung leicht bestimmt werden. Wir nehmen dafür an, dass die Iteration konvergiert und xn → α für n → ∞. Ausserdem sei die Iterationsvorschrift ϕ p-mal stetig differenzierbar und es gelte ϕ0 (α) = ϕ00 (α) = · · · = ϕ(p−1) (α) = 0,

44

ϕ(p) (α) 6= 0

(4.14)

für ein p ≥ 1. Eine Taylorentwicklung um den Fixpunkt α ergibt ϕ(xn ) = ϕ(α) + (xn − α)ϕ0 (α) + · · · + = ϕ(α) +

(xn − α)(p−1) (p−1) (xn − α)p (p) ϕ (α) + ϕ (ξn ) (p − 1)! p!

(xn − α)p (p) ϕ (ξn ), p!

wobei ξn zwischen α und xn liegt. Da ϕ(α) = α und xn+1 = ϕ(xn ) gilt 1 xn+1 − α = ϕ(p) (ξn ). (xn − α)p p! Da xn → α, gilt aufgrund der Stetigkeit von ϕ(p) lim

n→∞

xn+1 − α 1 = ϕ(p) (α) 6= 0. p (xn − α) p!

Somit ist die Konvergenz der Ordnung p gezeigt. Bisher wurde immer angenommen, dass die Fixpunktiteration konvergiert. Der folgenden Satz besagt, unter welchen Voraussetzungen an ϕ dies tatsächlich der Fall ist. Satz 4.3 Sei α ein Fixpunkt von ϕ und I = {x ∈ R : |x − α| ≤ }. Angenommen ϕ ∈ C p (I ) erfüllt die Bedingung (4.14) und M () := max |ϕ0 (t)| < 1. t∈I

Dann konvergiert die Fixpunktiteration (4.13) gegen α für jeden Startwert x0 ∈ I . Die Konvergenz ist von der Ordnung p. Beispiel 4.6 Es soll die Gleichung sin(0.5 · x) − x + 1 = 0 gelöst werden. Diese lässt sich in der Form (4.11) schreiben mit ϕ(x) = sin(0.5 · x) + 1. Die Bedingungen von Satz 4.3 sind erfüllt und die Iteration xn+1 = ϕ(xn ), n = 0, 1, . . . konvergiert gegen einen Fixpunkt und somit gegen die Lösung der Gleichung α ≈ 1.77572. Abbildung 4.2 zeigt die Konvergenz der Fixpunktiteration für den Startwert x0 = 0 und x0 = 40.

4.5 Die ∆2 -Methode von Aitken Die ∆2 -Methode von Aitken gehört zu den Methoden zur Konvergenzbeschleunigung einer gegebenen konvergenten Folge (xn )n mit xn → α für n → ∞. Sie basiert auf der Idee die gegebene Folge in einer andere Folge (˜ xn )n zu transformieren, die schneller gegen α konvergiert als die ursprüngliche Folge. Die Methode kann zur Nullstellenbestimmung von Funktionen verwendet werden, wenn (xn )n eine Folge aus einem anderen Verfahren zur Nullstellenbestimmung ist. Um die Methode herzuleiten, nehmen wir an, dass die ursprüngliche Folge (xn )n linear gegen α konvergiert, d.h. es gilt xn+1 − α = c(xn − α),

für n = 0, 1, . . . ,

für ein c mit |c| ≤ 1. Aus den Iterierten xn , xn+1 und xn+2 kann man mit Hilfe der Gleichungen xn+1 − α = c(xn − α),

xn+2 − α = c(xn+1 − α)

45

Abbildung 4.2: x-Achse: Anzahl Iterationen, y-Achse: Fehler |xn − α| der Fixpunktiteration für ϕ(x) = sin(0.5 · x) + 1.

c und α berechnen. Durch Subtraktion der Gleichungen findet man c= und daraus α=

xn+2 − xn+1 xn+1 − xn

xn xn+2 − x2n+1 . xn+2 − 2xn+1 + xn

(4.15)

Wir führen die Notationen ∆xn := xn+1 − xn und ∆2 xn := ∆xn+1 − ∆xn = xn+2 − 2xn+1 + xn ein und können (4.15) in der Form α = xn −

(∆xn )2 ∆2 xn

(4.16)

schreiben. Es liegt nun die Vermutung nahe, dass (4.16) einen guten Näherungswert für den gesuchten Grenzwert α der Folge (xn )n liefert, auch wenn die Voraussetzung, dass (xn )n eine linear konvergente Folge ist, nicht zutrifft. Die ∆2 -Methode von Aitken besteht somit darin, zu einer gegebenen Folge (xn )n die transformierte Folge (∆xn )2 (xn+1 − xn )2 x ˜ n = xn − = xn − (4.17) 2 ∆ xn xn+2 − 2xn+1 + xn zu berechnen. Der folgende Satz zeigt, dass die transformierte Folge schneller gegen α konvergiert, wenn die ursprüngliche Folge asymptotisch linear gegen α konvergiert. Satz 4.4 Sei |c| < 1 und es gelte für die Folge (xn )n , xn 6= α und xn+1 − α = (c + δn )(xn − α),

lim = 0.

n→∞

Dann existiert die Folge (4.17) für hinreichend grosses n und es gilt x ˜n − α = 0. n→∞ xn − α lim

46

Beweis. Für den Fehler en := xn − α gilt nach Voraussetzung en+1 = (c + δn )en . Es folgt xn+2 − 2xn+1 + xn = en+2 − 2en+1 + en = en ((c + δn+1 )(c + δn ) − 2(c + δn ) + 1) = en ((c − 1)2 + µn ),

wobei µn → 0.

xn+1 − xn = en+1 − en = en ((c − 1) − δn ). Die Folge x ˜n ist durch (4.17) wohldefiniert, da en 6= 0, c 6= 1 und µn → 0 und somit xn+2 − 2xn+1 + xn 6= 0 für grosse n. Weiterhin folgt für grosses n (xn+1 − xn )2 (∆xn )2 = xn − −α 2 ∆ xn xn+2 − 2xn+1 + xn ((c − 1) − δn )2 = en − en . (c − 1)2 + µn

x ˜ n − α = xn −

Somit gilt   ((c − 1) − δn )2 x ˜n − α = lim 1 − =0 lim n→∞ n→∞ xn − α (c − 1)2 + µn 

47