Fallstudie: 1D-Wellengleichung - Wissenschaftliches Rechnen (i788)

Betrachten das Lösen er Wellengleichung in einer Raumdimension mittels Finite-Differenzen-Verfahren genauer. Sei dazu $B=[0,a]$ das betrachtete Gebiet, sei $[0,T]$ das betrachtete Zeitintervall und sei $c>0$ die Ausbreitungsgeschwindigkeit der Wellen. Gesucht ist eine Funktion $u:[0,a]\times[0,T]\to\bbR$ , die

c^2\,u_{xx}(x,t)-u_{tt}(x,t)=0

(1)

erfüllt. Eine Anfangsbedingung sei durch

u(x,0)=f(x),\qquad u_t(x,0)=0

(2)

mit $f:[0,a]\to\bbR$ gegeben.

8.1Randbedingungen¶

Wir betrachten verschiedene Randbedingungen:

beidseitig feste Einspannung (Dirichlet),
beidseitig keine Übertragung von Vertikalkräften auf den Rand (Neumann),
linkes Ende ohne vertikale Kraftübertragung (Neumann), rechtes Ende frei.

Im ersten Fall soll also

u(0,t)=u(a,t)=0

(3)

gelten.

Im zweiten Fall sollen keine Vertikalkräfte auf den Rand übertragen werden. Stellt man sich die Funktion $u$ als Verlauf einer Gitarrensaite vor, soll die Zugkraft also nur horizontal auf die Befestigung am Rand wirken. In vertikale Richtung ist die Befestigung frei beweglich. Erhaltene aus diesen Überlegungen

u_x(0,t)=u_x(a,t)=0.

(4)

Im dritten zu betrachtenden Fall soll die Welle am rechten Rand frei auslaufen können. Das Verhalten am Rand soll also so sein, als würde der betrachtete Bereich $B$ nach rechts beliebig groß sein, wir aber nur den Ausschnitt $[0,a]$ sehen. Zusammen mit der Neumann-Bedingung für den linken Rand erhält man

u_x(0,t)=0,\quad u_t(a,t)+c\,u_x(a,t)=0.

(5)

Dass diese Randbedingung für unsere Zweck die richtige ist, sieht man, indem man $u(x,t)=g(x-c\,t)$ für eine beliebige Funktion $g$ einsetzt. Diese Wahl von $u$ löst die Wellengleichung einschließlich Randbedinung am rechten Rand und entspricht einer nach rechts laufenden Welle (IDVID 810). Nach rechts laufende Wellen können also ungehindert den betrachteten Bereich verlassen. Verzichtet man auf diese Randbedingung, so ist die Lösung der Wellengleichung nicht mehr eindeutig bestimmt.

8.2Diskretisierung¶

8.2.1Gitterpunkte¶

Setzen

\Delta x:=\frac{a}{n_x},\qquad\Delta t:=\frac{T}{n_t}

(6)

und erhalten so Gitterpunkte $(x_i,t_j)$ mit

x_i:=i\,\Delta x,\qquad t_j:=j\,\Delta t

(7)

für $i=0,1,\ldots,n_x$ und $j=0,1,\ldots,n_t$ .

Gesucht sind die Werte

u_{i,j}:=u(x_i,t_j)

(8)

von $u$ an den Gitterpunkten.

Die Gitterwerte der Funktion $f$ in der Anfangsbedingung bezeichnen wir entsprechend mit

f_i:=f(x_i).

(9)

8.2.2Innere Punkte¶

Für Punkte im Inneren des betrachteten Gebietes soll die Wellengleichung gelten, wobei wir die zweiten Ableitungen durch entsprechende Differenzenquotienten ersetzen:

c^2\,\frac{u_{i-1,j}-2\,u_{i,j}+u_{i+1,j}}{(\Delta x)^2}-\frac{u_{i,j-1}-2\,u_{i,j}+u_{i,j+1}}{(\Delta t)^2}=0

(10)

für $i=1,\ldots,n_x-1$ und $j=1,\ldots,n_t-1$ .

Umstellen nach $u_{i,j+1}$ liefert

u_{i,j+1}=\left(\frac{c\,\Delta t}{\Delta x}\right)^2\,(u_{i-1,j}-2\,u_{i,j}+u_{i+1,j})-(u_{i,j-1}-2\,u_{i,j}).

(11)

