Mathematik f¨ur Informatiker III

(1)

Andreas Griewank

(griewank@math.hu-berlin.de) Wiss. Mitarbeiter:

Jan Riehme (riehme@math.hu-berlin.de) Stefan K¨orkel (skoerkel@math.hu-berlin.de)

Julia Sternberg (jstern@math.hu-berlin.de)

Institut f¨ur Angewandte Mathematik Humboldt Universit¨at zu Berlin

9. M¨arz 2006

– 1–

Teil D

Differentialgleichungen mit Numerik

Vorl¨aufige Gliederung 1.Numerik im ¨Uberblick

2.Gleitkommadarstellung und -arithmetik 3.Lösung (nicht-)linearer Gleichungssysteme 4.Gewöhnliche Differentialgleichungen (=ODE) 5.Euler Verfahren für Systeme von ODEs 6.Interpolation mit Polynomen und Splines 7.Quadraturen = Numerische Integration 8.Randwertprobleme und Schwingende Seite

I

Ubung zu 1-3 abzugeben am 8.11 ¨

I

Ubung zu 4-6 abzugeben am 22.11. ¨

– 2–

Literaturhinweise I

Peter Hartmann,

Mathematik f¨ur Informatiker. 3. ¨uberarbeitete Auflage, 2004, Vieweg.

Bei Lehmann’s vorhanden, ca. 30e.

Gute Grundlage, ¨ausserst lesbar, ISBN: 3-528-23181-5 Guerino Mazzola, G´erard Milmeister, Jody Weissmann, Comprehensive Mathematics for Computer Scientists 1, 2004, Springer.

Ziemlich axiomatisch und knapp geschrieben. Zweiter Band in Vorbereitung. Definitiv für höhere Ansprüche. Begleitender Kurs im Internet verfügbar. ca 30e, ISBN: 3-540-20835-6

Gerhard Opfer,

Numerische Mathematik für Anfänger. Eine Einführung für Mathematiker, Ingenieure und Informatiker. 4. durchgesehene Auflage, 2002, Vieweg

– 3–

Literaturhinweise II

Hans-G¨org Roos, Hubert Schwetlick,

Numerische Mathematik. Das Grundwissen f¨ur jedermann.

Mathematik f¨ur Ingenieure und Naturwissenschaftler. 1999, Teubner Friedrich Stummel, Karl Hainer,

Praktische Mathematik. 1982, Teubner J.M. Ortega, W.C. Rheinboldt,

Iterative solution of nonlinear equations in several variables. 1970 Academic Press, Inc.

Josef Stoer,

Numerische Mathematik 1. Eine Einf¨uhrung - unter

Ber¨ucksichtigung von Vorlesungen von F.L. Bauer. 7. neubearbeitete und erweiterte Auflage, 1994, Springer.

– 4–

(2)

Numerik im ¨Uberblick – Was ist, was will ’Numerik’

D - 1 Numerik im ¨ Uberblick – Was ist, was will ’Numerik’

Ausgangsdilemma

Die Modellierung natur- oder sozialwissenschaftlicher Zusammenhänge bzw ’Systeme’ führt zu mathematischen ’Gleichungen’, die nur in ganz einfachen Fällen per Hand oder sonstwie ’exakt’ gelöst werden können.

Zum Beispiel können schon bei der unbestimmten IntegrationMapleund Mathematicanur in speziellen Ausnahmefällen eine Lösung als Formel angeben.

Es l¨asst sich sogar zeigen, dass eine solche ’symbolische’ L¨osung im Regelfall garnicht existiert.

– 5–

Praktischer Ausweg

Die mathematischen Gleichungen werden in Computerprogramme umgesetzt und, wenn es sich dabei um Differentialgleichungen handelt

’diskretisiert’.

Die resultierenden Systeme linearer oder nichtlinearer algebraischer Gleichungen werden dann annäherungsweise über dem Raster(=Screen) der Gleitkommazahlen gelöst

Die Ergebnisse werden ausgedruckt oder besser graphisch dargstellt.

– 6–

Stufen des ’Wissenschaftlichen Rechnens’

(i)Modellierung ( des Anwendungssystems )

(ii)Diskretisierung ( von Differentialgleichungen )

(iii)Dateneingabe ( f¨ur aktuelle Situation )

(iv)L¨osung ( durch Gleitkomma-Algorithmen )

(v)Datenausgabe ( in geeigneter Form )

Eventuell k¨onnen (iii) - (v) auch innerhalb einer Wiederholungsanweisung (Schleife, Schlaufe) ausgef¨uhrt werden (z.B. wenn die Ausgabe zur Echtzeitsteuerung eines System dient).

– 7–

Numerische Grundaufgaben und ihre L¨osbarkeit

Lineare algebraische Gleichungssysteme

Im Prinzip völlig im Griff. Variablenzahl jeweils durch Speichergrösse und Prozessorzahl und -geschwindigkeit beschränkt.

Nichtlineare algebraische Gleichungssysteme

Lokal, d.h. bei vorhandener guter Anfangsn¨aherung: wie linearer Fall.

Global: beliebig schwierig und eventuell unl¨osbar.

Anfangswertaufgaben f¨ur ODEs

Im Prinzip völlig im Griff unabhängig von Linearität.

Randwertaufgaben f¨ur ODEs

Standarddiskretisierung f¨uhrt auf lineare bzw nichtlineare algebraische Gleichungen und ist entsprechend l¨osbar.

Partielle Differentialgleichungen PDE

Nur im elliptischen Fall schnell lösbar, alles andere ist Forschungsgebiet und stösst jeweils an die Grenzen vorhandener Rechnerkapazitäten.

– 8–

(3)

Numerische Grundaufgaben und ihre L¨osbarkeit

Warnung

Alles wird beliebig viel schwieriger wenn

Ieinige Variablen ganzzahlig sein m¨ussen und / oder

Idie Lösung gegebenen Ungleichungen genügen muss wie in der Optimierung üblich.

– 9–

Mathematik f¨ur Informatiker III Gleitkommadarstellung und -arithmetik

D - 2 Gleitkommadarstellung und -arithmetik

Ein System von Gleitkommazahlen wird definiert durch:

IBasis (oder Radix)b(= ¨ublicherweise 2)

IMantissenl¨angel

IMinimaler Exponentemin IMaximaler Exponentemax

Teilmenge der reellen Zahlen

R

mit Darstellung

x= −1s

0.m1m2· · ·ml

| {z }

Mantissem

