Numerisches Lösen von Differentialgleichungen

Behandelt wird eine DGL der Form mit der Anfangsbedingung y(x0) = y0.

Die rechte Seite F(x,y) kann als "Steigungslieferant" für die Lösung y aufgefaßt werden in dem Sinn: Wenn eine Lösung durch P(x/y) verläuft, so gibt F(x,y) die Steigung der Lösung an. Trägt man für Punkte der Ebene diese "potentiellen Steigungen" als kleine Tangentenstückchen an, so erhält man ein sogenannte Richtungsfeld der DGL. Als Beispiel sei das Richtungsfeld der Gleichung y’ = –x× y angegeben.

(Programm: RICHTFEL.PAS)

Mit etwas Phantasie ist die Gaußsche Glockenkurve schon jetzt erkennbar.

Für die numerische Lösung ist es typisch, daß die gesuchte Funktion y nur an diskret liegenden aufeinander folgenden Stellen xi berechnet wird. Das heißt, es werden Folgen x0, x1, x2, x3, ... von Stellen und y0, y1, y2, y3,... von Funktionswerten rekursiv berechnet.

Es sollen nun sogenannte Einschritt-Verfahren angegeben werden, bei denen ein Wert yi+1 nur von yi , xi und der Schrittweite h = xi+1 – xi abhängt.

Grundlegende Rechnung:

Die Aufgabe, aus dem Funktionswert y(x) einen Funktionswert y(x+h) zu berechnen, wird also auf die Berechnung eines Integrals zurückgeführt.

1. Das Euler-Cauchy-Verfahren

Die einfachste Lösung ist, das Integral durch den Inhalt eines Rechtecks der Höhe F(x,y(x)) anzunähern:

Die Gleichung zum Euler-Cauchy-Verfahren
(1 Auswertung von F, 1 Multiplikation, 1 Addition):

Obige integrale Veranschaulichung für F(x,y) kann durch eine differentiale für y(x) ergänzt werden:

2. Das Verfahren von Heun

Viel besser als die Rechtecknäherung ist die Näherung durch ein Trapez. Dafür ist aber im Fall der Sekanten-Trapez-Regel der Funktionswert F(x+h,y(x+h)) vonnöten. Da y(x+h) erst berechnet werden soll, wird die lineare Näherung

benutzt, in der die Ableitung durch die rechte Seite der DGL ersetzt werden kann. Als Vorhersage (Prädiktor) wird also genau die Euler-Cauchy-Näherung benutzt.

(Programm: HEUN.PAS)

Die Skizze gibt maßstabsgerecht das Arbeiten des Heun-Verfahrens für die Gleichung y’ = –x× y , x = – 0,9 und h = 0,8 wieder.
Zuerst wird mit Hilfe von F(x,y(x)) die Steigung y’(x) ermittelt und damit der Prädiktor berechnet. Hier hat (genauer: die Lösung, die durch P(x+h/) verläuft) an der Stelle x+h eine so ähnliche Steigung wie y, daß die Funktionswerte F(x+h,y(x+h)) und F(x+h,) fast zusammenfallen. Das kleine Quadrat gibt die Lage der mit dem Heunverfahren berechneten Näherung für y(x+h) wieder.
Es fällt auf, daß diese Näherung trotz der "riesigen" Schrittweite schon sehr gut ist.

Die Gleichungen des Heun-Verfahrens
(3 Auswertungen von F, 2 Multiplikationen, 4 Additionen):

3. Das Runge-Kutta-Verfahren

Während das Verfahren von Heun aus der Sekanten-Trapez-Regel hervorgeht, erhält man aus der Simpson-Regel

das Runge-Kutta-Verfahren. Die Simpsonregel hat die Stützstellen , so daß ähnlich wie beim Verfahren von Heun Schätzwerte für benötigt werden. Diese Schätzwerte werden folgendermaßen gewonnen:

In einem ersten Schritt wird die Steigung m0 = F(x,y(x)) berechnet. Die Tangente g0 führt zu einem Schätzwert für. Mit Hilfe von y1 wird nun die Steigung m1 = F(x+h,y1) berechnet. m1 ist eine Näherung für . Bei etwa gleichbleibender Krümmung führt g0 zu einem zu großen Funktionswert y1 und g1, das ist die Gerade durch (x/y(x)) mit der Steigung m1, zu einem zu kleinen Funktionswert oder umgekehrt. 2× m1 + 2 × m2 ()ist also eine sehr gute Näherung für . Die Steigung m2 schließlich wird zur Konstruktion von g2 und diese zur Ermittlung von benutzt.

Gleichungen des Runge-Kutta-Verfahrens
(4 Auswertungen von F, 4 Multiplikationen (/2,× 2 nicht gezählt), 10 Additionen):

Graphisches Beispiel:

(Programm: RKUTTA.PAS)

Die Skizze gibt maßstabsgerecht das Arbeiten des Runge-Kutta-Verfahrens für die Gleichung y’ = –x× y , x = – 0,9 und h = 0,8 wieder. Trotz der großen Schrittweite ist graphisch kein Fehler erkennbar.

4. Vergleich von Genauigkeit und Aufwand

Eine Bewertung der Verfahren kann z.B. durch Vergleich mit der exakten Lösung geschehen. Im Folgenden wird die Gleichung y’ = –x × y mit den Anfangswerten x0=0, y0=1 zur Berechnung von (Gaußsche Glockenkurve) benutzt. Untersucht man die Abhängigkeit des relativen Fehlers r von der Schrittweite h, so vermutet man (z.B. Experimente mit DGLNUM1.PAS) einen Zusammenhang r(h) = c× hp, (p heißt Ordnung des Verfahrens). Eine doppeltlogarithmische Darstellung ( ) ermöglicht eine einfache Bestimmung von p. Die Diagramme (DGLNUM3.PAS) zeigen auch den Einfluß der Rundungsfehler, die eine Verkleinerung der Schrittweite unter einen - von Verfahren und Genauigkeit der verwendeten "Rechnerzahlen" abhängigen - optimalen Wert sinnlos machen

mit 11 Stellen (Datentyp real):

Den Steigungsdreiecken entnimmt man, daß das Euler-Verfahren die Ordnung p=1 (~h1), das Heun-Verfahren die Ordnung p=2 (~h2) und das Runge-Kutta-Verfahren die Ordnung p=4 (~h4) haben.

 

mit 19 Stellen (Datentyp extended):