Sind für ein festes $j$ die $u_{i,j-1}$ und die $u_{i,j}$ für alle $i=0,\ldots,n_x$ bekannt, so erhält man aus dieser Formel die $u_{i,j+1}$ für $i=1,\ldots,n_x$ . Wir können $u$ also Zeitschritt für Zeitschritt berechnen ohne ein Gleichungssystem lösen zu müssen.

Die Obergrenze $T$ für die Zeit spielt hier keine Rolle. Theoretisch kann die Iteration beliebig lang laufen.

8.2.3Anfangsbedingung¶

Wir haben

u_{i,0}=f_i

(12)

für $i=0,\ldots,n_x$ .

Die Ableitungen $u_t(x,0)$ ersetzen wir durch die zentrale Differenz

u_t(x_i,0)\approx\frac{u_{i,1}-u(x_i,-\Delta t)}{2\,\Delta t}.

(13)

Im Vergleich zu einem einseitigen Differenzenquotient erscheint dies sinnvoller, da der Zeitbereich $[0,T]$ nur einen Ausschnitt aus einem größeren Zeitbereich darstellt, der beobachtete Prozess also auch schon vor $t=0$ läuft. Andererseits bringt diese Wahl das (scheinbare) Problem mit sich, dass wir $u(x_i,-\Delta t)$ nicht kennen. Aus der kontinuierlichen Anfangsbedingung $u_t(x,0)=0$ wird nun die diskrete Anfangsbedingung

u_{i,1}=u(x_i,-\Delta t),\quad i=0,\ldots,n_x.

(14)

In Formel (11) für $j=0$ eingesetzt erhalten wir für den ersten Zeitschritt die Formel

u_{i,1}=\left(\frac{c\,\Delta t}{\Delta x}\right)^2(u_{i-1,0}-2\,u_{i,0}+u_{i+1,0})-({\color{red}u_{i,1}}-2\,u_{i,0}).

(15)

Auflösen nach $u_{i,1}$ liefert nun

u_{i,1}=\frac{1}{2}\,\left(\frac{c\,\Delta t}{\Delta x}\right)^2\,(u_{i-1,0}-2\,u_{i,0}+u_{i+1,0})+\,u_{i,0}.

(16)

sodass die unbekannten Werte $u(x_i,-\Delta t)$ gar keine Rolle spielen.

8.2.4Randbedingungen¶

8.2.4.1Dirichlet¶

Die diskretisierten Dirichlet-Randbedingungen sind

u_{0,j}=u_{n_x,j}=0,\qquad j=0,\ldots,n_t.

(17)

8.2.4.2Neumann links¶

Für die Neumann-Randbedingung links verwenden wir analog zur Zeitableitung in der Anfangsbedingung die zentrale Differenz

u_x(0,t_j)\approx\frac{u_{1,j}-u(-\Delta x,t_j)}{2\,\Delta x},

(18)

welche zur Bedingung

u_{1,j}=u(-\Delta x,t_j),\qquad j=0,\ldots,n_t

(19)

führt. Betrachten wir nun (11) für $i=0$ und setzen dort diese Beziehung ein, so erhalten wir mit

u_{0,j+1}=\left(\frac{c\,\Delta t}{\Delta x}\right)^2\,({\color{red}u_{1,j}}-2\,u_{0,j}+u_{1,j})-(u_{0,j-1}-2\,u_{0,j}).

(20)

eine Formel zur Berechnung des Randwertes $u_{0,j+1}$ , welche für $j\geq 1$ funktioniert. Für $j=0$ kann mittels Einsetzen von (14) für $i=0$ und Auflösen nach $u_{0,1}$ daraus widerum

u_{0,1}=\frac{1}{2}\,\left(\frac{c\,\Delta t}{\Delta x}\right)^2\,(u_{1,0}-2\,u_{0,0}+u_{1,0})+\,u_{0,0}

(21)

gewonnen werden.

8.2.4.3Neumann rechts¶

Analog zum linken Rand erhält man für den rechten Rand die Berechnungsvorschriften

u_{n_x,j+1}=\left(\frac{c\,\Delta t}{\Delta x}\right)^2\,(u_{n_x-1,j}-2\,u_{n_x,j}+u_{n_x-1,j})-(u_{n_x,j-1}-2\,u_{n_x,j}).

(22)

für $j\geq 1$ und

u_{n_x,1}=\frac{1}{2}\,\left(\frac{c\,\Delta t}{\Delta x}\right)^2\,(u_{n_x-1,0}-2\,u_{n_x,0}+u_{n_x-1,0})+\,u_{n_x,0}.

