Numerische Lösung der 1D-Wellengleichung

Als Beispiel zum Anfangswertproblem bei partiellen Differentialgleichungen behandeln wir die 1D-Wellengleichung

¶2 u(x, t) ¶t2    =    c2    ¶2 u(x, t) ¶x2

(1)

für 0 < x < L und t > 0. Das ist eine lineare partielle Differentialgleichung 2. Ordnung (vom hyperbolischen Typ) für die unbekannte Fun ktion u(x, t). Um eine eindeutige Lösung von (1) zu erhalten, müssen noch zusätzliche Bedingungen vorgegeben sein, z.B.:

ì ï í ï î u(x, 0) = f(x)  ¶u(x, t) ¶t ê ê t=0  = g(x) ü ï ý ï þ        0 £ x £ L

(2)

ì í î u(0, t) = 0 u(L, t) = 0 ü ý þ        t ³ 0

(3)

Mit den Anfangsbedingungen (2) und den Randbedingungen (3) beschreibt die PDGL (1) die Bewegung einer idealisierten Saite der Länge L, die an den Rändern bei x = 0 und x = L eingespannt ist und zur Zeit t = 0 die Anfangsauslenkung f(x) und die Anfangsgeschwindigkeit g(x) hat, Abb. 1. Die Lösungsfunktion u(x, t) ist dann die vertikale Auslenkung der Saite am Ort x zur Zeit t und der Parameter c stellt die Ausbreitungsgeschwindigkeit einer Störung (``Welle'') dar.

Zur numerischen Lösung der Wellengleichung (1) müssen sowohl die Zeitvariable t als auch die Ortsvariable x diskretisiert werden, d.h. die Zeit- und Ortsableitungen müssen aus Funktionswerten an vorgegebenen Stützstellen berechnet werden. Diese Vorgangsweise ist natürlich nur eine Approximation, und die Genauigkeit einer Näherungslösung von (1) hängt sowohl vom gewählten Zeitschritt Dt als auch von der gewählten Ortsauflösung Dx ab. Wir betrachten die Diskretisierung von (1) auf einem äquidistanten raum-zeitlichen Gitter mit den Stützstellen

ì í î xj = j ·Dx        j = 0, 1, ¼, N tn = n ·Dt        n = 0, 1, 2, ¼

(4)

Für die Approximation der Lösung von (1) am Gitter, u(x_j, t_n), schreiben wir kurz u_jⁿ, wobei sich der untere Index j auf den Ort und der obere Index n auf die Zeit bezieht.

Um eine Näherung für die Ortsableitung ¶² u / ¶x² zu erhalten, betrachten wir folgende Taylor-Entwic klungen von u um einen inneren Gitterpunkt x_j:

u(xj + Dx) = u(xj) + Dx  u¢(xj) +  1 2 (Dx)2 u¢¢(xj) +  1 6 (Dx)3 u¢¢¢(xj) + O((Dx)4)

u(xj - Dx) = u(xj) - Dx  u¢(xj) +  1 2 (Dx)2 u¢¢(xj) -  1 6 (Dx)3 u¢¢¢(xj) + O((Dx)4)

Addition der beiden letzten Gleichungen ergibt

u(xj + Dx) + u(xj - Dx) = 2 u(xj) + (Dx)2 u¢¢(xj) + O((Dx)4)

Daraus erhält man als Approximation von ¶² u / ¶x² eine Zentraldifferenz 2. Ordnung in Dx:

æ è  ¶2 u ¶x2 ö ø n j    =    uj+1n - 2 ujn + uj-1n (Dx)2  +  O((Dx)2)

(5)

Analog bekommt man als Näherung von ¶² u / ¶t² eine Zentraldifferenz 2.& nbsp;Ordnung in Dt:

æ è  ¶2 u ¶t2 ö ø n j    =    ujn+1 - 2 ujn + ujn-1 (Dt)2  +  O((Dt)2)

(6)

Einsetzen von (5) und (6) in die Wellengleichung (1) liefert die Diskretisierung

1 (Dt)2 ( ujn+1 - 2 ujn + ujn-1 )    =     c2 (Dx)2 ( uj+1n - 2 ujn + uj-1n )

