Mathematik f¨ ur Informatiker III

Mathematik f¨ ur Informatiker III

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)

9. M¨arz 2006

Inhaltsverzeichnis

D Differentialgleichungen mit Numerik 3

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

D - 1.1 Numerische Grundaufgaben und ihre L¨osbarkeit . . . 5

D - 2 Gleitkommadarstellung und -arithmetik . . . 6

D - 2.1 Gleitpunktoperationen . . . 7

D - 2.2 Zu Grundidee (i) – Rundung von Zwischenergebnissen . . . 7

D - 2.3 Zu Grundidee (ii) – Fortsetzung der Berechnung trotz Fehlers . . . 8

D - 3 Summation numerischer Reihen . . . 9

D - 3.1 Fehlerfortpflanzung. . . 9

D - 3.2 Rundungsfehlerabsch¨atzung bei Riemann . . . 11

D - 3.3 Konvergenzbeschleunigung (1. Stufe nach Wijngaard) . . . 13

D - 3.4 Schlussfolgerungen aus dem Summationsbeispiel . . . 14

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

D - 4.1 Linearisierung durch Jacobimatrix . . . 15

D - 4.2 Newton’s Methode im Vektorfall . . . 16

D - 4.3 Lokale Konvergenz von Newton . . . 16

D - 5 Gew¨ohnliche Differentialgleichungen (ODE) . . . 17

D - 5.1 Separable Differentialgleichungen . . . 17

D - 5.2 Lineare Differentialgleichungen erster Ordnung . . . 18

D - 5.3 Lineare Differentialgleichungen n-ter Ordnung. . . 18

D - 5.4 Lineare Differentialgleichungen mit konstanten Koeffizienten. . . 19

D - 6 Euler Verfahren f¨ur Systeme von ODEs . . . 20

D - 6.1 Systeme von ODEs und ihre numerische L¨osung . . . 20

D - 6.2 Eulers Methode und andere explizite ODE-L¨oser . . . 20

D - 6.3 Runge-Kutta Verfahren der Ordnung 2 und 4 . . . 22

D - 6.4 Visualisierung der Verfahrensordnung . . . 23

D - 6.5 Numerische Integration von Systemen . . . 24

D - 6.6 Langzeitverhalten von ODE – L¨osungen . . . 24

D - 7 Interpolation mit Polynomen und Splines . . . 26

D - 7.1 Interpolation mit Polynomen (Whd. 1.Semester) . . . 26

D - 7.2 Interpolation durch kubische Splines . . . 27

D - 8 Numerische Integration – Quadratur . . . 29

D - 8.1 Interpolatorische Quadraturformeln . . . 29

D - 8.2 Quadratur mit Extrapolation – Romberg’s Verfahren. . . 30

E Lineare und nichtlineare Optimierung 32 E - 1 Grundlagen der Optimierung . . . 34

E - 1.1 (Nicht)lineare Ausgleichsprobleme . . . 34

E - 1.2 Allgemeine lineare Funktionenapproximation . . . 35

E - 2 Lineare Optimierung . . . 38

E - 2.1 Einf¨uhrendes Beispiel: Barkeeper . . . 38

E - 2.2 Lineare Optimierungsprobleme . . . 39

E - 2.3 Transformationen. . . 40

E - 2.4 Geometrische Untersuchung . . . 40

E - 2.5 Berechnung der optimalen Ecke – Simplex-Algorithmus . . . 42

E - 2.6 Dualit¨at . . . 44

E - 2.7 Komplementarit¨at . . . 46

E - 2.8 Lineare ganzzahlige Optimierung . . . 47

E - 3 Grundlagen der Optimierung . . . 53

E - 3.1 Nichtlineare Optimierungsprobleme mit Komplexit¨at . . . 53

E - 3.2 Nichtlineare Ausgleichsprobleme . . . 54

E - 3.3 Klassen von Optimierungsverfahren . . . 55

E - 3.4 Unrestringierte nichtlineare Optimierung. . . 57

E - 3.5 Restringierte Nichtlineare Optimierung . . . 62

F Kombinatorik und Wahrscheinlichkeitsrechnung 65 F - 1 Endliche Wahrscheinlichkeitsr¨aume . . . 67

F - 1.1 Elementare Definitionen . . . 67

F - 1.2 Bedingte Wahrscheinlichkeit. . . 70

F - 1.3 Unabh¨angigkeit von Ereignissen . . . 72

F - 1.4 Produktexperimente . . . 73

F - 1.5 Zufallsvariablen. . . 74

F - 1.6 Erwartungswert, Varianz, Kovarianz . . . 77

F - 1.7 Das schwache Gesetz der großen Zahlen . . . 82

F - 2 Unendliche Wahrscheinlichkeitsr¨aume. . . 85

F - 2.1 Diskrete Wahrscheinlichkeitsr¨aume . . . 85

F - 2.2 Kontinuierliche Wahrscheinlichkeitsr¨aume . . . 86

2

Teil D

Differentialgleichungen mit Numerik

Inhaltsangabe

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

D - 1.1 Numerische Grundaufgaben und ihre L¨osbarkeit . . . 5

D - 2 Gleitkommadarstellung und -arithmetik . . . 6

D - 2.1 Gleitpunktoperationen . . . 7

D - 2.2 Zu Grundidee (i) – Rundung von Zwischenergebnissen . . . 7