b^e∼ −1s

m1bê⁻¹+m2bê⁻²+m3bê⁻³+. . .+mlbê⁻^l

Vorzeichenbit

s, Mantissem, Exponente

s∈

0,1 mi ∈ {0,1, . . . ,b−1} e∈ {emin,emin+ 1, . . . ,emax}

– 10–

Bin¨ardarstellung, d.h. Basis

b

= 2

ist die am h¨aufigsten verwendete Basis von Gleitkommazahlen Auchb= 10 wird zuweilen in Hardware verwendet.

Arten von Gleitkommazahlen

Inormalisierte Gleitpunktzahl:

m1 >0 =⇒ 1 b ≤m ≤

x b⁻^e <1

x=±0.m1m2m3· · ·ml·b^ewithm1>0 =⇒eindeutige Darstellung

Iunnormalisiert:m1= 0 zugelassen =⇒ keine Eindeutigkeit

Idenormalisiert:m1 = 0,e =emin

Vorsicht:

Rechnen mit denormalisierten Zahlen f¨uhrt zu verst¨arkten

Rundungseffekten. _{– 11}_–

Betragsm¨assig kleinste normalisierte Zahl

TINY TINY= 0.1·b^e^min=b^e^min⁻¹

Betragsm¨assig gr¨oßte normalisierte Zahl

HUGE

HUGE= 0.(b−1)(b−1)(b−1). . .(b−1). . .b^e^max=b^e^max(1−b^−l)

Epsilon (relative Maschinengenauigkeit)

ε

ist die kleinste Zahlεf¨ur die 1 +εin Gleitkommaarithmetik nicht 1 ergibt, d.h.ε≈b⁻^l

Merke:

IMantissenl¨angelbestimmt die Rechengenauigkeit.

IExponentenbereichemax−eminbestimmt den Wertebereich.

– 12–

(4)

Beispiel D.1 (Gleitpunktzahlsystem mit Basis 2 und Mantissenl¨ange 3)

PSfrag replacements

x= 0.m1m2m32^e Exponentenbereich−1≤e≤1 Normalisierte positive Zahlen: m1= 1, m2∈ {0,1} 3m3

Denormalisierte positive Zahlen: m1= 0, e=−1, m2∈ {0,1} 3m3

v5=v3∗v4

0

denormalisiert TINY=¹₄, HUGE=⁷₄, EPSILON=¹₈ 1

e −1−1−1−1−1−1−1−1 m1

m2

m3

1 16

1 8 3

16 1 4 5 16 3

8 7 16 1

2 5 8

3

4 7

8 1 ⁵₄ ³₂ ⁷₄

0 0 0 0

0 0 0 0 0

0 0 0

0 0 0 0 0

0 0 0

0

0 0 0 0

1 1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1

1 1 1

1 1 1 1 1

1 1 1

1 1 1 1

-1

– 13–

Beispiel D.2 (Einfache genaue Gleitkommazahlen im Salford Fortran 95 Compiler)

b= 2, l= 24, emin=−125, emax= 128 HUGE ≈ 2¹²⁸ = 2¹⁰12.8

≈ 10³12.8

≈ 10³⁸ TINY ≈ 2⁻¹²⁵⁻¹ = 2¹⁰_−12.6

≈ (10³)^−12.6 ≈ 10⁻³⁸ Epsilon ≈ 2⁻²⁴ = 2¹⁰_−2.4

≈ 10³_−2.4

≈ 10⁻⁷

Folgerung D.3

Bei Verwendung der Gleitkommazahlen des Salford Fortran 95 Compilers in Standardgenauigkeit wird mit etwasieben signifikanten

Dezimalstellengerechnet.

– 14–

Gleitpunktoperationen

Bemerkenswert

( 1.0 / 8.0 ) * 8.0=1.0 ( 1.0 / 5.0 ) * 5.06=1.0

Konsequenz

Gleitpunktoperationen st¨oren normale algebraische Rechenregeln, insbesondere Distributivit¨at:

Im Allgemeinen gilt (a+b)∗c6=a∗c+b∗c.

Man muss sich also ¨uber die Reihenfolge der Anwendung von Operationen Gedanken machen.

– 15–

Gleitpunktoperationen

Allgemein g¨ultiger Standard: ANSI - IEEE 754

(ANSI→American National Standards Institute und IEEE→Institute of Electrical and Electronics Egineering.)

Grundideen:

(i)Alle Zwischenergebnisse werden zur n¨achsten Gleitpunktzahl gerundet.

(ii)The show must go on. Auch bei Fehlern wird weiter gerechnet.

– 16–

(5)

Zu Grundidee (i) – Rundung von Zwischenergebnissen

Auch wennxundyim Gleitpunktbereich liegen, gilt dies im Allgemeinen nicht für das Ergebnisx◦y, wobei◦ ∈ {−,+,·, /}. Dann wirdx◦y zunächst mit erhöhter Genauigkeit berechnet und anschließend zur nächstliegenden Gleitpunktzahl gerundet.

Rundungsarten

∇(x◦y) nach unten gerundet

(gr¨oßte untere Schranke im Gleitpunktbereich)

∆(x◦y) nach oben gerundet

(kleinste obere Schranke im Gleitpunktbereich)

Verh¨altnis der Rundung nach oben und unten

Fallsegemeinsamer Exponent von ∆(x◦y) und∇(x◦y) ist, dann gilt

∆(x◦y) − ∇(x◦y) ≤2⁻^l2^e≤2⁻^l2· |x◦y|, da|x◦y| ≥¹22^e

q q

0.m·2^e 0.˜m·2^e

– 17–

Bezeichnet man also mit(x◦y)∈ {∇(x◦y),∆(x◦y)}die Gleitpunktzahl, die am n¨achsten zux◦yliegt, so gilt

|(x◦y)−x◦y| ≤ ¹2|∆(x◦y)− ∇(x◦y)| ≤2⁻^l|x◦y| ≤ eps· |x◦y| wobeieps= 2^−ldie relative Maschinengenauigkeit ist.

Alternative Schreibweise:

fl(x◦y) = (x◦y)∗(1 +ε), wobei |ε| ≤eps.

fl(x◦y) bezeichnet das in Gleitpunktarithmetik erzielte Ergebnis f¨ur x◦y.

Konsequenz f¨ur relativen Fehler:

fl(x◦y)−(x◦y) x◦y

≤ |ε| ≤eps

– 18–

Warnung:

Rundungsfehler entstehen in (fast) jeder einzelnen Operation und pflanzen sich fort.

Algorithmen (z.B. zur Matrixfaktorisierung) müssen deswegen auf ihre Stabilität, d.h. die Verstärkung oder Abdämpfung von Rundungsfehlern, untersucht werden.

Beispiel D.4

Gausssche Elimination ohne Pivotierung ist extrem instabil.

Gauss mit Pivotierung ist dagegen recht stabil.

– 19–

Frage

Was passiert, wennx◦yaußerhalb des Wertebereichs[-HUGE, HUGE]

liegt, d.h. entweder∇(x◦y) oder ∆(x◦y) nicht existiert?

Beispiel D.5 (Programm) REAL u,s,t

s = TINY(u)**2 ! ergibt 0

t = HUGE(u)*8 ! ergibt INF, signalisiert OVERFLOW

– 20–

(6)

Zu Grundidee (ii) – Fortsetzung der Berechnung trotz Fehlers

MitINFund-INFkann (soweit es geht)normalweiter gerechnet werden, ohne dass sich je wieder normale Zahlen ergeben.

(Einige) Rechenregeln

x + INF == INF für allex6=-INF x * INF == sign(x) * INF fürx6=0 x / 0 == sign(x) * INF fürx6=0 wobeisign(x)das Vorzeichen von x liefert.

Undefinierte Operationen wie0/0,INF/INF,INF-INFund0*INF ergeben den sehr speziellen WertNaN≈Not a Number.

Da einNaNnicht mit sich selbst oder etwas anderem verglichen werden kann, gilt

x 6=x .EQUIV. .TRUE.

genau dann wennxeinNaNist.

– 21–

Zu Grundidee (ii) – Fortsetzung der Berechnung trotz Fehlers

Infektionsprinzip:

Wenn immer einNaNals Argument oder Operator einer Operation auftritt sind die Ergebnisse wiederumNaNs.

Auf diese Weise wird der gesamte Berechnungszweig als ung¨ultig ausgewiesen.

– 22–

Mathematik f¨ur Informatiker III Summation numerischer Reihen

Fehlerfortpflanzung

D - 3 Summation numerischer Reihen Fehlerfortpflanzung

Erinnerung:

fl(x◦y) =x◦y∗(1+ε) mit −eps≤ε≤eps wobei ◦ ∈ {+,−,∗, /}

Prinzip Hoffnung f¨ur komplexe Berechnungen

Da Auf- oder Abrunden mehr oder minder zuf¨allig auftreten hebt sich deren Wirkung (hoffentlich) im Großen und Ganzen auf.

– 23–

Fehlerfortpflanzung

Positives Beispiel: Geometrische Reihe:

s= Xn

i=0

xⁱ=1−xⁿ⁺¹

1−x falls x6= 1 .

Einfach genaues Auswertungsprogramm in Fortran 95

INTEGER i,n

REAL(KIND=1) x,y,s REAL(KIND=2) check

s = 0 ! Partialsumme

y = 1 !jeweils Potenz von x

DO i = 0, n

s = s+y ; y = y*x END DO

check = x ; eps = EPSILON(x) check = (1-check**(n+1))/(1-check) WRITE(*,*) s,check,s/check-1,n*eps

– 24–

(7)

Fehlerfortpflanzung

Programm ergibt f¨ur

n

= 100 und

x

= 2.0/3.0

s check s/check - 1 n * eps 3.0000002 3.00000019 2·10⁻⁸ 1.2·10⁻⁵

Beobachtungen

IGleitpunktwert vonxist offenbar gr¨oßer als²₃(durch Rundung), da beide Summen gr¨oßer als

1 +2 3+

2 3

2

+· · ·+ 2

3 n

= 3 1− 2

3 n+1!

| {z }

≤1

≤3

IDer beobachtete relative Fehler zwischen einfach und doppelt genauer Lösung ist lediglich 2·10⁻⁸, d.h. von der Größenordnung der Maschinengenauigkeit, obwohl wir 100 Operationen durchgeführt haben. Die Rundungen scheinen sich partiell aufgehoben zu haben.

IEine exakte Absch¨atzung f¨ur denworst case(d.h. schlimmster Fall) ergibt den Wert (1 +eps)¹⁰⁰≈100·epsals relativen Fehler. Das

l¨asst sich wie folgt herleiten. – 25–

Fehlerfortpflanzung

Theoretische Schranke des Fehlers im obigen Programm

F¨uryi+1=fl(yi∗x) als berechneter Wert vonyimi-ten Schritt gilt:

y0= 1 y1=x

y2=fl(y1·x) =x²(1 +ε2)

y3=fl(y2·x) =x³(1 +ε2)(1 +ε3) =x³(1 + ˜ε3)² wobei|˜ε3| ≤eps y4=fl(y3·x) =x⁴(1 + ˜ε2)²(1 +ε4) =x⁴(1 + ˜ε4)³

...

yi=xⁱ(1 + ˜εi)ⁱ⁻¹ ...

yn=xⁿ(1 + ˜εn)ⁿ⁻¹

– 26–

Fehlerfortpflanzung

Entsprechend erh¨alt man f¨ur die Partialsummensi+1=fl(si+yi) als berechnete Werte von 1 +x. . .+xⁱ⁺¹

s1=fl(y0+y1) =fl(1 +x) = (1 +x)(1 +εn+1) s2=fl(s1+y2) =fl(s1+y2)(1 +εn+2)

= (1 +x)(1 +εn+1) +x²(1 +ε2) (1 +εn+2)

= (1 +x+x²)(1 + ˜εn+2)² f¨ur |˜εn+2| ≤eps sn= (1 +x+x²+· · ·+xⁿ)(1 + ˜ε2n)ⁿ≤s(1 +ε)ⁿ

so dass fallseps¹n ⇐⇒ n·eps1

|(sn/s−1)|=|(1 +ε)ⁿ|−1 = 1+n·ε+n·(n−1)

2 ε². . .−1≈n·|ε|≤n·eps

Ergebnis:

Worst case error

- Absch¨atzung:

|sn/s−1| ≈n·eps

– 27–

Fehlerfortpflanzung

Negatives Beispiel (d.h. Prinzip Hoffnung versagt) : Harmonische Reihe

X∞

i=1

1 i =











∞ (mathematisch, in exakter Arithmetik)

15.403 auf Griewank’s Laptop, in einfacher Genauigkeit (f¨ur alle hinreichend großen Summations-Schranken

= Zahl der Terme)

