[zurück]

Numerische Lösung von Differentialgleichungen

(Lisa Huber, Christoph Pöppe)
[vor]

Bisher konnten wir unsere Differentialgleichungen noch analytisch lösen. Die weitaus meisten Probleme sind dafür jedoch zu kompliziert. Deshalb stellen wir im Folgenden einige numerische Verfahren vor.

Wir betrachten ein Anfangswertproblem:

y¢(t) = dy(t)
dt
= f(t, y(t))
y(t0) = y0

bei dem wir davon ausgehen, dass es eine eindeutige Lösung besitzt (siehe den Beitrag von Karin Heitzmann). Wir nennen die eindeutige Lösung F(t). Mit yn bezeichnen wir die angenäherten Werte von F(tn) an den Stellen tn = t0 + nh, h = tn + 1 - tn, n = 0,1,2,... .

Hier gilt es höllisch aufzupassen! Diese Differentialgleichung, so wie sie dasteht, ist zunächst mal eine Wunschvorstellung, eine Forderung an das y. Das y haben wir aber erstmal noch nicht. Da steht zwar y¢(t) = f(t, y(t)), aber deswegen ist f noch lange nicht die Ableitung von y. Das gilt erst, wenn y die Lösung ist, und dann sind wir sowieso fertig. Solange wir die Lösung nicht haben, ist f - na ja, eben gar nix mit y. Es ist die rechte Seite der Differentialgleichung; einen schöneren Namen gibt es nicht dafür. Physikalisch gesprochen: f ist das (bekannte) Naturgesetz, und y ist das (vorläufig unbekannte) Verhalten eines Systems, das diesem Naturgesetz folgt und zum Zeitpunkt t0 im Zustand y0 ist. Um Lösung und Weiß-noch-nicht-ob-Lösung säuberlich auseinanderzuhalten, schreiben wir F für ersteres und y für letzteres. Die meisten Lehrbücher machen diese Unterscheidung (und die Praktiker verschlampern sie, weil denen klar ist, was jeweils gemeint ist).

Übrigens muss in der Differentialgleichung das y (bzw. die Lösung F) nicht unbedingt eine zahlenwertige Funktion sein. Alle jetzt folgenden Überlegungen funktionieren genauso, wenn die unbekannte Funktion ein Vektor ist, das heißt aus mehreren Komponenten besteht. Die rechte Seite f bildet dann einen Vektor y samt einem Skalar t auf einen Vektor f(t, y) ab. Physikalisch bedeutet das, dass der Zustand des Systems nicht nur durch eine Zahlengröße beschrieben wird, sondern durch mehrere - was der bei weitem interessantere Fall ist. Diese mehreren Größen - die Komponenten von y - können zum Beispiel die Koordinaten eines oder mehrerer Massenpunkte sein, die durch Kräfte aufeinander einwirken. Anders ausgedrückt: Man hat mehrere unbekannte Funktionen, die auch noch voneinander abhängig sind: ein System von Differentialgleichungen. Eine Gleichung höherer Ordnung kann auf eine einfache Weise in ein System von Gleichungen 1. Ordnung umgewandelt werden. Das angegebene Anfangswertproblem ist also keineswegs so speziell, wie es aussieht, sondern bereits so ziemlich das allgemeinste, was es gibt. Für die bildliche Veranschaulichung, die jetzt kommt, muss man sich allerdings das y doch wieder als einen Skalar vorstellen.

RichtungsfeldAn einem geeigneten Richtungsfeld (hier abgebildet ist das Richtungsfeld für y¢ = e - t - 2y) können wir uns die Vielfalt der Lösungen veranschaulichen und die Problematik der Näherung betrachten.

Die rechte Seite f der Differentialgleichung sagt uns zu jedem Punkt (t, y) in der (t, y)-Ebene, welche Steigung y¢(t) die Lösungskurve haben müsste (nämlich f(t, y)), wenn sie durch den Punkt (t, y) verliefe. (Ob sie wirklich durch diesen Punkt verläuft, wissen wir noch nicht.) Anders ausgedrückt: Wir suchen eine Kurve y(t), die in jedem ihrer Punkte die für diesen Punkt vorgeschriebene Tangente f(t, y(t)) hat.