D - 2.3 Zu Grundidee (ii) – Fortsetzung der Berechnung trotz Fehlers . . . 8

D - 3 Summation numerischer Reihen. . . 9

D - 3.1 Fehlerfortpflanzung. . . 9

D - 3.2 Rundungsfehlerabsch¨atzung bei Riemann . . . 11

D - 3.3 Konvergenzbeschleunigung (1. Stufe nach Wijngaard) . . . 13

D - 3.4 Schlussfolgerungen aus dem Summationsbeispiel . . . 14

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

D - 4.1 Linearisierung durch Jacobimatrix . . . 15

D - 4.2 Newton’s Methode im Vektorfall . . . 16

D - 4.3 Lokale Konvergenz von Newton . . . 16

D - 5 Gew¨ohnliche Differentialgleichungen (ODE) . . . 17

D - 5.1 Separable Differentialgleichungen . . . 17

D - 5.2 Lineare Differentialgleichungen erster Ordnung . . . 18

D - 5.3 Lineare Differentialgleichungen n-ter Ordnung. . . 18

D - 5.4 Lineare Differentialgleichungen mit konstanten Koeffizienten. . . 19

D - 6 Euler Verfahren f¨ur Systeme von ODEs . . . 20

D - 6.1 Systeme von ODEs und ihre numerische L¨osung . . . 20

D - 6.2 Eulers Methode und andere explizite ODE-L¨oser . . . 20

D - 6.3 Runge-Kutta Verfahren der Ordnung 2 und 4 . . . 22

D - 6.4 Visualisierung der Verfahrensordnung . . . 23

D - 6.5 Numerische Integration von Systemen . . . 24

D - 6.6 Langzeitverhalten von ODE – L¨osungen . . . 24

D - 7 Interpolation mit Polynomen und Splines . . . 26

D - 7.1 Interpolation mit Polynomen (Whd. 1.Semester) . . . 26

D - 7.2 Interpolation durch kubische Splines . . . 27

D - 8 Numerische Integration – Quadratur. . . 29

D - 8.1 Interpolatorische Quadraturformeln . . . 29

D - 8.2 Quadratur mit Extrapolation – Romberg’s Verfahren . . . 30

Literaturverzeichnis

[Hartmann] Peter Hartmann, Mathematik f¨ur Informatiker. 3. ¨uberarbeitete Auflage, 2004, View- eg. Bei Lehmann’s vorhanden, ca. 30e.

Gute Grundlage, ¨ausserst lesbar, ISBN: 3-528-23181-5

[Mazzola] 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¨ur h¨ohere Anspr¨uche. Begleitender Kurs im Internet verf¨ugbar. ca 30e, ISBN: 3- 540-20835-6

[Opfer] Gerhard Opfer, Numerische Mathematik f¨ur Anf¨anger. Eine Einf¨uhrung f¨ur Mathe- matiker, Ingenieure und Informatiker. 4. durchgesehene Auflage, 2002, Vieweg [Roos+Schwetlick] Hans-G¨org Roos, Hubert Schwetlick, Numerische Mathematik. Das Grund-

wissen f¨ur jedermann. Mathematik f¨ur Ingenieure und Naturwissenschaftler. 1999, Teubner

[Stummel+Hainer] Friedrich Stummel, Karl Hainer, Praktische Mathematik. 1982, Teubner [Ortega+Rheinboldt] J.M. Ortega, W.C. Rheinboldt, Iterative solution of nonlinear equations in

several variables. 1970 Academic Press, Inc.

[Stoer] 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

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

Ausgangsdilemma

Die Modellierung natur- oder sozialwissenschaftlicher Zusammenh¨ange bzw ’Systeme’ f¨uhrt zu mathematischen ’Gleichungen’, die nur in ganz einfachen F¨allen per Hand oder sonstwie ’exakt’

gel¨ost werden k¨onnen.

Zum Beispiel k¨onnen schon bei der unbestimmten IntegrationMapleundMathematicanur in speziellen Ausnahmef¨allen eine L¨osung als Formel angeben.

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

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¨aherungsweise ¨uber dem Raster(=Screen) der Gleitkommazahlen gel¨ost

Die Ergebnisse werden ausgedruckt oder besser graphisch dargstellt.

Stufen des ’Wissenschaftlichen Rechnens’

1. Modellierung ( des Anwendungssystems )

2. Diskretisierung ( von Differentialgleichungen )

3. Dateneingabe ( f¨ur aktuelle Situation )

4. L¨osung ( durch Gleitkomma-Algorithmen )

5. Datenausgabe ( in geeigneter Form )

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

D - 1.1 Numerische Grundaufgaben und ihre L¨ osbarkeit

Lineare algebraische Gleichungssysteme

Im Prinzip v¨ollig im Griff. Variablenzahl jeweils durch Speichergr¨osse und Prozessorzahl und - geschwindigkeit beschr¨ankt.

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¨ollig im Griff unabh¨angig von Linearit¨at.

Randwertaufgaben f¨ur ODEs

Standarddiskretisierung f¨uhrt auf lineare bzw nichtlineare algebraische Gleichungen und ist ent- sprechend l¨osbar.

Partielle Differentialgleichungen PDE