Frage:

Was passiert?

Antwort:

Die Summation bleibt irgendwannliegen, da die zus¨atzlichen Terme im Vergleich zur berechneten Teilsumme zu klein werden.

– 28–

(8)

Fehlerfortpflanzung

Erkl¨arung:

BetrachtekleinenSummandenyundgroßenSummanden x= 0.m1m2. . .ml·2êso dassx=x+ 2⁻^l+edie nächst größere Gleitpunktzahl zuxist undx=x−2⁻^l+eist die nächst kleinere Gleitpunktzahl zux.

PSfrag replacements

2^e⁻¹ x

2^−l+e 2^−l+e

2^e

x x

Konsequenz:

Falls|y|<¹₂2⁻^l+e= 2⁻^l^−1+e gilt immer fl(x+y) =x.

Eine hinreichende Bedingung ist:|y| ≤ |x| ·eps.

– 29–

Fehlerfortpflanzung

Am Beispiel derharmonischen Reihegilt nach (n−1) Termen:

x=

n−1

X

i=1

1 i &

Zn 1

1 zdz= ln(n).

Also bleibt die Summationliegen(d.h. die Partialsummen wachsen nicht mehr weiter) wenn

|y|=1

n≈ln(n)·eps was auf jeden Fall gilt wenn

n& 1 eps·ln(n)

– 30–

Fehlerfortpflanzung

Beispiel D.6 (Programm, das die harmonische Reihe summiert, bis die Partialsummen konstant bleiben:)

REAL(KIND=1) salt,sneu,one

salt = -1 ; sneu = 0 ; one = 1.0 ; n = 1 DO WHILE (sneu6=salt)

salt = sneu sneu = sneu+one/n n = n+1

END DO

WRITE(*,*) sneu,n

Ergebnis auf Griewank’s Laptop

sneu = 15.403. . . n = 2097152≈2·10⁶ Laufzeit ≈ ¹6Sekunde

D.h. obiger Schleifenk¨orper wird in etwa 10⁷mal pro Sekunden ausgef¨uhrt (entspricht ca. 10 Megaflops, d.h. 10 Millionen Operationen/Sekunde.)

– 31–

Fehlerfortpflanzung

Vergleich zur theoretischen Herleitung

n= 2097152 ergibt ln(n)∗n∗EPSILON(x) = 3.6

Frage:

Was passiert bei Ausf¨uhrung des obigen Programms, wenn statt mit einfacher Genauigkeit (d.h.KIND=1) nun mit doppelt genauen Gleitkommazahlen (d.h.KIND=2) gerechnet wird?

Antwort:

Das Programm l¨auftewig, daeps⁻¹und damit dann auchnum Faktor 2⁵³/2²⁴≈2²⁹≈¹210⁹gewachsen ist.

In Sekunden:

1 6·1

2·10⁹s = 10⁸

36·10³h = 25·10⁴h = 25.000 Stunden ≈1000 Tage.

– 32–

(9)

Rundungsfehlerabsch¨atzung bei Riemann

Verallgemeinerung der harmonischen Reihe:

Riemannsche Zetafunktion

ζ(x) =

X∞

k=1

1

k^x f¨ur x>1 Konvergenzbeweis mittels Integralschranke

PSfrag replacements X∞

k=1

1 k^x ≤1 +

Z∞

1

dy

y^x = 1 +y⁻^x+1 1−x

∞ 1

= 1− 1 1−x = x

x−1 1

k⁻^x

– 33–

∆ζn(x) = ζ(x)−ζn(x) = X∞

k−1

k⁻^x− Xn

k=1

k⁻^x

= X∞

k=n+1

k⁻^x ≤ Z_∞

n

k⁻^xdk = k¹⁻^x 1−x

∞ k=n

= 1

k¹⁻^x(1−x)

∞ k=n

= 0− 1

n^x⁻¹(1−x) = 1

n^x⁻¹(x−1) ≤tol

⇒ n ≥ ^x−1 s 1

tol(x−1)

– 34–

Partialsummen:

ζn(x) = Pn k=1 1

kx wachsen monoton mitnund sind nach oben durch_x₋₁^x beschr¨ankt, haben also einen eindeutigen Grenzwertζ(x).

Praktische Notwendigkeit: Diskretisierung

Hier, wie häufig in numerischer Mathematik muss mathematisches Problem durch Ausführung endlich vieler Operationen auf endlich vielen Variablen annäherungsweise gelöst werden. Hier einfach Annäherung von ζ(x) durchζn(x). Der entsprechende Abbruchfehler|ζ(x)−ζn(x)|kann hier einfach mit Hilfe einer Integralschranke abgeschätzt werden.

Unabh¨angig vom in der Numerischen Analysis betrachteten Diskretisierungsfehlerist der Rundungsfehler zu ber¨ucksichtigen.

– 35–

Rundungsfehlerabsch¨atzung

F¨urbi>0

fl . . . b1+b2

+b3

+b4

. . .+bn+1

+bn

= . . . b1+b2

1 +ε1

+b3

1 +ε2

+b4

1 +ε3

. . .+bn

1 +εn−1

= b11 + ˜ε1n−2

+b21 + ˜ε1n−1

+b31 + ˜ε2n−2

+. . .+bn1 + ˜εn−11

=⇒

fl b1+. . .+bn

− b1+b2+. . .+bn

≤ b1

h

1 +epsn−1

−1i +b2

h

1 +epsn−1

−1i

+. . .+bn1 +eps

≈

b1+b2

(n−1) + (n−2)b3+ (n−3)b4+. . .+bn

eps

Mit anderen Worten:

Der an derj+ 1-ten Stelle eingebrachte Summand wird (n−j) -mal in den Operationen von einer Rundung betroffen und tr¨agt entsprechend zur Gesamtfehlerschranke bei.

– 36–

(10)

Schlussfolgerung:

Um Rundungsfehler zu minimieren sollten Summen möglichst vom kleinsten zum größten Summanden gebildet werden. Bei konvergenten (hoffentlich monoton fallenden) Reihen sollte von hinten, d.h. rückwärts summiert werden.

Beispiel D.7 (ζ(2) auf G’s Laptop in einfacher Genauigkeit:)

ζ(2) = X∞

k=1

1 k²≡











π²/6 = 1.6449340. . . exakt

1.6447253 vorwärts bis. liegen bleiben n= 4097 1.6446900 rückwärts vom gleichen n= 4097 1.6449339 rückwärts mit n= 2²³= 8388608

