[zurück]

Finite-Differenzen-Verfahren zum Lösen partieller Differentialgleichungen

(Jörn-Thorsten Paßmann)
[vor]

Elliptische Gleichungen

Gegeben sei eine elliptische Differentialgleichung, die sogenannte Poisson-Gleichung. Sie hat die Form

2u
x2
(x, y) + 2u
y2
(x, y) = f(x, y)

und sei definiert auf einem rechteckigen Gebiet

R = {(x, y) | a < x < b, c < y < d}

mit der Randbedingung (S bezeichnet den Rand des Gebietes)

u(x, y) = g(x, y)   für (x, y) Î S  .

Die Differentialgleichung allein ist nicht eindeutig lösbar; für ein komplettes Problem müssen Rand- und/oder Anfangswertbedingungen hinzukommen. Das hier gegebene Problem ist eindeutig lösbar, sofern f und g auf ihrem Definitionsbereich stetig sind.

Ein Beispiel für ein physikalisches Problem, das durch diese Gleichung beschrieben wird, ist die stationäre Temperaturverteilung in einer dünnen rechteckigen Metallplatte. Die Platte wird durch die x-Koordinaten a und b nach links und rechts und die y-Koordinaten c und d nach unten und oben begrenzt. In diesem Fall gilt f(x, y) = 0 (keine äußere Wärmequelle), was die Gleichung zur sogenannten Laplace-Gleichung vereinfacht. g ist eine bekannte Funktion, die die Temperaturen beschreibt, die am Rand herrschen. u(x, y) ist die gesuchte Temperatur an einer bestimmten Stelle der Platte.

Ziel soll sein, eine Poisson-Gleichung numerisch zu lösen, d. h. die Werte der Funktion u an möglichst dicht verteilten diskreten Punkten näherungsweise zu berechnen.

Diese Punkte werden z. B. ausgewählt, indem das Gebiet r mit einem Netz überzogen wird. Man legt eine Anzahl n von Zerlegungen in x-Richtung und m in y-Richtung fest und definiert die Schrittweiten h := (b - a) / n und k := (d - c) / m. Sodann zieht man horizontale und vertikale Gitternetzlinien durch die Gitterpunkte (xi, yj) mit xi := a + ih und yj := c + jk.

Es sollen nun die Werte u(xi, yj) durch die Werte von u in umliegenden Punkten näherungsweise ausgedrückt werden. Dadurch entsteht ein Gleichungssystem mit den Werten in den Gitterpunkten als Unbekannten. Setzt man die Rand- beziehungsweise Anfangswerte ein, kann man es lösen und erhält Näherungen der gesuchten Funktionswerte in allen Gitterpunkten.

Dazu werden zunächst die zweiten partiellen Ableitungen an der Stelle (xi, yj) mit Hilfe der Taylor-Reihe durch die Werte und Ableitungen von u in (xi, yj) und umliegenden Punkten ausgedrückt. Die zweite partielle Ableitung nach x ergibt sich folgendermaßen:

u(xi + 1, yj) = u(xi, yj) + h
1!
u
x
(xi, yj) + h2
2!
2u
x2
(xi, yj) + h3
3!
3u
x3
(xi, yj) + h4
4!
4u
x4
(x(1)i, yj)   mit x(1)i Î [xi, xi + 1]

und

u(xi - 1, yj) = u(xi, yj) - h
1!
u
x
(xi, yj) + h2
2!
2u
x2
(xi, yj) - h3
3!
3u
x3
(xi, yj) + h4
4!
4u
x4
(x(2)i, yj)   mit x(2)i Î [xi - 1, xi]

Die Gleichungen addiert man und fasst die Restglieder mit O (h4) zusammen:

u(xi + 1, yj) + u(xi - 1, yj) = 2u(xi, yj) + 2h2
2!
2u
x2
(xi, yj) + O(h4) Û 2u
x2
(xi, yj) = u(xi - 1, yj) - 2u(xi, yj) + u(xi + 1, yj)
h2
- O(h2) (1)

Entsprechend ergibt sich für die zweite partielle Ableitung nach y

2u
y2
(xi, yj) = u(xi, yj - 1) - 2u(xi, yj) + u(xi, yj + 1)
k2
- O(k2)  .

Einsetzen in die Poisson-Gleichung führt zu

