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.:
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
Für die Approximation der Lösung von (1) am Gitter,
u(xj, tn), schreiben wir kurz ujn, 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 ¶2 u / ¶x2
zu erhalten, betrachten wir folgende Taylor-Entwic
klungen von u um einen
inneren Gitterpunkt xj:
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 ¶2 u / ¶x2 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 ¶2 u / ¶t2 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 xj des
Ortsgitters zu den Zeiten tn und tn-1, so liefert das
Schema (7) eine Näherungslösung ujn+1 für alle inneren
Punkte xj zum nächsten Zeitpunkt tn+1, vgl. Abb. 2.
Zur Vervollständigung des Differenzenschemas (7) fehlen aber noch
die Informationen aus den Anfangs- und Randbedingungen:
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(xj, 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
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
cA 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. cA = Dx / Dt.
Die CFL-Bedingung, cA = 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.
File translated from
TEX
by
TTH,
version 3.06.
On 26 Apr 2002, 17:29.