Bemerkung:

Durch Rückwärtssummation können deutlich mehr Summanden der Form 1/k⁻xmitn>4097 ihren Beitrag zur Gesamtsumme leisten. Mehr Summanden zu benutzten bedeutet aber, denDiskretisierungsfehlerzu verringern und damit den exakten Wertζ(x) besser zu approximieren.

– 37–

Absch¨atzung des Rundungsfehlers Vorw¨arts:

eps Xn

k=1

1

k²(n−k) =eps Xn

k=1

n k²

−1 k

≈eps

nπ² 6 −ln(n)

≈eps·n·π² 6

R¨uckw¨arts:

eps Xn

k=1

1 k²k=eps

Xn

k=1

1

k≈eps·ln(n)

Vergleich:

eps·n·π²

6 eps·ln(n)

– 38–

Konvergenzbeschleunigung (1. Stufe nach Wijngaard)

Beobachtung bei Riemann:

ζ(x) = 1 + 1

2^x+· · ·+ 1 100^x+ 1

101^x+ 1 102^x+· · ·

| {z }

sp¨atere Terme ¨andern sich nur langsam

Idee:

Erste grobe Ann¨aherung mitbk= 1 k^x

a1=b1+b2·2 +b4·4 +· · ·+ (b2ⁱ)·2ⁱ > ζ=b1+b2+b3+b4. . . Reihe der 2ⁱb2ⁱkonvergiert viel schneller alsP

bk. Die Korrektur erfolgt durch transformierte Terme

aj= X∞

i=1

bj2ⁱ

2ⁱ.

– 39–

Satz D.8

Satz: F¨ur bk=k⁻^xoder andere monoton konvergierende Reihen gilt im Grenzwert

X∞

k=1

bk= X∞

j=1

(−1)^j⁻¹aj .

Bemerkung

Bemerkung: Die neue Reihe ist alternierend, wobeiaj≥bj, d.h. die einzelnen Terme gehen nicht schneller gegen Null als die der Ursprungsreihe.

– 40–

(11)

Idee des Beweises:

Betrachte, wie oftbkinajauftritt

Vorz j\k 1 2 3 4 5 6 7 8 9 10 11 12

+ 1 1 2 − 4 − − − 8 − − − −

− 2 − 1 − 2 − − − 4 − − − −

+ 3 − − 1 − − 2 − − − − − 4

− 4 − − − 1 − − − 2 − − − −

+ 5 − − − − 1 − − − − 2 − −

− 6 − − − − − 1 − − − − − 2

+ 7 − − − − − − 1 − − 2 − −

P

mit Vorzeichen

1 1 1 1 1 1 1 . . . .

Bemerkung

Bei Riemann k¨onnen dieai=ai(x) sogar explizit berechnet werden.

– 41–

Schlussfolgerungen aus dem Summationsbeispiel

IDie Behandlung mathematischer und anderer Modellierungsprobleme bedingt das Auftreten vonAbbruchs-≡Diskretisierungsfehlernsowie Rundungsfehlern. Beide sollten abgesch¨atzt und m¨oglichst minimiert werden.

IGleitpunktarithmetik ist weder kommutativ noch assoziativ, distributiv usw.

Spezielle Konsequenz:Betragsm¨aßig fallende Reihen von hinten summieren!

IEs ist erstaunlich einfach, an die Grenzen der Gleitpunkt- und Ganzzahlarithmetik zu stoßen.

IViele Jobs (≡Programme, Daten) laufen entweder im Sekunden- oder Stundenbereich. Beobachtung der Abarbeitung im Minutenbereich ist relativ selten.

IMathematisch endlichist nicht gleichrechentechnisch endlich.

– 42–

Mathematik f¨ur Informatiker III L¨osung (nicht-)linearer Gleichungssysteme

D - 4 L¨osung (nicht-)linearer Gleichungssysteme

Methoden zur L¨osung des linearen Problemes

Ax

=

b

mit dim(x) = dim(b) =

n

ICramersche Regelxi= (−1)ⁱdet(Ai)/det(A) f¨uri= 1..n ( InAiwird diei−te Spalte vonAdurchbersetzt )

IGauss-Elimination≈P A = LU Faktorisierung

(PPermutation,Lunterhalb undUoberhalb dreiecksf¨ormig )

ISchmidt-Ortogonalisierung≈A = QR Faktorisierung (Qorthogonal,Roberhalb dreiecksf¨ormig )

IFixpunkt Iterationx←x−M F(x) mitF(x)≡Ax−b (M∈Rⁿ×nangen¨aherte Inverse so dassM A≈I)

Hinweise:

IFür (eindeutige) Lösbarkeit ist überalldet(A)6= 0 vorrauszusetzen.

IL¨oseLUx=bbzwQRx=bdurch Substitution/Transponierung.

IDie letzte Methode l¨asst sich auch auf nichtlinearesF(x) anwenden.

– 43–

Linearisierung des ’Freistoss’ Beispieles

Das nichtlineare System von 3 Gleichungen in 3 Unbekannten F1(x1,x2,x3) = x1∗x2−4.9∗x₁²−2 = 0 F2(x1,x2,x3) = 10∗ln(1 + 0.1∗x3∗x1)−25 = 0 F3(x1,x2,x3) = (x2−9.8∗x1)∗(_x¹₃+ 0.1∗x1) +√¹

3 = 0

hat dieJacobimatrix F⁰(x)≡ ∂

∂xF(x)≡

∂Fi

∂xj

i=1,2,3

j=1,2,3

≡







x2−9.8∗x1 x1 0

x3

1+0.1∗x₁∗x₃ 0 _1+0.1∗^x¹x₁∗x₃

z(x) x¹₃+ 0.1∗x1 −^x²^−9.8∗x²₃ ^x¹







mitz(x)≡ −9.8∗(_x¹₃+₁₀^x¹) +₁₀¹(x2−9.8∗x1) =^x₁₀²−9.8

1 x3+¹₅x1

– 44–

(12)

Linearisierung durch Jacobimatrix

Falls f¨urF:Rⁿ→Rⁿdien²Komponenten der Jacobimatrix F⁰(x) ≡ ∂

∂xF(x)≡

∂Fi

∂xj

i=1,...,n

j=1,...,n

bezüglich jeder der Variablenx1, . . . ,xnLipschitz-stetig sind, so lässt sich aus dem Hauptsatz der Differential- und Integralrechnung herleiten, dass für jedenSchritt s∈Rⁿgilt