u(xi - 1, yj) - 2u(xi, yj) + u(xi + 1, yj)
h2
+ u(xi, yj - 1) - 2u(xi, yj) + u(xi, yj + 1)
k2
= f(xi, yj) + O(h2 + k2)

Anschließend wird die Formel diskretisiert, d. h. für die korrekten Werte u(xi, yj) werden die approximierten (genäherten) Werte wij eingesetzt, und das Restglied wird vernachlässigt. Das führt zur zentralen Differenzenmethode mit der Formel

2 é
ê
ë
æ
ç
è
h
k
ö

ø
2

+ 1 ù
ú
û
wij - (wi + 1, j + wi - 1, j) - æ
ç
è
h
k
ö

ø
2

(wi, j + 1 + wi, j - 1) = - h2f(xi, yj)  .

Der durch das Vernachlässigen des Restgliedes entstandene Fehler liegt in der Größenordnung von h2 und k2 und kann durch kleine Schrittweiten klein gehalten werden.

Es empfiehlt sich, jetzt die inneren Gitterpunkte (xi, yj) von links nach rechts und von oben nach unten durchzunummerieren, ebenso wie ihre genäherten Funktionswerte: P1 = (x1, ym - 1) mit dem Näherungswert w1, P2 = (x2, ym - 1) mit w2 und so weiter. Nach der zentralen Differenzenmethode lässt sich für jeden Punkt eine Gleichung aufstellen, die den Wert w durch die umliegenden Werte ausdrückt. Durch Einsetzen der Randbedingungen erhält man ein Gleichungssystem, in dem nur die inneren Gitterpunkte als Unbekannte vorkommen. In diese Gleichungen werden auch alle dem Wert nicht benachbarten Werte eingefügt und mit dem Koeffizienten 0 versehen, so dass das Gleichungssystem so viele Gleichungen und so viele Variablen enthält, wie es innere Gitterpunkte gibt. Das so entstandene Gleichungssystem lässt sich als Matrix schreiben und kann schließlich z. B. mit einem der Verfahren, die Christian Moldenhauer in seinem Beitrag darstellt, aufgelöst werden.

Parabolische Gleichungen

Aufsteigendes Differenzenverfahren

Bei der sogenannten Wärme- oder Diffusionsgleichung handelt es sich um eine parabolische Gleichung. Sie hat die Form

u
t
(x, t) = a2 2u
x2
(x, t)   für 0 < x < l und t > 0

mit den Anfangs- und Randbedingungen:

u(0, t) = u(l, t) = 0 für t > 0 und u(x,0) = f(x) für 0 £ x £ l

Zur numerischen Lösung dieses Problems wird prinzipiell genauso wie im vorigen Abschnitt beschrieben vorgegangen. Die Näherungslösung für (2u) / (x2) (Gleichung (1)) kann (ersetze y durch t) übernommen und in die Wärmegleichung eingesetzt werden.

Die Taylorreihe von u um (xi, tj)

u(xi, tj + 1) = u(xi, tj) + k
1!
u
t
(xi, tj) + k2
2!
2u
t2
(xi, mj)   mit mj Î [tj, tj + 1]     (2)

lässt sich nach der ersten partiellen Ableitung nach t auflösen:

u
t
(xi, tj) = u(xi, tj + 1) - u(xi, tj)
k
- k
2
2u
t2
(xi, mj)     (3)

Einsetzen dieses Terms in die Wärmegleichung und Diskretisieren führen zum aufsteigenden Differenzenverfahren

wi, j + 1 - wij
k
- a2 wi + 1, j - 2wij + wi - 1, j
h2
= 0
Û wi, j + 1 = æ
ç
è
1 - 2a2k
h2
ö

ø
wij + a2 k
h2
(wi + 1, j + wi-1, j)  .

Es handelt sich hierbei um ein explizites Verfahren, da aus den für den ersten Zeitschritt durch die Anfangswertbedingung gegebenen Werten direkt die Werte für den nächsten Zeitschritt berechnet werden können und so weiter. Daher erfordert das Verfahren kein rechenintensives Lösen eines Gleichungssystems.

