h
3.5 Imagin¨arschrittableitung
Tats¨achlich kann man Differenzen (und damit Rundungsfehler bzw. Weghebeph¨anomene) zumindest bei ersten Ableitungen mit einem Trick v¨ollig vermeiden: F¨ur einen rein imagin¨aren Schritt ih lautet die Taylorreihe:
f(x+ih) =f(x) +ihf0(x)− h2
2 f00(x)− ih3
3! f000(x) +· · · (28)
Division durch h und Betrachtung nur des Imagin¨arteils auf beiden Seiten liefert:
f0(x) = Im[f(x+ ih)]
h + O(h2) (29)
Vorteile der Ableitungsformel Gl. 29:
• Fehlerordnung h2 (statt h1 wie bei simplen Vorw¨artsdifferenzen), bzw. selbe Fehlerordnung wie bei Zentraldifferenzen, aber mit nur einer Funktionsberechnung;
• keine differenzbasierten Rundungsfehler/Weghebung, damit sehr viel kleinere h-Werte m¨oglich;
• bei Sprachen mit “eingebauten” komplexen Variablen und Funktionen (wie z.B. Fortran) ist die Implementation sehr einfach.
3.6 Schlußbemerkungen
• H¨ohere Ordnung bewirkt nicht immer h¨ohere Genauigkeit.
• H¨ohere Ordnung bedeutet nicht unbedingt bessere Effizienz.
• Funktionsberechnungen k¨onnen teuer sein.
Beispiel: Gradienten/Frequenzen in der ab-initio-Quantenchemie.
Goldene Praxisregel: keine finiten Differenzen verwenden, sonderen analytische Ableitungen. Diese m¨ussen jedoch zus¨atzlich einprogrammiert werden. Das wird trotzdem gemacht ⇒ das lohnt sich! Beispiel: Anilin:
– eine SCF-Energieberechnung dauere 5 min
– bei 14∗3=42 Freiheitsgraden und Zentraldifferenzen (2 Punkte pro Freiheitsgrad) braucht der komplette Gradientenvek-tor
42∗2∗5 min = 420 min = 7 h
– ⇒ eine konvergierte Geometrieoptimierung braucht 1 Tag!
– mit “analytischen Gradienten” braucht der komplette Gradientenvektor nur ca. 5–10 min und die konvergierte Geome-trieoptimierung ca. 1 h.
• Es existieren weitere Alternativen zur Ableitungsberechnung, z.B. via Fouriertransformation:
dn
4 Integration
4.1 Vorbemerkungen:
Wann brauchen wir numerische Integration?
• nicht bei analytisch integrierbaren Funktionen, sondern
• bei analytisch gegebenen, aber nicht analytisch integrierbaren Funktionen,
• bei nur numerisch gegebenen Funktionen.
Einzige Voraussetzung an zu integrierende Funktion f(x):
f(x) muß berechenbar sein, f¨ur einige diskrete, vorgegebene Werte von x innerhalb des Integrationsintervalls [a, b].
Grundprinzip:
Zb
a
f(x)dx =
N
X
i=1
cif(xi) + Fehler(N, . . .) (31)
• unabh¨angig von der Form von f(),
• Lage der St¨utzstellen xi bestimmt Werte der Koeffizienten ci,
• N m¨oglichst klein,
• Fehler f¨ur gegebenes N m¨oglichst klein und absch¨atzbar, soll f¨ur wachsendes N schnell gegen Null gehen.
Alle Integrationsverfahren dieses Typs machen implizite Annahmen ¨uber den Verlauf von f(x) zwischen den Integrationsst¨ utz-stellen xi. Das ist
• einerseits ein Grundproblem bei allen numerischen Algorithmen, die mit Diskretisierungen arbeiten;
• andererseits theoretisch und praktisch nicht ganz so gef¨ahrlich wie es scheint.
4.2 Trapezregel (Euler-Verfahren)
(bei ¨aquidistanten St¨utzstellen: Newton-Cotes-Verfahren)
a x b
In jedem Intervall der Breite h wird f(x) durch Gerade (Polynom 1.Ordnung) approximiert. Daraus folgt:
• die Formel ist exakt nur dann, wenn f(x) eine Gerade ist;
• Wenn wir f(xi+1) durch eine Taylorreihe um f(xi) approximieren:
f(xi+1) = f(xi) + hf0(xi) + h2
2!f00(xi) +· · · (36) brechen wir hier also nach dem linearen Glied ab
⇒ der Fehler ist von der Ordnung O 1/N2
= O(h2) Kontrolle der Konvergenz in der Praxis:
• einmalige Berechnung des Integrals (mit einem geratenen Wert f¨ur h) ist gef¨ahrlich; nur mit weiteren Informationen zu rechtfertigen (z.B. hinreichend genaue Resultate mit diesem h-Wert bei sehr ¨ahnlichen Funktionen); sonst immer n¨otig:
• mehrmalige Berechnung des Integrals mit jeweils kleiner werdenden h-Werten, bis ¨Anderung der Integralapproximation kleiner als vorgegebene Schranke .
4.3 Programmiertips zur Integration
Ein typisches Integrationsprogramm hat drei wesentliche Teile:
1. Funktion, die f(x) f¨ur einen beliebigen input-Wert x berechnet;
2. Subroutine f¨ur die eigentliche Integration (bei festem N):
• Parameter:
– input: Integrationsgrenzen a, b
– input: (feste) Anzahl von St¨utzstellen N – output: Integralwert s
• Schleife mit Z¨ahlvariable i:
– berechne Integral si (Aufruf der subroutine) – Konvergenztest: ist |si −si−1| < ϑsi ?
wenn ja, Ergebnisausgabe und Stop; sonst weiter – (zum Testen si ausgeben)
– erh¨ohe N
Fehleranf¨alligste Stelle: Konvergenztest:
• Absolutwert der Differenz beachten!
• relativer Test besser als absoluter |si −si−1| < ϑ (dann ist ϑ abh¨angig von f(x) und Integrationsin-tervall)
• die Computergenauigkeit ist begrenzt
⇒ man kann ϑ auch zu klein w¨ahlen!
4.4 Simpsonregel und andere Verfahren
• Simpson/Keplersche Fassregel: f(x) wird im Intervall h nicht durch Gerade, sondern durch Parabel angen¨ahert; dies ergibt:
Zb
a
f(x)dx ≈ h
3(f1 + 4f2 + 2f3 + 4f4 +· · ·+ 2fN−2 + 4fN−1 +fN) (37)
https://de.wikipedia.org/wiki/Datei:Composite_Simpsons_rule2.png
• Implementation sehr ¨ahnlich zur Trapezregel. Achtung: Zahl der Punkte muss ungerade sein!
• Diese Formel ist exakt f¨ur Polynome 3. Ordnung; der Fehler ist von h¨oherer Ordnung als bei der Trapezregel:
O 1
N4
= O(h4) (38)
• Zahlreiche weitere Formeln dieser Art verf¨ugbar: Bis zur 7. Ordnung (https://de.wikipedia.org/wiki/Newton-Cotes-Formeln).
Danach aufgrund von negativen Gewichten nicht mehr sinnvoll (Rundungsfehler bei Differenzbildung ¨ahnlich großer Zahlen!)
• In manchen Situationen wichtig (s.u.): Es gibt auch sog.
”offene“ Formeln, bei denen die Funktionswerte f(a) und/oder f(b) an den Intervallgrenzen nicht berechnet werden m¨ussen, sondern nur Funktionswerte an Punkten streng innerhalb des Intervalls.
4.5 Adaptive Quadratur
• Manche Bereiche von f(x) ben¨otigen mehr Integrationspunkte als andere.
• Mehrgitter-Verfahren:
1. Starte bei einem groben Gitter (per Trapezintegration)
2. Verfeinere das Gitter durch Verfahren h¨oherer Ordnung (Trapez→Simpson)
3. Bestimme aus dem Genauigkeitsgewinn vom groberen zum feineren Gitter, wo feinere Intervalle n¨otigt sind.
4. Iteriere bis zur Konvergenz.
⇒ Die Gitterpunkte sammeln sich dort, wo
”viel passiert“.
• Probleme: Fehlerabsch¨atzung, Robustheit
• Beispiele von resultierenden Integrationsgittern:
0
4.6 Gauß-Integration
Bei allen obigen Newton-Cotes-Formeln werden nur Anzahl der ¨aquidistanten St¨utzstellen und Werte der Koeffizienten variiert, aber f¨ur Erzielung noch h¨oherer Genauigkeit kann man noch einen weiteren Parameter variieren: Positionierung von
nicht-¨aquidistanten St¨utzstellen. Man kann zeigen (z.B. Stoer/Bulirsch), daß
b
exakt ist f¨ur f(x) = Polynom vom Grad ≤ (2N −1), wenn die xi Nullstellen eines Orthogonalen Polynoms vom Grad N sind und die wi geeignet gew¨ahlte Gewichte.
• Exaktheit f¨ur 2N −1 ist viel mehr als bei Newton-Cotes (N −1 bzw. N)
• Theorie dahinter nicht ganz trivial; hier unn¨otig; Anwendung auch ohne m¨oglich.
• {xi, wi} aus Tabellen (Abramowicz/Stegun) oder mit speziellen Programmen (Numerical Recipes) f¨ur beliebige N berechen-bar (Hintergrund z.T. in Henrik Larssons Numerikskript WS15/16)
• bei bekannten {xi, wi} nicht langsamer als Newton-Cotes
• Fehlerformeln existieren (z.B. Stoer/Bulirsch), sind aber schwierig anwendbar (gehen nicht so einfach mit Nk oder hk).
• Verfeinerung des Integrationsgitters/Mehrgitterverfahren schwieriger als bei Newton-Cotes
• bei Implementation zu beachten:
– wg. Definitionsbereich der Orthogonalen Polynome gelten Gauß-Integralformeln direkt nur im Intervall [−1,+1] ⇒ die St¨utzpunkte m¨ussen auf das allgemeine Intervall [a, b] umskaliert werden, und das Integralergebnis ebenso (mit 1/2×(b−a), analog zur Formel hfi(b−a)).
– Legendre-Polynome sind (un)gerade bzgl. Ursprung, daher sind in Tabellen oft nur halb so viele Punkte und Gewichte wie scheinbar n¨otig angegeben.
Einschub: Orthogonale Polynome
• Wichtige Klasse von Funktionen, welche ein vollst¨andiges Orthogonalsystem bilden: (klassische) Orthogonale Polynome
• Unterschiedliche Arten je nach Definitionsbereich und Art des Skalarprodukts:
hm|ni ≡ Z b
a
pm(x)pn(x)W(x)dx = hnδm,n = hn×
(1, m = n,
0, m 6= n. (40)
• W(x) nennt man Wichtungsfunktion.
• Alle orthogonalen Polynome erf¨ullen die folgenden rekursiven Beziehungen:
p−1(x) = 0, (41)
p0(x) = 1, (42)
pj+1(x) = (x−aj)pj(x)−bjpj−1(x), j = 0,1,2, . . . (43)
• Die Polynome sind zudem L¨osungen von bestimmten Differentialgleichungen (DGL) und tauchen daher in der Physik sehr oft auf.
• Wichtige Eigenschaft: Ein orth. Polynom vom Grad n hat nverschiedene Nullstellen, welche jeweils zwischen den Nullstellen vom Polynom vom Grade n−1 zu finden sind.
Die wichtigsten orthogonalen Polynome:
4.7 Schwierige Integranden
Unstetige Integranden Bei bekannten Sprungstellen: Aufteilung des Integrals, sodass die Sprungstellen jeweils die Integrati-onsgrenzen sind. Unbekannte Sprungstellen auffinden durch adaptive Quadratur.
Nadelimpulse/Singularit¨aten Aufteilung des Integrals wie bei Sprungstellen, dann: Variablentransformation (s.u.), oder: tanh-sinh-Quadratur:
• Trapezregel nach Variablensubstitution x = g(t) = tanh π2 sinh(t)
, dadurch f¨allt der Integrand doppelt-exponentiell ab.
• R1
• Neben Singularit¨aten an den Grenzen auch gut geeignet f¨ur hochgenaue Integration.
Uneigentliche Integrale
• Benutzung von entsprechender Gauß-Quadratur (Laguerre, Hermite), sofern der Integrand daf¨ur geeignet ist (analytisch).
• |a| = ∞ oder |b| = ∞: Variablentransformation mit x = 1/t, dx = dt/t2, also
”offene“ Formel anwenden (ohne Berechnung von f(a) bzw. f(b)), um Singularit¨at bei 1/t2 zu vermeiden.
Achtung: Der Integrand muss schneller abfallen als x−2!
• Beide Grenzen ∞: Aufteilung in zwei Integrale.
• Oder: Integration ¨uber gen¨ugend großen Bereich falls bekannt ist, dass der Integrand schnell abklingt.
Stark oszillierende Integranden Gefahr von Rundungsfehlern (Subtraktion ¨ahnlich großer Zahlen), evtl. analytische Mittelung ¨uber Teilintervalle
4.8 Mehrdimensionale Integration
Bei einfachen Grenzen, schwach ver¨anderlichem Integranden und hohem Genauigkeitsanspruch: Wiederholte 1D-Integrationen:
Dabei sind Grenzen der 2., 3., . . . Integration durch Funktionen der vorhergehenden Variablen auszudr¨ucken, z.B. Integral ¨uber einen Kreis mit Radius 1 um den Ursprung:
1
Problem dieses Zugangs: Anzahl der St¨utzstellen skaliert bei naivem Vorgehen exponentiell mit der Anzahl der Dimensionen:
x x x x
Alternative: D¨unnbesetzte Gitter, z.B. Smolyak- oder Clenshaw-Curtis-Quadratur:
http://people.sc.fsu.edu/∼jburkardt/datasets/sparse_grid_cce/345.png
Hat der Integrand
”lokale Spitzen“ in unbekannten Stellen, ist hochdimensionale Integration hoffnungslos.
4.9 Einschub: Zufallszahlen
• “echte” Zufallszahlen: Detektion von externen Signalen (Mausbewegung, Spannungs-/Frequenzunterschiede, . . . );
dauert i.d.R. zu lange
• berechnete Pseudozufallszahlen: schell viele Zahlen erzeugbar; aber deterministisch und sich periodisch wiederholend.
• Starte von einer Zahl I0 (random seed) als und berechne daraus iterativ neue Zufallszahlen. Einfaches Beispiel:
Ij+1 = mod(a×Ij + c, m), (46)
mit Parametern a, c (mod liefert den Rest einer Division zur¨uck). Liefert ramdom integers zwischen 0 und m.
Erzeugung von random reals im Intervall [0,1] ¨uber
r = Ij
m. (47)
• In der Praxis deutlich kompliziertere Formeln. De-facto-Standard: Mersenne-Twister,1 wiederholt sich erst nach 219937 −1≈ 4,3·106001 Berechnungen.
• Einfacher Test auf Regelm¨aßigkeiten: Graphische Auftragung einer Zahlensequenz als Punkte in einem 2D-plot:
Schlechter versus guter Zufallszahlengenerator:2
1https://en.wikipedia.org/wiki/Mersenne_twister
2http://www.reedbeta.com/blog/2013/01/12/quick-and-easy-gpu-random-numbers-in-d3d11/
Wichtige Praxistricks:
• f¨ur Fehlersuche: immer wieder denselben random seed verwenden.
• f¨ur Produktionsrechnungen quasi-zuf¨alliger
random seed: aktuelle Zeit, inkl. (Milli)sekunden
Zerst¨orung von Korrelationen und damit Verbesserung eines existierenden Zufallszahlengenerators:
• Initialisierung: generiere N Zufallszahlen und speichere sie in array r(i)
• bei allen folgenden Aufrufen:
– skaliere vorherige Zufallszahl auf integer-Wert j im Intervall [1, N]
– r(j) wird als Zufallszahl-output ausgegeben – generiere eine neue Zufallszahl und speichere sie
in r(j)
Einfaches Programmieren mit Zufallszahlen
Verschiedene Verteilungen
• Gleichverteilung von random_number(r) in [0,1[ direkt verwenden
• gleichverteiltes ˜r im Intervall [x1, x2[ durch ˜r = r(x2 −x1) +x1
• andere Verteilungen (exponentiell, Gaußf¨ormig, Poisson,. . .) leicht aus Gleichverteilung konstruierbar (Numerical Recipes)
M¨unzwurf-Simulation (direkte Erzeugung zuf¨alliger bits w¨are effizienter, siehe Numerical Recipes)
Zufallsereignisse bekannter Wahrscheinlichkeit
z.B.: Ereignis a: 60%, Ereignis b: 30%, Ereignis c: 10%
0 1 Im allgemeinen Fall mitnEreignissen mit Wahrscheinlichkeitenwi, i =
1, . . . , n, 0 ≤ wi ≤ 1 sowie P
iwi = 1, wird diese Implementation wegen vieler if-Abfragen unn¨otig langsam.
Alternative:
Dieses Programmfragment w¨ahlt integer-Wert i mit Wahrscheinlichkeit w(i), unter obigen Voraussetzungen.
4.10 Monte-Carlo-Integration
Anschauliche Vor¨ubung: MC-Berechnung von π:
r r
0
Fl¨ache 4
Fl¨ache = πr2/4 r2 = 1
4π (50)
⇒ Algorithmus zur Berechnung von π:
1. Generiere zwei Zufallszahlen xi, yi (gleichverteilt in [0,1); implizite Annahme: r = 1) 2. Teste, ob Punkt (xi, yi) im Viertelkreis liegt (trifft zu, wenn x2i +yi2 ≤ 1)
3. π = 4× Anzahl Punkte im Viertelkreis Anzahl aller Punkte
4. Wenn noch nicht konvergiert, gehe zu (1) Konvergenz: nur ca. wie √
N, d.h. viermal so viele Punkte n¨otig, um den Fehler zu halbieren.
MC-Integration in 1D
zwei m¨ogliche Sichtweisen / Verfahren:
mittlere Fl¨ache: wie MC-Berechnung von π:
<f>
y(b−a)x,y-Punkte unterf(x) alle x,y-Punkte (53)
Offensichtlich ist Verfahren 1 einfacher: Nur eine statt zwei Zufallszahlen pro Funktionsberechnung.
Konvergenz ebenso schlecht wie bei MC-π-Berechnung ⇒ nicht konkurrenzf¨ahig zu Newton-Cotes, Gauß,. . . , in 1D, aber √ N -Konvergenz ist dimensionsunabh¨angig ⇒ MC-Integration gewinnt in h¨oheren Dimensionen immer:
MC-Integration in nD
In mehr als einer Dimension ist MC-Integration eine gute Wahl, wenn
• die Integrationsgrenzen kompliziert sind und/oder
• der Integrand schwach ver¨anderlich ist (keine
”lokalen Spitzen“) und/oder
• niedrige Genauigkeit ausreicht.
Weitere Praxisvorteile:
• trivial-parallel;
• Genauigkeit durch weitere Punkte nachtr¨aglich sehr einfach und beliebig steigerbar (wenn auch nur mit √ N)
4.11 Schlußbemerkungen zur Praxis
• meist wird nichts von alldem selbst programmiert, sondern man benutzt Bibliotheksroutinen; einmal selbst programmieren hilft aber extrem, mit Bibliotheksroutinen umzugehen.
• Quantenchemie: Billionen von Integralen schon bei kleinen Molek¨ulen
– wenige Gruppen weltweit haben extrem hochgez¨uchtete Spezialverfahren entwickelt; Mix von analytischer und numeri-scher Integration
– DFT: spezielle numerische Integrationen mit Kombinationen von atomzentrierten 3D-St¨utzstellengittern (¨ahnlich zur 1D-Gaußintegration)
• MC-artige Verfahren weit verbreitet in Statistischer Thermodynamik: MC-Berechnung hochdimensionaler Integrale in mo-lekularen Freiheitsgraden; oft verbunden mit Zusatztechniken wie “importance sampling” (Plazierung von Punkten nicht gleichverteilt-zuf¨allig, sondern bevorzugt dort, wo der Integrand groß ist) u.v.a.m.
• Elementarteilchenphysik: diskretisierte Quantenchromodynamik (Gitter-QCD) mit MC-Verfahren am HLRN:
– verwenden sehr viele Knoten/CPUcores mit exzellenter Skalierung – gr¨oßter Rechenzeitanteil irgendeiner Subdisziplin
– nehmen gerne noch mehr Rechenzeit (Konvergenz wie √
N ;-) )
5 Differentialgleichungen
5.1 Vergleich zur Integration
Zwischen der Integration und der L¨osung von Differentialgleichungen (DGL) besteht ein enger Zusammenhang: Das Aufsuchen der L¨osung y(x1) des Integrals
x1
Z
x0
f(x)dx = y(x1) (54)
ist ¨aquivalent zum Aufsuchen der L¨osung von
dy
dx = y0 = f(x) (55)
mit Anfangsbedingung y(x0) = 0 und an einem Punkt x = x1. Differentialgleichungen vom Typ Gl. 55 k¨onnen wir also mit den bekannten Integrationsverfahren bereits l¨osen.
Die allgemeine Formel f¨ur Differentialgleichungen dieses Typs ist jedoch:
y0 = f(x, y) (56)
Problem: Um y analog zu Gl. 54 berechnen zu k¨onnen, m¨ußte y bereits bekannt sein. ⇒ anderes, extrapolatives Verfahren n¨otig!
5.2 Anfangswert- und Randwertprobleme
Allgemeine DGL-L¨osungen (mit unbestimmten Konstanten) sind numerisch nicht erzeugbar. Zur Ermittlung einer speziellen L¨osung einer DGL 2.Ordnung y00 = f(x, y, y0) sind zwei Angaben n¨otig (zur Bestimmung der beiden Integrationskonstanten):
Anfangswertproblem: Zahlenwerte f¨ur y(x0) und y0(x0) Randwertproblem: Zahlenwerte f¨ur y(x0) und y(x1)
Der gr¨oßte Teil dieses Kapitels behandelt die viel einfacheren Anfangswertprobleme; ein Unterkapitel geht kurz auf Randwert-probleme ein.
5.3 Explizites Euler-Cauchysches Polygonzugverfahren F¨ur eine Differentialgleichung
y0 = dy
dx = f(x, y) (57)
mit Anfangsbedingung (x0, y0) sind weitere Punkte (x, y) der Funktion y(x) gesucht, die die Differentialgleichung Gl. 57 erf¨ullt.
F¨ur jeden Punkt (x, y) gibt das bekannte f(x, y) exakt die Steigung der Funktion y(x) an. Es liegt nahe, diese Information f¨ur eine linear-extrapolierende Konstruktion von y(x) zu verwenden:
x1
x0 x2 x
y
y0
Wir machen also die Approximation
dy
dx ≈ ∆y
∆x (58)
woraus sich ergibt
∆y = yi+1 −yi = f(xi, yi) ∆x = hf(xi, yi) (59) oder folgende iterativ-punktweise Erzeugung (Propagation) von y(x), ausgehend vom Anfangswert (x0, y0):
xi+1 = xi+ h , yi+1 = yi +hf(xi, yi) , i = 0,1,2, . . . (60) Nach dieser Herleitung k¨onnte bei jedem Schritt ein anderer Wert f¨ur h = ∆x gew¨ahlt werden. Das ist ohne weitere Kriterien jedoch schwierig (s.u.), daher w¨ahlt man hier f¨ur alle Schritte denselben h-Wert.
Gl. 59 entspricht einer Taylorreihe bis inkl. dem Term 1.Ordnung, daher ist der Fehler dieses Euler-Verfahrens O(h2). Bei der Euler-Integration wird interpoliert, hier wird jedoch extrapoliert. Daher ist das Euler-Verfahren bei DGLs i.d.R. nicht effizient genug.
5.4 Runge-Kutta-Verfahren
Man kann eine Verringerung des Fehlers erreichen, wenn man f¨ur die Steigung im Intervall nicht die Steigung am Anfang, sondern die Steigung in der Intervallmitte ansetzt, nach folgender Strategie:
• konstruiere einen Hilfspunkt in der Intervallmitte: x1/2 = x0 + 12h
• berechne den dortigen y-Wert: y1/2 = y0 + 12hf(x0, y0) (bis hierher ist das ein Euler-artiger Schritt mit h/2)
• die Steigung in (x1/2, y1/2) ergibt sich durch Einsetzen in f(): f(x1/2, y1/2)
• mit dieser Steigung geht man dann von x0 nach x1: y1 = y0 +hf(x1/2, y1/2)
(das ist wieder ein Euler-artiger Schritt, diesmal mit Schrittweite h) x0 x1 x2 x y
y0
x1/2
h
In der in diesem Gebiet ¨ublichen Notation f¨uhrt man zwei Hilfsgr¨oßen k1 und k2 ein und schreibt dies als:
k1 = hf(xi, yi) (61)
k2 = hf(xi + 1
2h, yi + 1
2k1) (62)
yi+1 = yi +k2 (63)
Dieses Schema heißt Runge-Kutta-Verfahren 2.Ordnung (modifiziertes Euler-Verfahren, midpoint-Verfahren). Die Erh¨ohung der Fehlerordnung auf O(h3) erkauft man sich dabei durch 2 Funktionsberechnungen pro Schritt (und zwei Euler-artige Schritte pro eigentlicher Schrittweite h).
Verallgemeinerung zu h¨oherer (m-ter) Ordnung:
Besonders oft verwendet: Runge-Kutta 4. Ordnung (verwandt mit Simpson-Regel bei Integration):
k1 = f(xi, yi) (66)
Auflistung der Koeffizienten oft per Butcher-Schema: α~ βββ
~γT
Euler RK 2.Ordng. RK 4.Ordng.
0 Preis f¨ur die Ordnung O(h5): 4 (aufwendige!) Funktionsberechnungen pro Schritt.
Aber oftmals gilt in der Praxis:
h bei Runge-Kutta 4. Ordnung ≥ 2h bei midpoint ≥ 4h bei Euler
4 Funktionsberechnungen werden durch l¨angere Schrittweite ¨uberkompensiert.
Wie bei der Integration gilt allerdings auch hier:
• h¨ohere Ordnung 6= h¨ohere Genauigkeit (bei vergleichbarem Aufwand)/Effizienz; Gleichheit nur bei hinreichend gutartigen Differentialgleichungen m¨oglich.
• Konvergenztests (verschiedene Schrittweiten h) extrem wichtig.
• Zu kleine Schrittweite h wegen begrenzter Computergenauigkeit auch nicht gut (siehe ¨Ubungen zur numerischen Differen-tiation).
Implementation
Algorithmen zur L¨osung von Differentialgleichungen mit fester Schrittweite sollten in drei Modulen implementiert werden, sehr
¨ahnlich zur Integration:
1. Funktion C zur Berechnung von f(x, y):
• input: x, y
• berechne output y0 = f(x, y)
2. subroutine B, die einen Schritt der L¨ange h macht:
• input: h, xi, yi
• Formel zur Berechnung von yi+1 aus yi, xi, h (und ggf. k1, k2, . . .) und einem (oder mehreren) Aufrufen der Funktion C
• output: yi+1, xi+1 = xi +h 3. Hauptprogramm A:
• user-interface: Eingabe von – Anfangsbedingungen: x0, y0
– xf, sowie entweder Schrittweite h oder Schrittanzahl n
(die jeweils andere Gr¨oße kann/muß intern berechnet werden) – Genauigkeitstoleranz
• X: Schleife ¨uber n Aufrufe der subroutine B
• Ausgabe von yf = y(xf) und ggf. analytischer Vergleichswert
• wenn Genauigkeit nicht ausreichend, h verkleinern (bzw. n erh¨ohen) und zur¨uck zu X
5.5 Adaptive Schrittweitenkontrolle
5.5.1 Grundideen
Die obige Technik des iterativen
”Nochmal-Propagierens“ mit verkleinertem h bis zur Konvergenz kann (fast) v¨ollig vermieden werden, wenn der Programmteil B an Teil A zur¨uckmeldet, ob h beibehalten werden kann oder nicht; daraufhin kann A ent-scheiden, ob und wie h ggf. vergr¨oßert oder verkleinert werden sollte. Dazu muß B eine Fehlerabsch¨atzung machen. Dazu gibt es verschiedene M¨oglichkeiten, z.B.:
• vergleiche einen Schritt mit gleichem h von zwei Verfahren der Ordnung N und N + 1 ⇒ z.B. Runge-Kutta-Fehlberg; oder:
• ”Schrittverdopplung“: vergleiche einen regul¨aren Schritt xi → xi+1 von Schrittweite 2h mit einem Doppelschritt xi →→ xi+1
mit Schrittweite h.
5.5.2 Schrittverdopplung
Schrittverdopplung f¨ur Runge-Kutta 4.Ordnung:
(2h-Schritt) : y(x+ 2h) = yRK,1 + (2h)5φ+O(h6) (71) 2×(h-Schritt) : y(x+ 2h) = yRK,2 + 2(h5)φ+O(h6) (72)
⇒ ∆ = |yRK,2 −yRK,1| ∝ h5 (73)
∆ soll kleiner als sein: ∆ ≤ = ATOL +|yRK,1| ×RTOL. Mit einem gew¨unschten, relativen Fehler k¨onnten wir dann sofort die f¨ur diese Zielgenauigkeit n¨otige Schrittweite hneu ausrechnen:
hneu = h Allerdings: hneu wurde f¨ur xi → xi+1 berechnet, soll aber auf xi+1 → xi+2 angewandt werden. Daher besser: Konservativere Absch¨atzung hneu = Ch
∆
1/5, C - 1 und zus¨atzliches Abfangen von zu großen/kleinen h-Werten (z.B. maximal Schrittweiten-verdopplung/halbierung)
5.5.3 Runge-Kutta-Fehlberg
Es k¨onnen Koeffizienten gefunden werden, die sowohl f¨ur ein RK-Verfahren nter Ordnung als auch f¨ur ein Verfahren (n−1)ter Ordnung gelten. Beispiel: Runge-Kutta-Fehlberg 5(4):
0
1/4 1/4
3/8 3/32 9/32
12/13 1932/2197 −7200/2197 7296/2197
1 439/216 −8 3680/513 −845/4104
1/2 −8/27 2 −3544/2565 1859/4104 −11/40
16/135 0 6656/12825 28561/56430 −9/50 2/55 (RKF5)
25/216 0 1408/2565 2197/4104 −1/5 0 (RKF4)
Der Fehler ist damit bis zur 5ten Ordnung direkt absch¨atzbar, ohne Schrittverdopplung. Die neue Schrittweite ergibt sich wieder nach Gl. 74.
Hier ist das RKF4-Resultat f¨ur die Fortsetzung der Rechnung und das RKF5-Resultat zur Schrittweitenkontrolle gedacht. Bei der Dormand-Prince-Methode (RKDP) ist es umgekehrt (und die Koeffizientenwerte sind etwas anders).
5.5.4 Einsch¨atzung
Generell ist die Fehlerabsch¨atzung mit signifikantem Mehraufwand verbunden:
• Derselbe Schritt muß mind. zweimal gerechnet werden.
• h-Verkleinerung/Vergr¨oßerung muß sinnvoll kontrolliert werden (Vermeidung von zu großen ¨Anderungen und underflow) Dennoch kann die adaptive Schrittweitenkontrolle sehr lohnend sein:
• Einfache Kontrolle der Zielgenauigkeit (⇒ Black-Box-Routinen)
• Jede Region bekommt w¨ahrend der Propagation automatisch ein an sie angepasstes h;
• Bei guter -Wahl werden Iterationen des Gesamtverfahrens ggf. v¨ollig ¨uberfl¨ussig.
5.6 Pr¨adiktor-Korrektor/Mehrschrittverfahren
Grundidee: Bei bisherigen Verfahren ist jeder ermittelte Punkt ein neuer Anfangspunkt; alte y-Werte werden
”vergessen“. Bei Pr¨adiktor-Korrektor-Verfahren werden stattdessen die folgenden beiden Schritte alternierend iteriert:
• Pr¨adiktor-Schritt: N −1 alte y-Werte werden benutzt, um einen neuen y-Wert yN zu extrapolieren;
• Korrektor-Schritt: F¨ur die jetzt N bekannten y-Werte wird dydx = f(x, y) in
b
R
a
f(x, y)dx = y(b) zur¨uckverwandelt und mit einer geeigneten N-Punkte-Integrationsformel ausgerechnet ⇒ verbesserter Wert von yN.
H¨aufig verwendet:
”Adams-Bashforth-Moulton“ (ABM); z.B. 3. Ordnung:
yn+1Pr¨adiktor = yn+ h
12(23y0n−16yn−10 + 5y0n−2) +O(h4), (Adams-Bashforth-Schritt) (75) ynKorrektor+1 = yn+ h
12(5yn+10
+ 8yn0 −yn−10 ) +O(h4), (Adams-Moulton-Schritt). (76)
• Anfangswerte werden oft per RK2 oder ¨ahnliche Verfahren generiert.
• Sehr gutes Verfahren, wenn die L¨osung der ODE nicht stark oszillativ und die Berechnung von f sehr aufwendig ist.
• Differenz Pr¨adiktor-Korrektor → Genauigkeitsabsch¨atzung → Schrittweiten- oder Ordnungskontrolle; aber:
• Programmierung sehr aufwendig (Buchhaltung alter Werte n¨otig);
• Programmieraufwand steigt mit Schrittweitenkontrolle noch einmal erheblich.
• Allerdings: Altes, daher ausgereiftes Verfahren ⇒ moderne Implementationen enthalten Jahrzehnte an Entwicklungsarbeit
5.7 Weitere Verfahren und Aspekte
Bulirsch-Stoer-Verfahren: große Schrittweite, mit Verfahren steigender Ordnung, Richardson-extrapoliert auf ∞-Ordnung.
Steife DGLs: genaue Definition schwierig; pragmatisch:
• gewisse Komponenten der L¨osung klingen sehr viel schneller ab als andere;
• nicht Genauigkeits-, sondern Stabilit¨atsanforderungen begrenzen die Schrittweite;
• explizite L¨osungsverfahren haben Probleme; funktionieren wenn ¨uberhaupt dann nur mit extrem kleiner Schrittweite.
Implizite Verfahren: haben deutlich weniger Probleme mit steifen DGLs.
Euler explizit: yi+1 = yi +hf(xi, yi) (77)
Euler implizit: yi+1 = yi +hf(xi+1, yi+1) (78)
yi+1 kann aus Gl. 78 nicht direkt berechnet werden, da es links und rechts vorkommt! ⇒ Iterative Nullstellensuche n¨otig, um yi+1 n¨aherungsweise zu bestimmen. Kann bei steifen DGLs trotzdem lohnend sein.
(Bsp. f¨ur steife DGL und explizite/implizite L¨osung: Henriks Skript WS15/16)
5.8 Randwertprobleme
Einfaches Schießverfahren: Formuliere Randwertproblem
y00 = f(x, y, y0), mit y(x0) = a, y(x1) =b (79) um in ein Anfangswertproblem
y00 = f(x, y, y0), mit y(x0) = a, y0(x0) =s (80) Rate Versuchswert s0 und propagiere von x0 zu x1.
Je nach y(x1) wird s0 in geeigneter Weise zu s1 ver¨andert. Iteration bis zur Konvergenz.
x0 x
y
y0
x1 s0
s1
s2
sn y1
Nachteile:
• Nicht jede gegebene Randbedingung (x1, y1) muss notwendigerweise erf¨ullbar sein ⇒ keine theoretische Garantie f¨ur eine L¨osung;
• Die Wahl von s0 kann so falsch sein, dass dazu gar kein y(x1) existiert;
• y1 kann so sensitiv vons abh¨angen, daß das Auffinden des korrektens praktisch unm¨oglich wird⇒keine praktische Garantie f¨ur eine L¨osung.
5.9 DGL-Systeme, h¨ohere Ordnung
Gekoppelte DGL-Systeme 1.Ordnung kommen h¨aufig vor, z.B. chemische Kinetik (Phantasiebeispiel):
dA
dt = −k1A−k2AB (81)
dB
dt = k3C2 −k4BD (82)
...
Sie sind mit allen obigen Verfahren behandelbar. Einzige ¨Anderung in der Implementation:y undf(x, y) sind keine Einzelvariablen mehr, sondern Vektoren.
DGLs h¨oherer Ordnung sind immer auf DGL-Systeme 1.Ordnung r¨uckf¨uhrbar, durch Einf¨uhrung geeigneter Hilfsgr¨oßen.
Beispiel:
Nach Newton gilt f¨ur 1 Teilchen in 1 Dimension:
ma = mdv Das ist eine DGL 2.Ordnung.
Mit dem Impuls p = mv ergeben sich die Hamiltonschen Bewegungsgleichungen:
dx Das ist ein System aus zwei gekoppelten DGLs 1.Ordnung.
5.10 Partielle Differentialgleichungen
Weitere Schwierigkeitsstufe, weil Variabilit¨at noch einmal h¨oher:
Integration → gew. DGL-System → partielles DGL-System
dy
dx = f(x) dydxi = f(x, y1, y2, . . .) ∂x∂yi
j = f(x1, x2, . . . , y1, y2, . . .) (86)
⇒ Weites, neues Feld mit viel weniger Systematik und Theorie, und mit vielen speziellen Methoden.
Zwei Standardstrategien:
• Ersatz von Differentialquotienten durch finite Differenzen, oder
• Entwicklung unbekannter L¨osungsfunktionen in Reihen aus bekannten Basisfunktionen.
⇒ Umwandlung der partiellen DGL in Probleme der linearen Algebra, L¨osung mit den dort etablierten Methoden (lineare Gleichungssysteme, Matrix-Eigenwertproblem).
5.10.1 Beispiel Finite Differenzen I
2. Fouriersches Gesetz der W¨armeleitung/2. Ficksches Gesetz (in 1D):
∂T
∂t = c ∂2T
∂x2 (87)
Approximation der zweiten Ableitung durch finite Differenzen
yi00 = yi+1 −2yi +yi−1
h2 +O(h3) (88)
ergibt auf der rechten Seite
c ∂2T
Definiert man die Temperatur am Punktxi alsTi, erh¨alt man auf diesem 1D-Ortsraumgitter also (inkl. zweier Randbedingungen):
Das ist ein System gew¨ohnlicher DGLs, also mit den ¨ublichen Methoden direkt l¨osbar.
5.10.2 Beispiel Finite Differenzen II
Die Schr¨odingergleichung f¨ur die Bewegung eines Teilchens der Masse m in einer Raumdimension x lautet:
HˆΨ = EΨ, (95)
Die Diskretisierung der 2. Ableitung von Ψ nach x d2Ψ
dx2 ≈ Ψ(xi+1) + Ψ(xi−1)−2Ψ(xi)
∆x2 (99)
f¨uhrt nach Einsetzen in Gleichung (98) zum Gleichungssystem
[Ψ(xi+1) + Ψ(xi−1) +aiΨ(xi)] = λΨ(xi) (100)
Somit lassen sich die Gleichungen (100) in Matrixform schreiben als:
was ein einfaches Eigenwertproblem einer tridiagonalen Matrix ist.
Nachteile dieser sehr einfachen Strategie sind die relativ ungenaue Diskretisierungsformel (⇒ große Matrizen) und Probleme bei der Verallgemeinerung auf mehr als eine Dimension.
5.10.3 Spezielle Verfahren in der theoretischen Chemie
An das physikalische Problem angepasste DGL-L¨oser sind deutlich besser als Standardverfahren. Beispiele:
• Klassische Molekulardynamiksimulationen: simple Verlet, leapfrog, velocity Verlet, symplektische Integratoren;
siehe chem0503, chem1004D
• Zeitunabh¨angige Schr¨odingergleichung (partielle DGL!): Quantenchemie-Verfahren (Coupled Cluster etc.), Galerkin/Finite Basisdarstellung, Numerov-Methode, Kollokation;
siehe chem1004D
• Zeitabh¨angige Schr¨odingergleichung: Short Iterative Lanczos, symplektische Integratoren, Chebycheff, Split-Operator, Taylor-Reihe, . . . ;
siehe chem1004D
Beispiel: Basisfunktionsentwicklung (Galerkin/Finite Basisdarstellung) Wesentlich allgemeiner und sehr leicht auf beliebige
Beispiel: Basisfunktionsentwicklung (Galerkin/Finite Basisdarstellung) Wesentlich allgemeiner und sehr leicht auf beliebige