(23)

8.2.4.4Freier Rand rechts¶

Im Prinzip können die beiden Ableitungen in der Bedingung $u_t(a,t)+c\,u_x(a,t)=0$ wie oben durch zentrale Differenzenquotienten ersetzt werden. Die auftretenden Funktionswerte außerhalb der Orts- und Zeitintervalle $[0,a]$ und $[0,T]$ können dann wieder eliminiert werden, was allerdings recht aufwendig ist. Der Einfacheit halber beschränken wir uns hier auf einseitige Differenzenquotienten (vorwärts in der Zeit, rückwärts im Ort):

u_t(a,t_j)+c\,u_x(a,t_j)\approx\frac{u_{n_x,j+1}-u_{n_x,j}}{\Delta t}+c\,\frac{u_{n_x,j}-u_{n_x-1,j}}{\Delta x}.

(24)

Nullsetzen und Umstellen nach $u_{n_x,j+1}$ liefert

u_{n_x,j+1}=u_{n_x,j}-\frac{c\,\Delta t}{\Delta x}\,(u_{n_x,j}-u_{n_x-1,j})

(25)

für $j\geq 0$ .

8.2.4.5Zusammenfassung¶

Zur Berechnung aller $u_{i,j}$ bietet sich nun folgendes Vorgehen an:

Für $j=0$ die Werte $f_{i}$ verwenden.
Für $j=1$ und $i=1,\ldots,n_x-1$ die aus der Anfangsbedingung abgeleitete Formel für innere Punkte verwenden.
Für $j=1$ und $i\in\{0,n_x\}$ die Formeln aus den Randbedingungen nutzen.
Für $j=2,3,\ldots$ :
1. Für $i=1,\ldots,n_x-1$ die Formel für innere Punkte nutzen.
2. Für $i\in\{0,n_x\}$ die Formeln aus den Randbedingungen nutzen.

8.3Kondition¶

Man kann zeigen, dass die Berechnung der Funktion $u$ auf dem Gitter gut konditioniert ist, wenn

\frac{c\,\Delta t}{\Delta x}\leq 1

(26)

gilt. Der Wert auf der linken Seite heißt auch Courant-Zahl.

8.4Konvergenz¶

Man kann zeigen, dass für

\frac{c\,\Delta t}{\Delta x}=1

(27)

die auf dem Gitter (ohne Rundungsfehler) berechneten Werte für $u$ exakt sind. Werden $\Delta x$ und $\Delta t$ gleichmäßig verkleinert, so konvergiert die Lösung des diskretisierten Problems also gegen die exakte Lösung.

Ist die Courant-Zahl kleiner als 1, so treten im Laufe der Zeit immer größere Abweichungen zwischen diskretisierter und exakter Lösung auf.

8.5Implementierung¶

Die Implementierung kann entsprechend den Formeln oben erfolgen.

import numpy as np
import matplotlib.pyplot as plt

def wave1d(a, T, c, f, nx, nt, left='d', right='d'):

    delta_x = a / nx
    delta_t = T / nt
    C = c * delta_t / delta_x
    print(f'Courant-Zahl: {C:.2f}')
    
    x = np.linspace(0, a, nx + 1)
    u = np.empty((nt + 1, nx + 1), dtype=float)

    # Anfangsbedingungen
    u[0, :] = f(x)
    u[1, 1:-1] = 0.5 * C ** 2 * (u[0, :-2] - 2 * u[0, 1:-1] + u[0, 2:]) + u[0, 1:-1]
    if left == 'd':
        u[1, 0] = 0
    else:
        u[1, 0] = 0.5 * C ** 2 * (2 * u[0, 1] - 2 * u[0, 0]) + u[0, 0]
    if right == 'd':
        u[1, -1] = 0
    elif right == 'n':
        u[1, -1] = 0.5 * C ** 2 * (2 * u[0, -2] - 2 * u[0, -1]) + u[0, -1]
    else:
        u[1, -1] = u[0, -1] - C * (u[0, -1] - u[0, -2])

    # Zeitschritte
    for j in range(1, nt):

        # innere Punkte
        u[j + 1, 1:-1] = C ** 2 * (u[j, :-2] - 2 * u[j, 1:-1] + u[j, 2:]) \
                         - (u[j - 1, 1:-1] - 2 * u[j, 1:-1])

        # Randbedingungen
        if left == 'd':
            u[j + 1, 0] = 0
        else:
            u[j + 1, 0] = 0.5 * C ** 2 * (2 * u[j, 1] - 2 * u[j, 0]) \
                          - (u[j - 1, 0] - 2 * u[j, 0])
        if right == 'd':
            u[j + 1, -1] = 0
        elif right == 'n':
            u[j + 1, -1] = 0.5 * C ** 2 * (2 * u[j, -2] - 2 * u[j, -1]) \
                           - (u[j - 1, -1] - 2 * u[j, -1])
        else:
            u[j + 1, -1] = u[j, -1] - C * (u[j, -1] - u[j, -2])
    
    return u