Die Nachteile liegen zum einen in der Fehlerordnung: Wegen des k im Restglied hängt die Größenordnung von h2 und k ab, was einen sehr kleinen Zeitschritt erfordert, will man hinreichend genaue Ergebnisse erhalten. Zum anderen gilt es, das Stabilitätskriterium a2(k / h2) £ 1/2 zu erfüllen, da sich ansonsten völlig fehlerhafte und unbrauchbare Werte ergeben.

<>Absteigendes Differenzenverfahren

In Gleichung (3) wird (u) / (t)(xi, tj) mit Hilfe von u(xi, tj + 1), dem Funktionswert im folgenden Zeitschritt, ausgedrückt. Berechnet man die erste partielle Ableitung nach der Zeit jedoch unter Verwendung des Wertes im vorangehenden Zeitschritt, u(xi, tj), so erhält man statt der Gleichung (3) folgenden Ausdruck:

u
t
(xi, tj) = u(xi, tj) - u(xi, tj - 1)
k
+ k
2
2u
t2
(xi, mj)   mit mj Î (tj - 1, tj)     (4)

Setzt man dies statt (3) in die Wärmegleichung ein, so erhält man nach Diskretisieren das absteigende Differenzenverfahren

wij - wi, j - 1
k
- a2 wi + 1, j-2wij + wi - 1, j
h2
= 0   .

Es hat ebenfalls den Fehler O(h2 + k), weist jedoch nicht die Stabilitätsprobleme des aufsteigenden Verfahrens auf. Da es ein implizites Verfahren ist, führt es zu einem Gleichungssystem, das gelöst werden muss.

Methode von Crank-Nicolson

Durch Addieren der Gleichungen des aufsteigenden Verfahrens im Schritt j,

wi, j + 1 - wij
k
- a2 wi + 1, j-2wij + wi - 1, j
h2
= 0   ,

und des absteigenden Verfahrens im Schritt j + 1,

wi, j + 1 - wij
k
- a2 wi + 1, j + 1-2wi, j + 1 + wi - 1, j + 1
h2
= 0   ,

entsteht die Gleichung der Methode von Crank-Nicolson:

wi, j + 1 - wij
k
- a2
2
é
ê
ë
wi + 1, j-2wij + wi - 1, j
h2
+ wi + 1, j + 1-2wi, j + 1 + wi - 1, j + 1
h2
ù
ú
û
= 0

Da die Restglieder der Gleichungen (3) und (4) in etwa den gleichen Betrag, aber unterschiedliche Vorzeichen besitzen, heben sie sich gegenseitig auf. Folglich weist dieses Verfahren den günstigen Fehler O(h2 + k2) auf.

Wie ist das? Das aufsteigende Verfahren setzt das (diskretisierte) uxx von heute (Zeitschritt j) in Beziehung mit der Differenz zwischen morgen (Zeitschritt j + 1) und heute. Das absteigende Verfahren setzt das uxx von morgen in Beziehung mit der Differenz zwischen morgen und heute. Gehört die Differenz zwischen morgen und heute eher zu morgen oder zu heute? Na ja, sie gehört zu beiden Zeitpunkten gleich schlecht, nämlich mit dem Fehler O(k). Nur: Daran können wir nichts tun. Wenn wir ut diskretisieren wollen und haben nur die Werte von u in den Zeitpunkten tj und tj + 1 zur Verfügung, gibt es keine bessere Näherung als den Differenzenquotienten. Interessant ist es, die Frage andersrum zu stellen: Für welchen Zeitpunkt t liefert der Differenzenquotient (u(tj + 1) - u(tj)) / k die beste Näherung an ut(t)? Antwort: für den Zeitpunkt genau in der Mitte, t = (tj + tj + 1) / 2, heute um Mitternacht. Dann nämlich läppern sich in der Taylorreihe die Glieder erster Ordnung genau weg, und es bleibt ein Fehler zweiter Ordnung übrig. Wenn das so ist, müssen wir natürlich auch das genäherte uxx heute um Mitternacht auswerten. Wie macht man das, mit den Gitterpunkten, die wir zur Verfügung haben? Man nimmt den Mittelwert zwischen heute und morgen - was sonst? Und siehe da - es entsteht das Verfahren von Crank-Nicolson.

Literatur

J. Douglas Faires, Richard L. Burden, Numerische Methoden, Spektrum Akademischer Verlag 1994


[zurück] [Inhaltsverzeichnis] [vor]