Durch jeden Punkt der (t, y)-Ebene verläuft genau eine solche Lösungskurve, vorausgesetzt, die rechte Seite f ist nicht zu exotisch. Das sagt uns der Existenz- und Eindeutigkeitssatz (siehe unten). Lösungskurven können sich nicht schneiden; sonst gäbe es im Schnittpunkt zwei verschiedene Steigungen. Es gibt aber nur einen Wert von f(t, y). Das hindert nicht, dass Lösungskurven auseinander- oder zusammenlaufen. Aus der großen Schar der Lösungskurven (der möglichen Verhaltensweisen des Systems) suchen wir diejenige, die durch den Punkt (t0, y0) verläuft (zum Zeitpunkt t0 im Zustand y0 ist).

Es ist ein bisschen wie Autofahren auf einer unendlich vielspurigen Autobahn: Immer schön auf der Spur bleiben, das heißt die Richtung halten, die an der Stelle, wo du gerade bist, angesagt ist. Wenn du nicht aufpasst, machst du einen unwillkürlichen Spurwechsel. Auf einer echten Autobahn kracht es dann meistens. Auf einem Richtungsfeld ist es - na ja, eben ein Fehler.

Somit stellt sich uns vor allem die Aufgabe, einen solchen Fehler möglichst klein zu halten.

Das einfache Eulersche Verfahren

Das erste numerische Verfahren wurde von Euler entwickelt. Der erste Schritt besteht darin, dass wir mit unserem bekannten, exakten Wert y0 durch Einsetzen in die rechte Seite der Differentialgleichung die Steigung in diesem Punkt bekommen. Mit deren Hilfe konstruieren wir die Tangente an die Kurve im Punkt y0 = F(t0). Entlang dieser Tangente f(t0, y0) bewegen wir uns um eine konstante Schrittweite h weiter bis zu dem Punkt y1. Das soll ein Näherungswert für F(t1) sein: y1 = y0 + hf(t0, y0). Jetzt setzen wir wieder den errechneten Wert in die rechte Seite der Differentialgleichung ein, bekommen einen Steigungswert, konstruieren die Tangente und erlangen den nächsten Punkt y2. Allgemein können wir schreiben:

yn + 1 = yn + hf(tn, yn)

mit n = 0, 1, 2, 3, ... . Je kleiner hierbei die Schrittweite, desto geringer ist die Abweichung der Näherung von der exakten Lösung, desto höher allerdings der Rechenaufwand.

Eulerverfahren ist wie ein total übermüdeter Autofahrer: Fährt los, nickt ein, schreckt hoch, stellt das Steuer richtig für die Spur, auf der er gerade ist, pennt wieder eine Sekunde, schreckt hoch, und so weiter. Der hält natürlich um so besser die Spur, je öfter er die Augen aufmacht. Einmal Augen aufmachen ist in diesem Bild aber soviel wie ein kompletter Rechenschritt auf dem Computer - egal wieviel Zeit zwischen zwei "Augenblicken" verstreicht. Deswegen ist für das Rechnen mit dem Computer unpraktikabel, was auf der Autobahn eine gute Idee ist: immer die Augen aufzuhalten.

IntegralEs ist auch möglich, das Eulersche Verfahren mathematisch über das Integral über (t) herzuleiten:

F(tn + 1) - F(tn) = òtntn +1  (t) dt = òtntn +1  f[t, F(t)] dt

Þ F(tn + 1) = F(tn) + òtntn + 1  f[t, F(t)] dt

Beim Euler-Verfahren ersetzen wir den Integranden durch seinen Wert an der Stelle t = tn, eine Konstante, und können ihn somit aus dem Integral herausziehen und über die Konstante 1 integrieren:

F(tn + 1) » F(tn) + f[tn, F(tn)](tn + 1 - tn)

Ersetzt man jetzt F(tn) durch unsere Näherungswerte yn, so ergibt sich wiederum das Euler-Verfahren.

