Skript zur Vorlesung:
”Einfuhrung in die numerische Mathematik “
Numerische Analysis
Prof. Dr. P.E. Kloeden
Institut fur Mathematik
Johann Wolfgang Goethe Universitat
Zimmer 101, Robert-Mayer-Straße 10
Telefon: (069) 798 28622 — Sekretariat (069) 798 22422
email: [email protected]
6. Oktober 2009
Inhaltsverzeichnis
1 Approximationsformeln 3
1.1 Funktionenberechnung . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Polynomberechnungen . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Nullstellen und Fixpunkte 9
2.1 Fixpunkte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Konvergenzordnung 21
3.1 Computerarithmetische Iterationen . . . . . . . . . . . . . . . . . 24
3.2 Nullstellen (nochmal) . . . . . . . . . . . . . . . . . . . . . . . . 25
4 Nullstellen zum letzten Mal 30
4.1 Vereinfachungen der Newton–Methode . . . . . . . . . . . . . . . 31
4.2 Nullstellen von Polynomen . . . . . . . . . . . . . . . . . . . . . . 35
4.3 Komplexwertige Nullstellen eines Polynoms . . . . . . . . . . . . 36
5 Interpolationspolynome 40
5.1 Lagrange’sche Interpolationspolynome . . . . . . . . . . . . . . . 41
5.2 Newton’sche Interpolationspolynome . . . . . . . . . . . . . . . . 44
6 Fehlerabschatzung 49
6.1 Differenzenoperatoren . . . . . . . . . . . . . . . . . . . . . . . . 52
6.2 Aquidistante Stutzstellen . . . . . . . . . . . . . . . . . . . . . . 53
7 Spline–Interpolation 56
7.1 Lineare Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.1.1 Fehlerabschatzung . . . . . . . . . . . . . . . . . . . . . . 58
7.2 Quadratische Splines . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.3 Kubische Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
7.3.1 Die Minimaleigenschaft kubischer Splines . . . . . . . . . 62
7.3.2 Fehlerabschatzung . . . . . . . . . . . . . . . . . . . . . . 62
1
2
8 Orthogonale Polynome 65
8.1 Ein bißchen Lineare Algebra . . . . . . . . . . . . . . . . . . . . . 65
8.2 Eigenschaften von Orthogonalpolynomen . . . . . . . . . . . . . . 69
8.3 Die”beste“ Polynomapproximation . . . . . . . . . . . . . . . . . 71
8.4 Die”besten“ Stutzstellen . . . . . . . . . . . . . . . . . . . . . . 73
9 Numerische Differentiation 76
9.1 Rundung in numerischer Differentiation . . . . . . . . . . . . . . 79
10 Numerische Integration 82
10.1 Interpolatorische Quadraturformeln . . . . . . . . . . . . . . . . . 85
10.2 Rundung in numerischer Integration . . . . . . . . . . . . . . . . 89
11 Summierte Integrationsformeln 91
11.1 Extrapolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
11.2 Romberg-Integration . . . . . . . . . . . . . . . . . . . . . . . . . 97
12 Gauß-Quadratur 102
Kapitel 1
Approximationsformeln
Wir konnen einfache Approximationsformeln fur numerische Berechnungen aus
den Grunddefinitionen/Ergebnisse der mathematischen Analysis herleiten, z.B.
Ableitung / Differenzenquotient
f ′(x) = limh→0
f(x + h)− f(x)
h≃ f(x + h)− f(x)
h,
Integral / Riemannsche Summe
∫ 1
0
f(x) dx = limN→∞
N∑
j=0
(xj+1 − xj)f(xj) ≃N∑
j=0
(xj+1 − xj)f(xj),
Taylor-Entwicklung / Polynom
ex =∞∑
j=0
1
j!xj ≃
N∑
j=0
1
j!xj .
Eine Berechnung mit einer solchen Formel ist nicht exakt, d.h. wir haben
einen Approximationsfehler, z.B.
∣∣∣∣e1 −
N∑
j=0
1
j!
∣∣∣∣ =
∞∑
j=N+1
1
j!
Wie groß sind solche Approximationsfehler?
Konnen wir”bessere“ Formeln finden?
=⇒ Grundlagen der Numerik!
ABER die Approximationsfehler hier sind nur”theoretische“ Fehler — in der
Praxis konnen wir die Zahlen selbst nur ungefahr berechnen: der Zahlenkorper
eines Computers ist nur endlich.
3
KAPITEL 1. APPROXIMATIONSFORMELN 4
1.1 Funktionenberechnung
Literatur Stummel/Hainer, Kap. 1
Sei f : R → R eine Funktion mit der Taylor-Entwicklung von x0 :
f(x) =
∞∑
j=0
1
j!f (j)(x0) · (x− x0)
j
z.B. ex =∑∞
j=01j!x
j und arctanx =∑∞
j=0(−1)j
2j+1 x2j+1
Wir konnen das Polynom
fN(x) =
N∑
j=0
1
j!f (j)(x0) (x− x0)
j
benutzen, um f(x) zu berechnen. Wie groß soll N sein, um eine erwunschte
Genauigkeit zu versichern ? d.h.
FehlerN =∣∣∣
∞∑
j=N+1
1
j!f (j)(x0) · (x− x0)
∣∣∣ ≤ ε.
Es gibt viele Tricks/Ergebnisse aus der Analysis, die wir benutzen konnen,
um diesen Fehler zu schatzen.
1. Beispiel das Leibniz-Kriterium fur eine alternierende Reihe S. Sei S =∑∞j=0(−1)jaj mit aj > 0 fur alle j:
⇒∣∣∣S −
N∑
j=0
(−1)jaj
∣∣∣ ≤ aN+1
z.B. fur arctan(12 ) =
∑∞j=0
(−1)j
2j+1 (12 )2j+1 gilt
∣∣∣ arctan
(1
2
)−
N∑
j=0
(−1)j
2j + 1
(1
2
)2j+1 ∣∣∣ ≤ 1
2N + 32−2N−3.
Wahle N , so dass
1
2N + 32−2N−3 ≤ ε aber
1
2N + 12−2N−1 > ε ⇒ |FehlerN | ≤ ε
2. Beispiel die Taylor-Entwicklung mit Restterm. Sei f : R → R mindestens
(N + 1)-mal stetig differenzierbar in einer Umgebung von x0.
f(x) =
N∑
j=0
1
j!f (j)(x0) · (x− x0)
j+
1
(N + 1)!f (N+1)(ξx0,x) ·
(x− xN+1
0
),
KAPITEL 1. APPROXIMATIONSFORMELN 5
wo ξx0,x ∈ (x0, x) meistens unbekannt ist.
z.B. fur e1 =∑N
j=01j! + 1
(N+1)!eξ, wo ξ ∈ (0, 1).
aber 0 < eξ ≤ e1 < 3, weil ex monoton ist.
⇒
∣∣∣∣∣∣e1 −
N∑
j=1
1
j!
∣∣∣∣∣∣<
3
(M + 1)!.
Sei ε = 10−3. Wir haben 36! = 0.004 . . . und 3
7! = 0.0005 . . .. Nimm N = 6:
⇒
∣∣∣∣∣∣e1 −
6∑
j=0
1
j!
∣∣∣∣∣∣< 10−3
Ein so gewahltes N ist oft viel zu groß — overkill!
1.2 Polynomberechnungen
Literatur Stummel/Hainer, Kap. 1
Trivial! ? Aber pass auf! Rundungsfehler konnen sich akkumulieren.
Betrachte das Polynom
pN(x) = aNxN + aN−1xN−1 + . . . + a1x + a0 .
Aufwand:
Additionen = N ,
Multiplikationen = N +∑N
j=1 j = 12N(N + 3).
Obiger Ansatz ist nicht nur aufwendig sondern auch numerisch nicht sinnvoll,
speziell fur kleine wie auch grosse x. Es geht besser
pN (x) =[. . .[ [
aNx + aN−1︸ ︷︷ ︸1. Schritt
]x + aN−2
︸ ︷︷ ︸2. Schritt
]x + . . . + a1
]x + a0
⇒ Additionen = N , Multiplikationen = N .
KAPITEL 1. APPROXIMATIONSFORMELN 6
Zusatzlich: keine Schwierigkeiten fur kleine oder grosse x
Algorithmus:
aN = aN
ak = ak + x0ak+1 for k = N − 1, . . . , 1, 0
⇒ pN (x0) = a0
Betrachte jetzt das neue Polynom
pN−1(x) := aNxN−1 + aN−1xN−2 + . . . + a2x + a1.
Dann haben wir
pN (x) = r0 + (x− x0) pN−1(x)
wobei r0 = a0 = pN (x0).
Insbesondere
pN (x)− pN (x0)
x− x0= pN−1(x)
Nimm den Limes als x→ x0 und benutze die Stetigkeit von pN−1
⇒ p′N (x0) = pN−1 (x0)
Aber wir konnen pN−1(x0) wie oben berechnen:
KAPITEL 1. APPROXIMATIONSFORMELN 7
pN−1(x0) = a′′1 = r1
mit pN−1(x) = r1 + (x− x0) pN−2(x), usw. Zusammengefasst:
Horners-Methode
Fur j = 0, 1, . . ., N definiere
pN−j(x) = a(j)N xN−j + . . . + a
(j)j+1x + a
(j)j
wobei
a(j+1)N = a
(j)N
a(j+1)k = a
(j)k + x0 a
(j+1)k+1 , k = N − 1, . . . , j.
⇒ rj =1
j!
dj
dxjpN (x0) = a
(j+1)j
und
p(j)N (x0) = j! rj
Alles wird verstandlich mit folgender Tabelle
a(0)N a
(0)N−1 . . . a
(0)2 a
(0)1 a
(0)0
↓ ↓ ↓ ↓ ↓
a(1)N → a
(1)N−1 . . . → a
(1)2 → a
(1)1 → a
(1)0
↓ ↓ ↓ ↓
a(2)N → a
(2)N−1 . . . → a
(2)2 a
(2)1
↓ ↓
... Berechnung: ↓ ·1
a(N)N
−→·x0
+
Horners-Methode ergibt pN (x0), p′N(x0), . . . , p
(N)N (x0)
KAPITEL 1. APPROXIMATIONSFORMELN 8
Beispiel Stummel/Hainer, Seite 14
p(x) = x4 = 1.x4 + 0.x3 + 0.x2 + 0x + 0.1
j = 0 1: 0 0 0 0
: ↓ ↓ ↓ ↓
j = 1 1:→ x0 → x20 → x3
0 → x40
: ↓ ↓ ↓ տ
j = 2 1:→ 2x0 → 3x20 → 4x3
0 p(x0)/0!
: ↓ ↓ տ
j = 3 1:→ 3x0 → 6x20 p′(x0)/1!
: ↓ տ
j = 4 1:→ 4x0 p′′(x0)/2!
տ
p′′′(x0)/3!
Bemerkung
pN−j(x) = rj + (x− x0) pN−j−1(x)
⇒ pN(x) =
N∑
j=0
rN−j(x− x0)N−j
⇒ p(j)N (x0) = j! rj .
Kapitel 2
Nullstellen und Fixpunkte
Literatur Oevel, Kap. 4; Schwarz, Kap. 5; Stummel/Hainer, Kap. 2
maple–Arbeitsblatt Polynome
Sei f : R → R stetig, x∗ ∈ R und gelte f(x∗) = 0 dann heißt x∗
Nullstelle von f
f(x∗) = 0
Wann/wo ?
Zwischenwertsatz aus der Analysis ergibt eine hinreichende Bedingung
f stetig auf einem Intervall [a, b] mit f(a)f(b) < 0
⇒ es existiert mindestens eine Nullstelle x∗ von f mit x∗ ∈ (a, b)
9
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 10
Wie konnen wir ein solches x∗ berechnen ?
Die einfachste Methode ist die Intervallhalbierungsmethode oder
Bisektionsmethode .
Sei f : [a, b] → R stetig mit f(a)f(b) < 0. Definiere c = a+b2 , als den
Mittelpunkt des Intervalls [a, b]
f(c) = 0 ⇒ x∗ = c Ende!
f(c) 6= 0
f(a)f(c) < 0 ⇒ nimm [a, c] statt [a, b]
f(c)f(b) < 0 ⇒ nimm [c, b] statt [a, b]
Definiere I0 = [a0, b0] = [a, b], dann wiederhole mit I1 = [a1, b1], dann mit
I2 = [a2, b2], usw,
a0
a1
b0
b1
a2 b2
x
f(x)
Abbildung 2.1: Bisektionsverfahren
Algorithmus
• Fur n = 0 setze I0 = [a0, b0] mit a0 = a und b0 = b und
fur n > 0 setze In = [an, bn] und definiere cn =bn + an
2
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 11
• Falls f(cn) = 0 dann setze x∗ = cn Ende!
• Falls f(an)f(cn) < 0 dann setze In+1 = [an, cn] und wiederhole
• Falls f(cn)f(bn) < 0 dann setze In+1 = [cn, bn] und wiederhole.
Bemerkungen
1). f(a)f(b) < 0 ⇒ ∃ x∗ ∈ In ∀ n ≥ 0
2). Lange (In) = 12 · Lange (In−1) = . . . = 1
2n Lange (I0) → 0 fur n → ∞
⇒ konvergiert immer !
3). a priori Abschatzung wieviele Iterationen?
|x∗ − cn| ≤ Lange (In) =1
2n(b − a) < ε
⇒ erwunschte Genauigkeit, wenn n ≥ log(
b−aε
)/ log 2 .
Beispiel Schwarz: Seite 246 + ..
f(x) = cosx · coshx + 1
Fur die Bisektionmethode haben wir (Schwarz Tab. 5.6).
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 12
n an bn cn f(cn)
0 1.8 1.9 1.85 > 0
1 1.85 1.9 1.875 > 0
2 1.875 1.9 1.8875 < 0
3 1.875 1.8875 1.88125 < 0
4 1.875 1.88125 1.878125 < 0
5 1.875 1.878125 1.8765625 < 0
......
......
...
21 1.875104048 1.875104096 1.875104072 < 0
22 1.875104048 1.875104072 1.875104060 > 0
23 1.875104060 1.875104072 1.875104066 > 0
24 1.875104066 1.875104072 1.875104069 ≈ x∗ > 0
Konvergiert sehr langsam
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 13
Regula falsi Methode
Betrachte ein Intervall In = [an, bn] mit f(an)f(bn) < 0
Statt dem Mittelpunkt von In, nimm den Schnittpunkt mit der x-Achse der
Gerade zwischen (an, f(an)) und (bn, f(bn)).
y − f(an)
x− an=
f(bn)− f(an)
bn − an
Schnittpunkt y = 0 ⇒ x = cn wobei
cn =anf(bn)− bnf(an)
f(bn)− f(an)
Naturlich: f(bn) 6= f(an)
Dann wahle In+1 = [an+1, bn+1] wie oben.
Geht es schneller? Wir haben nur Lange (In+1) < Lange (In)
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 14
Regula falsi Methode: hier ist 12 Lange(In) ≤ Lange(In+1) ≤ Lange(In)
Beispiel : Siehe das obige Beispiel von Schwarz mit f(x) = cosx ·cosh x+1
Fur die Regula falsi-Methode haben wir (Schwarz Tab. 5.7).
n an bn cn f(cn)
0 1.80 1.90 1.873697942 5.8127 · 10−3
1 1.873697942 1.90 1.875078665 1.0512 · 10−4
2 1.875078665 1.90 1.875103609 1.9021 · 10−6
3 1.875103609 1.90 1.875104061 3.19 · 10−8
4 1.875104061 1.90 1.875104068 2.9 · 10−9
2.1 Fixpunkte
Die obigen Algorithmen sind etwas kompliziert — geht es etwas einfacher ?
Bemerkung Definiere g(x) := x + µf(x) mit µ 6= 0. Dann gilt
f(x∗) = 0 ⇔ g(x∗) = x∗
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 15
Fixpunkte sind auch fur sich genommen sehr interessant und haben viele
Anwendungen in der Mathematik.
Wir konnen oft einen Fixpunkt durch sukzessive Iterationen berechnen.
xn+1 = g(xn), n = 0, 1, 2, . . .
dh, mit x1 = g(x0), x2 = g(x1), x3 = g(x2), . . ..
⇒ Rekursionsformel,
auch eine Differenzengleichung oder ein zeitdiskretes dynamisches System mit
einer Ruhelage x∗
Angenommen xn → x∗ dann gilt wegen der Stetigkeit von g
xn+1 = g(xn)
↓ ↓
x∗ = g(x∗)
d.h., der Limeswert muß ein Fixpunkt sein.
Beispiel Schwarz, Seite 233 (Beispiel 5.1)
Betrachte g(x) = e−x fur x ≥ 0 und definiere
f(x) = g(x)− x = e−x − x,
die eine stetige Funktion ist. Dann
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 16
f(0) = g(0)− 0 = 1 > 0, f(1) = g(1)− 1 = e−1 − 1 < 0
⇒ f(x∗) = 0 oder g(x∗) = x∗ fur mindestens ein x∗ in (0, 1)
Wahle x0 ∈ [0, 1]
Fur die Fixpunktiterationen haben wir (Schwarz Tab. 5.1).
n xn n xn n xn
0 0.55000000 10 0.56708394 20 0.56714309
1 0.57694981 11 0.56717695 21 0.56714340
2 0.56160877 12 0.56712420 22 0.56714323
3 0.57029086 13 0.56715412 23 0.56714332
4 0.56536097 14 0.56713715 24 0.56714327
konvergiert sehr langsam!
Geht es immer? NEIN !
Kein Fixpunkt
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 17
Zu steil
Eine hinreichende Bedingung ist der Banachsche Fixpunktsatz. Wir betrach-
ten die skalare Version hier!
Eine Abbildung g : [a, b] → R heißt Kontraktion oder kontrahierende Ab-
bildung, wenn eine Konstante K ∈ [0, 1) existiert, so dass fur alle x, y ∈ [a, b] gilt:
|g(x)− g(y)| ≤ K|x− y|
K heißt Kontraktionskonstante
Beispiel g(x) = −e−x, x ≥ a > 0
Mittelwertsatz ∃ ξx,y ∈ (x, y) (oder (y, x), wenn y < x) mit
g(x)− g(y) = g′(ξx,y)(x− y) ∀x, y ≥ a
In diesem Fall mit g′(x) = e−x gilt |g′(x)| ≤ e−a fur x ≥ a > 0
⇒ |g(x)− g(y)| ≤ e−a|x− y| ∀ x, y ≥ a
Fixpunktsatz (Siehe Oevel, Satz 4.2, S. 44–45)
Es sei g : [a, b] → [a, b] eine Kontraktion uber dem abgeschlossenen Intervall
[a, b] mit der Kontraktionskonstanten K < 1. Dann
1). existiert ein eindeutiger Fixpunkt g(x∗) = x∗ ∈ [a, b],
2). konvergiert die durch xn+1 = g(xn) definierte Folge fur jeden Startwert
x0 ∈ [a, b] gegen den Fixpunkt x∗
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 18
3). gelten die Abschatzungen
|xn − x∗| ≤ K
1−K|xn − xn−1|
︸ ︷︷ ︸a posteriori
≤ Kn
1−K|x1 − x0|
︸ ︷︷ ︸a priori
fur jede so konstruierte Folge.
Beweisskizze
1. Schritt Fur n ≥ 0 und j ≥ 1, wegen der Dreiecksungleichung, gilt
|xn+j − xn| ≤ |xn+j − xn+j−1|+ . . . + |xn+2 − xn+1|+ |xn+1 − xn|
aber
|xn+2 − xn+1| = |g(xn+1)− g(xn)| ≤ K |xn+1 − xn|
|xn+3 − xn+2| = |g(xn+2)− g(xn+1)| ≤ K2 |xn+1 − xn|
usw.
⇒ |xn+j − xn| ≤(Kj−1 + . . . + K + 1
)|xn+1 − xn|
=1−Kj
1−K|xn+1 − xn|
≤ 1
1−K|xn+1 − xn|
Aber
|xn+1 − xn| = |g(xn)− g(xn−1)|
= K |xn − xn−1| ≤ . . . ≤ Kn |x1 − x0|
deshalb gilt:
|xn+j − xn| ≤Kn
1−K|x1 − x0|
⇒ {xn} ist Cauchy-Folge: ∀ ε > 0 ∃ N(ε) ≥ 1 mit
|xn+j − xn| < ε ∀j ≥ 1, n ≥ N(ε).
⇒ xn → x∗ wegen der Vollstandigkeit des metrischen Raums (R, | · |)
xn+1 = g(xn) ⇒ g(x∗) = x∗ Existenz!
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 19
2. Schritt (Eindeutigkeit)
Sei x∗ = g(x∗) 6= x∗∗ = g(x∗∗). Dann gilt
|x∗ − x∗∗| = |g(x∗)− g(x∗∗)| ≤ K |x ∗ −x∗∗| < |x∗ − x∗∗|
Widerspruch!
3. Schritt (Abschatzungen)
|xn+j − xn| ≤ K1−K |xn−1 − xn| ≤ Kn
I−K |x1 − x0| wie oben.
Und |xn+j − xn| → |x∗ − xn| als j →∞ ∀n
weil limj→∞
xn+j = limj→∞
xj = x∗.
Bemerkungen
1). Wir konnen die a priori Abschatzung
|xn − x∗| ≤ Kn
1−K|x1 − x0|
benutzen, wieviele Iterationen gebraucht werden (fur eine erwunschte Genauig-
keit) um im Voraus zu schatzen;
n ≥ N(ε) ⇒ |xn − x∗| ≤ . . . ≤ ε
N(ε) ist oft zu groß hier.
2). Die a posteriori Abschatzung ist oft besser oder genauer und gibt uns
einen Stoppbefehl mechanismus:
|xn − x∗| ≤ K
1−K|xn−1 − xn|
Sei ε > 0 die gewunschte Genauigkeit und sei N(ε) das erste n mit
|xn−1 − xn| ≤1−K
Kε
Dann gilt: |xn − x∗| ≤ ε fur n ≥ N(ε).
Beispiel g(x) = e−x fur x ≥ 0
KAPITEL 2. NULLSTELLEN UND FIXPUNKTE 20
g(x) ist Kontraktion uber x ≥ a, wo a > 0:
|g(x)− g(y)| ≤ e−a|x− y| ∀x, y ≥ a > 0 .
Wir mussen ein abgeschlossenes Intervall [a, b] mit a > 0 finden, fur welches
g([a, b]) ⊆ [a, b]
d.h. g bildet [a, b] in sich ab!
Bemerkung g(x) ist monoton abfallend
a ≤ x ≤ b ⇒ g(b) ≤ g(x) ≤ g(a)
Deshalb brauchen wir ein b > a mit
a ≤ g(b) ≤ g(a) ≤ b
Probiere a = 1/4 mit g(a) = e−1/4 ≃ 0.78: Dann fur b = g(a) gilt g(b) ≃ 0.46.
Alles ist OK hier
⇒ g ist Kontraktion und bildet [1/4, e−1/4] in sich ab.
⇒ Fixpunktsatz gilt hier! ∃ ! Fixpunkt x∗ ∈ [1/4, e−1/4].
Kapitel 3
Konvergenzordnung
Betrachte eine Differenzengleichung erster Ordnung oder 1-Schritt–Iterations-
verfahren
xn+1 = g(xn), n = 0, 1, 2, . . . ,
wobei g : [a, b] → [a, b] eine Kontraktion mit Kontraktionskonstante K < 1 ist.
Setze voraus, dass xn → x∗ fur n → ∞. Dann — g ist stetig! — gilt
x∗ = limn→∞
xn+1 = limn→∞
g(xn) = g(x∗)
d.h. x∗ = g(x∗) Fixpunkt von g.
Wie schnell konvergiert xn gegen x∗ ?
Schreibe en fur den Fehler der n-ten Iteration, d.h.
en = xn − x∗ .
Dann ist en+1 = xn+1 − x∗ = g(xn)− g(x∗) und deshalb gilt
|en+1| = |g(xn)− g(x∗)| ≤ K |xn − x∗| = K |en|
Definiere En = |en|. Dann haben wir
En+1 ≤ KEn n = 0,1,2,. . . , Differenzen-UN-gleichung
und deshalb gilt
limn→∞
En+1
En≤ K (falls En 6= 0 ∀n)
oder 0 ≤ En ≤ KnE0 → 0 fur n → ∞
21
KAPITEL 3. KONVERGENZORDNUNG 22
d.h. En → 0 als n→∞ mindestens linear oder mit erster Ordnung.
En+1
En
Kann es schneller sein?
Nicht immer: betrachte g(x) = 12x + 1
2 mit Fixpunkt x∗ = g(x∗) = 1
Hier gilt En+1 ≡ 12En fur n = 0, 1, 2 und deshalb konvergiert En = 1
2n E0
→ 0 genau linear
kann es besser sein, wenn g nicht linear ist?
Leider, auch nicht immer!
Sei g stetig differenzierbar mit g′(x∗) 6= 0. Von der Taylor-Entwicklung haben
wir
En+1 = |g(xn)− g(x∗)| = |g′(ξn)| |xn − x∗| = |g′(ξn)| En
Daher gilt
limn→∞
En+1
En= lim
n→∞|g′(ξn)| = |g′(x∗)| 6= 0,
weil ξn →x∗ und g′ stetig ist. (Tatsachlich gilt 0 < |g′(x∗)| < 1 hier)
Man sagt: die Konvergenz ist asymptotisch linear!
Im Allgemeinen konnen wir nur behaupten, dass eine Folge sukzessiver Itera-
tionen mindestens linear konvergiert, und in Sonderfallen nur linear konvergiert!
Aber setze jetzt voraus, dass g p-mal stetig differenzierbar ist mit
g′(x∗) = · · · = g(p−1)(x∗) = 0
und
g(p)(x∗) 6= 0.
KAPITEL 3. KONVERGENZORDNUNG 23
Dann ergibt die Taylor-Entwicklung:
g(xn)− g(x∗) =1
p!g(p)(ξn) (xn − x∗)p fur ξn ∈ x∗xn.
Davon haben wir
En+1 = |xn+1 − x∗| = |g(xn)− g(x∗)|
=1
p!
∣∣∣g(p)(ξn)∣∣∣ |xn − x∗|p
=1
p!
∣∣∣g(p)(ξn)∣∣∣ Ep
n
und deshalb
limn→∞
En+1
Epn
= limn→∞
1
p!
∣∣∣g(p)(εn)∣∣∣ =
1
p!
∣∣∣g(p)(x∗)∣∣∣ 6= 0
oder
En+1 ≃1
p!
∣∣∣g(p)(x∗)∣∣∣
︸ ︷︷ ︸Kp
Epn
In diesem Fall sagen wir, dass die Konvergenz p ter–Ordnung hat
p = 1 ⇒ linear, p > 1 ⇒ superlinear.
Konvergenzordnungen: linear, superlinear und quadratisch (von oben nach
unten)
Anschaulicher ist aber die logarithmische Darstellung
KAPITEL 3. KONVERGENZORDNUNG 24
Konvergenzordnungen: linear, superlinear und quadratisch (von oben nach
unten)
Im Allgemeinen: xn → x∗ mit Ordnung p > 0, wenn eine Konstante Kp ∈(0,∞) existiert, so dass
limn→∞
En+1
Epn
= Kp
mit
0 < q < p ⇒ limn→∞
En+1
Eqn
= 0, q > p ⇒ limn→∞
En+1
Eqn
=∞
Wir haben En+1 ≃ KpEpn, falls En ≪ 1
Bemerkungen
1). Die Konvergenzordnung p muss nicht immer eine Ganzzahl sein — spater
sehen wir, dass p ∼= 1.618 fur die Sekantenmethode.
2). Die obige Konvergenz ist”nur“ theoretisch, d.h. ohne Rundungsfehler oder
andere Berechnungsfehler.
3.1 Computerarithmetische Iterationen
Betrachte ein Iterationsverfahren
xn+1 = g(xn),
wobei g eine Kontraktion mit Konstante K < 1 ist. Durch Rundung (und an-
derer Art Approximationen, z.B. von Funktionenauswertung) haben wir
x0 = x0 + r0, . . . , xn+1 = g(xn) + rn+1, . . .
statt x0, . . ., xn+1 = g(xn), . . ., wobei
|rn| ≤ ε
KAPITEL 3. KONVERGENZORDNUNG 25
fur n = 0, 1, 2, . . .
Deshalb mussen wir En = |xn − x∗| statt En = |xn − x∗| abschatzen. Mit
der Dreiecksungleichung erhalten wir
En+1 = |xn+1 − x∗| = |g(xn)− x∗ + rn+1|
≤ |g(xn)− g(x∗)|+ |rn+1| ≤ K |xn − x∗|+ ε
d.h., En+1 ≤ KEn + ε statt En+1 ≤ K En.
Durch Induktion folgt, dass
En ≤ KnE0 +(1 + K + . . . + Kn−1
)ε
Aber E0 = |x0 − x∗| = |x0 + r0 − x∗| ≤ |x0 − x∗| + ε
⇒ En ≤ Kn |x0 − x∗|+ (1 + K + . . . + Kn) ε
d.h.,
En ≤ Kn |x0 − x∗|+ 1−Kn+1
1−Kε
Aus dieser Ungleichung konnen wir nicht schliessen, dass En = |xn − x∗| → 0
fur n→∞.
Aber fur n groß genug, z.B. fur
n ≥ n(ε) = log
(ε
|x0 − x∗|
)/ log K
haben wir Kn |x0 − x∗| ≤ ε und deshalb
En ≤ ε +1−Kn+1
1−Kε ≤ ε +
1
1−Kε =
2−K
1−K· ε sehr grob!
d.h. En = |xn − xn| ≤ 2−K1−K ε fur n ≫ 1
oder |xn − x∗| ≪ 1 fur n≫ 1, falls ε≪ 1
⇒ Rundung hat hier keine bosartige Wirkung. Auf einem Bildschirm
wurde es aussehen, als ob die Folge konvergiert!
3.2 Nullstellen (nochmal)
Sei g(x) = x + µf(x) mit µ = Konstante 6= 0.
KAPITEL 3. KONVERGENZORDNUNG 26
f(x∗) = 0 ⇔ g(x∗) = x∗
Nullstelle Fixpunkt
Hier gilt g′(x) = 1 + µf ′(x). Fur f ′(x∗) 6= 0 konnen wir µ 6= 0 wahlen, so
dass
|g′(x∗)| = |µ| ·∣∣∣∣1
µ+ f ′(x∗)
∣∣∣∣ < 1 .
Daher wird g eine Kontraktion uber einer Umgebung von x∗ sein.
Betrachte eine Folge sukzessiver Iterationen
xn+1 = g(xn) = xn + µf(xn), n = 0, 1, 2, . . .
Manchmal ist die Konvergenz (falls xn → x∗ = g(x∗)) nicht besser als linear.
Vielleicht konnen wir µ = µn sukzessiv verandern, um die Konvergenz zu
beschleunigen. Aber wie ?
Newton hat µn = −1/f ′(xn) gewahlt.
⇒ xn+1 = xn −f(xn)
f ′(xn)
Diese Vorschrift hat eine geometrische Bedeutung!
xx xii+1
f(x)
Newton Methode
KAPITEL 3. KONVERGENZORDNUNG 27
xn+1 ist der Schnittpunkt mit der x-Achse der Tangentengerade
y = f(xn) + (x − xn)f ′(xn)
d.h., y = 0 fur x = xn+1. Wir brauchen f ′(xn) 6= 0!
⇒ Newton-Methode
Bemerkungen
1). xn+1 = g(xn) mit g(x) = x− f(x)/f ′(x)
2). f stetig differenzierbar und f ′(x∗) 6= 0 ⇒ f ′(x) 6= 0 auf einer Umge-
bung von x∗
⇒ g stetig auf dieser Umgebung
Beispiel f(x) = cosx coshx + 1
k x(k) f(x(k)) f ′(x(k)) f(x(k))/f ′(x(k))
0 1.800000000 2.939755852 · 10−1 −3.694673552 7.95674946 · 10−2
1 1.879567405 −1.8530467 · 10−2 −4.165295668 −4.44877590 · 10−3
2 1.875118629 −6.0253 · 10−5 −4.138222447 −1.4560116 · 10−5
3 1.875104069 −1.0 · 10−9 −4.138133987 −2.4165 · 10−10
4 1.875104069
Die Tabelle ist aus Schwarz entnommen.
Satz Sei f 3-mal stetig differenzierbar mit f(x∗) = 0 und f ′(x∗) 6= 0. Dann
besitzt die Newton–Methode die Konvergenzordnung mindestens p = 2.
Beweis xn+1 = g(xn)
a). g(x) = x− f(x)
f ′(x)⇒ g(x∗) = x∗
b). g′(x) =f(x)f ′′(x)
[f ′(x)]2⇒ g′(x∗) = 0
⇒ mindestens p > 1 (superlinear!)
c). g′′(x) =[f ′(x)f ′′(x) + f(x)f ′′′(x)]f ′(x)− 2f(x)[f ′′(x)]2
[f ′(x)]3existiert und ist
stetig mit g′′(x∗) =f ′′(x∗)
f ′(x∗)
KAPITEL 3. KONVERGENZORDNUNG 28
⇒ mindestens p ≥ 2.
Im Allgemeinen gilt: g′′(x∗) 6= 0. Daher konnen wir nicht mehr als Ordnung
zwei erwarten.
Konvergiert die Newton-Methode immer? NEIN!
Hier ist |x0 − x∗| zu groß. Aber
g′(x) =f(x)f ′′(x)
[f ′(x)]2
ist stetig mit g′(x∗) = 0 ⇒ ∃ RK > 0 (vielleicht sehr klein), so dass
|g′(x)| ≤ K < 1 fur |x− x∗| ≤ RK
⇒ g ist Kontraktion auf [x∗−Rk, x∗ +RK ] mit Kontraktionskonstante K,
d.h.,
|g(x)− g(y)| ≤ K|x− y| ∀ x, y ∈ [x∗ −RK , x∗ + RK ]
Wir haben auch
|g(x)− x∗| = |g(x)− g(x∗)| ≤ K|x− x∗|
≤ KRK , falls x ∈ [x∗ −RK , x∗ + RK ],
< RK , weil K < 1 ist.
KAPITEL 3. KONVERGENZORDNUNG 29
⇒ g bildet |x∗ −RK , x∗ + RK ] in sich ab.
Daher konnen wir den Fixpunktsatz verwenden: die Newton-Methode kon-
vergiert, wenn |x0 − x∗| klein genug ist.
Praktische Probleme: Wie klein? Wo ist x∗?
Kapitel 4
Nullstellen zum letzten Mal
maple–Arbeitsblatt Beispiel zu Verfahren von Mueller
Wir wollen x∗ berechnen, wobei f(x∗) = 0.
Wir werden verschiedene Variationen und Verallgemeinerungen der Newton–
Methode untersuchen
Die Newton–Methode lautet
xn+1 = g(xn) mit g(x) = x− f(x)
f ′(x)
• einfach, 1-Schritt, mindestens quadratische Konvergenzordnung (falls
f ′(x∗) 6= 0)
• aber f ′(x) kann oft kompliziert sein — aufwendig herzuleiten oder zu
berechnen — oder unbekannt, z.B. durch Daten interpoliert. Verschiedene Me-
thoden versuchen die Ableitung zu vermeiden oder nicht so oft zu berechnen.
Beispiel : Schwarz Beispiel 5.5 mit f(x) = cosx · coshx + 1
Fur die Newton-Methode haben wir (Schwarz Tab. 5.9).
n xn f(xn) f ′(xn) f(xn)/f ′(xn)
0 1.800000000 2.939755852 · 10−1 −3.694673552 7.95674046 · 10−2
1 1.879567405 −1.8530467 · 10−2 −4.165295668 −4.44877590 · 10−3
2 1.875118629 −6.0253 · 10−5 −4.138222447 −1.4560116 · 10−5
3 1.875104069 −1.0 · 10−9 −4.138133987 −2.4165 · 10−10
4 1.875104069 ≈ x∗
30
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 31
4.1 Vereinfachungen der Newton–Methode
Die vereinfachte Newton-Methode lautet
xn+1 = xn −f(xn)
f ′(x0)
Hier ist g(x) = x− f(x)
f ′(x0)mit g(x∗) = x∗. Im Allgemeinen haben wir
g′(x∗) = 1− f ′(x∗)
f ′(x0)6= 0
Vereinfachte Newton-Methode
alle Geraden haben die Neigung = f ′(x0)
Beispiel Schwarz Beispiel 5.5 mit f(x) = cosx cosh x + 1
Fur die vereinfachte Newton-Methode haben wir (Schwarz Tab. 5.10).
n xn f(xn) f ′(x0) f(xn)/f ′(x0)
0 1.880000000 −2.0332923 · 10−2 −4.167932940 −4.878419 · 10−3
1 1.875121581 −7.2469 · 10−5 −1.738728 · 10−5
2 1.875104194 −5.19 · 10−7 −1.245222 · 10−7
3 1.875104069 −1.0 · 10−9 −2.399 · 10−10
4 1.875104069 ≈ x∗
• einfach und ziemlich schnell fur ein”gutes“ x0, aber ohne die quadrati-
sche Konvergenz der Newton-Methode.
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 32
Die Sekantenmethode vermeidet f ′(x) und die Tangentengerade — benutzt
eine Gerade durch 2 Punkte (x0, f(x0)) und (x1, f(x1)), d.h.,
y − f(x0)
x− x0=
f(x1)− f(x0)
x1 − x0.
Man berechnet die nachste Iteration durch
y = 0 ⇒ x = x2 =x0f(x1)− x1f(x0)
f(x1)− f(x0).
xi x
f(x)
xi+1xi-1
Sekantenmethode
Wiederhole mit x1 und x2 . . . mit xn−1 und xn, um x3, . . ., xn+1 zu berech-
nen.
xn+1 =xn−1f(xn)− xnf(xn−1)
f(xn)− f(xn−1)
• 2-Schritt-Iterationsverfahren, aber braucht nur eine neue Funktionsberech-
nung pro Schritt.
xn+1 = G(xn, xn−1)
mit
G(x, y) =
yf(x)− xf(y)
f(x)− f(y)y 6= x und f(x) 6= f(y)
x− f(x)
f ′(x)y = x
⇒ G stetig (L’Hospital-Regel!)
Beispiel Schwarz Beispiel 5.5 mit f(x) = cosx cosh x + 1
Fur die Sekantenmethode haben wir (Schwarz Tab. 5.8).
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 33
n xn f(xn)
0 1.80 2.39755852 · 10−1
1 1.90 −1.049169460 · 10−1
2 1.873697942 5.8127364 · 10−3
3 1.875078664 1.051126 · 10−4
4 1.875104095 −1.09 · 10−7
5 1.875104069 ≈ x∗ −1.0 · 10−9
Bemerkung 1.) Wir konnen dieses 2-Schrittverfahren in R1 als ein 1-Schrittverfahren
in R2 umschreiben.
Definiere x =
(x
y
)und g(x) =
(y
G(y, x)
)
Dann gilt: xn+1 = g(xn) fur n = 0, 1, . . . mit x0 =
(x0
x1
), d.h.,
(xn
xn+1
)=
(xn
G(xn, xn−1)
)
Wir konnen einen Fixpunktsatz in R2 verwenden, um Konvergenz zu un-
tersuchen oder zu versichern. Aber ohne Beschrankungen auf x0 und x1 ist die
Konvergenz nicht versichert. Solche Beschrankungen sind im Voraus nicht of-
fensichtlich.
Bemerkung 2.) Die Regula-falsi-Methode sieht sehr ahnlich aus, aber ist ein
bisschen komplizierter
Betrachte In = [an, bn] mit f(an)f(bn) < 0 fur n = 0, 1, . . . und berechne
cn =anf(bn)− bnf(an)
f(bn)− f(an)
Dann wahle
In+1 = [an+1, bn+1] =
[an, cn], falls f(an)f(cn) < 0
[cn, bn], falls f(cn)f(bn) < 0
Konvergenz cn → x∗ ist versichert — aber die Konvergenzordnung ist nur p =
1, d.h., linear (falls f ′(x∗) 6= 0, usw.).
Im Vergleich konvergiert die Sekantenmethode superlinear mit Konvergenz-
ordnung
p =1 +√
5
2≃ 1.618
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 34
Dies ist die positive Losung von
p2 − p− 1 = 0 goldener Schnitt!
Beweis Schreibe en = xn − x∗, wo f(x∗) = 0. Dann gilt
en+1 =en−1f(xn)− enf(xn−1)
f(xn)− f(xn−1)(∗)
Betrachte die Taylor-Entwicklung
f(xj) = f(x∗) + f ′(x∗)(xj − x∗) +1
2f ′′(x∗)(xj − x∗)2 + · · ·
≃ f ′(x∗) ej +1
2f ′′(x∗) e2
j + · · ·
Im obigen Zahler:
en−1f(xn)− enf(xn−1) ≃ en−1
[f ′(x∗) en +
1
2f ′′(x∗) e2
n
]
− en
[f ′(x∗) en−1 +
1
2f ′′(x∗) e2
n−1
]
=1
2f ′′(x∗) enen−1 (en − en−1)
Im obigen Nenner:
f(xn)− f(xn−1) ≃ f ′(x∗) · (en − en−1).
Deshalb erhalten wir aus (∗) die Approximation en+1 ≃ 12
f ′′(x∗)f ′(x∗) enen−1 oder
En+1 ≃ K EnEn−1,
wobei En = |en| und K = 12
∣∣∣ f′′(x∗)
f ′(x∗)
∣∣∣.Fur eine Konvergenzordnung p > 0 brauchen wir En+1 ≃ C Ep
n, wo die Kon-
stante C ∈ (0,∞).
Betrachte die Gleichungen
(a) En+1 = K EnEn−1
und
(b) En+1 = C Epn, (c) En = C Ep
n−1
Benutze (b) in (a) ⇒ C Epn = K EnEn−1 oder C Ep−1
n = K En−1.
Schreibe die letzte Gleichung hoch p und benutze (c) an der rechten Seite,
d.h.,
Cp E(p−1)pn = Kp Ep
n−1 = Kp C−1 En
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 35
⇒ E(p−1)p−1n = Kp C−1−p
d.h. Ep2−p−1n ≡ Konstante ∀n
daher muß p2 − p− 1 = 0 mit p > 0.
In der Steffensen-Methode ersetzen wir die Ableitung in der Newton-Methode
durch
f ′(x) ≃ f(x + f(x))− f(x)
f(x)Taylor-Entwicklung!
Wir erhalten das 1-Schrittverfahren
xn+1 = xn −f(xn)2
f(xn + f(xn))− f(xn)
Hier ist
g(x) = x− f(x)2
f(x + f(x))− f(x)fur x 6= x∗
mit
g(x∗) = x∗ − 0
0!!!!
Aber g(x)→ x∗ fur x→ x∗ (L’Hospital) ⇒ definiere g(x∗) = x∗
Diese Methode hat eine superlineare, aber nicht quadratische Konvergenz-
ordnung p ∈ (1, 2). Obwohl sie ein”einfaches“ 1-Schrittverfahren ist, braucht
sie 2 Funktionsberechnungen pro Schritt — etwas aufwendig!
4.2 Nullstellen von Polynomen
Betrachte ein Polynom
p(x) = aNxN + aN−1xN−1 + . . . + a1x + a0
mit aN 6= 0 und zudem seien alle Koeffizienten reellwertig. Dann haben wir
p(x) = aNxN
[1 +
aN−1
aNx+ . . . +
a1
aNxN−1+
a0
aNxN
]
≃ aNxN fur |x| ≫ 0
Deshalb genugen die Nullstellen von p (reellwertige und komplexwertige) der
Ungleichung
|z| ≤ R
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 36
fur ein geeignetes R <∞.
Benutze die Newton-Methode
xn+1 = xn −p(xn)
p′(xn)
mit (in jedem Schritt) der Horner-Methode, um p(xn) und p′(xn) zu berechnen.
⇒ quadratische Konvergenz fur eine einfache Nullstelle, aber superlineare
Konvergenz fur eine vielfache Nullstelle.
Sei x∗1 eine “numerisch berechnete“ Nullstelle, d.h.,
p(x∗1) ≃ 0.
Um eine neue Nullstelle von p zu berechnen — besser mit einem neuen
geeigneten Startwert x0 anzufangen als die Deflation-Methode zu benutzen, d.h.,
1. Definiere q(x) = p(x)/(x− x∗)
2. Berechne eine Nullstelle von q
Die Deflation-Methode kann schief gehen wegen Rundungsfehler. (siehe Schwarz,
Kap. 6.1 fur ein Beispiel).
4.3 Komplexwertige Nullstellen eines Polynoms
Wie konnen wir eine komplexwertige Nullstelle eines Polynoms berechnen ?
Fur die Newton Methode gilt
x0 ∈ R ⇒ xn ∈ R ∀n.
Vielleicht konnen wir x0 ∈ C \R wahlen — OK, aber wie? (es ist manchmal
schwierig f ′(x) fur x ∈ C zu berechnen!)
Eine andere Moglichkeit ist die Methode von Muller.
Die Sekantenmethode benutzt eine Gerade durch 2 sukzessive Stellen, um
die nachste Iterationsstelle zu finden. Die Muller-Methode benutzt eine Parabel
durch 3 sukzessive Stellen.
(x0, f(x0))
(x1, f(x1))
(x2, f(x2))
bekannt
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 37
Abbildung 4.1: Mueller-Methode
Betrachte p(x) = A(x− x2)2 + B(x− x2) + C durch die 3 obigen Stellen
f(x0) = A(x0 − x2)2 + B(x0 − x2) + C
f(x1) = A(x1 − x2)2 + B(x1 − x2) + C
f(x2) = C
ergibt ein lineares System
(x0 − x2)2 (x0 − x2) 1
(x1 − x2)2 (x1 − x2) 1
0 0 1
A
B
C
=
f(x0)
f(x1)
f(x2)
Die Matrix hier hat die Determinante
det = (x0 − x1)(x0 − x2)(x1 − x2)
6= 0 falls x0, x1, x2 paarweise verschieden sind
⇒ eindeutig losbar fur A, B, C. Wahle x3 mit p(x3) = 0, d.h.
x3 = x2 +
−C/B falls A = 0, B 6= 0√−C/A falls A 6= 0, B = 0
−2sgn(B) · C|B|+
√B2 − 4AC
sonst
d.h., mit (x0, x1, x2) berechnen wir x3
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 38
Wiederhole, also mit (x1, x2, x3) berechne x4 . . . und mit (xn−2, xn−1, xn)
berechne xn+1
⇒ 3-Schritt-Iterationsverfahren
Fehleranalyse ⇒ En+1 ≃∣∣∣∣f (3)(x∗)
6f ′(x∗)
∣∣∣∣ EnEn−1En−2
Ansatz En+1 = CEpn ⇒ p3 − p2 − p− 1 = 0
⇒ superlineare Konvergenz p ≃ 1.839
Wichtig xn+1 kann komplexwertig sein
⇒ die nachsten A, B, C sind dann komplexwertig usw.
Wir konnen komplexwertige Nullstellen mit der Muller-Methode berechnen,
auch mit reellwertigen Anfangswerten x0, x1, x2. (Die Programmierung ist et-
was kompliziert - wie schreibt man sgn(B) fur ein komplexwertiges B ?)
Beispiel (Aus Burden & Faires) Betrachte das Polynom 4tes Grades
p(x) = 16x4 − 40x3 + 5x2 + 20x + 6
besitzt 2 reellwertigen Nullstellen und ein Paar komplexkonjugierter Nullstellen.
Wir konnen alle Nullstellen mit der Muller-Methode berechnen
KAPITEL 4. NULLSTELLEN ZUM LETZTEN MAL 39
1. Fall Anfangswerte x0 = 0.5, x1 = 1.0 und x2 = 1.5
n xn f(xn)
3 1.28785 −1.37624
4 1.23746 0.126941
5 1.24160 0.219440 · 10−2
6 1.24168 0.257492 · 10−4
7 1.24168 0.257492 · 10−4
2. Fall Anfangswerte x0 = 2.5, x1 = 2.0 und x2 = 2.25
n xn f(xn)
3 1.96059 −0.611255
4 1.97056 0.748825 · 10−2
5 1.97044 −0.295639 · 10−4
6 1.97044 −0.295639 · 10−4
3. Fall Anfangswerte x0 = 0.5, x1 = −0.5 und x2 = 0
hier benutzen wir ı =√−1
n xn f(xn)
3 −0.555556 + 0.598352ı −29.4007− 3.89872ı
4 −0.435450 + 0.102101ı 1.33223− 1.19309ı
5 −0.390631 + 0.141852ı 0.375057− 0.670164ı
6 −0.357699 + 0.169926ı −0.146746− 0.00744629ı
7 −0.356051 + 0.162856ı −0.183868 · 10−2 + 0.539780 · 10−3ı
8 −0.356062 + 0.162758ı 0.286102 · 10−5 + 0.953674 · 10−6ı
Siehe auch maple–Arbeitsblatt”Beispiel zum Verfahren von Mueller“
Kapitel 5
Interpolationspolynome
Literatur Oevel: Kap. 8.1; Schwarz, Kap. 3.1–3.2; Stummel/Hainer, Kap. 3.1
Betrachte
• N + 1 paarweise verschiedene Zahlen x0, x1, . . . , xN in R (Stutzstellen)
• N + 1 Zahlen y0, y1, . . ., yN in R (Funktionswerte)
• gesucht sei eine Funktion f : R → R mit
yj = f(xj), j = 0, 1, . . . , N.
ENTWEDER ist f unbekannt, z.B. die y0, y1, . . . , yN sind Versuchsdaten
und wir wollen eine geeignete Funktion f mit f(xj) = yj finden, um fur x 6= xj
zu interpolieren, d.h. wir benutzen f(x) statt dem unbekannten f(x);
ODER f ist bekannt, aber zu kompliziert und wir wollen eine einfachere
Funktion f mit f(xj) = yj finden, die gunstiger zu benutzen ist, z.B. wir haben
eine lineare Funktion oder ein quadratisches Polynom in der Sekantenmethode-/
Muller-Methode benutzt, um eine Approximation der Nullstellen zu finden.
Diese Aufgaben sind Interpolationsaufgaben.
In den beiden Fallen konnen wir ein Polynom fur die Funktion f finden.
Dieses Polynom ist eindeutig, falls der Grad ≤ N ist.
Satz (Polynominterpolation)
Zu je N + 1 paarweise verschiedenen Zahlen x0, x1, . . ., xN in R und je
N +1 beliebigen Zahlen y0, y1, . . ., yN in R gibt es stets ein eindeutig bestimmtes
Interpolationspolynom P hochstens N -ten Grades mit der Eigenschaft
P (xj) = yj, j = 0, 1, . . . , N
40
KAPITEL 5. INTERPOLATIONSPOLYNOME 41
Beweis Wir fangen mit der Voraussetzung an, dass
P (x) = γ0 + γ1x + . . . + γNxN
ein solches Polynom ist, d.h. mit der Eigenschaft
P (xj) = yj , j = 0, 1, . . . , N.
Diese Eigenschaft bedeutet, dass die Koeffizienten γ0, γ1, . . ., γN dem fol-
genden linearen Gleichungssystem genugen mussen:
γ0 + γ1x0 + . . . + γNxN0 = y0
γ0 + γ1x1 + . . . + γNxN1 = y1
. . . . . . . . . . . . . . . . . . . . . . . . . . .
γ0 + γ1xN + . . . + γNxNN = yN
oder aquivalent:
1 x0 . . . xN0
1 x1 . . . xN1
......
......
1 xN . . . xNN
γ0
γ1
...
γN
=
y0
y1
...
yN
Die obige Matrix[xj
i
]ist eine Vandermond’sche Matrix mit der Determi-
nante
det[xj
i
]=
N∏
0≤i<j≤N
(xj − xi)
Fur paarweise verschiedene Stutzstellen gilt det[xj
i
]6= 0. Daher ist das
lineare Gleichungssystem eindeutig losbar fur γ0, γ1, . . ., γN .
Aber Vandermond’sche Matrizen sind schwierig algebraisch und numerisch
zu behandeln, insbesondere fur N ≫ 1. Es gibt andere Strategien, wobei wir
das obige lineare Gleichungssystem losen konnen.
5.1 Lagrange’sche Interpolationspolynome
Der Ansatz von Lagrange ist, dass wir das Polynom P durch
P (x) =N∑
j=0
f(xj)LN,j(x) mit f(xj) = yj fur j = 0, 1, . . . , N
KAPITEL 5. INTERPOLATIONSPOLYNOME 42
darstellen konnen, wobei LN,1(x), . . ., LN,N(x) Polynome des N -ten Grades
sind mit der Eigenschaft
LN,j(xi) = δi,j =
{1 falls i = j
0 sonst
(δi,j ist das Kronecker δ-Symbol)
1.) P ist Polynom hochstens N -ten Grades (Ausloschung moglich!)
2.) P (xi) =N∑
j=0
f(xj)LN,j(xi) =N∑
j=0
f(xj)δij = f(xi) for i = 0, 1, . . . , N
Wie sehen die Polynome LN,j aus? Es ist leicht zu zeigen, dass die folgenden
Polynome die erwunschten Eigenschaften besitzen:
LN,j(x) ≡N∏
i=0i6=j
x− xi
xj − xi
LN,j heißt das j-Lagrange’sche Interpolationspolynom des N -ten Grades fur die
Stutzstellen x0, x1, . . ., xN . (Man schreibt oft Lj statt LN,j fur festes N).
Beispiel Betrachte die Funktion f(x) = ex mit Stutzstellen 0, 1, 2. Dann gilt
x0 = 0, f(x0) = 1, x1 = 1, f(x1) = e, x2 = 2, f(x2) = e2
Die entsprechenden Lagrange’schen Interpolationspolynome (N = 2 hier)
lauten
L2,0(x) =(x− x1)(x − x2)
(x0 − x1)(x0 − x2)=
(x − 1)(x− 2)
(0 − 1)(0− 2)=
1
2x2 − 3
2x + 1
L2,1(x) =(x− x0)(x − x2)
(x1 − x0)(x1 − x2)=
(x − 0)(x− 2)
(1 − 0)(1− 2)= −x2 + 2x
L2,2(x) =(x− x0)(x − x1)
(x2 − x0)(x2 − x1)=
(x − 0)(x− 1)
(2 − 0)(2− 1)=
1
2x2 − 1
2x
Dann ist das gesuchte Interpolationspolynom
P (x) =
2∑
j=0
f(xj)L2,j(x)
= 1 · L2,0(x) + e · L2,1(x) + e2 · L2,2(x)
KAPITEL 5. INTERPOLATIONSPOLYNOME 43
= 1 ·(
1
2x2 − 3
2x + 1
)+ e ·
(−x2 + 2x
)+ e2 ·
(1
2x2 − 1
2x
)
=
(1
2e2 − e +
1
2
)x2 +
(−1
2e2 + 2e− 3
2
)x + 1
=1
2(e− 1)2x2 − 1
2(e− 1)(e− 3)x + 1
Siehe Vergleich des Interpolationspolynoms und des Taylor-Polynoms um x0
= 0, d.h.,2∑
j=0
1
j!f (j)(x0) (x− x0)
j = 1 + x +1
2x2
KAPITEL 5. INTERPOLATIONSPOLYNOME 44
> restart;> f1:=exp(x);
f1 := ex
f2:=convert(taylor(f1,x,3),polynom);
f2 := 1 + x +1
2x2
f3:=1-1/2*(a-1)*(a-3)*x+1/2*(a-1)ˆ2*xˆ2;
f3 := 1− 1
2(a− 1) (a− 3)x +
1
2(a− 1)2 x2
a:=evalf(exp(1));
a := 2.718281828
plot({f1,f2,f3},x=0..3,y=0..10,color=black);
0
2
4
6
8
10
y
0.5 1 1.5 2 2.5 3x
Nachteil der Lagrange’schen Darstellung eines Interpolationspolynoms: wir
mussen fast alles vom Anfang wieder herleiten, wenn wir eine zusatzliche Stutz-
stelle xN+1 einfuhren wollen.
Vorteil: Wir konnen die f(xi) leicht verandern.
5.2 Newton’sche Interpolationspolynome
Wir fangen mit dem Ansatz von Newton an, dass wir das gesuchte Interpolati-
onspolynom P durch
KAPITEL 5. INTERPOLATIONSPOLYNOME 45
P (x) = c0 + c1(x− x0) + c2(x− x0)(x− x1) + . . .
+cN (x− x0)(x− x1) . . . (x− xN−1)
darstellen konnen. Dann werden wir versuchen, die geeigneten Koeffizienten c0,
c1, . . ., cN herzuleiten, um die Eigenschaft
P (xj) = f(xj), j = 0, 1, . . . , N
zu haben.
j = 0 : P (x0) = c0 = f(x0)
j = 1 : P (xi) = c0 + c1(x1 − x0) = f(x1)
j = 2 : P (x2) = c0 + c1(x2 − x0) + c2(x2 − x0)(x2 − x1) = f(x2)
usw. ⇒
c0 = f(x0), c1 =f(x1)− f(x0)
x1 − x0
c2 =f(x2)− c0 − c1(x2 − x0)
(x2 − x0)(x2 − x1)
=
f(x2)− f(x0)−[f(x1)− f(x0)
x1 − x0
](x2 − x0)
(x2 − x0)(x2 − x1)
=
f(x2)− f(x1)
x2 − x1− f(x1)− f(x0)
x1 − x0
x2 − x0
Schreibe x2 − x0 oben als x2 − x0 = (x2 − x1) + (x1 − x0)
Wichtige Beobachtung: fur jedes k ≥ 0 hangen die c0, c1, . . ., ck nur von
(x0, f(x0)), (x1, f(x1)), . . ., (xk, f(xk)) ab.
⇒ Wir konnen durch Rekursion fortsetzen. Dafur brauchen wir”dividierte
Differenzen“. Definiere fur k = 0
f[xj] := f(xj) fur j = 0, 1, . . . , N ,
und dann rekursiv fur k ≥ 1
f[xj0 ,xj1 ,...,xjk] :=
f[xj1 ,...,xjk] − f[xj0 ,...,xjk−1
]
xjk− xj0
KAPITEL 5. INTERPOLATIONSPOLYNOME 46
fur jede Permutation {xj0 , . . . , xjk} von k + 1 der Stutzstellen {x0, . . ., xN}
Diese Terme heißen dividierte Differenzen. Wir konnen sie systematisch durch
Rekursion berechnen
f[x0]
ցf[x0,x1]
ր ցf[x1] f[x0,x1,x2]
ց ր ցf[x1,x2] f[x0,x1,x2,x3]
ր ց ր . . .
f[x2] f[x1,x2,x3]
ց ր . . .
f[x2,x3]
ր . . .
f[x3]
. . .
Satz ck = f[x0,x1,...,xk] fur k = 0, 1, . . ., N .
Der Beweis folgt aus dem Ansatz, der Definition und der Konstruktion.
Beispiel : Betrachte f(x) = ex mit x0 = 0, x1 = 1 und x2 = 2.
f[x0] = 1
ցf[x0,x1] = e−1
1−0 = e− 1
ր ցf[x1] = e f[x0,x1,x2] = (e2−e)−(e−1)
2−0 = 12 (e− 1)2
ց րf[x1,x2] = e2−e
2−1 = e2 − e
f[x2] = e2 ր
Daher sind die Koeffizienten
c0 = f[x0] = 1, c1 = f[x0,x1] = e− 1, c2 = f[x0,x1,x2] =1
2(e− 1)2
und das gesuchte Interpolationspolynom ist
P (x) = c0 + c1(x− x0) + c2(x − x0)(x − x1)
KAPITEL 5. INTERPOLATIONSPOLYNOME 47
= 1 + (e− 1)(x− 0) +1
2(e− 1)2(x− 0)(x− 1)
= 1 + (e− 1)
[1− 1
2(e− 1)
]x +
1
2(e− 1)2x2
= 1− 1
2(e− 1)(e− 3)x +
1
2(e− 1)2x2
genau wie die Lagrange’sche Darstellung des eindeutigen Interpolationspoly-
noms.
Das Newton’sche Interpolationspolynom fur die Daten (xj , f(xj)), mit j =
0, 1, . . . , N , ist definiert durch
P[x0,x1,...,xN ](x) = f[x0] + f[x0,x1] (x− x0) + . . .
+f[x0,x1,...,xN ] (x− x0)(x− x1) . . . (x− xN−1)
⇒ P[x0,x1,...,xN ](x) hat hochstens den Grad N und
P[x0,x1,...,xn](xj) = f[x0] + f[x0,x1](xj − x0) + . . .
+f[x0,...,xj ](xj − x0) . . . (xj − xj−1)
= . . . (Rest ≡ 0)
= f(xj)
Wir konnen solche Polynome durch Rekursion aufbauen: fur k = 0 gilt
P[x0](x) ≡ f[x0] = f(x0)
und dann fur k ≥ 1 haben wir
P[x0,...,xk](x) = P[x0,...,xk−1](x) + f[x0,...,xk] (x− x0) . . . (x− xk−1)
d.h. wir mussen nur einen zusatzlichen Term addieren!
⇒ sehr gunstig, falls wir eine zusatzliche Stutzstelle einfuhren wollen.
Nutzliche Eigenschaften
1). f[x0,...,xk] = 1k!
dk
dxk P[x0,...,xn](x)
KAPITEL 5. INTERPOLATIONSPOLYNOME 48
2). P[xj0 ,...,xjN ](x) ≡ P[x0,...,xN ](x)
fur alle Permutationen {xj0 , . . ., xjN} von {0, 1, . . . , N}
Beweis von 2). Definiere
Q(x) := P[xj0 ,...,xjn ](x) − P[x0,...,xN ](x)
⇒ Q ist ein Polynom mit Grad ≤ N und
Q(xj) = f(xj)− f(xj) = 0, j = 0, 1, . . .N
d.h. Q hat (N +1)-Nullstellen, aber am hochsten Grad N . Daher haben wir
Q(x) ≡ 0 ∀x.
Kapitel 6
Fehlerabschatzung
maple–Arbeitsblatt Interpolationspolynome: das Beispiel von Runge
Sei f : R → R (N + 1)-mal stetig differenzierbar. Die Taylor-Entwicklung
fur f um x0 lautet:
f(x) =
N∑
j=0
1
j!f (j)(x0) (x − x0)
j Taylor-Polynom N -ten Grad = fN(x)
+1
(N + 1)!f (N+1)(ξ) (x− x0)
N+1 Restterm mit ξ ∈ (x0, x)
Dafur haben wir die Fehlerabschatzung
|f(x)− fN (x)| ≤ MN+1
(N + 1)!|x− x0|N+1
wobei MN+1 := maxa≤ξ≤b
∣∣fN+1(ξ)∣∣ mit x0, x ∈ [a, b].
Wir wollen eine ahnliche Fehlerabschatzung fur Interpolationspolynome her-
leiten. Dafur werden wir die Newton’sche Darstellung benutzen.
Daten (xj , f(xj)) j = 0, 1, . . . , N
Wir definieren das Newton’sche Interpolationspolynom P[x0,x1,...,xN ] rekur-
siv durch.
P[x0](x) ≡ f[x0] = f(x0)
fur k = 0 und
P[x0,...,xk](x) := P[x0,...,xk−1](x) + f[x0,...,xk] (x− x0) . . . (x− xk−1)
fur k ≥ 1.
49
KAPITEL 6. FEHLERABSCHATZUNG 50
Dann gilt insbesondere
(i) P[x0,...,xN ](x) hat hochstens Grad N ;
(ii) P[x0,...,xN ](xj) = f(xj), j = 0, 1, . . . , N .
Wir haben nie gesagt, dass die Stutzstellen x0, x1, . . . , xN aufsteigend sein
mussen. Daher definiere das Intervall
I[x0,...,xN ] =
[min
i=0,...,Nxi, max
i=0,...,Nxi
]
Satz (Oevel: Satz 8.4b, S. 338)
Sei P[x0,...,xN ](x) das Interpolationspolynom einer (N+1)-mal stetig differen-
zierbaren Funktion f mit den paarweise verschiedenen Stutzstellen x0,x1,. . .,xN .
Dann gibt es einen Zwischenwert ξ ∈ I[x0,...,xN ] mit
f(x) = P[x0,...,xN ](x) +1
(N + 1)!f (N+1)(ξ) · (x− x0) . . . (x − xN )
fur x ∈ I[x0,...,xN ].
Vergleiche Taylor-Entwicklung!
Beweis
1). Sei x = xj . Von der Konstruktion gilt
f(xj) = P[x0,...,xN ](xj).
Daher gilt trivialerweise die Fehlerformel.
2). Sei x 6= xj fur j = 0, 1, . . . , N . Definiere (als Funktion von y mit festem x)
△(y) := f(y)− P[x0,...,xN ](y)− C[x0,...,xN ,x] · (y − x0) . . . (y − xN )
wobei
C[x0,...,xN,x] :=f(x)− P[x0,...,xN ](x)
(x− x0) . . . (x− xN )
Es ist offensichtlich, dass
1). △(y) (N + 1)-mal stetig differenzierbar ist.
2). △(y) mindestens N + 2 Nullstellen hat, namlich x0, . . ., xN und x.
Ordne diese Nullstellen aufsteigend an:
y(0)0 < y
(0)1 < . . . < y
(0)N+2
KAPITEL 6. FEHLERABSCHATZUNG 51
Satz von Rolle ⇒ ∃ y(1)j ∈
(y(0)j , y
(0)j+1
)mit △′
(y(1)j
)= 0 fur j = 0, . . .,
N + 1
Setze induktiv fort: k = 1, 2, . . . , N
∃ y(k+1)j ∈
(y(k)j , y
(k)j+1
)mit △(k+1)
(y(k+1)j
)= 0 fur j = 0, . . . , N + 1− k.
Insbesondere fur k = N haben wir
∃ ξ = y(N+1)0 ∈
(y(N)0 , y
(N)1
)mit △(N+1)(ξ) = 0
Aber △(N+1)(y) = f (N+1)(y)− (N + 1)! C[x0,...,xN ,x]
Deshalb gilt:
C[x0,...,xN ,x] =1
(N + 1)!f (n+1)(ξ) mit ξ ∈ I[x0,...,xN ]
Korollar (Oevel: Satz 8.4a, S. 338)
Sei f k-mal stetig differenzierbar und x0, . . ., xk paarweise verschieden.
Dann existiert ξk ∈ I[x0,...,xk] mit
f[x0,...,xk] =1
k!f (k)(ξk).
Beispiel f(x) = ex mit x0 = 0, x1 = 1 und x2 = 2. Dann ist
f(x) = P[0,1,2](x) +1
3!f (3)(ξ) · (x− 0)(x− 1)(x− 2)
KAPITEL 6. FEHLERABSCHATZUNG 52
mit
P[0,1,2](x) = 1− 1
2(e− 1)(e− 3)x +
1
2(e− 1)2x2
Es gilt
(i) f (3)(ξ) = eξ ≤ e2 < 9 fur ξ ∈ [0, 2]
(ii) max0≤x≤2
|x(x − 1)(x− 2)| = 23√
3
Daher lautet die Fehlerabschatzung
∣∣f(x)− P[0,1,1](x)∣∣ ≤ 1
3!
∣∣∣f (3)(e)∣∣∣ |x(x − 1)(x− 2)|
≤ 1
69
2
3√
3∀ x ∈ [0, 2]
=1√3≃ 0.577
ziemlich grob, aber gleichmaßig auf [0, 2]
6.1 Differenzenoperatoren
Schreibe fj = f(xj) fur j = 0, 1, . . ., N und paarweise verschiedenen Stutzstellen
x0, . . ., xN .
Wir definieren zwei Arten Differenzenoperatoren (1 ≤ j ≤ N)
(1) Vorwartsdifferenz/ aufsteigend
∆fj := fj+1 − fj , ∆ = Delta
(2) Ruckwartsdifferenz/ absteigend
∇fj := fj − fj−1 ∇ = Nabla
Wir werden diese Operationen sukzessive verwenden. Deswegen definieren
wir
∆0fj := fj, ∇0fj := fj , k = 0
∆k+1fj := ∆kfj+1 −∆kfj, ∇k+1fj := ∇kfj −∇kfj−1 k > 0
Dann kann man zeigen, dass
KAPITEL 6. FEHLERABSCHATZUNG 53
∆kfj =
k∑
l=0
(−1)k−l
(k
l
)fj+l
und
∇kfj =
k∑
l=0
(−1)l
(k
l
)fj−l
wobei (k
l
):=
k!
l!(k − l)!
Bemerkung: Die ∆/∇-Notation kommt aus Stummel/Hainer. Die Bucher
von Oevel und Schwarz benutzen andere Notationen.
6.2 Aquidistante Stutzstellen
Betrachte die aquidistanten Stutzstellen
xj = x0 + jh, j = 0, 1, . . . , N
mit aquidistanter Schrittweite h > 0 und schreibe
fj = f(xj), j = 0, 1, . . . , N.
Wir konnen die obigen Differenzenoperatoren benutzen, um eine kompakte
Darstellung des Interpolationspolynoms in diesem Fall herzuleiten. Wir fangen
mit dem Newton’schen Interpolationspolynom an:
P[x0,...,xN ](x) = f[x0] + f[x0,x1] (x− x0) + . . .
. . . + f[x0,...,xN ] (x− x0) . . . (x− xN−1)
Fur aquidistante Stutzstellen haben wir:
f[xj] = fj = ∆0fj
f[xj,xj+1] =fj+1 − fj
xj+1 − xj=
1
h∆fj
f[xj,xj+1,xj+2] =f[xj+1,xj+2] − f [xj , xj+1]
xj+2 − xj=
1h∆fj+1 − 1
h∆fj
2h=
1
2h2∆2fj
Durch Induktion gilt dann im Allgemeinen:
KAPITEL 6. FEHLERABSCHATZUNG 54
f[xj,...,xj+k] =1
k!hk∆kfj
Deshalb konnen wir das Newton-Polynom umschreiben:
P[x0,...,xN ](x) = ∆0f0 +(x− x0)
1!h1∆1f0 + . . . +
(x− x0) . . . (x− xN−1)
N !hN∆Nf0
Schreibe x = x0 + th mit t = x−x0
h ∈ R (t muss nicht ganzzahlig sein).
Dann haben wir
P[x0,...,xN ](x0 + th) = ∆0f0 + 11! t∆
1f0 + . . .
+ 1N ! t(t− 1) . . . (t−N + 1) ·∆Nf0
=N∑
k=0
(tk
)∆kf0
wobei(
t
k
)=
t!
k!(t− k)!=
t(t− 1) . . . (t− k + 1) . . .
k!(t− k)(t− k − 1) . . .=
t(t− 1) . . . (t− k + 1)
k!
mit (0
0
)=
0!
0!0!= 1
Die obige Formel heißt
”Newton-Gregory-Interpolationsformel“
Es gibt eine entsprechende Formel mit den Ruckwartsdifferenzenoperatoren:
P[x0,...,xN ](xN − th) =
N∑
k=0
(−1)k
(t
k
)∇kfN
Solche Interpolationsformeln sind sehr wichtig fur theoretische Entwicklun-
gen.
ABER Es gibt Schwierigkeiten mit aquidistanten Stutzstellen — der Fehler
kann sehr groß in der Nahe der Randstellen sein.
Fehler (x) = f(x)− P[x0,...,xN ](x)
= 1(N+1)! f (N+1)(ξ) (x− x0)(x − x1) . . . (x− xN )
= 1(N+1)! hN+1 f (N+1)(ξ) t(t− 1) . . . (t−N)
mit t = (x− x0)/h ∈ [0, N ] (nicht zwingend ganzzahlig!)
KAPITEL 6. FEHLERABSCHATZUNG 55
|Fehler (t)| ≤ MN+1 hN+1
(N + 1)!|t(t− 1) . . . (t−N)|
wobei MN+1 := maxx0≤ξ≤xN
|f (N+1)(ξ)| <∞.
Die Ursache der Schwierigkeiten ist, dass
EN := max0≤t≤N
|t(t− 1) . . . (t−N)|
sehr groß sein kann mit EN →∞ fur N →∞
Eine wichtige Folge ist, dass das Interpolationspolynom P[x0,...,xN ](x) sehr
große Schwankungen in der Nahe der Randstellen x0, xN besitzen kann. Wie
gut ist dann die Approximation ?
Runges Gegenbeispiel
f(x) =1
1 + x2mit x ∈ [−5, 5]
der Fehler hier ist nicht gleichmaßig beschrankt ist
Siehe maple–Arbeitsblatt”Interpolationspolynome: das Beispiel von Run-
ge“
Fragen
1.) Konnen wir (und wie) die Stutzstellen wahlen, um diesen Fehler zu
minimieren?
2.) Sollen wir N ≫ 1 benutzen? −→ Spater:”SPLINES“
Kapitel 7
Spline–Interpolation
Literatur Oevel, Kap. 8.2.1-4; Schwarz, Kap. 3.7.1-3.
Aufgabe: Interpolation einer nur durch Stutzwerte gegebenen Funktion f(x):
Stutzstellen x0 < x1 < . . . < xn oder”Knoten“
Funktionswerte fi = f(xi), i = 0 . . . n
Interpolationspolynome hoheren Grades weisen in der Regel unerwunschte
Oszillationen auf — Runges Beispiel!
Zweckmaßiger: Interpolation auf Teilintervallen [xi, xi+1] mit Polynomen
niedrigen Grades.
Anschaulich:
1). Verbinde die Punkte (xi, fi) durch Geradenstucke
Das englische Wort spline bedeutet”dunner Stab“
”halte Stab an den Punkten (xi, fi) fest“
2).”Ziehe eine glatte Kurve“ durch die Punkte = Kurvenlineal
56
KAPITEL 7. SPLINE–INTERPOLATION 57
Beispiele:
(1) S(x) = ai + bix auf [xi, xi+1] — linearer Spline
(2) S(x) = ai + bix + cix2 auf [xi, xi+1] — quadratischer Spline
(3) S(x) = ai + bix + cix2 + dix
3 auf [xi, xi+1] — kubischer Spline
S(x) heißt Splinefunktion vom Grad k
mit
i) S ist auf [xi, xi+1] Polynom vom Grad ≤ k
ii) S(xi) = fi, i = 0, . . . , n
iii) S(x) ist stetig differenzierbar bis zur Ordnung k − 1
Bemerkung: Polynome immer stetig und differenzierbar.
⇒{
I) S(xi − 0) = S(xi + 0) = fi, i = 1, . . . , n− 1
II) S(m)(xi − 0) = S(m)(xi + 0), i = 1, . . . , n− 1, m = 1, . . . , k − 1
(Limes von der linken Seite und Limes von der rechtenen Seite)
Allgemeine Darstellung: S(x) auf [xi, xi+1] Polynom k-ten Grades
S(x) = a[k]i xk + a
[k−1]i xk−1 + . . . + a
[1]i x + a
[0]i
mit (k + 1) · n freien Parametern und
⇒
Interpolationsbedingung : (n + 1) Gleichungen
Glattheitsbedingung : k · (n− 1) Gleichungen
d.h. es verbleiben k−1 Freiheitsgrade. Fur k > 1 mussen zusatzliche Bedingun-
gen gestellt werden um die Eindeutigkeit zu sichern.
S [k][x0,...,xn] = { Spline vom Grad k zu den Stutzstellen x0 < . . . < xn }
Es zeigt sich, dass S [k][x0,...,xn] einen linearen Unterraum von Ck−1[x0, xn] bil-
det.
Dim(S [k]
[x0,...,xn]
)= k + n
durch n + 1 Punkte und k − 1 Freiheitsgrade
KAPITEL 7. SPLINE–INTERPOLATION 58
7.1 Lineare Splines
Aufgabe: Bestimme einen/den linearen Spline S(x) ∈ S [1][x0,...,xn] zu den gege-
benen Stutzstellen x0 < . . . < xn mit den Stutzwerten f0, . . ., fn.
Auf [xi, xi+1] ist S(x) eine lineare Funktion
S(x) = ai + bi(x− xi), x ∈ [xi, xi+1]
• Interpolationsbedingung: S(xi) = fi, i = 0, . . . , n
• Stetigkeitsbedingung: S(xi − 0) = S(xi + 0), : i = 1, . . . , n− 1
⇒ 2n Unbekannte zu den 2n Gleichungen ⇒ eindeutige Losung!
Klar: S(xi) = ai = fi
• Ferner gilt
ai−1 + bi−1 ((xi − 0)− xi−1) = ai + bi ((xi + 0)− xi)
ai−1 + bi−1(xi − xi−1) = ai
Eindeutige Losung
S(x) =(xi+1 − x)fi + (x− xi) fi+1
xi+1 − xiauf [xi, xi+1]
7.1.1 Fehlerabschatzung
Sei f eine zweifach stetig differenzierbare Funktion. Aus der linearen Interpola-
tion auf [xi, xi+1] folgt:
f(x)− S(x) =1
2f ′′(ξ) (x− xi)(x − xi+1), ξ ∈ [xi, xi+1].
Da |x− xi| · |x− xi+1| ≤(xi+1 − xi)
2
4gilt
maxx∈[x0,xn]
|f(x)− S(x)| ≤ 1
8h2 max
ξ∈[x0,xn]|f ′′(ξ)|
wobei h = maxi
(xi+1 − xi).
⇒ mit h→ 0 konvergiert S(x) gleichmaßig gegen f(x).
KAPITEL 7. SPLINE–INTERPOLATION 59
Wie sieht eine Basis fur den linearen Raum S[1][x0,...,xn] aus?
Betrachte die sogenannten”Dach“-Funktionen
gi(x) =
(x− xi−1)/hi−1, x ∈ [xi−1, xi]
(xi+1 − x)/hi, x ∈ [xi, xi+1]
0, sonst,
wobei hi = xi+1 − xi. Dann gilt
S(x) =n∑
k=0
fk gk(x)
7.2 Quadratische Splines
Da Dim(S [2]
[x0,...,xn]
)= n + 2 muß eine zusatzliche Randbedingung gestellt wer-
den (n + 1 Stutzstellen), z.B.
S′(x0) = f ′0 oder S′(xn) = f ′
n
mit gegebenen f ′0 und f ′
n.
Die allgemeine Darstellung auf [xi, xi+1] lautet
S(x) = ai + bi(x− xi) + ci(x − xi)2
Interpolationsbedingung S(xi) = ai = fi, i = 0 . . . n
Stetigkeitsbedingung S(xi − 0) = S(xi + 0), i = 1 . . . n− 1
Differenzierbarkeitsbedignungen S′(xi − 0) = S′(xi + 0), i = 1 . . . n− 1
⇒ (n + 1) + 2(n− 1) = 3n− 1 Gleichungen fur 3n Unbekannte
Die Losung des linearen gestalteten Gleichungssystems fur f ′i
f ′i + f ′
i+1 = 2fi+1 − fi
xi+1 − xi, i = 0, . . . , n− 1
lautet
ai = fi, bi = f ′i , ci = 2
f ′i+1 + f ′
i
xi+1 − xi, i = 0, . . . , n− 1
KAPITEL 7. SPLINE–INTERPOLATION 60
7.3 Kubische Splines
In der Praxis wichtigster Fall
Dim(S [3]
[x0...xn]
)= n + 3
Fur die Eindeutigkeit sind zwei zusatzliche Bedingungen aufzustellen!
Randbedingungen
a) naturliche S′′(x0) = S′′(xn) = 0
b) vollstandige S′(x0) = f ′0 und S′(xn) = f ′
n
c) periodische︸ ︷︷ ︸S(x0)=S(xn)
S′(xn) = S′(x0) und S′′(xn) = S′′(x0)
d) not-a-knot Stetigkeit von S′′ an x1 und xn−1
Satz: Unter a),b),c) oder d) ist der kubische Spline eindeutig bestimmt.
Bemerkungen: Aus d) folgt, dass S(x) auf [x0, x2] und auf [xn−2, xn] ein Po-
lynom dritten Grades ist. (x1 und xn−1 sind keine eigentlichen Knoten.)
Allgemeine Form: auf [xi, xi+1]
S(x) = ai + bi(x− xi) + ci(x− xi)2 + di(x − xi)
3
Interpolationsbedingung S(xi) = fi
⇒ ai = fi
Glattheitsbedingungen
(I) : S(xi − 0) = S(xi + 0)
(II) : S′(xi − 0) = S′(xi + 0)
(III) : S′′(xi − 0) = S′′(xi + 0)
Die Bedingungen (I) und (III) liefern
KAPITEL 7. SPLINE–INTERPOLATION 61
bi =fi+1 − fi
hi− hi
6(f ′′
i+1 + 2f ′′i ), ci =
f ′′i
2, di =
f ′′i+1 − f ′′
i
6hi,
Zusatzlich (II)
bi−1 + 2ci−1hi−1 + 3di−1h2i−1 = bi, i = 1, . . . , n− 1
Nach dem Einsetzen erhalten wir ein lineares Gleichungssystem zur Bestim-
mung der f ′′0 , . . ., f ′′
n
hi−1f′′i−1 + 2(hi−1 + hi)f
′′i + hif
′′i+1 = 6
fi+1 − fi
hi− 6
fi − fi−1
hi−1︸ ︷︷ ︸=:δi
Zur Bestimmung der n + 1 Unbekannten stehen nur n − 1 Gleichungen zur
Verfugung. Die fehlenden Gleichungen werden durch die Randbedingungen ge-
liefert.
• Fall a) Naturliche Randbedingungen S′′(x0) = S′′(xn) = 0
f ′′0 = 0 und f ′′
n = 0
Aufstellen des linearen Gleichungssystems in Matrixform
2(h0 + h1) h1 0 ©h1 2(h1 + h2) h2
0 h2 2(h2 + h3) h3
. . .. . .
. . .. . .
© . . .. . .
. . . hn−2
hn−2 2(hn−2 + hn−1)
·
f ′′1
····
f ′′n
=: ~δ
⇒ symmetrische, tridiagonale Matrix, die zusatzlich diagonal dominant ist.
Dieses Gleichungssystem kann schnell (0(n)-Berechnungen) und numerisch
stabil gelost werden.
• Falle b), c) und d) fuhren auf ahnlicher Weise auf”fast“ tridiagonale Glei-
chungssysteme die schnell und stabil gelost werden konnen.
Bemerkungen: Sind h1 = . . . hn−1 = h, Konstante, so kann durch h dividiert
werden und es ergibt sich
KAPITEL 7. SPLINE–INTERPOLATION 62
4 1
1 4 1
1 4 1
1 4 1. . .
. . .
1 4
~f ′′ =1
h~δ
7.3.1 Die Minimaleigenschaft kubischer Splines
Kubische Splines sind durch die Eigenschaft ausgezeichnet, dass sie die”mittlere
Krummung“ minimieren:
E(f) =
xn∫
x0
(f ′′(x))2
dx
Bemerkungen: Physikalisches Prinzip (Elastizitatstheorie)
Deformationsenergie K(f) =
xn∫
x0
(f ′′(x))2
1 + (f ′(x))2dx minimal
Fur |f ′(x)| ≪ 1 geht K(f) in E(f) uber.
Interpolation:
Wird ein dunner (biegsamer) Stab an den Stutzwerten (xi, fi) festgehalten, so
nimmt der Stab zwischen den Stutzstellen, die Kurve, die sich aus der Losung
des Minimierungsproblems fur K(f) ergibt, an.
E(f) ist einfach mit Hilfe von Standardmethoden der Variationsrechnung
(siehe Oevel S. 354) auswertbar und liefert
f ′′′′(x) = 0.
d.h. die minimierende Funktion ist stuckweise ein Polynom maximal dritten
Grades.
Bemerkungen: Die naturliche Randbedingung S′′(x0) = S′′(xn) = 0 besagt,
dass in x0 und xn keine Krummung und damit keine Krafte auftreten.
7.3.2 Fehlerabschatzung
Wie gut ist die Naherung zwischen den Stutzstellen?
Hier betrachten wir nur vollstandige Randbedingung (andere Falle sind ahnlich).
KAPITEL 7. SPLINE–INTERPOLATION 63
Satz (Interpolationsfehler kubischer Splines.)
Sei f vierfach stetig differenzierbar und S der zugehorige kubische Spline mit
S(x0) = f(x0), . . . . . . , S(xn) = f(xn
mit den vollstandigen Randbedingungen
S′(x0) = f ′(x0) und S′(xn) = f ′(xn) .
Mit h := maxi=1...n
(xi − xi−1) und
K =maxi(xi − xi−1)
mini(xi − xi−1)· max
ξ∈[x0,xn]
∣∣∣f (4)(ξ)∣∣∣
gilt fur alle x ∈ [x0, xn]
|S(x) − f(x)| ≤ h4K
|S′(x)− f ′(x)| ≤ 2h3K
|S′′(x)− f ′′(x)| ≤ 2h2K
|S′′′(x)− f ′′′(x)| ≤ 2hK
Beweis: Taylorentwicklung und Auswertung der Matrix, siehe Oevel S. 350.
Bemerkungen: Damit nicht nur gleichmaßige Konvergenz
|S(x)− f(x)| = 0(h4)
sondern auch gleichmaßiger Konvergenz fur die Ableitung
Bemerkungen: Falls die Funktion f(x) den Randbedingungen nicht genugt,
ist der Fehler in der Regel nur quadratisch |S(x)− f(x)| = 0(h2).
Bemerkungen: Ein kubischer Spline kann wieder mittels Basisfunktionen
(B-Splines) dargestellt werden:
S(x) =∑
j
αjBj(x)
Die B-Splines Bj sind von den Stutzstellen x0, . . ., xn und den Randbedingun-
gen abhangig.
KAPITEL 7. SPLINE–INTERPOLATION 64
Die Interpolationsbedingungen
∑
j
αjBj(xi) = fi
liefern ein (einfaches) lineares Gleichungssystem fur αj ! (siehe Oevel S. 362)
Anwendung: numerische Integration, Losungsdarstellung von Differential-
gleichungen
Weitere Ansatze: Approximation, Least squares–Approximation
Kapitel 8
Orthogonale Polynome
Literatur Schwarz, Kap.4.3; Stummel/Hainer, Kap. 7.4
maple–Arbeitsblatt Tschebyscheff-Entwicklung
Wir haben schon verschiedene Arten von Polynome mit Grad hochstens N
betrachtet, die eine (N + 1)-mal stetig differenzierbare Funktion f : R → R
irgendwie darstellen, d.h. Taylor-Polynome und Interpolationspolynome.
Gibt es ein”bestes“ Polynom? Was bedeutet
”bestes“ hier?
Wir konnen diese Fragen mit orthogonalen Polynomen beantworten. Wir
konnen auch orthogonale Polynome benutzen, um die”besten“ Stutzstellen fur
Interpolationspolynome zu finden, d.h., im Sinne der Term
maxx∈I|(x− x0) · · · (x− xN )|
in der Fehlerabschatzung
maxx∈I
∣∣f(x)− P[x0,...,xN ](x)∣∣ ≤ MN+1
(N + 1)maxx∈I|(x − x0) . . . (x− xN )|
uber alle paarweise verschiedene x0, . . . , xN zu minimieren.
8.1 Ein bißchen Lineare Algebra
Vektoren in Rd x = (x0, . . . , xd)
Skalarprodukt < x, y > :=d∑
i=1
xiyi
65
KAPITEL 8. ORTHOGONALE POLYNOME 66
⇒ < x, x > =d∑
i=1
x2i ≥ 0 fur alle x ∈ Rd mit < x, x > = 0 genau
dann, wenn x = 0
x, y heißen orthogonal, wenn < x, y >= 0, d.h. senkrecht in R2 oder R3
Norm ‖x‖ :=√
< x, x > =
√d∑
i=1
x2i
Wir wollen diese Begriffe zu dem linearen Funktionenraum C[a, b] aller steti-
gen Funktionen f : [a, b]→ R verallgemeinern —”linear“ bezuglich punktweiser
Addition und skalarer Multiplikation, d.h.
(f + g)(x) := f(x) + g(x) ∀x ∈ [a, b]
(αf)(x) := αf(x) ∀α ∈ R .
Das einfachste Skalarprodukt auf C[a, b] lautet:
< f, g > :=
b∫
a
f(x)g(x) dx
Es ist oft gunstig, eine Gewichtsfunktion in das Skalarprodukt aufzunehmen,
d.h. eine positive stetige Funktion w : [a, b] → R — d.h. mit w(x) > 0 fur alle
x ∈ [a, b].
Dann definieren wir das Skalarprodukt < ·, · >w mit Gewicht w durch
< f, g >w :=
b∫
a
w(x)f(x)g(x) dx
Es gilt
< f, f >w =
b∫
a
w(x)f(x)2 dx ≥ 0 ∀f ∈ C[a, b]
mit < f, f >w = 0 genau dann, wenn f(x) ≡ 0 (wegen der Stetigkeit von f und
Positivitat von w).
Wir konnen auch die entsprechende Norm auf C[a, b] definieren:
‖f‖w :=√
< f, f >w =
√∫ b
a
w(x)f(x)2 dx
Nun konnen wir den Begriff der”Orthogonaliltat“ (bezuglich eines Skalar-
produkts) in C[a, b] einfuhren: f , g ∈ C[a, b] heißen orthogonal bzg. des Ska-
larprodukts < ·, · >w, wenn
KAPITEL 8. ORTHOGONALE POLYNOME 67
< f, g >w= 0
Beispiel a = −1, b = 1 und w(x) ≡ 1
< f, g > =
1∫
−1
f(x)g(x) dx
Betrachte f(x) = sin(πx) und g(x) = cos(πx).
< f, g > =
∫ 1
−1
sin(πx) cos(πx) dx
=1
2πsin2(πx)
∣∣∣1
−1= 0 orthogonal!
Eine Familie {ϕ0, ϕ1, . . . , ϕk, . . .} von Polynomen mit den Eigenschaften
1.) ϕk hat Grad k, k = 0, 1, 2, . . .
2.) < ϕi, ϕj >w = 0 fur alle i 6= j paarweise orthogonal
heißt Orthogonalpolynomsystem (bzg. des Skalarprodukts < ., . >w).
Wir konnen den Gram-Schmidt-Algorithmus verwenden, um ein Orthogo-
nalpolynomsystem fur jedes Skalarprodukt zu konstruieren.
k = 0 ϕ0(x) ≡ 1
k > 0 ϕk(x) = xk −k−1∑
j=0
< xk, ϕj >w
< ϕj , ϕj >wϕj(x)
genau wie fur Vektoren im Rd!
Beispiel (1): Betrachte < f, g > =∫ 1
−1f(x)g(x) dx
k = 0 ϕ0(x) ≡ 1
k = 1 ϕ1(x) = x− < x, ϕ0 >
< ϕ0, ϕ0 >· ϕ0(x)
= x− 0∫ 1
−112dx
· 1 = x
KAPITEL 8. ORTHOGONALE POLYNOME 68
k = 2 ϕ2(x) = x2 − < x2, ϕ0 >
< ϕ0, ϕ0 >· ϕ0(x) − < x2, ϕ1 >
< ϕ1,ϕ1 >· ϕ1(x)
= x2 −∫ 1
−1 x2 · 1dx∫ 1
−112dx
· 1− 0∫ 1
−1x2dx
= x2 − 1
3.
Die Polynome ϕ0(x) = 1, ϕ1(x) = x, ϕ2(x) = x2−1/3, . . . heißen Legendre Polynome.
Sie genugen der Rekursionsformel (schneller!)
ϕk+1(x) = x · ϕk(x)− k2
4k2 − 1· ϕk−1(x) k ≥ 1
mit ϕ0(x) ≡ 1 und ϕ1(x) ≡ x.
Beispiel(2): Betrachte < f, g >w =1∫
−1
1√1− x2
f(x)g(x) dx mit der Gewicht-
funktion w(x) =1√
1− x2(nur fur x ∈ (−1, 1) definiert, aber OK !).
k = 0 ϕ0(x) = 1
k = 1 ϕ1(x) = x− < x, ϕ0 > w
< ϕ0, ϕ0w· ϕ0(x)
= x− 0∫ 1
−1
1√1− x2
dx
· 1 = x
k = 2 ϕ2(x) = x2 − < x2, ϕ0 >w
< ϕ0, ϕ0 >w· ϕ0(x)− < x2, ϕ1 >w
< ϕ1, ϕ1 >w· ϕ1(x)
= x2 −
∫ 1
−1
x2
√1− x2
dx
∫ 1
−1
1√1− x2
dx
· 1− 0∫ 1
−1
x2
∫1− x2
dx
· x
benutze x = cosϕ fur welches
dx = − sinϕ · dϕ = −√
1− cos2 ϕ · dϕ
oder
dx√1− x2
= −dϕ⇒ ϕ2(x) = x2 − 2∫ π/2
0cos2 ϕdϕ
2∫ π/2
0dϕ
· 1 = x2 − 1
2
KAPITEL 8. ORTHOGONALE POLYNOME 69
usw.
Bemerkung: Wir konnen jedes ϕk mit einer Konstante multiplizieren, oh-
ne die Orthogonalitatseigenschaft zu verletzen. In diesem Fall betrachten wir
2ϕ2(x) = 2x2−1 statt ϕ2(x). Dann haben wir Tschebyscheff Polynome, die der
folgenden Rekursionsformel genugen
Tk+1(x) = 2xTk(x)− Tk−1(x) k ≥ 1
mit T0(x) = 1 und T1(x) = x
Ein sehr wichtiger Vorteil ist:
Tk(cosϕ) = cos(kϕ) k ≥ 0
Der Beweis folgt durch die Rekursionsformel und trigonometrische Iden-
titaten (siehe maple-Worksheet).
8.2 Eigenschaften von Orthogonalpolynomen
Orthogonalpolynomsysteme haben einige wichtige Eigenschaften. Sei PN der li-
neare Raum aller Polynome vom Grade hochstens N und sei {ϕ0, ϕ1, . . . , ϕk, . . .}ein Orthogonalpolynomsystem bzg. eines Innterprodukts < ., . >w
Eigenschaft 1: {ϕ0, ϕ1, . . . , ϕN} ist eine Basis des linearen Raums PN .
Beweis Hausarbeit!
Eigenschaft 2: Jedes Polynom p ∈ PN hat die Darstellung
p(x) =
N∑
k=0
ckϕk(x) mit ck =< p, ϕk >w
< ϕk, ϕk >w, k = 0, 1, . . . , N
Beweis Betrachte ein festes j ∈ {0, 1, . . . , N}. Dann gilt
⟨p−
N∑k=0
ckϕk, ϕj
⟩
w
= 〈p, ϕj〉w −N∑
k=0
ck 〈ϕk, ϕj〉w
= < p, ϕj >w −cj < ϕj , ϕj >w
= 0 (Definition von cj)
⇒ p−N∑
k=0
ckϕk ≡ 0
KAPITEL 8. ORTHOGONALE POLYNOME 70
Eigenschaft 3: < p, ϕN+1 >w = 0 ∀ p ∈ PN .
Beweis p =N∑
k=0
ckϕk mit ck wie oben.
< p, ϕN+1 >w =
⟨N∑
k=0
ckϕk, ϕN+1
⟩
w
=
N∑
k=0
ck < ϕk, ϕN+1 >w
=
N∑
k=0
ck · 0 = 0 (Orthogonalitat)
Eigenschaft 4: ϕk hat genau k einfache reelle Nullstellen x(k)1 , . . . , x
(k)k auf dem
offenen Intervall (a, b) fur k ≥ 1.
Beweis Als Polynom k-ten Grades hat ϕk bis k reelle Nullstellen, angenommen
ϕk hat weniger als k Nullstellen,
x1, . . . , xj , j < k,
mit Vielfachheit
n1, . . . , nj ⇒ n1 + . . . + nj ≤ k.
Daher konnen wir schreiben
ϕk(x) = Qk(x)
j∏
i=1
(x− xi)ni
wobei Qk(x) Polynom des Grades k − n1 − . . . − nj ≥ 0 mit keinen reellen
Nullstellen
⇒ Qk(x) > 0 ∀ x ∈ (a, b) oder Qk(x) < 0 ∀x ∈ (a, b).
Betrachte jetzt das Polynom
p(x) =
j∏
i=1
(x− xi)mi
wobei
mi =
0, falls ni gerade
1, falls ni ungerade
p hat Grad m1 + . . . + mj ≤ j < k und
KAPITEL 8. ORTHOGONALE POLYNOME 71
< p, ϕk >w =
b∫
a
w(x)Q(x)
j∏
i=1
(x− xj)ni+mi dx
wobei der erste Term und der letzte Term nichtnegativ sind (ni +mi ist gerade!)
⇒ < p, ϕk >w 6= 0.
Aber p ∈ Pj mit j < k, also mit der 3ten Eigenschaft gilt < p, ϕk > = 0.
Widerspruch! (zu ϕk vom Grade k)
Beispiel Betrachte die Tschebyscheff-Polynome mit x = cosϕ:
Tk(x) = Tk(cosϕ) = cos(kϕ) = 0
mit kϕ = (2j − 1)π2 fur j = 1, . . . , k (dann wiederholt sich periodisch!)
⇒ x(j) = cos
(2j − 1
k
π
2
)j=1,. . . ,k
8.3 Die”beste“ Polynomapproximation
Sei f : [a, b]→ R stetig und sei {ϕ0, ϕ1, . . .} ein Orthogonalpolynomsystem bzg.
eines Linearprodukts < ., . >w.
Betrachte PN mit der Basis {ϕ0, . . . , ϕN}. Jedes Polynom pN ∈ PN hat bzg.
dieser Basis eine Darstellung
pN =
N∑
k=0
ckϕk
Gibt es ein Polynom p∗N ∈ PN mit der Eigenschaft
||f − p∗N ||w = minpN∈PN
||f − pN ||w
JA! p∗N =N∑
k=0
c∗kϕk mit
c∗k =< f, ϕk >w
< ϕk, ϕk >w
Bemerkung Diese Aufgabe ist eine Art verallagemeinerte”kleinste Quadratmit-
tel“–Aufgabe: finde das Minimum der Funktion.
F (c0, c1, . . . , cN ) : = ‖f − pN‖2w =
⟨f −
N∑
k=0
ckϕk, f −N∑
k=0
ckϕk
⟩
w
KAPITEL 8. ORTHOGONALE POLYNOME 72
d.h. lose∂F
∂cj= 0 j = 0, 1, . . . , N
Wegen der Orthogonalitat haben wir
F (c0, c1, . . . , cN ) = < f, f >w −2N∑
k=0
ck < f, ϕk >w
+
N∑
k=0
c2k < ϕk, ϕk >w
Deshalb gilt
∂F
∂cj= −2 < f, ϕj >w +2cj < ϕj , ϕj >w
und∂2F
∂c2j
= 2 < ϕj , ϕj >w > 0
Geometrisch < f − p∗N , p∗N >w = 0 d.h. p∗N ist die orthogonale Projektion
der Funktion f ∈ C[a, b] auf den Teilraum PN
PN
f
f ⊥ p∗N
p∗N
Beispiel Betrachte die Funktion f(x) = ex fur −1 ≤ x ≤ 1 und das Skalarpro-
dukt < f, g > =1∫
−1
f(x)g(x) dx.
Bezuglich obiger Norm verwenden wir die Legendre Polynome ϕj , also
ϕ0(x) = 1, ϕ1(x) = x, ϕ2(x) = x2 − 1/3, . . . .
Betrachte P1, d.h. alle Polynome von grad 0 oder 1. Wir wollen
p∗1(x) = c∗0ϕ0(x) + c∗1ϕ1(x) = c∗0 + c∗1x
wobei
c∗0 =< f, ϕ0 >
< ϕ0, ϕ0 >=
∫ 1
−1ex, 1dx
∫ 1
−1 12dx=
1
2(e− e−1)
und
c∗1 =
∫ 1
−1 ex, xdx∫ 1
−1x2dx
=2e−1
2/3= 3e−1
KAPITEL 8. ORTHOGONALE POLYNOME 73
⇒ die”beste“ lineare Approximation in P1 der Funktion f(x) = ex auf −1 ≤
x ≤ 1 (bzg. < ., . > oben) lautet
p∗1(x) =1
2(e− e−1) + 3e−1x
Mit einem anderen Skalarprodukt wurden wir im Allgemeinen ein anderes
Polynom bekommen.
8.4 Die”besten“ Stutzstellen
Wir wollen die Stutzstellen x0, . . ., xN ∈ [a, b] wahlen, um die Fehlerabschatzung
des Interpolationspolynoms P[x0,...,xN ](x) auf [a, b] zu minimieren:
∣∣f(x)− P[x0,...,xN ](x)∣∣ ≤ MN+1
(N + 1)!max
a≤x≤b|(x− x0) . . . (x− xN )|
Insbesondere: Wir wollen
maxa≤x≤b
|(x− x0) . . . (x− xN )|
minimieren durch eine geeignete Wahl der Stutzstellen x0, . . . , xN
Satz (Schwarz: Satz 4.12, Seite 218)
Unter allen Polynomen Pn(x) von Grad n ≥ 1, deren Koeffizient von xn
gleich Eins ist, hat Tn(x)/2n−1 die kleinste Maximumnorm im Intervall [−1, 1],
d.h. es gilt
minpn
max−1≤x≤1
|pn(x)| = max−1≤x≤1
1
2n−1|Tn(x)| = 1
2n−1
Beweis Tn(x) = 2n−1x+. . . durch die Rekursionsformel fur Tschebyscheff’sche
Polynome. Mit x = cosϕ haben wir
Tn(x) = Tn(cosϕ) = cos(nϕ)
⇒ Tn hat (n + 1) Extremstellen x0, . . . , xn in [−1, 1] mit Tn(xj) = ±1 (al-
ternierendes Verzeichnis).
Sei Pn(x) ein Polynom n-ten Grades mit Hochstkoeffizient gleich Eins, so
dass
|Pn(x)| < 1
2n−1∀x ∈ [−1, 1].
Dann gilt
KAPITEL 8. ORTHOGONALE POLYNOME 74
Pn(x0) < Tn(x0)/2n−1 = +1
2n−1
Pn(x1) > Tn(x1)/2n−1 = − 1
2n−1
Pn(x2) < Tn(x2)/2n−1 = +1
2n−1
und so fort.
Folglich nimmt das Differenzpolynom
Φ(x) := Pn(x)− Tn(x)/2n−1
in den Extremalstellen x0 > x1 > . . . > xn Werte mit alternierenden Vorzei-
chen. Aus Stetigkeitsgrunden besitzt Φ (mindestens) n Nullstellen. Aber Φ hat
hochstens (n−1)-ten Grades, weil Pn und Tn/2n−1 Hochstterme xn haben. We-
gen des Hauptsatzes der Algebra ist dies ein Widerspruch.
⇒ ein solches Polynom Pn existiert nicht!
Seien x0, x1, . . ., xN paarweise verschiedene Stutzstellen in [−1, 1] und be-
trachte das Polynom
ϕ(x) = (x− x0) . . . (x− xN ).
Dieses Polynom hat Grad N +1 und Hochstkoeffizient gleich Eins. Von dem
obigen Satz (mit n = N + 1) haben wir
1
2N= max
−1≤x≤1
∣∣TN+1(x)/2N∣∣ ≤ max
−1≤x≤1|ϕ(x)|
⇒ mit den Nullstellen x(N+1)j , j = 0, 1, . . ., N , des Tschebyscheff Polynoms
TN+1(x) konnen wir den Term |(x− x0) . . . (x− xN )| in der Fehlerabschatzung
des Interpolationspolynoms minimieren
|f(x)− Ph
x(N+1)0 ,...,x
(N+1)N
i(x)| ≤ MN+1
(N + 1)!· 1
2N
fur x ∈ [−1, 1]
Bemerkung Fur ein Intervall [a, b] benutze die Transformation
x→ t =b− a
2· x +
b + a
2
KAPITEL 8. ORTHOGONALE POLYNOME 75
und die Stutzstellen
tj =b− a
2x
(N+1)j +
b + a
2
Beispiel Sei f(x) = ex fur −1 ≤ x ≤ 1.
Betrachte 2 Stutzstellen x0, x1 ∈ [−1, 1] und das Interpolationspolynom
P[x0,x1](x) = f[x0] + f[x0,x1] (x − x0) = xx0 +ex1 − ex0
x1 − x0(x− x0)
Das Tschebyscheff Polynom fur n = N + 1 = 2 (N = 1 hier) lautet
T2(x) = 2x2 − 1
mit Nullstellen x0 = −1/√
2 und x1 = +1/√
2.
Daher besitzt das Interpolationspolynom
e−1/√
2 +e1/
√2 − e−1/
√2
√2
(x + 1/√
2)
die kleinste Fehlerabschatzung aller Interpolationspolynome P[x0,x1](x) auf dem
Intervall [−1, 1], d.h.
|Fehler| ≤ M2
2!· 1
21≤ e1
2· 12
=e
4≃ 0.67
weil
f(x) = ex, f (2)(x) = ex, ⇒ max−1≤x≤1
≤ |f (2)(x)| ≤ e1
Bemerkung Das Polynom hier ist im Allgemeinen nicht dasselbe wie im
obigen least-squares Fall.
Siehe maple–Arbeitsblatt Tschebyscheff-Entwicklung
Kapitel 9
Numerische Differentiation
Literatur Schwarz, Kap. 3.2.2; Stummel/Hainer, Kap. 3.3
maple–Arbeitsblatt Differenzenquotienten: Diskretisierungsfehler versus Run-
dungsfehler
In der Analysis ist es oft leichter, die Ableitung einer Funktion herzuleiten
als das Integral der Funktion, z.B.
f(x) = e−x2
,d
dx
{e−x2
}= −2x e−x2
,
∫e−x2
dx =???
In der Numerik ist alles umgekehrt — obwohl es einfache Formeln fur nu-
merische Differentiation gibt, ist die numerische Differentiation im Allgemeinen
ein gefahrlicher Prozess wegen Abrundung (Ausloschung)
f ′(x0) = limh→0
f(x0 + h)− f(x0)
h≃ f(x0 + h)− f(x0)
h
d.h. Quotient von 2 kleinen Zahlen von denen eine aus einer Subtraktion
hervorgeht
f(x0 + h)− f(x0) ≃ 0, und h ≃ 0
Aber mehr spater — zunachst werden wir die Differenzenquotientenformel
herleiten. Dafur (insbesondere, um die entsprechende Fehlerabschatzung herzu-
leiten) benutzen wir eine Taylor-Entwicklung. Sei f 2-mal stetig differenzierbar:
f(x0 + h) = f(x0) + f ′(x0) h +1
2f ′′(ξh) h2 mit ξh ∈ [x0, x0 + h]
⇒ f ′(x0) =f(x0 + h)− f(x0)
h︸ ︷︷ ︸Differenzenquotient
− 1
2f ′′(ξh) · h︸ ︷︷ ︸
Fehler
Vorwartsdifferenzenquotient
76
KAPITEL 9. NUMERISCHE DIFFERENTIATION 77
D+h f(x0) : =
f(x0 + h)− f(x0)
h=
1
h∆f(x0) h > 0
mit Fehlerabschatzung
∣∣f ′(x0)−D+h f(x0)
∣∣ ≤ K · h
wobei
maxx∈I|f ′′(x)| .
= 2K.
fur ein geeignetes Intervall I.
Mit dem Landau-Symbol”großes Oh“ [siehe Oevel, Kap. 2], schreiben wir
D+h f(x0) = f ′(x0) + O(h)
d.h. eine Approximation erster Ordnung.
Ahnlicherweise gibt es einen Ruckwartsdifferenzenquotient
D−h f(x0) =
f(x0)− f(x0 − h)
h=
1
h∇f(x0) h > 0
fur welches
D−h f(x0) = f ′(x0) + O(h)
d.h., auch von erster Ordnung.
Konnen wir eine hohere Ordnung erreichen? JA!
Dafur verwenden wir einen der grundsatzlichen Tricks der Numerik.
Sei f 3-mal stetig differenzierbar. Betrachte die Taylor-Entwicklungen
f(x0 + h) = f(x0) + f ′(x0)h +1
2f ′′(x0)h
2 +1
3!f ′′(ξ+
h ) · h3
mit ξ+h ∈ [x0, x0 + h] und
f(x0 − h) = f(x0)− f ′(x0)h +1
2f ′′(x0)h
2 − 1
3!f ′′′(ξ−h ) · h3
mit ξ−h ∈ [x0 − h, x0].
Dann subtrahiere die zweite Gleichung von der ersten Gleichung :
f(x0 + h)− f(x0 − h) = 2f ′(x0) · h + 0.h2 +h3
6
[f ′′′(ξ+
h ) + f ′′′(ξ−h )]
KAPITEL 9. NUMERISCHE DIFFERENTIATION 78
Aber von dem Zwischenwertsatz (fur f ′′′) existiert ein
ξh ∈[ξ−h , ξ+
h
]
mit
f ′′′(ξh) =1
2
[f ′′′(ξ+
h ) + f ′′′(ξ−h )].
Deswegen haben wir
f(x0 + h)− f(x0 − h) = 2f ′(x0) h +h3
3f ′′′(ξh)
oder
f ′(x0) =f(x0 + h)− f(x0 − h)
2h︸ ︷︷ ︸Zentraldifferenzenquotient
−h2
3f ′′′(ξh).
Der Zentraldifferenzenquotient hier ist gleich
1
2
[1
h∆f(x0) +
1
h∇f(x0)
]
d.h. das arithmetisches Mittel der beiden Differenzenquotienten erster Ordnung.
Wir haben die Fehlerabschatzung∣∣∣∣f
′(x0)−f(x0 + h)− f(x0 − h)
2h
∣∣∣∣ ≤ K · h2
und schreiben daher
f(x0 + h)− f(x0 − h)
2h= f ′(x0) + 0(h2)
d.h. eine Approximation zweiter Ordnung.
Der Trick war, die beiden Terme zweiter Ordnung in den Taylor-Entwicklungen
aus zu loschen.
Fur hohere Ableitungen geht es ahnlich. z.B.
f(x0 ± h) = f(x0)± f ′(x0)h +1
2f ′′(x0)h
2 ± 1
3f ′′′(ξ±h )h3
Addiere!
f(x0 + h) + f(x0 − h) = 2f(x0) + 0 · h + f ′′(x0) · h2 +h3
3
[f ′′′(ξ+
h )− f ′′′(ξ−h )]
⇒ f ′′(x0) =f(x0 + h)− 2f(x0) + f(x0 − h)
h2− h
3
[f ′′′(ξ+
h )− f ′′′(ξ−h )]
KAPITEL 9. NUMERISCHE DIFFERENTIATION 79
eine Approximation erster Ordnung?
Aber f ′′′(ξ+h )− f ′′′(ξ−h ) konnte sehr klein sein. Tatsachlich hat diese Formel
eine bessere Ordnung, falls f 4-mal stetig differenzierbar ist. Dies folgt wegen
des obigen Tricks.
f(x0 ± h) = f(x0)± f ′(x0) h +1
2f ′′(x0) h2
± 1
3!f (3)(x0) h3 +
1
4!f (4)(ξ±h ) h4
Addiere: Die f ′ undf (3) Terme loschen sich aus. ⇒
f ′′(x0) =f(x0 + h)− 2f(x0) + f(x0 − h)
h2− 1
4!h2[f (4)(ξ+
h ) + f (4)(ξ−h )]
Zwischenwertsatz fur f (4) ⇒ ∃ξh ∈ [ξ−h , ξ+h ] mit
f (4)(ξh) =1
2[f (4)(ξ−h ) + f (4)(ξ+
h )].
Deshalb haben wir
f ′′(x0) =f(x0 + h)− 2f(x0) + f(x0 − h)
h2− h2
12f (4)(ξh)
Die Fehlerabschatzung hier lautet
∣∣∣f ′′(x0)−f(x0 + h)− 2f(x0) + f(x0 − h)
h2
∣∣∣ ≤ Kh2,
wobei maxx∈I|f (4)(x)| = 12K ist.
⇒ f(x0 + h)− 2f(x0) + f(x0 − h)
h2= f ′′(x0) + 0(h2)
d.h. eine Approximation zweiter Ordnung.
NBf(x0 + h)− 2f(x0) + f(x0 − h)
h2=
1
h2∆2f(x0 − h)
9.1 Rundung in numerischer Differentiation
Jetzt werden wir zeigen, warum Rundung in numerischer Differentiation gefahrlich
ist. Wir betrachten den Vorwartsdifferenzenquotient
D+h f(x0) =
f(x0 + h)− f(x0)
h= f ′(x0) + O(h)
oder |f ′(x0)−D+h f(x0)| ≤ Kh.
KAPITEL 9. NUMERISCHE DIFFERENTIATION 80
Wegen Rundung und anderer Art Fehler (z.B. Funktionenapproximation)
berechnen wir
f(x0) = f(x0) + e0, f(x0 + h) = f(x0 + h) + eh
mit |e0|, |eh| ≤ ε, und
D+h f(x0) =
f(x0 + h)− f(x0)
h
=[f(x0 + h)− f(x0)] + [eh − e0]
h= D+
h f(x0) +eh − e0
h
Der wesentliche Fehler hier lautet
∣∣∣f ′(x0)− D+h f(x0)
∣∣∣ ≤∣∣f ′(x0)−D+
h f(x0)∣∣+∣∣∣D+
h f(x0)− D+h f(x0)
∣∣∣
≤ Kh +|eh − e0|
h≤ Kh +
2ε
h
weil |e0|, |eh| ≤ ε ⇒ |eh − e0| ≤ |e0| + |eh| ≤ 2ε.
Definiere
F (h) = Kh +2ε
h
min F (h) = F (h∗) mit h∗ =
√2ε
K
Fur h≪ h∗ konnte alles schief gehen.
Aber
1). F (h) ist nicht der echte Fehler, nur eine Abschatzung von oben —”worst-
case“ (schlimmster Fall)
2). Gute Ergebnisse sind moglich, aber h soll nicht zu klein sein.
Beispiel (Siehe Folie)
f(x) = sinx, f ′(x) = cosx, x0 =π
3 · 2Siehe maple–Arbeitsblatt
”Differenzenquotienten: Diskretisierungsfehler ver-
sus Rundungsfehler“
KAPITEL 9. NUMERISCHE DIFFERENTIATION 81
0.02
0.03
0.04
0.05
0.06
0.07
0 0.02 0.04 0.06 0.08 0.1h
Rundung in der numerischen Differentiation
Kapitel 10
Numerische Integration
Literatur Oevel, Kap. 9.1; Schwarz, Kap. 8.3; Stummel/Hainer, Kap. 4.1
Betrachte das Integral
I =
b∫
a
f(x) dx
einer Funktion f : [a, b] → R, die mindestens stetig ist.
Sei F eine Stammfunktion von f , d.h. F ′ = f , dann gilt wegen des Haupt-
satzes der Integralrechnung:
I =
b∫
a
f(x) dx = F (b)− F (a)
Aber leider sind Stammfunktionen oft unbekannt. z.B. f(x) = e−x2
.
Wir mussen Integrale solcher Funktionen numerisch berechnen. Numerische
Integrationsformeln sind von Formeln fur die Flache einfacher geometrischer
Gestalten aufgebaut ⇒”Quadratur“
a bc
f(a)
f(b)
f(c)
Rechteck-Regel
b∫
a
f(x) dx ≃ (b− a) · f(ξ),
d.h. Lange × Hohe, wobei z.B.
ξ = a oder b (Randpunkte), oder
ξ = (a + b)/2 (Mittelpunkt)
82
KAPITEL 10. NUMERISCHE INTEGRATION 83
Diese Formel ist genau (d.h. ohne Fehler) fur f(x) ≡ Konstante (konstante
Funktionen)
Trapez-Regelb∫
a
f(x) dx ≃ (b− a)
(f(b) + f(a)
2
)
d.h. Lange × Hohenmittelwert — genau, falls f(x) = αx + β, ∀α, β (lineare
Funktionen)
Parabelregel (Simpsons-Regel)
Parabel: p(x) = A(x− a)2 + B(x − a) + C durch die Punkte
(a, f(a)),
(a + b
2, f
(a + b
2
)), (b, f(b))
Lose das lineare Gleichungssystem
⇒
f0 = f(a) = A · 02 + B · 0 + C
f1 = f(
a+b2
)= A ·
(b−a2
)2+ B ·
(b−a2
)+ C
f2 = f(b) = A(b − a)2 + B(b− a) + C
mit Losung
⇒
(b − a)2A = 2f2 − 4f1 + 2f0
(b − a) · B = 4f1 − f2 − 3f0
C = f0
Wir haben die Approximation
b∫
a
f(x) dx ≃b∫
a
p(x) dx
=1
3A(b − a)3 +
1
2B(b − a)2 + C(b− a)
=b− a
6
[2A(b− a)2 + 3B(b− a) + 6C
]
=b− a
6[f0 + 4f1 + f2]
d.h. die Approximationsformel
b∫
a
f(x) dx ≃ b− a
6
[f(a) + 4f
(a + b
2
)+ f(b)
]
KAPITEL 10. NUMERISCHE INTEGRATION 84
die genau ist, falls f quadratisches Polynom ist.
Im Allgemeinen leiten wir Formeln her, wie
b∫
a
f(x)dx ≃ (b− a)
N∑
j=0
cjf(xj) (∗)
mit Gewichten c0, c1, . . ., cN ∈ R und Stutztstellen x0, x1, . . ., xN ∈ [a, b], wo
a ≤ x0 < . . . < xN ≤ b
(die Stutztstellen x0 und xN mussen nicht immer Randpunkte sein).
Beispiele
1). Rechteckregeln N = 0, c0 = 1 und x0 = a, b oder (b + a)/2
2). Trapezregel N = 1, c0 = c1 = 12 und x0 = a, x1 = b.
3). Simpsons-Regel N = 2 mit
c0 =1
6, c1 =
4
6, c2 =
1
6, x0 = a, x1 =
a + b
2, x2 = b
Hier sind die Gewichte c0, . . ., cN ausgewahlt, so dass die Formel (∗) fur
Polynome bis zum N -ten Grad genau ist — z.B. fur die Polynome 1, x, . . ., xN .
⇒∫ b
a
xk dx ≡ (b − a)
N∑
j=0
cjxkj , k = 0, . . . , N
d.h.
N∑
j=0
cjxkj =
1
k + 1
(bk+1 − ak+1
b− a
)= dk, k = 0, . . . , N
oder in Matrix-Vektor-Form
1 1 1 . . . 1
x0 x1 x2 . . . xN
x20 x2
1 x22 . . . x2
N...
...
xN0 xN
1 xN2 . . . xN
N
c0
c1
· · ·cN
=
d0
d1...
dN
eine (N + 1)× (N + 1) Vandermond’sche Matrix mit Determinante
det[xi
j
]=∏
i<j
(xi − xj)
KAPITEL 10. NUMERISCHE INTEGRATION 85
⇒ invertierbar, im Prinzip eindeutig losbar,aber es gibt bessere Methoden!
Bemerkungen
1). k = 0 ⇒ x0j ≡ 1 ⇒
N∑j=0
cj = 1
2). In Kapitel 5 haben wir bereits Interpolationspolynome besprochen — fur
Interpolationspolynome gibt es Fehlerabschatzungen, die uns hier sehr nutzlich
sind.
10.1 Interpolatorische Quadraturformeln
(Vom Benutzer) vorgegebene Stutzstellen a ≤ x0 < . . . < xN ≤ b:
b∫
a
f(x)dx ≃ (b− a)
N∑
j=0
cjf(xj)
Solche Quadraturformeln heißen Newton-Cotes-Formeln
abgeschlossen falls x0 = a, xN = b
offen sonst
Die Lagrange’sche Darstellung des Interpolationspolynoms ist gunstiger hier,
weil sie die f(xj)-Werte explizit enthalt.
Sei f : [a, b] → R eine (N + 1)-mal stetig differenzierbare Funktion mit
f(x) = P[x0,...,xN ](x) + Fehler
d.h.
f(x) =
N∑
j=0
f(xj) · LN,j(x)
︸ ︷︷ ︸P[x0,...,xN ](x)
+1
(N + 1)!f (N+1)(ξ) ·
N∏
j=0
(x− xj)
︸ ︷︷ ︸Fehler RN (x)
⇒ RN (x) ≡ 0 fur f Polynom hochstens N -ten Grades.
Integriere:
b∫
a
f(x) dx =
b∫
a
P[x0,...,x](x) dx +
b∫
a
RN (x) dx
= (b− a) ·N∑
j=0
1
b− a
b∫
a
LN,j(x) dx
· f(xj) +
b∫
a
RN (x) dx
KAPITEL 10. NUMERISCHE INTEGRATION 86
⇒ Die geeigneten Gewichte lauten
cj =1
b− a
b∫
a
LN,j(x) dx j = 0, 1, . . . , N
Die Fehlerabschatzung folgt durch:
EN (f) :=
∣∣∣∣∣∣
b∫
a
f(x) dx −b∫
a
Px0,...,xN(x) dx
∣∣∣∣∣∣
=
∣∣∣∣∣∣
b∫
a
RN (x) dx
∣∣∣∣∣∣
=
∣∣∣∣∣∣
b∫
a
1
(N + 1)!f (N+1)(ξ)
N∏
j=0
(x − xj) dx
∣∣∣∣∣∣
≤ 1
(N + 1)!max
a≤ξ≤b
∣∣∣f (N+1)(ξ)∣∣∣ ·
b∫
a
N∏
j=0
|x− xj | dx.
Betrachte die Transformation:
x ∈ [a, b]→ t =x− a
b− a∈ [0, 1].
und definiere tj =xj − a
b− afur j = 0, 1, . . . , N
⇒ x− xj = (b − a)(t− tj), dx = (b − a) · dt
⇒N∏
j=0
(x− xj) = (b− a)N+1∏
(t− tj)
⇒b∫a
N∏j=0
|x− xj |dx = (b− a)N+21∫0
N∏j=0
|t− tj |dt
Wir haben daraus die Fehlerabschatzung
EN (f) ≤ (b − a)N+2
(N + 1)!max
a≤ξ≤b
∣∣∣f (N+1)(ξ)∣∣∣ ·
1∫
0
N∏
j=0
|t− tj | dt
KAPITEL 10. NUMERISCHE INTEGRATION 87
Beispiele
1). Rechteckregel (N = 0)
E0(f) ≤ (b− a)2 · maxa≤ξ≤b
|f ′(ξ)| ·1∫
0
|t− t0| dt
wobei
1∫
0
|t− t0| dt =
t0∫
0
−(t− t0) dt +
1∫
t0
(t− t0) dt
= −1
2(t− t0)
2∣∣∣t0
0+
1
2(t− t0)
2∣∣∣1
t0
=1
2t20 +
1
2(1− t0)
2
⇒ min∫· dt = 1/4 fur t0 = 1/2 und fur die Randpunkte t0 = 0 und 1
gilt∫· dt = 1/2.
2). Trapezregel (N = 1) t0 = 0, t1 = 1 hier
E1(f) ≤ (b− a)3
2!max
a≤ε≤b|f ′′(ε)| ·
1∫
0
|(t− 0)(t− 1)|dt
︸ ︷︷ ︸1∫
0
t(1− t)dt = 1/6
⇒ E1(f) ≤ (b− a)2
12max
a≤ξ≤b|f ′′(ξ)|
3). Simpsons-Regel (N = 2) t0 = 0, t1 = 1/2, t2 = 1 hier
E2(f) ≤ (b − a)4
3!max
a≤ξ≤b|f (3)(ξ)| ·
1∫
0
|(t− 0)(t− 1/2)(t− 1)|dt
wobei
1∫
0
|t(t− 1/2)(t− 1)|dt =
1/2∫
0
t(t− 1/2)(t− 1)dt +
1∫
1/2
t(t− 1/2)(1− t)dt
KAPITEL 10. NUMERISCHE INTEGRATION 88
=1
26+
1
26=
1
25=
1
32
⇒ E2(f) ≤ (b− a)4
192· max
a≤ξ≤b|f (3)(ξ)|
Aber diese Abschatzungen sind oft zu grob und scharfere sind moglich, ins-
besondere falls N gerade ist.
Satz (Schwarz, Satz 8.4, Seite 398)
Sei N gerade, f (N + 2)-mal stetig differenzierbar, und
(b− a)N∑
j=0
cjf(xj)
eine geschlossene Newton-Cotes-Formel auf [a, b]
Dann ist der Quadraturfehler gegeben durch
∣∣∣∣∣∣
b∫
a
f(x)dx− (b− a) ·N∑
j=0
cjf(xj)
∣∣∣∣∣∣
=
∣∣∣∣∣∣(b− a)N+3
(N + 2)!f (N+2)(ξ) ·
1∫
0
t
N∏
j=0
(t− tj)dt
∣∣∣∣∣∣
wobei a < ξ < b.
D.h. wir haben eine zusatzliche Potenzordnung gewonnen.
Beispiel Simpsons-Regel (N = 2) t0 = 0, t1 = 1/2, t2 = 1
1∫
0
t
N∏
j=0
(t− tj)dt =
1∫
0
t2(t− 1/2)(t− 1)dt
=1
2
1∫
0
t2(2t− 1)(t− 1)dt =1
120
⇒ E2(f) ≤ (b− a)5
2880· max
a≤ξ≤b|f (4)(ξ)|
statt
E2(f) ≤ (b − a)4
192max
a≤ξ≤b|f (3)(ξ)|
KAPITEL 10. NUMERISCHE INTEGRATION 89
Viel besser, falls b− a < 1 ist!
⇒ Wir sollen nur Formeln mit N gerade benutzen,
Aber oft haben wir b− a > 1. Was dann?
⇒ summierte Formeln! — spater
10.2 Rundung in numerischer Integration
Die Koeffizienten c0, . . ., cn einer Newton-Cotes-Formel
(b− a)
N∑
j=0
cj · f(xj)
eines Integralsb∫a
f(x)dx genugen der Eigenschaft
N∑
j=0
cj = 1
und, falls 0 ≤ N ≤ 6 (Hausarbeit!), auch
cj ≥ 0, j = 0, . . . , N.
(Tatsachlich sollen wir c(N)j schreiben, weil die Koeffizienten von N abhangen .)
Wir beschranken uns zum Fall: 0 ≤ N ≤ 6. Wegen Rundung und anderer
Fehlerarten (z.B. Funktionenapproximation) berechnen wir
fj statt f(xj)
mit dem Fehler ej = fj−f(xj), wobei |ej | ≤ ε≪ 1 ist. Deshalb ist der wesentliche
Fehler∣∣∣∣∣∣
b∫
a
f(x)dx− (b − a)N∑
j=0
cj fj
∣∣∣∣∣∣
=
∣∣∣∣∣∣
b∫
a
f(x)dx − (b− a)N∑
j=0
cj{f(xj) + ej}
∣∣∣∣∣∣
≤
∣∣∣∣∣∣
b∫
a
f(x)dx − (b− a)
N∑
j=0
cjf(xj)
∣∣∣∣∣∣+ (b− a)
N∑
j=0
|cj ej|
KAPITEL 10. NUMERISCHE INTEGRATION 90
≤ EN (f) + (b − a)
N∑
j=0
cj |ej | , cj > 0 hier
≤ EN (f) + (b − a)
N∑
j=0
cjε = EN (f) + (b− a) ε
d.h., der zusatzliche numerische Fehler (b − a) ε hangt nicht von N ab.
Nur ein Sonderfall, aber er zeigt, dass numerische Integration robust oder
glattend ist. Vgl. numerische Differentiation.
Kapitel 11
Summierte
Integrationsformeln
Literatur Oevel, Kap. 9.1–2; Schwarz, Kap. 8.1; Stummel/Hainer, Kap. 4.2.1
maple–Arbeitsblatt Adaptive numerische Integration
Der Fehler einer Newton-Cotes-Formel genugt einer Abschatzung der Form:
|EN (f)| :=
∣∣∣∣∣∣
b∫
a
f(x) dx− (b− a)
N∑
j=0
cj · f(xj
∣∣∣∣∣∣≤ KN(b − a)N+k (k = 2 oder 3)
wobei KN eine Konstante ist.
Im Allgemeinen betrachten wir aquidistante Stutzstellen: xj+1 = xj + h, j
= 0, . . ., N − 1 z.B. fur eine geschlossene Formel haben wir
xj = a + jh, h =b − a
N.
Dann haben wir
(b − a)N+k = hN+k ·NN+k
und die obige Fehlerabschatzung lautet
|EN (f)| ≤ KNNN+k · hN+k
z.B. Simpsons-Regel N = 2, k = 3 mit h = b−a2
E2(f) = − (b− a)5
2880f (4)(ξ) = −h5
90f (4)(ξ).
OK falls h = b−a2 < 1. Sonst?
91
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 92
Wir konnten immer N groß genug wahlen, so dass
h =b− a
N< 1
Aber
1). die Newton-Cotes-Formel ist dann sehr kompliziert
2). KN · hN+k konnte sehr groß sein — erinnere: die Interpolationspoly-
nome fur N ≫ 1 und aquidistante Stutzstellen haben sehr große Schwankungen.
Besser Benutze eine einfachere Formel (N -klein) auf Teilintervallen und
summiere z.B. schreibe
[a, b] =
N−1⋃
j=0
[xj , xj+1]
und benutze die Trapezregel auf dem Teilintervall [xj , xj+1] fur j = 0, . . . , N−1,
d.h.
xj+1∫
xj
f(x) dx ≃ (xj+1 − xj)
[f(xj+1) + f(xj)
2
]
dann summiere:
b∫
a
f(x) dx =
N−1∑
j=0
xj+1∫
xj
f(x) dx ≃N−1∑
j=0
(xj+1 − xj)
[1
2f(xj) +
1
2f(xj+1)
]
Alles ist einfacher mit aquidistanten Stutzstellen:
h =b− a
N, xj = a + jh, xj+1 − xj = h
Dann lautetxj+1∫
xj
f(x) dx ≃ h
2[f(xj) + f(xj+1)]
und daher gilt
b∫
a
f(x) dx =
N−1∑
j=0
xj+1∫
xj
f(x) dx
≃ h
2
N−1∑
j=0
[f(xj) + f(xj+1)]
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 93
=h
2[f(x0) + 2f(x1) + . . . + 2f(xN−1)︸ ︷︷ ︸+f(xN)]
rechte/ dann linke Randpunkt sukzessiver Teilintervalle
d.h. wir erhalten die summierte Trapez-Regel
b∫
a
f(x) dx ≃ h
2
f(x0) + f(xN ) + 2
N−1∑
j=1
f(xj)
Fur die Simpsons-Regel benutzen wir N = 2M (gerade) Teilintervalle und
betrachten jedes Mal ein Paar sukzessiver Teilintervalle [x2j , x2j+1] und [x2j+1, x2j+2]
fur j = 0, 1, . . ., M − 1. Es gilt
x2j+2∫
x2j
f(x) dx ≃ x2j+2 − x2j
6[f(x2j) + 4f(x2j+1) + f(xj+2)]
=2h
6[f(x2j) + 4f(x2j+1) + f(x2j+2)]
und dann:
b∫
a
f(x) dx =
M−1∑
j=0
x2j+2∫
x2j
f(x)dx
≃ h
3
M−1∑
j=0
[f(x2j) + 4f(x2j+1) + f(x2j+2)]
=h
3
f(x0) + f(x2M ) + 4
M−1∑
j=0
f(x2j+1) + 2
M−2∑
j=0
f(x2j+2)
d.h. wir erhalten die summierte Simpsons-Regel
b∫
a
f(x) dx ≃ h
3
f(x0) + f(x2M ) + 4
M−1∑
j=0
f(x2j+1) + 2
M−2∑
j=0
f(x2j+2)
wobei die letzte Summe die internen Randpunkte enthalt.
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 94
Wir brauchen Fehlerabschatzungen fur diese summierten Integrationsfor-
meln.
Trapez-Regel
xj+1∫
xj
f(x) dx =h
2[f(xj) + f(xj+1)]−
h3
12f ′′(ξj), ξj ∈ [xj , xj+1]
⇒b∫
a
f(x) dx =h
2
f(x0) + f(xN ) +
N−1∑
j=1
f(xj)
− h3
12
N−1∑
j=0
f ′′(ξj)
f ′′ stetig: Zwischenwertsatz ⇒ es existiert ξ∗ ∈ [a, b] mit
1
N
N−1∑
j=0
f ′′(ξj)
︸ ︷︷ ︸
= f ′′(ξ∗)
arithmetisches Mittel liegt zwischen min und max
⇒ summierter Fehler = −h3N
12· f ′′(ξ∗) = − h2
12(b − a)f ′′(ξ∗)
︸ ︷︷ ︸
eine Potenz weniger!
weil b− a = Nh.
Simpsons-Regel
summierter Fehler = −h5
90
M−1∑
j=0
f (4)(ξj) = −h5
90·Mf (4)(ξ∗)
f (4) stetig: Zwischenwertsatz ⇒ es existiert ξ∗ ∈ [a, b] mit
1
M
M−1∑
j=0
f (4)(ξj) = f (4)(ξ∗)
⇒ summierter Fehler = − h4
180· 2Mhf (4)(ξ∗) = − h4
180· (b− a)f (4)(ξ∗)
weil h =b− a
2Mist.
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 95
In beiden Fallen haben wir eine Potzenz verloren, aber jetzt haben wir hk
statt (b− a)k+1 mit h < 1, und es ist nicht wichtig, ob b− a > 1 hier ist.
Adaptive Integration: Wie oben, aber wir wahlen die Lange des nachsten
Teilintervalls im Laufe der Berechnung, um zu versichern, dass der Fehler auf
diesem Teilintervall nicht zu groß ist — nutzlich,
z.B. falls f sich sehr schnell steil verandert.
Aber wie konnen wir diese Idee automatisieren?
Siehe maple–Arbeitsblatt”Adaptive numerische Integration“
11.1 Extrapolation
Literatur Schwarz, Kap. 8.1; Stummel/Hainer, Kap. 4.2.2
Ein grundsatzlicher Trick der Numerik heißt Extrapolation.
Sei h > 0 ein Parameter (z.B. Schrittweite) und sei T (h) eine Approximation
fur T (0). Wir setzen voraus, dass T (h) der folgenden Fehlerentwicklung genugt:
T (h) = T (0) + K2 · h2 + 0(h4) (∗)
wobei K2 eine Konstante ist, die unabhangig von h ist.
Mit h/2 statt h erhalten wir
T (h/2) = T (0) + K2h2
4+ 0(h4)
weil 0((h2 )4) = 0( 1
16h4) = O(h4) — das O-Symbol hier absorbiert die Konstan-
ten.
Mulipliere mal 4
4T (h/2) = 4T (0) + K2 h2 + 0(h4) (∗∗)
Subtraktion der Gleichung (∗∗) von (∗) oben ergibt
4T (h/2)− T (h) = 3T (0) + 0.h2︸︷︷︸+0(h4)
Ausloschung!
und dann gilt:
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 96
4T (h/2)− T (h)
3︸ ︷︷ ︸= T (0) + 0(h4)
Approximation vierter Ordnung
Extrapolation: Die Herleitung einer Approximation hoherer Ordnung aus
Approximationen niedriger Ordnung durch Ausloschung von Fehlertermen.
Wir konnen diesen Trick wiederholen, falls wir die Fehlerentwicklung erwei-
tern konnen: z.B. sei
T (h) = T (0) + K2h2 + K4h
4 + 0(h6) (∗ ∗ ∗)
mit Konstanten K2 und K4. Mit h/2 hatt h haben wir
T (h/2) = T (0) + K2 ·h2
4+ K4
h4
16+ 0(h6) (∗ ∗ ∗∗)
Vier mal Gleichung (∗ ∗ ∗∗) minus Gleichung (∗ ∗ ∗) ergibt
4T (h/2)− T (h) = 3T (0)− 3
4K4h
4 + 0(h6)
oder
4T (h/2)− T (h)
3= T (0)− 1
4K4h
4 + 0(h6)
Definiere
T0(h) = T (h), T1(h/2) =4T0(h/2)− T0(h)
3.
Wir haben
T1(h/2) = T (0)− 1
4K4h
4 + 0(h6) (5∗)
und dann mit 12 · h/2 = h/4 statt h/2
T1(h/4) = T (0)− 1
4K4 ·
h4
16+ 0(h6) (6∗)
Von 16× (6∗) − 1× (5∗) erhaltenen wir:
16T1(h/4)− T1(h/2) = 15T (0) + 0(h6)
oder
T2(h/4) = T (0) + 0(h6)
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 97
wobei
T2(h/4) :=16T1(h/4)− T1(h/2)
15=
42T1(h/22)− T1(h/2)
42 − 1.
Im Allgemeinen: fangen wir an mit
T (h) = T (0) + K2 h2 + . . . + K2N h2N + 0(h2N+2)
dann konnen wir alles wie oben bis zum n = N wiederholen. Zu jedem Schritt
erhalten wir
Tn(h/2n) = T (0) + 0(h2n+2),
eine Approximation (2n + 2)-ter Ordnung, wobei
Tn
(h
2n
)=
4nTn−1
(h
2n
)− Tn−1
(h
2n−1
)
4n − 1
n = 1, . . . , N
mit T0(h) = T (h).
11.2 Romberg-Integration
Literatur Schwarz, Kap. 8.1; Stummel/Hainer, Kap. 4.2.2
Romberg-Integration baut auf die Trapez-Regel mit sukzessiver Halbierung
des Abstandes zwischen den Stutzstellen. Die Extrapolations-Methode liegt da-
hinter.
Wir werden diese Methode im Zusammenhang des Integrals
T (0) =
b∫
a
f(x)dx
und der summierten Trapezregel
T (h) =h
2
f(x0) + f(xN ) + 2
N−1∑
j=1
f(xj)
.
verwenden (h = (b− a)/N hier).
Dafur haben wir den Fehler
T (h) = T (0) +h2
12f (2)(ξ),
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 98
falls f 2-mal stetig differenzierbar ist.
Aber fur f (2N +2)-mal stetig differenzierbar durch die Euler-Maclaurinsche
Summenformel erhalten wir
T (h) = T (0) + K2h2 + Kn h4 + . . . + K2N h2N + 0(h2N+2)
mit Konstanten
K2j :=B2j
(2j)!
[f (2j−1)(b)− f (2j−1)(a)
]
wobei die Bj Bernoulli-Zahlen sind, d.h. die Koeffizienten in der Potenzreihe
x
ex − 1=
∞∑
j=0
1
j!Bj xj
z.B. B0 = 1, B1 = 0, im Allgemeinen B2j+1 = 0, B2j 6= 0 [siehe Oevel: Seite
433-437].
Wir haben tatsachlich
T (h) =h
2
f(a) + f(b) + 2
N−1∑
j=1
f(a + jh)
mit xj = a + jh und h = (b− a)/N
⇒ T
(h
2
)=
h/2
2
f(a) + f(b) + 22N−1∑
j=1
f
(a + j
h
2
)
summiere gerade und ungerade Indizes getrennt
=h
4[f(a) + f(b)] +
h
42
N−1∑
j=1
f(a + jh) +h
42
N∑
j=1
f
(a + (2j − 1)
h
2
)
d.h.
T
(h
2
)=
1
2T (h) +
h
2
N∑
j=1
f
(a + (2j − 1)
h
2
)
d.h. wir mussen f nur fur die neuen, d.h. ungeraden Stutzstellen berechnen.
Nach dem ersten Schritt der Extrapolations-Methode (T0(h) = T (h) hier)
erhalten wir
T1
(h
2
)=
4T0(h2 )− T0(h)
3
=1
3
T0(h) + 2h
N∑
j=1
f
(a + (2j − 1)
h
2
)
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 99
=1
3· h2
f(a) + f(b) + 2
N−1∑
j=1
f
(a + 2j
h
2
)+ 4
N∑
j=1
f
(a + (2j − 1)
h
2
)
d.h. die (summierte) Simpsons-Regel !
Weitere Extrapolationsschritte ergeben auch bekannte Newton-Cotes-Formeln
hoherer Ordnung, z.B. T2(h/4), . . .
Romberg-Integration fangt ganz einfach an, d.h. mit N = 1, h = b − a und
der einfachen Trapez-Regel.
Definiere T(k)n durch
T(k)0 = T0
(b− a
2k
)Trapez-Regel mit Schrittweite hk =
b − a
2k
und
T (k)n = Tn
(b− a
2k
)n-ter Extrapolationsschritt
fur k = 0, 1, 2, . . . und n = 0, 1, . . ..
Von oben haben wir
1). Trapezregel-Schritthalbierung hk+1 =1
2hk
T(k+1)0 =
1
2T
(k)0 + hk+1
2k∑
j=1
f(a + (2j − 1)hk+1)
2). Extrapolation
T (k+1)n =
4nT(k+1)n−1 − T
(k)n−1
4n − 1= T
(k)n−1 +
4n
4n − 1
[T
(k+1)n−1 − T
(k)n−1
]
Wir bauen systematisch eine Tabelle auf.
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 100
h0 : T(0)0
ց
h1 : T(1)0 −→ T
(1)1
ց ց
h2 : T(2)0 −→ T
(2)1 −→ T
(2)2
ց ց ց
h3 : T(3)0 −→ T
(3)1 −→ T
(3)2 −→ T
(3)3
......
......
Trapez-R. Simpson-R.
O(h2k) O(h4
k) O(h6k) O(h8
k)
Insbesondere genugen die diagonalen Terme
T (n)n =
b∫
n
f(x) dx + O(h2n+2n )
Romberg ⇒ benutze T(n)n als die erwunschte Approximation mit n, so dass
∣∣∣T (n+1)n+1 − T (n)
n
∣∣∣ < ε ←− erwunschte Genauigkeit
Vorteile: einfach, rekursiv — musse das engultige n nicht vom Anfang
wahlen — schnell konvergierend
Nachteil: 2n + 1 Funktionenauswertungen — steigt exponentiell mit n auf.
Romberg-Integration ist die meistbenutzte Quadraturmethode in Software-
Bibliotheken.
Variation: wir mussen h nicht nur halbieren.
Beispiel 1).: Betrachte das Integral∫ 1
0cos(
πx2
)dx — die Integrandfunktion
is in [0, 1] beliebig oft stetig differenzierbar
KAPITEL 11. SUMMIERTE INTEGRATIONSFORMELN 101
n T(n)0 T
(n)1 T
(n)2 T
(n)3
0 0.50000000
1 0.60355339 0.63807119
2 0.62841744 0.63670546 0.63661441
3 0.63457315 0.63662505 0.63661969 0.63661977
der wahre Wert:∫ 1
0 cos(
πx2
)dx = π
2 = 0.63661977
Beispiel 2).: Betrachte das Integral∫ 1
0 x3/2 dx
n T(n)0 T
(n)1 T
(n)2 T
(n)2 T
(n)3 T
(n)4
0 0.50000000
1 0.426777 0.402369
2 0.407018 0.400432 0.400303
3 0.401812 0.400077 0.400054 0.400050
4 0400463 0.400014 0.400009 0.400009 0.400009
5 0.400118 0.400002 0.400002 0.400002 0.400002 0.400002
der wahre Wert:∫ 1
0x3/2 dx = 2
5 = 0.40000000
Obwohl die Integrandfunktion in [0, 1] nur einmal stetig differenzierbar ist,
lauft das Romberg-Verfahren ziemlich gut, nur ein bißchen langsam!
Kapitel 12
Gauß-Quadratur
Literatur Oevel, Kap. 9.4; Schwarz, Kap. 8.3.2; Stummel/Hainer, Kap. 4.3
Betrachte eine Newton-Cotes-Formel
(b− a)N∑
j=0
cj · f(xj)
fur ein Integralb∫
a
f(x) dx, wobei die Stutzstellen x0, . . . , xN ∈ [a, b] vorgegeben
sind. Die entsprechenden Koeffizienten genugen
cj =1
b− a
b∫
a
LN,j(x)dx,
wobei
LN,j =∏
i6=j
(x− xi))
(xj − xifur j = 0, 1, . . . , N.
Diese Formel ist genau fur Polynome bis zum N -ten Grad ⇒ cj auch dem
linearen Gleichungssystem genugen:
1 1 . . . 1
x0 x1 . . . xN
...
xN0 xN
1 . . . xNN
c0
c1
...
cN
=
d0
d1
...
dN
wobei
dk =1
b− a
b∫
a
xkdx =bk+1 − ak+1
(b − a)(k + 1).
.
102
KAPITEL 12. GAUSS-QUADRATUR 103
Fur aquidistante Stutzstellen mit Abstand xj+1 − xj = h > 0 gilt es
|Fehler| = 0(hN+k
)k = 2 oder 3
Frage: Konnen wir die Genauigkeit verbessern durch eine”optimale“ Wahl er
Stutzstellen?
Dann werden wir 2N + 2 Unbekannten c0, . . ., cN und x0, . . ., xN haben.
Dafur brauchen wir 2N+2 Gleichungen, die wir finden konnen, mit dem Ansatz,
dass die Formel
(b− a)N∑
j=0
cj f(xi)
genau ist fur die Polynome 1, x, . . ., x2N+1, d.h.
N∑
j=0
cjxkj =
bk+1 − ak+1
(b − a)(k + 1)fur k = 0, 1, . . . , 2N + 1
⇒
c0 + c1 + . . . + cN = 1
c0x0 + c1x1 + . . . + cNxN =b + a
2
......
......
c0x2N+10 + c1x
2N+11 + . . . + cNx2N+1
N =b2N+2 − a2N+2
(b− a)(2N + 2)
⇒ ein nichtlineares Gleichungssystem
Beispiel 1). N = 0 ⇒ 2 Unbekannten und 2 Gleichungen ⇒
c0 = 1, c0x0 =b + a
2⇒ c0 = 1, x0 =
b + a
2
Mittelpunk-Rechteckregel
b∫
a
f(x)dx ≃ (b− a) · f(
b + a
2
)
Beispiel 2). N = 1 ⇒ 4 Unbekannten und 4 Gleichungen
c0 + c1 = 1
c0x0 + c1x1 =b + a
2
c0x20 + c1x
21 =
1
3
b3 − a3
b + a
c0x30 + c1x
31 =
1
4
b4 − a4
b + a
KAPITEL 12. GAUSS-QUADRATUR 104
⇒ offensichtlich!
c0 =1
2, x0 = −b− a
2· 1√
3+
b + a
2
c1 =1
2, x1 = +
b− a
2· 1√
3+
b + a
2
Gauß-Legendre-Regel
b∫
a
f(x)dx ≃ (b− a) ·[1
2f
(−b− a
2· 1√
3+
b + a
2
)+
1
2f
(b− a
2· 1√
3+
b + a
2
)]
Es geht weiter: N = 2, 3, usw
Erstaunlich! Das Gleichungssystem ist immer losbar fur
(1) N + 1 reelle x0, . . . , xN ∈ (a, b)
(2) N + 1 positive c0, . . . , cN
Der Beweis ist leichter und allgemeiner mit geeigneter orthogonalen Polyno-
men ϕ0(x), . . ., ϕ2N+1(x) statt den einfachen Polynomen 1, x, . . ., x2N+1.
Um unendliche Intervalle wie [0,∞) und (−∞,∞) sowie [a, b] zu betrachten,
werden wir ein Integral mit einer Gewichtsfunktion w(x) untersuchen und die
Quadratur-Formeln etwas anders aufschreiben, namlich
b∫
a
w(x)f(x) dx ≃N∑
j=1
αj f(xj)
d.h. mit einer Gewichtfunktion w(x) ≥ 0; j = 1,. . . , N statt j = 0,1, . . . , N und
αj statt (b− a)cj
Satz (Siehe Oevel, Satz 9.8, Seite 451).
Die Gauß-Quadratur-Formel
b∫
a
w(x)f(x) dx ≃N∑
j=1
αj f(xj)
ist genau fur alle Polynome bis zum (2N − 1)-ten Grad, falls
1. Die Stutzstellen x1, . . ., xN die Nullstellen des orthogonalen Polynoms
φN (x) bzg. des Skalarprodukts < ., . >w sind, und
KAPITEL 12. GAUSS-QUADRATUR 105
2. die Gewichte α1, . . ., αN durch
αj =
b∫
a
w(x)
N∏
i=1i6=j
(x− xi)
(xj − xi), j = 1, . . . , N
definiert sind.
Erinnere, dass
LN−1,j(x) =
N∏
i=1i6=j
(x− xi)
(xj − xi), j = 1, . . . , N
(wir haben N − 1 hier statt N , weil j fangt mit 1 statt 0 an)
⇒ Definition wie im Newton-Cotes-Fall fur die Funktion w(x)f(x) statt
f(x) und die Koeffizienten αj statt (b− a)cj .
Beweis Sei Pk der lineare Raum aller Polynome hochstens k-ten Grades.
1). Sei f ∈ PN−1. Dann gilt es
f(x) =
N∑
j=1
f(xj)LN−1,j(x)
und deshalb
b∫
a
w(x)f(x) dx =
b∫
a
w(x)N∑
j=1
f(xj)LN−1,j(x) dx
=
N∑
j=1
f(xj)
b∫
a
w(x)LN−1,j(x) dx
︸ ︷︷ ︸αj
⇒ die Formel ist genau
2). Sei f ∈ P2N−1 \ PN−1, d.h. f ist Polynom mit Grad ∈ (N − 1, 2N − 1].
Dann gilt es
f(x) = g(x)ϕN (x) + h(x)
wobei g(x) und h(x) Polynome hochsten (N − 1)-ten Grades sind. Dann gilt es
b∫
a
w(x)f(x) dx −N∑
j=1
αjf(xj) =
KAPITEL 12. GAUSS-QUADRATUR 106
=
b∫
a
w(x)g(x)ϕN (x) dx +
b∫
a
w(x)h(x) dx −N∑
j=1
αjg(xj)ϕN (xj)−N∑
j=1
αjh(xj)
= 0 +
b∫
a
w(x)h(x) dx − 0−N∑
j=1
αjh(xj) = 0
weil die Formel genau fur Polynome in PN−1 z.B. fur h :
b∫
a
w(x)h(x) dx =
N∑
j=1
αjh(xj),
undb∫
a
w(x)g(x)ϕN (x) dx = 〈g, ϕN 〉w = 0,
weil g ∈ PN−1 ⇒ g⊥ϕN , sowie
N∑
j−1
αjg(xj)ϕN (xj) = 0
weil xi,. . .,xN Nullstellen von ϕN sind.
Beweis der Bemerkung, dass die Koeffizienten α1, . . . , αN positiv sind.
Betrachte die Funktion
f(x) = (LN−1,k(x))2 (k fest)
d.h. f ∈ P2N−2, weil LN−1,j ∈ PN−1
⇒ 0 <
b∫
a
w(x)(LN−1,k(x))2 dx =N∑
j=1
αj(LN−1,k(xj))2 = αk
weil LN−1,k(xj) =
{1 falls j = k
0 sonst
Auch gilt es
N∑
j=1
αj =
b∫
a
w(x) dx.
Beweis: Nimm f(x) ≡ 1!
Es gibt auch einen expliziten Ausdruck fur den Fehler (siehe z.B. Oevel, Seite
458)b∫
a
w(x)f(x) dx −N∑
j=1
αj · f(xj) =< ϕN , ϕN >w
(2N)!· f (2N)(ξ)
KAPITEL 12. GAUSS-QUADRATUR 107
mit ξ ∈ [a, b].
Deshalb fur f ∈ P2N−1 haben wir f (2N)(x) ≡ 0 und die Formel ist exakt.
Beispiel 1). Betrachte das Integralb∫
a
f(x) dx
Transformiere x ∈ [a, b] ↔ t ∈ [−1, 1] mit
x =b− a
2· t +
b + a
2
Dann gilt
b∫
a
f(x) dx =b− a
2
1∫
−1
f
(b− a
2t +
b + a
2
)dt =
b− a
2
1∫
−1
f(t) dt
wobei
f(t) ≡ f
(b− a
2t +
b + a
2
).
Wir haben Integrale wie1∫
−1
g(t)dt, d.h. mit Gewichtfunktion w(t) ≡ 1. Die
entsprechenden orthogonalen Polynome sind die Legendre’schen Polynome.
N = 1 ϕ1(t) = t auf dem Intervall [−1, 1] mit
Nullstelle t1 = t(1)1 = 0
Koeffizient α1 =1∑
j=1
αj =1∫
−1
w(t) dt =1∫
−1
1 dt = 2
⇒1∫
−1
f(t) dt ≃ 2 · f(0)
Im Allgemeinen: Mittelpunkt-Rechteckregel
b∫
a
f(x) dx ≃ (b − a) f
(b + a
2
)
N = 2 ϕ2(t) = t2 − 1/3
Nullstellen t(2)1 = −1/
√3 und t
(2)2 = 1/
√3
KAPITEL 12. GAUSS-QUADRATUR 108
Koeffizienten
α1 =
∫ 1
−1
1 · t− 1/√
3
−1/√
3− 1/√
3dt
= −√
3
2
∫ 1
−1
(t− 1/√
3) dt
= −√
3
2
(1
2t2 − t√
3
) ∣∣∣1
−1=−√
3
2· −2√
3= +1
ahnlicherweise α2 = +1 (=∫ 1
−1dt − α1!)
⇒1∫
−1
f(t) dt ≃ f
(− 1√
3
)+ f
(1√3
)
und dann in den ursprunglichen Koordinaten
b∫
a
f(x)dx ≃ (b− a) ·[1
2f
(−b− a
2· 1√
3+
b + a
2
)+
1
2f
(b− a
2· 1√
3+
n + a
2
)]
Dieser Fall heißt”Gauß-Legendre“-Quadratur.
Beispiel 2). Betrachte das Integral∞∫0
f(x) dx
∞∫
0
f(x) dx =
∞∫
0
e−x [exf(x)] dx =
∞∫
0
e−xF (x) dx
wobei F (x) ≡ exf(x).
Hier lautet das Innerprodukt
〈F, G〉w =
∞∫
0
e−x F (x)G(x) dx
mit der Gewichtsfunktion w(x) = e−x. (Wir mussen uns zu geeigneten Funktio-
nen F, G beschranken, so dass die Integrale existieren, d.h. konvergieren).
Die entsprechenden orthogonalen Polynome heißen Laguerre’sche Polynome:
ϕ0(x) = 1, ϕ1(x) = x− 1, ϕ2(x) = x2 − 4x + 2
N = 1 ϕ1(x) = 0 fur x1 = 1
KAPITEL 12. GAUSS-QUADRATUR 109
α1 =
1∑
j=1
αj =
∞∫
0
w(x) dx =
∞∫
0
e−x dx = 1
⇒∞∫
0
e−xF (x) dx ≃ α1F (x1) = 1 · F (1)
oder mit F (x) = ex · f(x)
∞∫
0
f(x) dx ≃ α1ex1f(x1) = 1e1f(1) = ef(1)
N = 2 ϕ2(x) = 0 fur x(2)1 = 2−
√2 und x
(2)2 = 2 +
√2
Die Koeffizienten lauten
α1 =
∞∫
0
e−x
(x− 2−
√2
+2−√
2− 2−√
2
)dx
=−1
2√
2
∞∫
0
e−x(x− 2−√
2)dx =1 +√
2
2√
2
und
α2 =
∞∫
0
w(x) dx − α1 =
∞∫
0
e−x dx− 1 +√
2
2√
2=−1 +
√2
2√
2
⇒∞∫
0
e−xF (x) dx ≃ 1 +√
2
2√
2F(2−√
2)
+−1 +
√2
2√
2F(2 +√
2)
und deshalb mit f(x) = e−xF (x) oder F (x) = exf(x)
∞∫
0
f(x)dx ≃ 1 +√
2
2√
2e2−
√2 f(2−√
2)
+−1 +
√2
2√
2e2+
√2 f(2 +√
2)
⇒ heißt Gauß-Laguerre-Quadratur
Bemerkung: Wir benutzen nur endlich viele Stutzstellen hier, obwohl die
Lange des Intervalls unendlich ist.
Literaturverzeichnis
[1] M. Hanke-Bourgeois, Grundlagen der Numerischen Mathematik und des
Wissenschaftlichen Rechnens. Stuttgart: Teubner 2002.
[2] W. Oevel, Einfuhrung in die Numerische Mathematik, Heidelberg: Spek-
trum 1996.
[3] R.Plato Numerische Mathematik kompakt, Wiesbaden: Vieweg 2000.
[4] A. Quarteroni, R. Sacco, F. Saleri, Numerische Mathematik 1, Berlin:
Springer 2002.
[5] H.R. Schwarz, Numerische Mathematik, Stuttgart: Teubner 1997.
[6] F. Stummel, K. Hainer, Praktische Mathematik, Stuttgart: Teubner 1982.
110