kF(x+s)−[F(x) +F⁰(x)s]k ≤ γksk² Hierbei istF⁰(x)sein Matrix-Vektor Produkt undk · kist eine Vektor- bzw. Matrixnorm (siehe Abschnitt B-3) mit

kF⁰(x)−F⁰(y)k ≤ γkx−yk

Fx(s)≡F(x) +F⁰(x)sist als Funktion des variablen Vektorssdie Linearisierung ( verallgemeinerte Tangente ) vonFan der Stellex.

– 45–

Newton’s Methode im Vektorfall

Setzt man die LinearisierungFx(s) =F(x) +F⁰(x)szu null so erh¨alt man das lineare Gleichungssystem

As=b mit A=F⁰(x) und b=−F(x)

Die Lösung lässt sich ausdrücken als

s=A⁻¹b=−F⁰(x)⁻¹F(x) und heisstNewtonschritt.

Wiederholte Berechnung vonsund anschliessende Inkrementierung x←x+sergibt Newton’s Methode

x^(k+1)≡x^(k)+s^(k) mit F⁰(x^(k))s^(k)=−F(x^(k)) f¨ur k= 0,1, . . . Hierbei z¨ahlt der hochgestellte Index (k) die Iterationen.

– 46–

Newton’s Methode im Vektorfall

Warnung:

IDas Verfahren muss abgebrochen werden wenndet(F⁰(x^(k))) null oder sehr klein ist.

IIm letzteren Falle werden die Schrittes^(k)typischerweise sehr gross und f¨uhren h¨aufig zu Argumentenx^(k+1)woFgarnicht mehr ausgewertet werden kann.

IZur Vermeidung dieses Problems wirds^(k)manchmal mit einem D¨ampfungsfaktorα^(k)<1 multipliziert, der dannSchrittweite genannt wird. Wir iterieren also effektiv

x^(k+1)=x^(k)−α^(k)F⁰(x^(k))⁻¹F(x^(k))

Die Bestimmung eines geeignetenα^(k)heisst auch Strahlsuche (engl:

Line Search).

– 47–

Lokale Konvergenz von Newton

Satz D.9 (Satz von Kantorovich)

Sei die Vektorfunktion F:Rⁿ→Rⁿeinmal differenzierbar und besitze ihre Jacobimatrix F⁰(x)∈Rⁿ×ndie Lipschitzkonstanteγ.

Weiterhin sei x⁽⁰⁾ein Punkt an dem F⁰(x⁽⁰⁾)regul¨ar ist und somit eine Inverse F⁰(x⁽⁰⁾)⁻¹existiert. Mitk · kals induzierte Matrix-Norm folgt dann aus

F⁰(x⁽⁰⁾)⁻¹

2 F(x⁽⁰⁾)

≤ 1 2γ

dass Newton’s Methode zu einer L¨osung x^(∗)mit F(x^(∗)) = 0konvergiert.

Die Konvergenzgeschwindigkeit ist quadratisch in dem Sinne dass f¨ur eine Konstante c und alle k gilt

x^(k+1)−x^(∗) ≤c

x^(k)−x^(∗)

2

Bemerkung:

Je nichtlinearer ein Problem umso grösser istγund desto stärker ist damit die Bedingung anx⁽⁰⁾. Wird praktisch nie überprüft !!!! – 48–

(13)

Mathematik f¨ur Informatiker III Gew¨ohnliche Differentialgleichungen (ODE)

D - 5 Gew¨ohnliche Differentialgleichungen (ODE)

(nachHartmann, Mathematik f¨ur Informatiker)

Definition D.10 (Gew¨ohnliche Differentialgleichungen (ODE))

Eine Gleichung, in der neben der unabh¨angigen Variablenxund einer gesuchten Funktiony=y(x) auch deren Ableitungen^d_dxⁿ^yn=y⁽ⁿ⁾(x) bis zur Ordnungnauftreten, heisstGew¨ohnliche Differentialgleichung n-ter Ordnung (ODE).

Sind ausserdem einx0aus dem Definitionsbereich vony(x) und zugeh¨orige Wertey(x0),y⁽¹⁾(x0), . . . ,y⁽ⁿ⁻¹⁾(x0) gegeben, so spricht man von einemAnfangswertproblem.

– 49–

Separable Differentialgleichungen

Definition D.11 (Separable Differentialgleichung)

Eine DifferentialgleichungF(x,y,y⁰) = 0 erster Ordnung heisst separabel, wenn sie sich in der Form

y⁰=f(x)g(y)

darstellen l¨asst, wobeif:I−→R,g:J−→Rstetige Funktionen auf den IntervallenI⊆R,J⊆Rsind.

Satz D.12 (L¨osbarkeit: Anfangswertproblem separabler ODE)

Eine separable Differentialgleichung erster Ordnung mit der

Anfangsbedingung y(x0) =y0f¨ur x0∈I , y0∈J, hat im Intervall J eine eindeutige L¨osung y(x) :I−→J, falls

g(y)6= 0 ∀y∈J.

– 50–

Seien

G(y) :=

Zy y0

1

g(y)dy, F(x) :=

Zx x0

f(x)dx die Stammfunktionen von_g(y)¹ bzw.f(x).

Dabei wurden f¨ur Integrationsvariable und Obergrenze der Integration das gleiche Symbol verwendet.

AufJistG⁰(y) =_g(y)¹ 6= 0 (Voraussetzung Satz D.12), daher istG streng monoton und besitzt eine UmkehrfunktionG⁻¹. Dann ist aber

y(x) :=G⁻¹(F(x))

die L¨osung des Anfangswertproblemsy⁰=f(x)g(y),y(x0) =y0.

– 51–

Probe:

G(y(x)) =F(x) =⇒ G⁰(y(x))y⁰(x) =F⁰(x) =_g(y(x))¹ y⁰(x) =f(x)

=⇒ y⁰(x) =f(x)g(y(x)) Anfangswert:y(x0) =y0

F(x0) = 0 =⇒ y(x0) =G⁻¹(F(x0)) =G⁻¹(0) G(y0) = 0 =⇒ G⁻¹(0) =y0

=⇒ G⁻¹(0) =y0=y(x0)

Satz D.13

Das Anfangswertproblem y⁰(x) =f(x)g(y), mit Funktionen f:I−→R, g:J−→R, und dem Anfangswert y(x0) =y0∈J, hat die eindeutige L¨osung y , die man erh¨alt, wenn man die folgende Gleichung nach y

