Explizites Euler-Verfahren

Das eulersche Polygonzugverfahren oder explizite Euler-Verfahren (auch Euler-Cauchy-Verfahren oder Euler-vorwärts-Verfahren) ist das einfachste Verfahren zur numerischen Lösung eines Anfangswertproblems.

Es wurde von Leonhard Euler 1768 in seinem Buch Institutiones Calculi Integralis vorgestellt. Cauchy benutzte es, um einige Eindeutigkeitsresultate für gewöhnliche Differentialgleichungen zu beweisen.

Das Verfahren wird manchmal in der Physik als Methode der kleinen Schritte bezeichnet.

Das Verfahren

Zwei Schritte des expliziten Euler-Verfahrens

Zur numerischen Lösung des Anfangswertproblems:

${\displaystyle {\dot {y}}(t)=f(t,y),\quad \quad y(t_{0})=y_{0}}$

für eine gewöhnliche Differentialgleichung wähle man eine Diskretisierungs-Schrittweite ${\displaystyle h>0}$, betrachte die diskreten Zeitpunkte

${\displaystyle t_{k}=t_{0}+kh,\quad \quad k=0,1,2,\dots }$

und berechne die Werte

${\displaystyle y_{k+1}=y_{k}+hf(t_{k},y_{k}),\quad k=0,1,2,\dots }$

Die berechneten Werte ${\displaystyle y_{k}}$ stellen Approximationen an die tatsächlichen Werte ${\displaystyle y(t_{k})}$ der exakten Lösung des Anfangswertproblems dar. Je kleiner die Schrittweite ${\displaystyle h}$ gewählt ist, desto mehr Rechenarbeit ist nötig, aber desto genauer werden auch die approximierten Werte.

Eine Modifikation des Verfahrens besteht hier darin, dass man die Schrittweite variabel wählt. Eine sinnvolle Veränderung der Schrittweite setzt einen Algorithmus zur Schrittweitensteuerung voraus, der den Fehler im aktuellen Schritt abschätzt und dann die Schrittweite für den nächsten Schritt dementsprechend wählt.

Wird ein Verfahren über ${\displaystyle y_{k+1}=y_{k}+hf(t_{k+1},y_{k+1})}$ definiert, erhält man das implizite Euler-Verfahren. Dieses ist A-stabil und daher für steife Anfangswertprobleme besser geeignet.

Herleitung

Für die Herleitung von Einschrittverfahren wird das Anfangswertproblem meist in die dazu äquivalente Integralgleichung umgeformt^[1]

${\displaystyle {\begin{aligned}{\dot {y}}&=f(t,y)\,,\qquad y(t_{0})=y_{0}\\\Longleftrightarrow \quad y(t)&=y_{0}+\int _{t_{0}}^{t}f(s,y(s))\,\mathrm {d} s\,.\end{aligned}}}$

Nun besteht die Idee, beim expliziten Euler-Verfahren eine simple Quadraturformel für das Integral zu benutzen: die linksseitige Boxregel. Man wählt in jedem ${\displaystyle k}$-ten Schritt den Integranden als konstanten Wert an der linken Intervallgrenze^[2]

${\displaystyle y(t_{k+1})=y(t_{k})+\int _{t_{k}}^{t_{k+1}}f(s,y(s))\,\mathrm {d} s\approx y_{k}+\int _{t_{k}}^{t_{k+1}}f(t_{k},y_{k})\,\mathrm {d} s=y_{k}+hf(t_{k},y_{k})=:y_{k+1}\,.}$

Eigenschaften

Stabilitätsgebiet des expliziten Euler-Verfahrens

Das explizite Euler-Verfahren hat Konsistenz- und Konvergenzordnung 1. Die Stabilitätsfunktion ist ${\displaystyle R(z)=1+z}$ und sein Stabilitätsgebiet daher das Innere des Kreises um −1 mit Radius 1 in der komplexen Zahlenebene.

Verbessertes explizites Euler-Verfahren

Anstatt die Boxregel für die numerische Integration zu verwenden, kann man auch die Mittelpunktsregel anwenden.