Auch mit der Taylor-Reihe um den Punkt tn (ausgewertet in tn + h = tn + 1), von der wir natürlich ausgehen, dass sie existiert, lässt sich das Eulersche Verfahren finden.

F(tn + h) = F(tn) + (tn)h + F¢¢(t*n)h2
Û F(tn + 1) = F(tn) + f[tn, F(tn)]h + F¢¢( t*n)h2
(t*n ist wieder eine unbekannte Zwischenstelle zwischen tn und tn + 1)
Wenn wir nämlich nach dem 2. Glied abbrechen, d. h. das Restglied nicht berücksichtigen, und für F(tn + 1), F(tn) wieder die Näherungswerte yn + 1, yn verwenden, so ergibt sich hieraus ebenfalls das Eulersche Verfahren. Uns muss jedoch bewusst sein, dass dieses Verfahren Fehler macht, die wir im Folgenden betrachten wollen.

Fehlerabschätzung:

Nennen wir den globalen Approximationsfehler En. Das ist der Unterschied zwischen der exakten Lösung und der angenäherten zum Zeitpunkt tn: En = F(tn) - yn.

Jetzt kommt die Sache mit der Sünde, die euch so erheitert hat. Und zwar: Die exakte Lösungskurve ist der Pfad der Tugend. Aber den kennt man nur im Himmel. Auf Erden sind wir allzumal Sünder, das heißt wir weichen vom Pfad der Tugend ab. Das ist unvermeidlich, wenn man in diskreten Zeitschritten rechnet, und anders können wir es nicht - wenn wir die Lösung nicht auf analytischem Wege finden. Beim Euler-Verfahren geht man halt ein Stück die Tangente entlang statt die Lösungskurve; das ist die Sünde, die man in jedem Schritt aufs Neue begeht. Allgemein verwendet man bei jedem Schritt nur eine Näherungsformel zur Lösungsbestimmung, denn eine exakte Formel gibt es nicht. Außerdem war man aber schon - abgesehen vom allerersten Schritt, als man sich noch im Stande der Unschuld befand - auf der falschen Lösungskurve. Das ist die Folge der früheren Fehltritte: die Erbsünde eben.

Versuchen wir, die Gesamtsünde En in Neusünde und Erbsünde zu zerlegen. Dazu bilden wir die Taylor-Reihe von F(t) um tn mit Restglied 2. Ordnung:

F(tn + h) = F(tn) + (tn)h + F¢¢(t*n)h2

mit einem gewissen t*n zwischen tn und tn + h. Das formen wir wie vorher um und subtrahieren die Formel des Eulerschen Verfahrens. Für die Gesamtsünde zur Zeit tn + 1 ergibt sich

En + 1 = F(tn + 1) - yn + 1 = F(tn) - yn + h{f[tn, F(tn)] - f[tn, yn]} + F¢¢(t*n)h2

An dieser ganzen Sünde ist neu nur der letzte Term; denn wenn wir im n-ten Schritt noch auf dem Pfad der Tugend wären, wäre yn = F(tn), und die beiden ersten Summanden würden wegfallen. Also ist die Neusünde - offizieller Name: der lokale Diskretisierungsfehler -

en + 1 = F(tn + 1) - yn + 1 = 1
2
F¢¢(t*n)h2,

das ist proportional zu h2 und der 2. Ableitung von F(t).

Was haben wir davon? F kennen wir nicht (sonst hätten wir die Lösung und könnten uns das Sündigen sparen) und t*n schon gar nicht. Wenn wir die Sünde nun schon nicht ausrechnen können, wollen wir sie wenigstens in Grenzen halten. Das gelingt auch häufig, denn meistens ist F¢¢(t) beschränkt: |F¢¢(tn)| £ M mit einem gewissen M, das wir zwar meistens auch nicht kennen, das aber wenigstens nicht von h abhängt. Dann ist auch der lokale Diskretisierungsfehler begrenzt: |en + 1| £ [(Mh2) / 2].