Nur im elliptischen Fall schnell l¨osbar, alles andere ist Forschungsgebiet und st¨osst jeweils an die Grenzen vorhandener Rechnerkapazit¨aten.

Warnung

Alles wird beliebig viel schwieriger wenn

• einige Variablen ganzzahlig sein m¨ussen und / oder

• die L¨osung gegebenen Ungleichungen gen¨ugen muss wie in der Optimierung ¨ublich.

D - 2 Gleitkommadarstellung und -arithmetik

Ein System von Gleitkommazahlen wird definiert durch:

• Basis (oder Radix)b(= ¨ublicherweise 2)

• Mantissenl¨angel

• Minimaler Exponentemin

• Maximaler Exponentemax

Teilmenge der reellen Zahlen Rmit Darstellung x= −1s

0.m1m2 · · · ml

| {z }

Mantissem

b^e∼ −1s

m1b^e−1+m2b^e−2+m3b^e−3+. . .+mlb^e−l

Vorzeichenbit s, Mantisse m, Exponente s∈

0,1 mi ∈ {0, 1, . . . , b−1} e∈ {emin, emin+ 1, . . . , emax} Bin¨ardarstellung, d.h. Basisb= 2

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

Arten von Gleitkommazahlen

• normalisierte Gleitpunktzahl:

m1 > 0 =⇒ 1

b ≤ m ≤ x b^−e

< 1

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

• unnormalisiert:m1 = 0 zugelassen =⇒ keine Eindeutigkeit

• denormalisiert:m1 = 0, e = emin

Vorsicht:

Rechnen mit denormalisierten Zahlen f¨uhrt zu verst¨arkten Rundungseffekten.

Betragsm¨assig kleinste normalisierte Zahl TINY TINY= 0.1·b^emin =b^emin−1

Betragsm¨assig gr¨oßte normalisierte ZahlHUGE

HUGE= 0.(b−1)(b−1)(b−1). . .(b−1). . . b^emax=b^emax(1−b^−l) Epsilon (relative Maschinengenauigkeit) ε

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

• Mantissenl¨angel bestimmt die Rechengenauigkeit.

• Exponentenbereichemax−emin bestimmt den Wertebereich.

Beispiel D.1(Gleitpunktzahlsystem mit Basis2und Mantissenl¨ange3).

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.

6

PSfrag replacements

x= 0.m1m2m32^e Exponentenbereich−1≤e≤1

Normalisierte positive Zahlen: m¹= 1, m²∈ {0,1} 3m³

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

v⁵ =v³∗v⁴ 0

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

e −1 −1 −1 −1 −1 −1 −1 −1 m¹

m2

m³

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

D - 2.1 Gleitpunktoperationen

Bemerkenswert

( 1.0 / 8.0 ) * 8.0 = 1.0[1em]( 1.0 / 5.0 ) * 5.0 6= 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.

Allgemein g¨ultiger Standard: ANSI - IEEE 754

(ANSI →American National Standards Institute und IEEE → Institute of Electrical and Elec- tronics Egineering.)

Grundideen:

1. Alle Zwischenergebnisse werden zur n¨achsten Gleitpunktzahl gerundet.

2. The show must go on. Auch bei Fehlern wird weiter gerechnet.

D - 2.2 Zu Grundidee (i) – Rundung von Zwischenergebnissen

Auch wennx undy im Gleitpunktbereich liegen, gilt dies im Allgemeinen nicht f¨ur das Ergebnis x◦y, wobei ◦ ∈ {−,+,·, /}. Dann wirdx◦y zun¨achst mit erh¨ohter Genauigkeit berechnet und anschließend zur n¨achstliegenden 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

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

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

Alternative Schreibweise:

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

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

Konsequenz f¨ur relativen Fehler:

f l(x◦y)−(x◦y) x◦y

≤ |ε| ≤ eps

Warnung:

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

Algorithmen (z.B. zur Matrixfaktorisierung) m¨ussen deswegen auf ihre Stabilit¨at, d.h. die Verst¨arkung oder Abd¨ampfung von Rundungsfehlern, untersucht werden.

Beispiel D.4. Gausssche Elimination ohne Pivotierung ist extrem instabil.

Gauss mit Pivotierung ist dagegen recht stabil.

Frage

Was passiert, wennx◦y außerhalb des Wertebereichs[-HUGE, HUGE]liegt, d.h. entweder∇(x◦y) oder ∆(x◦y) nicht existiert?

BeispielD.5(Programm). REAL u,s,t s = TINY(u)**2 ! ergibt 0 t = HUGE(u)*8

! ergibt INF, signalisiert OVERFLOW

D - 2.3 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¨ur allex 6= -INF x * INF == sign(x) * INF f¨urx 6= 0 x / 0 == sign(x) * INF f¨ur x 6= 0

wobeisign(x)das Vorzeichen von x liefert.

Undefinierte Operationen wie0/0,INF/INF,INF-INFund0*INFergeben den sehr speziellen Wert NaN≈Not a Number.

Da einNaNnicht mit sich selbst oder etwas anderem verglichen werden kann, gilt x 6= x .EQUIV. .TRUE.

genau dann wennxeinNaNist.

Infektionsprinzip:

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

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

8

D - 3 Summation numerischer Reihen

D - 3.1 Fehlerfortpflanzung

Erinnerung:

f l(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.

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 Programm ergibt f¨ur n= 100 undx= 2.0/3.0

s check s/check - 1 n * eps

3.0000002 3.00000019 2·10⁻⁸ 1.2·10⁻⁵ Beobachtungen

• Gleitpunktwert von xist 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

• Der beobachtete relative Fehler zwischen einfach und doppelt genauer L¨osung ist lediglich 2·10⁻⁸, d.h. von der Gr¨oßenordnung der Maschinengenauigkeit, obwohl wir 100 Operationen durchgef¨uhrt haben. Die Rundungen scheinen sich partiell aufgehoben zu haben.

• Eine exakte Absch¨atzung f¨ur den worst case (d.h. schlimmster Fall) ergibt den Wert (1 + eps)¹⁰⁰≈100·epsals relativen Fehler. Das l¨asst sich wie folgt herleiten.

Theoretische Schranke des Fehlers im obigen Programm F¨uryi+1=f l(yi∗x) als berechneter Wert vony imi-ten Schritt gilt:

y0 = 1 y1 =x

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

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

...

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

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

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

s1 =f l(y0+y1) = f l(1 +x) = (1 +x)(1 +εn+1) s2 =f l(s1+y2) = f l(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

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.

Erkl¨arung:

Betrachte kleinen Summanden y und großen Summanden x = 0.m1m2. . . ml·2^e so dass x = x+ 2^−l+e die n¨achst gr¨oßere Gleitpunktzahl zu x ist und x =x−2^−l+e ist die n¨achst kleinere Gleitpunktzahl zux.

PSfrag replacements

2^e−1 x

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

2^e

x x

Konsequenz:

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

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

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

x=

n−1

X

i=1

1 i &

Z n 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)

10

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 (sneu 6= 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 ≈ ¹6 Sekunde

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

Vergleich zur theoretischen Herleitung

n= 2097152 ergibt ln(n)∗n∗EP SILON(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¨auft ewig, da eps⁻¹ und damit dann auch n um 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.

D - 3.2 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

∆ζ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^1−x 1−x

∞ k=n

= 1

k^1−x(1−x)

∞ k=n

= 0− 1

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

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

⇒ n ≥ ^x−1 s 1

tol(x−1)

Partialsummen:

ζn(x) = Pn k=1

1

k^x wachsen monoton mit n und sind nach oben durch _x−^x₁ beschr¨ankt, haben also einen eindeutigen Grenzwertζ(x).

Praktische Notwendigkeit: Diskretisierung

Hier, wie h¨aufig in numerischer Mathematik muss mathematisches Problem durch Ausf¨uhrung endlich vieler Operationen auf endlich vielen Variablen ann¨aherungsweise gel¨ost werden. Hier ein- fach Ann¨aherung vonζ(x) durchζn(x). Der entsprechende Abbruchfehler|ζ(x)−ζn(x)|kann hier einfach mit Hilfe einer Integralschranke abgesch¨atzt werden. Unabh¨angig vom in der Numerischen Analysis betrachtetenDiskretisierungsfehlerist der Rundungsfehler zu ber¨ucksichtigen.

Rundungsfehlerabsch¨atzung F¨urbi>0

f l . . . b1+b2

+b3

+b4

. . .+bn+1

+bn

= . . . b1+b2

1 +ε1 +b3

1 +ε2 +b4

1 +ε3

. . .+bn

1 +εn−1

= b1 1 + ˜ε1

n−2

+b2 1 + ˜ε1

n−1

+b3 1 + ˜ε2

n−2

+. . .+bn 1 + ˜εn−1

1

=⇒

f l b1+. . .+bn

− b1+b2+. . .+bn

≤ b1

h 1 +epsn−1

−1i +b2

h 1 +epsn−1

−1i

+. . .+bn 1 +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.

Schlussfolgerung:

Um Rundungsfehler zu minimieren sollten Summen m¨oglichst vom kleinsten zum gr¨oßten Sum- manden gebildet werden. Bei konvergenten (hoffentlich monoton fallenden) Reihen sollte von hin- ten, d.h. r¨uckw¨arts 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¨arts bis. liegen bleiben n= 4097 1.6446900 r¨uckw¨arts vom gleichen n= 4097 1.6449339 r¨uckw¨arts mit n= 2²³= 8388608 Bemerkung:

Durch R¨uckw¨artssummation k¨onnen deutlich mehr Summanden der Form 1/k⁻x mit n > 4097 ihren Beitrag zur Gesamtsumme leisten. Mehr Summanden zu benutzten bedeutet aber, denDis- kretisierungsfehler zu verringern und damit den exakten Wertζ(x) besser zu approximieren.

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

Vergleich:

eps·n·π²

6 eps·ln(n)

D - 3.3 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 als P

bk. Die Korrektur erfolgt durch transformierte Terme

aj= X∞

i=1

bj2ⁱ

2ⁱ.

Satz D.8. Satz: F¨urbk=k^−xoder andere monoton konvergierende Reihen gilt im Grenzwert X∞

k=1

bk= X∞

j=1

(−1)^j−1aj .

Bemerkung

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

Idee des Beweises:

Betrachte, wie oftbk in aj auftritt

V orz 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.

D - 3.4 Schlussfolgerungen aus dem Summationsbeispiel

• Die Behandlung mathematischer und anderer Modellierungsprobleme bedingt das Auftreten von Abbruchs- ≡Diskretisierungsfehlern sowie Rundungsfehlern. Beide sollten abgesch¨atzt und m¨oglichst minimiert werden.

• Gleitpunktarithmetik ist weder kommutativ noch assoziativ, distributiv usw. Spezielle Konsequenz:Betragsm¨aßig fallende Reihen von hinten summieren!

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

• Viele Jobs (≡Programme, Daten) laufen entweder im Sekunden- oder Stundenbereich. Be- obachtung der Abarbeitung im Minutenbereich ist relativ selten.

• Mathematisch endlichist nicht gleich rechentechnisch endlich.

14

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

Methoden zur L¨osung des linearen Problemes Ax=b mitdim(x) = dim(b) =n

• Cramersche Regelxi= (−1)ⁱdet(Ai)/det(A) f¨uri= 1..n ( In Ai wird diei−te Spalte vonAdurchb ersetzt )

• Gauss-Elimination≈P A = LU Faktorisierung

(P Permutation,L unterhalb undU oberhalb dreiecksf¨ormig )

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

• Fixpunkt Iteration x←x−M F(x) mitF(x)≡Ax−b (M∈R^n×n angen¨aherte Inverse so dassM A≈I ) Hinweise:

• F¨ur (eindeutige) L¨osbarkeit ist ¨uberalldet(A)6= 0 vorrauszusetzen.

• L¨oseLU x=bbzwQRx=bdurch Substitution/Transponierung.

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

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) +^√1₃ = 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∗x1∗x3 0 _1+0.1^x_∗¹_x1∗x3 z(x) _x¹₃ + 0.1∗x1 −^x2−9.8_x²₃^∗x1





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

1 x3 +¹₅x1

D - 4.1 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¨uglich jeder der Variablenx1, . . . , xn Lipschitz-stetig sind, so l¨asst sich aus dem Hauptsatz der Differential- und Integralrechnung herleiten, dass f¨ur jeden Schritts∈Rⁿ gilt

kF(x+s) − [F(x) + F⁰(x)s]k ≤ γksk²

Hierbei istF⁰(x)s ein Matrix-Vektor Produkt und k · 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 ) vonF an der Stelle x.

D - 4.2 Newton’s Methode im Vektorfall

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

As=b mit A=F⁰(x) und b=−F(x) Die L¨osung l¨asst sich ausdr¨ucken als

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

Wiederholte Berechnung von s und anschliessende Inkrementierungx ← x+s ergibt 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.

Warnung:

• Das Verfahren muss abgebrochen werden wenn det(F⁰(x^(k))) null oder sehr klein ist.

• Im letzteren Falle werden die Schritte s^(k) typischerweise sehr gross und f¨uhren h¨aufig zu Argumentenx^(k+1)woF garnicht mehr ausgewertet werden kann.

• Zur Vermeidung dieses Problems wirds^(k) manchmal mit einem D¨ampfungsfaktorα^(k)<1 multipliziert, der dann Schrittweitegenannt 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).