aufl¨ost: Zy

y0

1 g(y)dy=

Zx x0

f(x)dx

– 52–

(14)

Lineare Differentialgleichungen erster Ordnung

Definition D.14 (Lineare Differentialgleichung)

Differentialgleichungen, bei denen die Funktiony=y(x) und ihre Ableitungen nur in linearem Zusammenhang auftreten heissenLineare Differentialgleichungen.

Lineare Differentialgleichungenerster Ordnunghaben die Form y⁰+a(x)y=f(x).

Ist die Funktionf(x)≡0 auf der rechten Seite identisch Null, so heisst die Gleichunghomogen, sonstinhomogen.

Die FunktionF(x) auf der rechten Seite heisstQuellfunktion.

– 53–

Lineare Differentialgleichungen erster Ordnung

Satz D.15 (L¨osung homogener linearer ODE)

Ist a(x)auf dem Intervall I stetig, so lautet die vollst¨andige L¨osung der linearen Differentialgleichung y⁰+a(x)y= 0

y(x) =c·e⁻^A(x) wobei c∈Rund A(x)eine Stammfunktion von a(x)ist.

Satz D.16 (L¨osung inhomogener linearer ODE)

Die inhomogen lineare Differentialgleichung y⁰+a(x)y=f(x), f,a:I−→Rstetig, x0∈I , besitzt die vollst¨andige L¨osung

y= Zx

x0

f(t)e^A(t)dt+c

·e⁻^A(x) wobei c∈Rund A(x)eine Stammfunktion von a(x)ist.

– 54–

Lineare Differentialgleichungen n-ter Ordnung

Definition D.17 (Lineare ODE n-ter Ordnung)

Eine Differentialgleichung der Form

y⁽ⁿ⁾+a1(x)y⁽ⁿ⁻¹⁾+· · ·+an−1(x)y⁰+an(x)y =f(x) heisstlineare Differentialgleichung n-ter Ordnung.

Dabei sind die Funktionenf,ai:I−→Rauf dem Intervall stetig.

Dieaiheissen Koeffizientenfunktionen,fheisst Quellfunktion.

Istf= 0, so heisst die Gleichunghomogen, sonstinhomogen.

– 55–

Satz D.18 (Existenz und Eindeutigkeit der L¨osung)

Sei

y⁽ⁿ⁾+a1(x)y⁽ⁿ⁻¹⁾+· · ·+an−1(x)y⁰+an(x)y =f(x) eine lineare Differentialgleichung n-ter Ordnung mit ai,f:I−→Rund x0∈I .

Dann gibt es zu den Anfangswerten

y(x0) =b0, y⁰(x0) =b1, . . . y⁽ⁿ⁻¹⁾(x0) =bn−1

genau eine L¨osung y=y(x)dieses Anfangswertproblems.

Diese L¨osung existiert auf dem ganzen Intervall I .

– 56–

(15)

Satz D.19 (L¨osungsstruktur linearer ODE n-ter Ordnung)

Die Menge H der L¨osungen y:I−→Rder homogenen linearen Differentialgleichung y⁽ⁿ⁾+a1(x)y⁽ⁿ⁻¹⁾+· · ·+an−1(x)y⁰+an(x)y = 0 mit ai:I−→Rbildet einen reellen Vektorraum der Dimension n.

Eine Basis des L¨osungsraumes H nennt manFundamentalsystem.

Jede L¨osung y der inhomogenen Gleichung

y⁽ⁿ⁾+a1(x)y⁽ⁿ⁻¹⁾+· · ·+an−1(x)y⁰+an(x)y =f(x)mit f:I−→R hat die Form

y=ys+yh

wobei xh∈H eine L¨osung der homogenen und yseine spezielle L¨osung der inhomogenen Differentialgleichung ist.

– 57–

Lineare Differentialgleichungen mit konstanten Koeffizienten

F¨ur inhomogene lineare Differentialgleichungen n-ter Ordnung (siehe Definition D.17) existiert kein allgemeines L¨osungsverfahren.

F¨ur den Fallkonstanter Koeffizientenfunktionen ai(x)∈Rkann jedoch ein Fundamentalsystem angegeben werden:

L¨osung des homogenen Systems

y⁽ⁿ⁾+a1y⁽ⁿ⁻¹⁾+· · ·+an−1y⁰+any = 0 L¨osungsansatz: Exponentialfunktion y(x) =e^λx und damit

y(x) =e^λx, y⁰(x) =λe^λx, y⁰⁰(x) =λ²e^λx, . . . , y⁽ⁿ⁾(x) =λⁿe^λx Einsetzen in die Differentialgleichung liefert

λⁿe^λx+a1λⁿ⁻¹e^λx+· · ·+an−1λe^λx+ane^λx = (λⁿ+a1λⁿ⁻¹+· · ·+an−1λ+an)e^λx = 0

– 58–

Definition D.20 (Charakteristisches Polynom)

Das Polynom

p(λ) :=λⁿ+a1λⁿ⁻¹+· · ·+an−1λ+an

heisst charakteristisches Polynom der homogenen linearen Differentialgleichung n-ter Ordnung mit konstanten Koeffizienten

y⁽ⁿ⁾+a1y⁽ⁿ⁻¹⁾+· · ·+an−1y⁰+any = 0.

Fortsetzung: L¨osung des homogenen Systems

Aus den Nullstellenλi,i= 1. . .nmitp(λi) = 0 des charakteristischen Polynoms kann ein Fundamentalsystem f¨ur die homogene

Differentialgleichung n-ter Ordnung konstruiert werden.

Dazu ist eine Fallunterscheidung nach derVielfachheit der Nullstellenλi

n¨otig:

– 59–

λ∈R

ist einfache Nullstelle

Dann ist e^λx eine L¨osung der Differentialgleichung.

λ

=

α

+

iβ∈C

ist einfache komplexe Nullstelle

e^α^xcosβx und e^α^xsinβx sind L¨osungen der Differentialgleichung.

λ∈R

ist

k-fache reelle Nullstelle xⁱe^λx, i= 0, . . . ,k−1 sindklinear unabh¨angige L¨osungen.

λ

=

α

+

iβ∈C

ist

k-fache komplexe Nullstelle xⁱe^α^xcosβx, xⁱe^α^xsinβx, i= 0, . . . ,k−1 sind die 2klinear unabh¨angige L¨osungsfunktionen.

Beispiel D.21

SieheHartmann, Mathematik f¨ur Informatiker, S.352 ff.