Die einzelnen Zeitschritt können als Animation visualisiert werden.

import matplotlib.animation as anim
from IPython.display import HTML

def show_anim(u, a, T):

    nt = u.shape[0] - 1
    nx = u.shape[1] - 1
    x = np.linspace(0, a, nx + 1)

    # Bilsrate auf 25 fps begrenzen
    fps_real = nt / T
    fps_show = 25
    skip_frames = int(np.floor(fps_real / fps_show))
    n_frames = (nt + 1) // (skip_frames + 1)
    
    # Figure/Axes erstellen
    fig, ax = plt.subplots()
    ax.set_xlim(0, a)
    ax.set_ylim(-1.1, 1.1)
    ax.set_xlabel('x')
    ax.set_ylabel('u(x,t)')
    ax.grid()
    
    # Linienobjekt erstellen (wird während Animation verändert)
    line = ax.plot(0, 0, '-b')[0]
    
    # Update-Funktion
    def update(frame):
        line.set_data(x, u[(skip_frames + 1) * frame, :])
    
    # Animation erstellen und anzeigen
    fa = anim.FuncAnimation(fig, update, frames=n_frames, interval=1000/fps_show)
    display(HTML(fa.to_jshtml()))
    plt.close()  # verhindert automatischen Aufruf von plt.show() durch Jupyter

8.6Experimente¶

Verwenden folgende Parameter:

a = 2
T = 4
c = 1
nx = 100
nt = int(np.ceil(c * nx * T / a))  # Courant = 1 oder wenigstens <= 1
print(f'{nt} Zeitschritte')

200 Zeitschritte

8.6.1Dirichlet-Randbedingungen¶

f = lambda x: np.sin(1/a * np.pi * x)
u = wave1d(a, T, c, f, nx, nt, 'd', 'd')
show_anim(u, a, T)

Courant-Zahl: 1.00

8.6.2Neumann-Randbedingungen¶

Möchten wir Neumann-Randbedingungen verwenden, so sollte die Anfangsbedingung diese auch erfüllen:

f = lambda x: np.cos(2/a * np.pi * x)
u = wave1d(a, T, c, f, nx, nt, 'n', 'n')
show_anim(u, a, T)

Courant-Zahl: 1.00

8.6.3Links Neumann-, rechts freie Randbedingung¶

f = lambda x: np.exp(-(x - 0.5 * a) ** 2 / 0.05)
u = wave1d(a, T, c, f, nx, nt, 'n', 'f')
show_anim(u, a, T)

Courant-Zahl: 1.00

8.6.4Links Dirichlet-, rechts Neumann-Randbedingung¶

Verlängern den Beobachtungszeitraum um eine volle Periode des Schwingungsvorgangs zu sehen:

T = 8
nt = int(np.ceil(c * nx * T / a))  # Courant = 1 oder wenigstens <= 1
print(f'{nt} Zeitschritte')

u = wave1d(a, T, c, f, nx, nt, 'd', 'n')
show_anim(u, a, T)

400 Zeitschritte
Courant-Zahl: 1.00

8.6.5Courant-Zahl zu groß¶

T = 4
nt = int(np.ceil(0.99 * c * nx * T / a))  # Courant > 1
print(f'{nt} Zeitschritte')

f = lambda x: np.sin(1/a * np.pi * x)
u = wave1d(a, T, c, f, nx, nt, 'd', 'd')
show_anim(u, a, T)

198 Zeitschritte
Courant-Zahl: 1.01

Obwohl die Courant-Zahl hier nur ganz wenig über der Eins liegt, wird die Berechnung schlecht konditioniert, sodass die unmeidbaren Rundungsfehler innerhalb weniger Rechenschritte extrem verstärkt werden und die numerische Lösung nichts mehr mit der exakten Lösungen des diskretisierten Problems zu tun hat.