sven.xn--kppel-jua.orgsven.köppel.org/uni/ws0910/Numerik/Scripte/Kloeden Numerische...

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

),


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 .


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:


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)


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


• 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).


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


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)


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∗


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


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


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∗


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!


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


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.


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


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| ≤ ε


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.


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


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∗)


⇒ 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.


⇒ 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.


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).


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


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


⇒ 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


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


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


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


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


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)


= 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


> 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


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


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)


= 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)


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


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)


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


∆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:


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!)


|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


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).


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


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


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


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).


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.


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


< 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


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


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


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


< 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


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


⇒ 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


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


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 )]


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 )]


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.


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“


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)

]


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)


⇒ 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


⇒ 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


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


=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)(ξ)|


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|


≤ 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)]


=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.


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.


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:


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)


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)(ξ),


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

)


=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.


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


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


⇒ 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


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) =


=

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)(ξ)


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


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


α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

sven.xn--kppel-jua.orgsven.köppel.org/uni/ws0910/Numerik/Scripte/Kloeden Numerische...

Documents

Transcript of sven.xn--kppel-jua.orgsven.köppel.org/uni/ws0910/Numerik/Scripte/Kloeden Numerische...