– 60–

(16)

Mathematik f¨ur Informatiker III Euler Verfahren f¨ur Systeme von ODEs

Systeme von ODEs und ihre numerische L¨osung

D - 6 Euler Verfahren f¨ur Systeme von ODEs Systeme von ODEs und ihre numerische L¨osung

In vielen Anwendungen wird der Zustand eines Systems zum Zeitpunktt durch einen Vektor

x(t) = [x1(t),x2(t), . . . ,xn(t)]^> mit n>0 beschrieben. Die Änderungsgeschwindigkeit ˙x≡dx(t)/dtdes Zustandes nach der Zeit ergibt sich häufig als FunktionF(x(t)) mitF:Rⁿ→Rⁿ eben dieses Zustandes. Also erhalten wir das System gewöhnlicher Differentialgleichungen

˙

x(t) =F(x(t)) kurz x˙ =F(x)

Das System heisst autonom, da die Zeittauf der rechten Seite nicht explizit, sondern nur mittelbar überx=x(t) vorkommt. Dieses ist keine Einschränkung da ein nichtautonomes System ˙x(t) =F(t,x(t)) sich autonom umschreiben lässt indem mantals nullte Zustandskomponente x0(t) hinzufügt und somit für ¯x≡(x0,x1, . . . ,xn)^Terhält

d dt¯x ≡

x˙0

˙ x

= ˙t

˙ x

= 1

F(¯x)

≡F(x)

– 61–

Auch ODEs h¨ohere Ordnungen lassen sich in Systeme von ODEs erster Ordnung umschreiben, indem man z.B. die erste Ableitungy⁰als neue abh¨angige Variablev≡y⁰definiert und danny⁰⁰durchv⁰ersetzt. So wird zum Beispiel aus einer nichtautonomen Differentialgleichung zweiter Ordnung

y⁰⁰ =f(t,y,y⁰)

das autonome System erster Ordnung in den drei Variableny0≡t, y1≡yundy2≡y⁰



 y0⁰

y₁⁰ y₂⁰



=



 1 y2

f(y0,y1,y2)





Entsprechend lassen sich Anfangsbedingungen umschreiben.

Die Umformulierung als System 1.Ordnung eröffnet die Möglichkeit numerische Standardmethoden und Software für die Lösung autonomer Systeme erster Ordnung mit Anfangsbedingungen zur Anwendung zu bringen.

– 62–

Satz D.22 (Existenz und Eindeutigkeit der L¨osung)

Sei F:D ⊂Rⁿ−→Rⁿin einem offenem GebietDlokal Lipschitz-stetig.

Dann existiert für jeden Punkt yo∈ Dein Intervall(a,b)30und eine eindeutige Lösung y(t)∈ Dder ODEy˙=F(y)für a<t<b mit y(0) =y0.

Bemerkung:

(i)F¨ur die Existenz einer L¨osung ist die Stetigkeit vonFhinreichend.

Vorraussetzung von Lipschitz - Stetigkeit ist f¨ur die Eindeutigkeit der L¨osung und die Konvergenz numerischer Verfahren erforderlich.

(ii)Das Intervall (a,b) kann so gross gew¨ahlt werden, dassy(b) den Rand vonDerreicht.

– 63–

Eulers Methode und andere explizite ODE-L¨oser

Die meisten ODEs haben keine geschlossen darstellbare L¨osung.

Die L¨osung kann aber durch numerische Methoden mit (mehr oder weniger) beliebiger Genauigkeit approximiert werden.

Numerische Approximationen sind auch alles, was zur Berechnung der mathematischen Standardfunktionene^x, sinxetc. zur Verf¨ugung steht, da diese Funktionen als L¨osung von ODEs definiert sind.

Die einfachste numerische Methode zur L¨osung von ODEs ist das Explizite (Vorw¨arts) Eulersche Polygonzugverfahren.

– 64–

(17)

Explizite (Vorw¨arts) Euler-Methode

Seiy(t)die exakte L¨osung vony(t) =˙ f(t,y(t))mity(0) =y0.

h 2h 3h tk=k·h T t

y

y(k·h)

y(T) exakter Wert

yk

yn=yt/h

imk-ten Schritt berechneter Wert

˙

y(k·h) =f(tk,yk)

≡Anstieg derTangen- tey(t) der L¨osung˙ y(t) intk y(0) =y0

Gesucht wird alsoyk≈y(tk)f¨urk= 0, . . . ,^T_h mittk=k·h:

yk+1≡yk+h f(tk,yk) ≈y(tk+1)

– 65–

Beispiel D.23 (Autonome lineare ODE)

˙

y=λy mit λ∈R und y0= 1 Anwendung von Eulers Methode:

y1 = y0+hλy0 = (1 +hλ)y0

y2 = y1+hλy1 = (1 +hλ)y1 = (1 +hλ)²y0

...

yk = (1 +λh)^ky0 = (1 +λh)^k ...

yn = (1 +λh)ⁿy0 = (1 +λh)^T^h

Vergleich mit exakter L¨osung:

y(t) = exp(λt) ergibt am EndpunktT y(T) =e^λT≡lim

h→0(1 +λh)^T^h= lim

n→∞

1 +λT

n n

– 66–

Erl¨auterung

Die angenäherte LösungyT/hkonvergiert gegen die exakte Lösungy(T) der ODE wenn die Schrittweiteh=T/ngegen Null geht. Das bedeutet aber dass die Anzahl der Eulerschritte und damit der

Berechnungsaufwand gegen∞gehen.

Frage:

Kann der ApproximationsfehlerkyT/h−y(T)kals Funktion der Schrittweiteh=T/ndargestellt und somit zur Bestimmung einer vern¨unftigen Schrittzahlngenutzt werden?

Antwort:

JA!

Im vorliegenden speziellen Fall gilt

hlim→0

yT/h

y(T)−1 1

h=−¹2Tλ² und somit erf¨ullt der Fehler

yT/h−y(T) =h(−¹2Tλ²) +O(h²)

– 67–

Beweis.

hlim→0

e^−λT(1 +λh)^T/h−1 h

= lim

h→0e^{−λT d}_dheT/hln(1+λh)

= lim

h→0e^−λT(1 +λh)^Tλ/λh

−T

h²ln(1 +λh) + Tλ h(1 +λh)

= lim

h→0

1 2hT

− − λ (1 +λh)+ λ

(1 +λh)+ λ²h (1 +λh)²

= −¹2Tλ²

– 68–