D - 4.3 Lokale Konvergenz von Newton

Satz D.9 (Satz von Kantorovich). Sei die VektorfunktionF :Rⁿ →Rⁿ einmal differenzierbar und besitze ihre Jacobimatrix F⁰(x)∈R^n×n die Lipschitzkonstante γ.

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

F⁰(x⁽⁰⁾)⁻¹

2

F(x⁽⁰⁾) ≤ 1

2γ

dass Newton’s Methode zu einer L¨osungx^(∗) mitF(x^(∗)) = 0 konvergiert.

Die Konvergenzgeschwindigkeit ist quadratisch in dem Sinne dass f¨ur eine Konstantecund allek gilt

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

x^(k)−x^(∗)

2

Bemerkung:

Je nichtlinearer ein Problem umso gr¨osser istγund desto st¨arker ist damit die Bedingung anx⁽⁰⁾. Wird praktisch nie ¨uberpr¨uft !!!!

16

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 Variablen x und einer gesuchten Funktion y =y(x) auch deren Ablei- tungen _dx^dnyn =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.

D - 5.1 Separable Differentialgleichungen

Definition D.11 (Separable Differentialgleichung).Eine DifferentialgleichungF(x, y, y⁰) = 0 erster Ordnung heisstseparabel, wenn sie sich in der Form

y⁰=f(x)g(y)

darstellen l¨asst, wobei f : I −→ R, g : J −→ R stetige Funktionen auf den Intervallen I ⊆R, J ⊆Rsind.

Satz D.12 (L¨osbarkeit: Anfangswertproblem separabler ODE). Eine separable Differen- tialgleichung erster Ordnung mit der Anfangsbedingung y(x0) = y0 f¨ur x0 ∈ I, y0 ∈ J, hat im IntervallJ eine eindeutige L¨osungy(x) :I −→J, falls

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

Seien

G(y) :=

Z y y0

1

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

Z x x0

f(x)dx die Stammfunktionen von _g(y)¹ bzw.f(x).

Dabei wurden f¨ur Integrationsvariable und Obergrenze der Integration das gleiche Symbol ver- wendet.