${\displaystyle \int _{t_{k}}^{t_{k+1}}f(s,y(s))\,\mathrm {d} s\approx hf\left(t_{k}+{\frac {h}{2}},y\left(t_{k}+{\frac {h}{2}}\right)\right)\,.}$

Nun wendet man wieder das explizite Euler-Verfahren zur Approximation von ${\displaystyle y\left(t_{k}+{\frac {h}{2}}\right)}$ an

${\displaystyle y\left(t_{k}+{\frac {h}{2}}\right)\approx y_{k}+{\frac {h}{2}}f(t_{k},y_{k})=:y_{k+{\frac {1}{2}}}\,.}$

Zusammen führt dies auf das verbesserte explizite Euler-Verfahren (oder Euler-Verfahren mit kleinerer Schrittweite)^[3]

${\displaystyle y_{k+1}=y_{k}+hf\left(t_{k}+{\frac {h}{2}},y_{k+{\frac {1}{2}}}\right)\quad {\text{mit}}\quad y_{k+{\frac {1}{2}}}=y_{k}+{\frac {h}{2}}f(t_{k},y_{k})\,.}$

Verallgemeinerungen

Es lässt sich im Wesentlichen durch zwei verschiedene Ideen auf effizientere Verfahren verallgemeinern.

• Die erste Idee ist, bei der Berechnung des nächsten Schrittes mehr als nur einen der zuvor berechneten Werte mit einzubeziehen. Auf diese Weise erhält man Verfahren höherer Ordnung in der Klasse der linearen Mehrschrittverfahren.

• Die zweite Idee ist, bei der Berechnung des nächsten Schrittes die Funktion ${\displaystyle f(t,y)}$ auf dem Intervall ${\displaystyle [t_{k},t_{k+1}]}$ an mehreren Stellen auszuwerten. Auf diese Weise erhält man die Klasse der Runge-Kutta-Verfahren.

Die Klasse der allgemeinen linearen Verfahren bezieht beide Ideen der Verallgemeinerung mit ein und enthält die Klasse der linearen Mehrschrittverfahren sowie die Klasse der Runge-Kutta-Verfahren als Spezialfall.

• Darüber hinaus gibt es auch eine Erweiterung des Euler-Verfahrens für stochastische Differentialgleichungen: das Euler-Maruyama-Verfahren.

Anwendungsbeispiel: Erdnaher freier Fall

In der Physik wendet man die Methode der kleinen Schritte beispielsweise bei der Bewegung eines Fallschirmspringers (im Fall) mit Luftwiderstand an.

Physikalischer Hintergrund

• Beim freien Fall in Erdnähe würde die Geschwindigkeit ${\displaystyle v}$ eines fallenden Körpers – bei Vernachlässigung des Luftwiderstandes – um 9,81 m/s pro Sekunde steigen. Dann wäre der freie Fall eine gleichmäßig beschleunigte Bewegung. Ein Fallschirmspringer, der sich aus einem stationären Ballon fallen lässt, wird zunächst immer schneller, seine Geschwindigkeit nimmt stetig zu. Seine Beschleunigung entspricht dabei der Fallbeschleunigung g=9,81 m/s² in der Nähe der Erdoberfläche und ist größer als die eines Autos: Nach einer Sekunde hat er theoretisch eine Geschwindigkeit von v = 9,81 m/s (ca. 35 km/h), nach zwei Sekunden 19,62 m/s (ca. 71 km/h), nach drei Sekunden 29,43 m/s (ca. 106 km/h). In einem echten freien Fall, d. h. im Vakuum, würde die Geschwindigkeit linear weiter entsprechend ansteigen.
• Tatsächlich wirkt auf den Fallschirmspringer jedoch auch der Luftwiderstand, welcher quadratisch mit der Geschwindigkeit zunimmt. Die resultierende Beschleunigung entspricht daher nur am Anfang der Erdbeschleunigung, danach nimmt sie ab, bis sie nach wenigen Sekunden null wird – der Fallschirmspringer fällt ab diesem Zeitpunkt (je nach Wahl der Anfangsbedingungen) mit seiner Fallgrenzgeschwindigkeit von z. B. ca. 55 m/s (ca. 198 km/h).