und mit der Courant-Zahl  a := c  (Dt / Dx)   das Differenzenschema

ujn+1    =   2  ( 1 - a2 )  ujn  +   a2  ( uj+1n + uj-1n )  -   ujn-1   .

(7)

(7) ist ein explizites Verfahren 2. Ordnung in Dt und Dx. Kennt man die Lösung u für alle inneren Punkte x_j des Ortsgitters zu den Zeiten t_n und t_n-1, so liefert das Schema (7) eine Näherungslösung u_jⁿ⁺¹ für alle inneren Punkte x_j zum nächsten Zeitpunkt t_n+1, vgl. Abb. 2. Zur Vervollständigung des Differenzenschemas (7) fehlen aber noch die Informationen aus den Anfangs- und Randbedingungen:

(a)       u(x, 0) = f(x)   , 0 £ x £ L :       uj0 = f(xj)   , j = 0, 1, ¼, N (b)       ut(x, 0) = g(x)   , 0 £ x £ L :       ut(xj, 0) = g(xj)   , j = 0, 1, ¼, N (c)       u(0, t) = 0   , t ³ 0 :       u0n = 0   , n = 0, 1, 2, ¼ (d)       u(L, t) = 0   , t ³ 0 :       uNn = 0   , n = 0, 1, 2, ¼

Da man beim expliziten Verfahren (7) zwei Startwerte benötigt, implementiert man die zweite Anfangsbedingung (b) durch Einführen eines fiktiven Zeitschrittes t_-1:

Verwendet man zur Approximation von ¶u(x_j, t) / ¶t eine Zentraldifferenz 2. Ordnung,

¶u(xj, t) ¶t    =     u(xj, t + Dt) - u(xj, t - Dt) 2 Dt  +   O((Dt)2)   ,

erhält man speziell für t = 0

¶u(xj, t) ¶t ê ê t=0     =     uj1 - uj-1 2 Dt    =   g(xj)   ,

oder

uj-1    =   uj1  -  2  Dt  g(xj)   .

(8)

Für n = 0 lautet das ursprüngliche Differenzenschema (7):

uj1    =   2  (1 - a2)  uj0 + a2  (uj+10 + uj-10) - uj-1   .

Einsetzen von (8) in die letzte Gleichung liefert eine Implementation der zweiten Anfangsbedingung (b) in Form eine s speziellen Differenzenschemas für den ersten Zeitschritt n = 1:

uj1    =   ( 1 - a2 )  uj0  +    1 2  a2  ( uj+10 + uj-10 )  +   Dt  g(xj)   .

(9)

Die Frage nach der Stabilität des expliziten Verfahrens (7) wird durch die Courant-Friedrichs-Lewy (CFL) Bedingung beantwortet: (7) ist bedingt stabil für

a   =   c    Dt Dx    £   1   .

(10)

Dx und Dt können also nicht unabhängig voneinander gewählt werden: Jede Verfeinerung der räumlichen Auflösung Dx erfordert auch eine Verkleinerung des Zeitschrittes Dt. Die Bedeutung der CFL-Bedingung (10) liegt darin, daß sie eine Beziehung herstellt zwischen der Informations-Ausbreitungsgeschwindigkeit c_A des Algorithmus und der physikalischen Ausbreitungsgeschwindigkeit c einer Lösung der Wellengleichung (1). Nach Abb. 2 wird Information durch das Differenzenverfahren um höchstens eine räumliche Gittermasche pro Zeitschritt transportiert, d.h. c_A = Dx / Dt. Die CFL-Bedingung, c_A = Dx / Dt ³ c besagt also, daß die Informations-Ausbreitungsgeschwindigkeit des Algorithmus mindestens so groß sein muß wie die physikalische Ausbreitungsgeschwindigkeit c.

Anmerkung:
Am Stabilitätslimit a = 1 bzw.  Dx = c  Dt liefert das explizite Verfahren (7) - abgesehen von Rundungsfehlern - die exakte Lösung von (1).

Abbildung 1: Anfangs- und Randbedingungen.

Abbildung 2: Informationsfluß beim expliziten Verfahren.