Auf J ist G⁰(y) = _g(y)¹ 6= 0 (Voraussetzung SatzD.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. 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 Anfangswertproblemy⁰(x) =f(x)g(y), mit Funktionenf :I −→R, g:J −→R, und dem Anfangswerty(x0) =y0∈J, hat die eindeutige L¨osungy, die man erh¨alt, wenn man die folgende Gleichung nach y aufl¨ost:

Z y y0

1 g(y)dy=

Z x x0

f(x)dx

D - 5.2 Lineare Differentialgleichungen erster Ordnung

Definition D.14 (Lineare Differentialgleichung). Differentialgleichungen, bei denen die Funk- tion y =y(x) und ihre Ableitungen nur in linearem Zusammenhang auftreten heissen Lineare 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 Gleichung homogen, sonstinhomogen.

Die FunktionF(x) auf der rechten Seite heisstQuellfunktion.

Satz D.15 (L¨osung homogener linearer ODE). Ista(x)auf dem IntervallI stetig, so lautet die vollst¨andige L¨osung der linearen Differentialgleichungy⁰+a(x)y= 0

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

Satz D.16 (L¨osung inhomogener linearer ODE). Die inhomogen lineare Differentialglei- chungy⁰+a(x)y=f(x), f, a:I −→Rstetig,x0∈I, besitzt die vollst¨andige L¨osung

y= Z x

x0

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

·e^−A(x) wobeic∈Rund A(x)eine Stammfunktion von a(x)ist.

D - 5.3 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.

Dieai heissen Koeffizientenfunktionen,f heisst Quellfunktion.

Istf = 0, so heisst die Gleichung homogen, sonstinhomogen.

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 mitai, f :I −→R undx0∈I. Dann gibt es zu den Anfangswerten

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

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

Diese L¨osung existiert auf dem ganzen IntervallI.

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 = 0mit ai:I −→Rbildet einen reellen Vektorraum der Dimension n.

Eine Basis des L¨osungsraumesH nennt manFundamentalsystem.

Jede L¨osungy der inhomogenen Gleichungy⁽ⁿ⁾+a1(x)y⁽ⁿ⁻¹⁾+· · ·+an−1(x)y⁰+an(x)y = f(x) mitf :I −→Rhat die Form

y=ys+yh

wobei xh∈ H eine L¨osung der homogenen und ys eine spezielle L¨osung der inhomogenen Diffe- rentialgleichung ist.

18

D - 5.4 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 Fall konstanter Koeffizientenfunktionen ai(x) ∈R kann 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 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. . . nmit p(λi) = 0 des charakteristischen Polynoms kann ein Funda- mentalsystem f¨ur die homogene Differentialgleichung n-ter Ordnung konstruiert werden.

Dazu ist eine Fallunterscheidung nach derVielfachheit der Nullstellenλi n¨otig:

λ∈R ist einfache Nullstelle

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

λ=α+iβ ∈Cist einfache komplexe Nullstelle

e^{α x}cosβ x und e^{α x}sinβ x sind L¨osungen der Differentialgleichung.

λ∈R istk-fache reelle Nullstelle

xⁱe^{λ x}, i= 0, . . . , k−1 sindklinear unabh¨angige L¨osungen.

λ=α+iβ ∈Cist k-fache komplexe Nullstelle

xⁱe^{α x}cosβ x, xⁱe^{α x}sinβ x, i= 0, . . . , k−1 sind die 2klinear unabh¨angige L¨osungsfunktionen.

Beispiel D.21. SieheHartmann, Mathematik f¨ur Informatiker, S.352 ff.

D - 6 Euler Verfahren f¨ ur Systeme von ODEs

D - 6.1 Systeme von ODEs und ihre numerische L¨ osung

In vielen Anwendungen wird der Zustand eines Systems zum Zeitpunkttdurch einen Vektor x(t) = [x1(t), x2(t), . . . , xn(t)]^> mit n >0

beschrieben. Die ¨Anderungsgeschwindigkeit ˙x≡dx(t)/dtdes Zustandes nach der Zeit ergibt sich h¨aufig als FunktionF(x(t)) mitF:Rⁿ→Rⁿeben dieses Zustandes. Also erhalten wir das System gew¨ohnlicher Differentialgleichungen

˙

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

Das System heisst autonom, da die Zeitt auf der rechten Seite nicht explizit, sondern nur mit- telbar ¨uber x = x(t) vorkommt. Dieses ist keine Einschr¨ankung da ein nichtautonomes System

˙

x(t) = F(t, x(t)) sich autonom umschreiben l¨asst indem man t als nullte Zustandskomponente x0(t) hinzuf¨ugt und somit f¨ur ¯x≡(x0, x1, . . . , xn)^T erh¨alt

d dtx¯ ≡

x˙0

˙ x

= t˙

˙ x

= 1

F(¯x)

≡ F(x)

Auch ODEs h¨ohere Ordnungen lassen sich in Systeme von ODEs erster Ordnung umschreiben, indem man z.B. die erste Ableitung y⁰ als neue abh¨angige Variable v ≡ y⁰ definiert und dann y⁰⁰ 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≡y undy2≡y⁰



 y₀⁰ y₁⁰ y₂⁰



 =



 1 y2

f(y0, y1, y2)





Entsprechend lassen sich Anfangsbedingungen umschreiben.

Die Umformulierung als System 1.Ordnung er¨offnet die M¨oglichkeit numerische Standardmethoden und Software f¨ur die L¨osung autonomer Systeme erster Ordnung mit Anfangsbedingungen zur Anwendung zu bringen.

Satz D.22 (Existenz und Eindeutigkeit der L¨osung). Sei F : D ⊂ Rⁿ −→ Rⁿ in einem offenem Gebiet Dlokal Lipschitz-stetig.

Dann existiert f¨ur jeden Punktyo∈ Dein Intervall(a, b)30und eine eindeutige L¨osungy(t)∈ D der ODEy˙=F(y) f¨ura < t < bmity(0) =y0.

Bemerkung:

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

2. Das Intervall (a, b) kann so gross gew¨ahlt werden, dass y(b) den Rand vonDerreicht.

D - 6.2 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 Standard- funktionene^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.

Explizite (Vorw¨arts) Euler-Methode

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

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 derTan- gentey(t) der˙ L¨osungy(t) in tk

y(0) =y⁰

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

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

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)^Th Vergleich mit exakter L¨osung:

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

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