Also: Wenn wir die Schrittweite halbieren, begehen wir nur noch ein Viertel der Neusünde - pro Schritt. Um zum selben Ziel zu gelangen, müssen wir jetzt aber doppelt so viele Schritte machen. Allgemein: Wenn wir n Schritte benötigen, um von t0 zu `t = t0 + nh zu gelangen, und bei jedem Schritt der Fehler höchstens Mh2 / 2 beträgt, ist der Fehler nach n Schritten höchstens nMh2 / 2. Da n = (`t - t0) / h ist, erhalten wir daraus als Schranke für die Summe aller Neusünden (`t - t0)Mh / 2. Der Fehler wird also nicht größer als eine Konstante mal h, d. h. er lässt sich durch Übergang zu kleineren Schrittweiten beliebig klein machen.

Diese Argumentation unterstellt, die Fehler jedes einzelnen Schrittes würden sich einfach aufaddieren, berücksichtigt also nicht die Fehlerfortpflanzungseffekte. (Wie würde man die nennen? Erbs-Erbsünde in Analogie zu Zinseszins?) Durch eine etwas kompliziertere Argumentation lässt sich jedoch zeigen, dass der Fehler En tatsächlich abschätzbar ist durch h mal eine gewisse Konstante; nur ist die Konstante etwas größer als die oben angegebene (`t - t0)M / 2.

Das verbesserte Eulersche Verfahren (Heun-Formel)

IntegralBeim einfachen Eulerschen Verfahren wird ein Integral durch eine Rechtecksfläche angenähert (siehe oben). Die Höhe des Rechtecks ist dabei der Funktionswert f[tn, F(tn)] am linken Rand des zu untersuchenden Intervalls. Das verbesserte Eulersche Verfahren approximiert den Integranden genauer, und zwar durch eine Trapezfläche, indem man den Mittelwert zwischen den Eckpunktwerten (f[tn, F(tn)] + f[tn + 1, F(tn + 1)]) / 2 nimmt. Wieder ersetzen wir die exakten durch die Näherungswerte: F(tn + 1), F(tn) durch yn + 1, yn. Daraufhin erhalten wir yn + 1 = yn + (h / 2) (f[tn, yn] + f[tn + 1, yn + 1]). Da die Unbekannte yn + 1 jedoch in der rechten Seite der Gleichung vorkommt, ersetzen wir sie hier durch den Wert den wir mittels der Euler-Formel erhalten. Folglich gilt:

yn + 1 = yn + h
2
( f[tn, yn] + f[tn + h, yn + hf(tn, yn)])

Das ist das Euler-Verfahren mit Beichte! Wir machen zuerst einen Fehltritt, das heißt, wir berechnen einen vorläufigen Wert für yn + 1 nach dem Euler-Verfahren. Im Lichte der Erfahrung, die wir dadurch gewonnen haben, bringen wir eine Korrektur an und machen dadurch unseren Fehler wenigstens zum Teil wieder gut. Allgemein heißen Verfahren dieses Typs (deren es viele gibt) Prädiktor-Korrektor-Verfahren.

Es lässt sich beweisen, dass der lokale Diskretisierungsfehler bei diesem Verfahren durch eine Konstante mal h3 und der globale auf einem endlichem Intervall durch eine Konstante mal h2 beschränkt ist, d. h. die Fehlergrenzen werden enger. Allerdings erfordert eine ausreichend komplizierte Funktion f einen beträchtlichen Rechenaufwand. Für das Heun-Verfahren müssen wir f doppelt so oft auswerten wie für das einfache Verfahren; oder: Das Heun-Verfahren ist so aufwendig wie ein Euler-Verfahren mit der halben Schrittweite. Es bringt allerdings mehr an Genauigkeit.

Überhaupt sollte kein schiefes Bild entstehen! Das Euler-Verfahren ist zwar das einzige, das wir uns etwas ausführlicher ankucken; aber verglichen mit anderen, tatsächlich eingesetzten Verfahren ist es so elend schlecht, dass es nur als abschreckendes Beispiel taugt. Eigentlich hat es den Namen des großen Leonhard Euler (eines der größten Mathematiker überhaupt) gar nicht verdient. Es ist nur nützlich, weil es so einfach zu verstehen ist und trotzdem bereits die wesentlichen Eigenschaften aller Verfahren für Anfangswertprobleme zeigt.