Anwendung einer Tabellenkalkulation

Mit Hilfe einer Tabellenkalkulation kann man derartige Probleme aber in viele einfache und vor allem lösbare Teilaufgaben zerlegen, deren Ergebnisse man durch das Computerprogramm zur Gesamtlösung zusammensetzen lässt. Die Vorteile liegen auf der Hand:

• Man benötigt keine Kenntnisse in höherer Mathematik
• Die Integration wird durch Summieren ersetzt. Das Ergebnis ist zwar nicht exakt, genügt aber den meisten praktischen Anforderungen.
• Anhand von Zwischenergebnissen erkennt man sofort kleine Irrtümer, die sich korrigieren lassen.
• Die vielen überprüfbaren Zwischenergebnisse steigern das Vertrauen in das Resultat.
• Durch Hinzufügen weiterer relevanter Formeln kann die Lösung schrittweise der Realität angepasst werden.

Die Vorgehensweise ist immer gleich: Mit elementaren Formeln werden relevante Größen wie Kraft, Beschleunigung oder Temperatur für einen gewissen Zeitpunkt berechnet – das sind die Anfangswerte für den nächsten Zeitpunkt. Die Ergebnisse sind nur dann korrekt, wenn sich von einem Zeitpunkt zum nächsten nur wenig ändert. Wie groß diese Änderungen und vor allem jeder Zeitschritt sein dürfen, kann man den Ergebnissen leicht entnehmen. Komplexe Formeln, wie sie beispielsweise bei der Wettervorhersage vorkommen, lassen sich gar nicht anders auswerten.

Einzelformeln des freien Falls mit Luftwiderstand

In der folgenden Berechnung wird angenommen, dass ein kugelförmiger Eisen-Meteor der Masse m = 4 g und der Querschnittsfläche A = 1 cm² mit der Geschwindigkeit v = 15 km/s in die Atmosphäre eindringt und abgebremst wird. Gesucht sind Geschwindigkeit und Bremsverzögerung als Funktion der Höhe. Diese Werte werden in bekannte Formeln eingesetzt und für jeden Zeitschritt neu berechnet. Die Einzelergebnisse werden in der Tabelle zu den gesuchten Größen kombiniert und zum Schluss graphisch ausgegeben. Man startet das Verfahren in ausreichend großer Höhe h, wo der Luftwiderstand noch vernachlässigbar ist.

• Die Gravitationsbeschleunigung der Erde wird mit zunehmendem Abstand h über der Erdoberfläche kleiner. Dafür gilt

${\displaystyle a_{\mathrm {gravi} }=9{,}81\,{\frac {\mathrm {m} }{\mathrm {s} ^{2}}}\cdot \left({\frac {6370000\,\mathrm {m} }{6370000\,\mathrm {m} +h}}\right)^{2}}$

• In dieser Höhe h beträgt die Dichte der Luft nur noch (Barometrische Höhenformel)

${\displaystyle \rho (h)=\rho ({\text{Boden}})\cdot e^{-{\frac {h}{8400\,\mathrm {m} }}}}$

• Der Strömungswiderstand in Luft F_Luft bei der Geschwindigkeit v hängt auch von dieser Dichte ab

${\displaystyle F_{\text{Luft}}=0{,}5\cdot \rho (h)\cdot C_{w}\cdot A\cdot v^{2}}$

• Bei Flugrichtung zum Erdmittelpunkt ist die effektive Beschleunigung auf den Meteor der Masse m die Differenz von Gravitationsbeschleunigung und Bremsbeschleunigung

${\displaystyle a_{\mathrm {gesamt} }=a_{\mathrm {gravi} }-{\frac {F_{\mathrm {Luft} }}{m}}}$

• Mit diesem Zwischenergebnis lässt sich einen Zeitschritt dt später die dann gültige Geschwindigkeit errechnen

${\displaystyle v_{\mathrm {neu} }=v_{\mathrm {alt} }+a_{\mathrm {gesamt} }\cdot dt}$