n→∞

1 +λT

n n

Erl¨auterung

Die angen¨aherte L¨osungyT /hkonvergiert gegen die exakte L¨osungy(T) der ODE wenn die Schritt- weiteh=T /n gegen 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²)

Beweis.

hlim→0

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

= lim

h→0e^{−λT d}_dhe^{T /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 λ²

Folgerung D.24 (Approximationsfehler der Euler-Methode). F¨ur alle Lipschitz-stetigen Probleme (d.h. die rechte SeiteF(t, y,y)˙ der ODE ist Lipschitz-stetig) liefert das Euler-Verfahren eine numerische L¨osung mit

yT /h−y(T) = c(T)h+O(h²).

Deshalb nennt man diese Methode auch Verfahren erster Ordnung:

Die Verdopplung der Approximationsgenauigkeit durch Halbierung der Schrittweitehverdoppelt den Berechnungsaufwand.

Frage:

Gibt es Verfahren der Fehlerordnungpso dass

kyn−y(T)k=c(T)h^p+O(h^p+1)

gilt und damit die Halbierung der Schrittweite hzu einer Reduktion des Fehlers um den Faktor (¹₂)^p f¨uhrt ?

Anwort: JA!

p=2 Mittelpunkt - Regel oderHeun’sches Verfahren p=4 Runge-Kutta 4. Ordnung

p=5 Runge-Kutta-Fehlberg

D - 6.3 Runge-Kutta Verfahren der Ordnung 2 und 4

Mittelpunkt-Regel

• tk+1/2 = tk+ 0.5hk; tk+1 = tk+hk

• yk+1/2 = yk+ 0.5hkf(tk, yk)

• yk+1 = yk+hkf(tk+1/2, yk+1/2) Runge-Kutta 4 (Standardwahl)

• tk+1/2 = tk+ 0.5hk; tk+1 = tk+hk

• yk+1/4 = yk+ 0.5hkf(tk, yk)

• yk+1/2 = yk+ 0.5hkf(tk+1/2, yk+1/4)

• yk+3/4 = yk+ hkf(tk+1/2, yk+1/2)

• yk+1 = yk+^h₆^k

f(tk, yk) + 2f(tk+1/2, yk+1/4) + 2f(tk+1/2, yk+1/2) +f(tk+1, yk+3/4)

22

D - 6.4 Visualisierung der Verfahrensordnung

F¨ur einen beliebigen numerischen Integrator folgt aus der vorrausgesetzten Beziehung kyT /h−y(T)k = c(T)h^p+O(h^p+1) ≈ c(T)h^p

durch Logarithmierung, dass

−log kyT /h−y(T)k

≈ p(−log(h))−log(c(T))

Die linke Seite ist ein Maß der korrekt berechneten Dezimalstellen in der L¨osung. Sie ist nun ann¨aherungsweise eine affine Funktion von−log(h) also eine Gerade, deren Steigung gerade die Ordnungpder Methode ist. [0.3cm]

Um die Ordnung eines Verfahrens zu pr¨ufen kann man die Schrittweite zum Beispiel wie hk = T /2^k f¨ur k = 1,2. . . variieren und die entsprechenden Fehler −log kyT /h_k−y(T)k

¨uber den Abzissenwerten−log(hk) =klog(2)−log(T) auftragen.

0 5 10 15 20 25 30

0 2 4 6 8 10 12

Euler Midpoint RK-4

Frage:

Wie kann die Schrittweite in Hinblick auf den gesch¨atzten Fehler gew¨ahlt werden?

Antwort:

Durch Vergleich der Ergebnisse f¨ur verschiedene Schrittweitenhoder verschiedener Methoden.

Beispiel D.25(Mittelpunkt - Regel).

yn = y(T) +c(T)h²+O(h³) y2n = y(T) +c(T)¹₄h²+O(h³)

=⇒ yn−y2n = c(T)³₄h²+O(h³)

=⇒ c(T) ≈ ⁴3

yn−y2n

h² ≡ c(T˜ )

=⇒ ky2n−y(T)k ≈ ⁴3ky2n−ynk ist eine Fehlerabsch¨atzung f¨ur die Mittelpunktregel.

Folgerung D.26 (Einfache Schrittweitensteuerung). Wenn die numerische L¨osung mit einer absoluten Genauigkeit vonτ >0gew¨unscht wird, dann w¨ahlt man bei der Mittelpunktsregel

h=p² τ /˜c(T) Allgemeiner empfiehlt sich f¨ur ein Verfahren der Ordnung p

h=p^p τ /˜c(T)

Hierbei ist die Fehlerkonstante˜c(T)STARK vom Verfahren abh¨angig.

Nimmt man dennoch an, dass f¨ur Euler, Mittelpunkt und Runge-Kutta 4 die c =c(T) ¨ahnlich gross sind, so ergeben sich Rechenaufw¨ande von

1·c/τ, 2·p

c/τ , 4·p⁴ c/τ

Auswertungen der rechten Seite. Bei gr¨osserer geforderter Genauigkeit, also kleineremτ sind Ver- fahren h¨oherer Ordnung zu bevorzugen,vorrausgesetzt die rechte Seite der ODE ist pmal diffe- renzierbar.

D - 6.5 Numerische Integration von Systemen

Runge-Kutta Methoden sind direkt auf Systeme

˙

y(t) =f(y(t))∈Rⁿ bzw y(t) =˙ f(t, y(t))∈Rⁿ

anwendbar. W¨ahrend die unabh¨angige Variabletund die entsprechenden SchrittweitenhSkalare bleiben, sind alle anderen Gr¨ossen jetzt Vektoren der L¨angen.

Die Euler Rekursion

yk+1=yk+hkF(tk, yk)∈Rⁿ

erfordert also das h-fache des RichtungsvektorsF(tk, yk)∈ Rⁿ zu dem alten Zustandsvektoryk

zu addieren, um den neuen Zustandsvektor yk+1 ∈ Rⁿ zu erhalten. Es ist davon auszugehen, dass diese Vektormultiplikation und -addition vom Aufwand her gegen¨uber der Auswertung der Rechten SeiteF(t, y) vernachl¨assigbar ist.

Die Konvergenzordnungen bleiben erhalten, wobei der Abstand zwischen der ann¨ahenden und der genauen L¨osung jetzt als eine VektornnormkyT /h−y(T)kder Differenz zwischenyT /h undy(T) zu bestimmen ist.

Lineares Beispiel f¨ur Euler

Das autonome System linearer Differentialgleichungen x(t)˙

˙ y(t)

=

−y(t) x(t)

mit

x(0) y(0)

= 1

0

hat die analytische L¨osung [x(t), y(t)] = [cos(t), sin(t)]. Die Anwendung der Eulermethode mit Schrittweite h ergibt

xn+1

yn+1

= xn

yn

+h

−yn

xn

=

xn−hyn

yn+hxn

=

1 −h

h 1

xn

yn

=ρ

cos(α) −sin(α) sin(α) cos(α)

xn

yn

=ρⁿ

cos(nα) −sin(nα) sin(nα) cos(nα)

x1

y1

wobeiρ≡√

1 +h² undα= arcsin(h/p

1 +h²) .

D - 6.6 Langzeitverhalten von ODE – L¨ osungen

Bemerkung zum Langzeitverhalten

H¨aufig ist von Interesse (z.B. in der Klimavorhersage), wie sich L¨osungeny(t) der ODE ˙y=F(y) f¨ur sehr grosset qualitativverhalten, und zwar unabh¨angig vom Anfangswerty(t0) =y0. D.h. man will wissen, ob das dynamische System sich einschwingt, einen Gleichgewichtszutand erreicht, zuf¨alliges (d.h. chaotisches) Verhalten o.¨a. zeigt.

Im folgenden machen wir Aussagen f¨ur autonome Systeme der Zustandsraumdimension n, die entspechend auch f¨ur nichtautonome Systeme der Dimension n−1 gelten.

(I) Fallsn= 1muss und sonst (n >1) kann einer der beiden folgenden F¨alle eintreten:

1. y(t) strebt einem station¨aren Grenzwert y_∞= lim

t→∞y(t) zu Beispiel: ˙y=λ(y−a), a∈R, λ <0, y0beliebig

24

t y

y_∞

y(t) =c e^λt+a, c <0 y(t) =c e^λt+a, c >0

(b) y(t) explodiert (blow up)

tlim→t^∗ky(t)k=∞ f¨ur endliche Zeitt^∗ (kritische Zeit)

Beispiel: y˙ =y² mit y(0) =y0>0

=⇒ dy

y² =dt =⇒ Z 1

y²dy= Z

dt =⇒ −1

y =t+c =⇒ y(t) =− 1 t+c AW: y0=−1

c >0

=⇒ c= −1 y0

<0

=⇒ y(t) = 1

1 y0 −t

t y

t^∗ y(t) = 1¹

y0−t

(II) Asymptotisch periodische L¨osung

Falls die Zustandsdimensionn= 2 ist muss, ansonsten kanny(t) sich asymptotisch einer periodi- schen L¨osungy_∗(t) n¨ahern, f¨ur die gilt

y_∗(t+T) =y_∗(t) f¨ur allet >0 und feste PeriodeT.

Beispiel: siehe obigesLineares Beispiel f¨ur Euler

(III) Chaotisches Verhalten

Falls Dimensionn > 2 (einschliesslich n= 2 im nichtautonomen Fall) kann die L¨osungy(t) der ODE sich chaotisch verhalten, d.h. auch nach sehr langer Zeit l¨asst sich keine periodische oder station¨are Struktur erkennen.

Beispiel: Lorenz - Attraktor ( ¨Ubung 2)