Eine wesentliche Eigenschaft ist: Wenn es um die Qualität eines Verfahrens geht, kommt es auf das genaue Ausmaß der Sünde gar nicht besonders an. Man macht sich in der Regel auch nicht die Mühe, die Konstante M (bzw. die etwas größere, die eigentlich korrekt ist) auszurechnen, selbst wenn man die Daten zur Verfügung hat. Das einzige, was bei der ganzen Fehlerabschätzerei am Ende interessiert, ist der Exponent an dem h. Das ist die Ordnung des Verfahrens. Das Euler-Verfahren hat Ordnung 1, das Heun-Verfahren Ordnung 2. Eine Halbierung der Schrittweite vermindert den zu befürchtenden Fehler bei einem Verfahren 1. Ordnung auf die Hälfte, bei 2. Ordnung auf ein Viertel, bei 3. Ordnung auf ein Achtel. Je höher die Ordnung des Verfahrens, desto besser das Ergebnis (desto höher im Allgemeinen auch der Aufwand).

Es ist eine gute Idee, dasselbe Problem mit zwei verschiedenen Schrittweiten zu rechnen. Der Unterschied in den Ergebnissen gibt einem - ohne dass man theoretischen Aufwand treiben müsste - einen Hinweis darauf, wie groß das M in diesem Fall ist. Daraus kann man wiederum ausrechnen, wie groß man das h wählen muss, damit der globale Fehler unter einer vorgegebenen Schranke bleibt. Und das ist die Standard-Forderung an jedes gute Verfahren.

Immer nur mit entsetzlich kleiner Schrittweite h die Zeitachse lang zu tippeln bringt nichts - nur Rechenaufwand: Durch die guten Werke allein wirst du nicht selig. Die Kunst besteht darin, die Schritte so groß zu nehmen, wie es der globalen Genauigkeitsforderung gerade noch vereinbar ist. Das kann zu verschiedenen Zeiten verschieden sein: Wenn das System in heftiger Bewegung ist (scharfe Kurve auf der gedachten Autobahn), muss man sehr häufig hinkucken, damit man alle Einzelheiten mitkriegt. Auf langen Geraden, wo nahezu nichts passiert, genügt ab und zu ein verschlafenes Blinzeln. Hochentwickelte Verfahren beherrschen die Kunst, die Schrittweite den Verhältnissen automatisch anzupassen.

Wie sehen andere Verfahren aus? Ein paar Stichworte müssen genügen. Man kann die Taylor-Entwicklung weiter treiben als nur bis zum zweiten Glied, muss dann aber auch Werte aus früheren Zeitschritten mitverwenden, weil sonst keine eindeutigen Formeln zustandekommen. Das läuft darauf hinaus, dass man die Lösungskurve nicht durch die Tangente (die in diesem Punkt anschmiegsamste Gerade) nähert, sondern durch die anschmiegsamste Parabel, kubische Parabel und so weiter. Das sind die Mehrschrittverfahren höherer Ordnung. Oder man treibt das Spiel aus Sünde und Beichte mehrmals hintereinander über gewisse Teil-Zeitschritte. Das läuft auf die sogenannten Runge-Kutta-Verfahren hinaus. Oder man rechnet von t0 bis `t mit verschiedenen Schrittweiten h1, h2, h3..., bekommt verschiedene Werte für F(`t) in Abhängigkeit vom jeweiligen Wert von h und extrapoliert daraus einen Wert für die (himmlische) Schrittweite h = 0. Das sind die Extrapolationsverfahren. Alle diese Verfahren kann man mit automatischer Schrittweitensteuerung ausstatten; das funktioniert am elegantesten bei den Extrapolationsverfahren.

Literatur

W. E. Boyce/R. C. DiPrima: Gewöhnliche Differentialgleichungen


[zurück] [Inhaltsverzeichnis] [vor]