• und daraus der Ort, an dem sich der Meteor dann befindet. Damit startet ein neuer Zyklus.

${\displaystyle h_{\mathrm {neu} }=h_{\mathrm {alt} }-v_{\mathrm {neu} }\cdot dt}$

Die Berechnung erfolgt schrittweise mit elementaren Mitteln und entspricht einer einfachen Integration, die bei ausreichend kleinem dt brauchbare Ergebnisse liefert. Speziell für die letzten beiden Schritte existieren bessere, aber auch aufwendigere Verfahren, die in Numerische Integration beschrieben sind. Oft ist deren Anwendung übertrieben, wenn nur ein schneller Überblick gewünscht wird oder – wie in diesem Beispiel – die Formel für den Strömungswiderstand für Überschallgeschwindigkeit nicht exakt gilt.

Numerische Lösung

Berechnungstabelle für freien Fall mit Luftwiderstand

Abbremsung eines Meteors in der Atmosphäre

Zunächst werden die Parameter in den Zellen J1 bis J5 und die Startwerte in A3, B3, C3 festgelegt, diese Werte werden fast überall in der Tabelle benötigt. In anderen Programmiersprachen würde man von „globalen Variablen“ sprechen. Die eben aufgezählten Formeln werden in benachbarten Spalten der Tabellenkalkulation programmiert, die Zwischenergebnisse werden im Regelfall in weiter rechts liegenden Spalten weiterverarbeitet. Die „Weiterschaltung“ in die folgende Zeile erfolgt dadurch, dass das Ergebnis der Zelle G3 verwendet wird, um den Inhalt der Zelle B4 nach dem folgenden Zeitschritt zu berechnen. Zum Schluss kopiert man die Formeln der 3. bzw. 4. Zeile in die nächsten 2000 Zeilen – gleichzeitig wird das Ergebnis berechnet.

Von ausschlaggebender Wichtigkeit für die physikalische Korrektheit der Ergebnisse ist die sinnvolle Wahl des Zeitschrittes dt, der möglichst klein sein soll und in der nebenstehenden Tabelle den – für diese Aufgabenstellung – recht hohen Wert 0,2 s hat. Das führt in der Umgebung der Zelle G20 zu gerade noch akzeptierbaren Wertesprüngen von etwa 40 %. Allerdings bewirkt auch eine Vergrößerung auf dt = 1 s noch keine gravierenden Änderungen, was die Robustheit dieses Lösungsverfahrens demonstriert.

Im nebenstehenden Bild wird neben der Tabelle die Gesamtbeschleunigung in Abhängigkeit von der Höhe dargestellt. Die überraschenden Ergebnisse:

• Die Meteore werden fast unabhängig von ihrer Masse in etwa 40 km Höhe am stärksten gebremst und können dabei in Bruchstücke zerlegt werden oder verglühen.
• Die Geschwindigkeiten in den letzten Kilometern über der Erdoberfläche betragen stets etwa 40 m/s – wenn die Bruchstücke bis dahin nicht verglüht sind. Der berechnete Geschwindigkeitsverlauf ist im unteren Bild dargestellt.

Weiterführende Untersuchungen

Das beschriebene Verfahren lädt dazu ein, Parameter wie Größe und Anfangsgeschwindigkeit zu variieren und deren Auswirkungen auf die berechneten Ergebnisse zu untersuchen. Diese Art von „experimenteller Mathematik“ kann zu größerem Verständnis der enthaltenen Physik führen als die Auswertung der komplexen Formeln im vorhergehenden Absatz.

Weblinks

Wikibooks: Abschnitt zur Methode der kleinen Schritte in „Informationstechnische Grundlagen (ITG)“

• Methode der kleinen Schritte bei konstant beschleunigter Bewegung (LEIFI)

Literatur

• E. Hairer, S.P. Norsett, G. Wanner: Solving Ordinary Differential Equations I, Springer Verlag
• M. Hermann: Numerik gewöhnlicher Differentialgleichungen, Anfangs- und Randwertprobleme, Oldenbourg Verlag, München und Wien, 2004, ISBN 3-486-27606-9