Die 2D-Wellengleichung

Simulation der 2D Wellengleichung mit der Methode der Finiten Differenzen in Python.

Die Grundlagen der Methode der Finite Differenzen

Eine Seite Python code zum Lösen der Wellengleichung mit absorbierenden Randbedingungen. (Klicken zum vergrößern)
Simulation der zweidimensionalen Wellengleichung in Python.

Heute lernen wir, wie man Wellenausbreitung in einem zweidimensionalen Raum mit der Methode der Finiten Differenzen simulieren kann. Mathematisch gesehen ist die Wellengleichung eine hyperbolische partielle Differentialgleichung zweiter Ordnung. Sie sieht so aus:

\begin{equation} \frac{\partial^2 u}{\partial t^2}(x,y,t) - \alpha^2 \left( \frac{\partial^2 u}{\partial x^2}(x,y,t) + \frac{\partial^2 u}{\partial y^2}(x,y,t) \right) = 0 \label{eq:pde1} \end{equation}

Die Gleichung beschreibt, wie sich eine Störung der Größe u zeitlich und räumlich in einem Medium bewegt, in dem sich Wellen mit der Geschwindigkeit α ausbreiten. Sie kombiniert dafür die zweiten Ableitungen des Feldes u nach den Ortsdimensionen mit der zweiten Ableitung nach der Zeit. Was genau u ist, spielt keine Rolle. Es könnte beispielsweise für eine Wassertiefe oder den Luftdruck stehen. Im Wasser gilt diese Gleichung strenggenommen nur für Flachwasserwellen, bei denen die Wellenlänge sehr viel größer als die Wassertiefe ist (z.B. Tsunamis).

Ich möchte hier nicht über partielle Ableitungen referieren, sondern lediglich das absolute Minimum dessen präsentieren, was zum Verständnis der Lösung dieser Differentialgleichung nötig ist. Gleichung \(\ref{eq:pde1}\) sieht zu kompliziert aus, wir schreiben sie deshalb um. Glücklicherweise mögen Physiker Operatoren. Ein Operator ist ein Symbol, das für eine, oder eine Serie von Operationen steht. Für unsere Formel benötigen wir den Laplace-Operator. Dieser steht für die Summe der zweiten partiellen Ableitungen nach allen Raumdimensionen einer Feldgröße. Weiterhin verwenden wir eine Konvention, die besagt, dass eine partielle Ableitung einfach durch Tiefstellen der Größe symbolisiert wird, nach der abgeleitet wird. Die Wellengleichung sieht dann so aus:

Wirkungsweise der diskretisierten Version des Differentialoperators für die Lösung der Wellengleichung Diskretisierung der Wellengleichung. Zellen, die für die Zeitableitung benötigt werden sind Grün eingefärbt, rote Zellen werden für den Laplaceoperator benötigt.
Der diskrete Laplace-Operator am Gitterknoten i,j berechnet sich aus einer gewichteten Summe der Werte der Nachbarzellen. Es gibt verschiedene diskrete Versionen dieses Operators, die sich durch die Anzahl der Nachbarzellen und ihre Gewichtung unterscheiden. \begin{equation} u_{tt} - \alpha^2 \Delta u = 0 \label{eq:wave2d} \end{equation}

Diese Gleichung ist lediglich eine andere Schreibweise von Gleichung \(\ref{eq:pde1}\). Hier wird aber deutlicher, dass wir für das Lösen der Differentialgleichung nur zwei Dinge benötigen: Die zweite Ableitung der Feldgröße nach der Zeit und das Ergebnis der Anwendung der Laplace-Operators.

Da wir numerisch vorgehen, müssen wir das Problem diskretisieren. Unsere Simulationsdomäne wird ein Gitter sein. Wir benötigen zwei Ortsdimensionen und eine Zeitdimension. Bei der Zeitdimension interessieren uns allerdings nur der aktuelle und die beiden vorangegangenen Zeitschritte. Wir benötigen drei Felder, die den Zustand der Feldgröße an den drei Zeitpunkten speichern: \(u^{(0)}, u^{(1)}, u^{(2)}\). Der Zeitindex 0 steht für die Zukunft, der Index 1 für den aktuellen Zustand und der Index 2 für den vorangegangenen Zustand des Feldes.

Die Knotenpunkte des Gitters haben in der Ortsdimension den Abstand \(h\) voneinander, in der Zeitdimension sei der Abstand \(t\). Jeder der Gitterknotenpunkte repräsentiert einen Näherungswert des Feldes \(u(x,y)\) an der Stelle \(u(i*h, j*h)\).

Beginnen wir mit der Bestimmung der Zeitableitung \(u_{tt}\) an einem Gitterknoten i,j zum Zeitindex 1. Für die numerische Berechnung dieser Ableitung verwenden wir den Differenzenquotient 2. Ordnung. Dieser benötigt neben dem aktuellen Wert der Feldgröße \(u^{(1)}_{(i,j)}\) einen Wert aus der Zukunft \(u^{(0)}_{(i,j)}\) und einen aus der Vergangenheit \(u^{(2)}_{(i,j)}\).

\begin{equation} u^{(1)}_{{tt}_{i,j}} = \frac{u^{(0)}_{i,j} - 2 u^{(1)}_{i,j} + u^{(2)}_{i,j}}{t^2} \label{eq:timediff} \end{equation}

Den Wert aus der Zukunft kennen wir natürlich noch nicht. Um dessen Berechnung wird es im folgenden aber gehen. Neben der zweiten Zeitableitung benötigen wir noch den Laplace-Operator für die 2. Ortsableitung in x- und y-Richtung. Die hier verwendete Version des zweidimensionalen, diskreten Laplace-Operators \(\mathbf{D}^2_{i,j}\) verwendet für die Berechnung die Werte in den 4 direkten Nachbarzellen in x und y-Richtung mit einfacher Wichtung und subtrahiert den Wert des Feldes im Knotenpunktes i,j mit vierfacher Wichtung:

\begin{equation} \quad \mathbf{D}^2_{i,j}=\begin{bmatrix}0 & 1 & 0\\1 & -4 & 1\\0 & 1 & 0\end{bmatrix} \end{equation}

Woraus sich für den Gitterknoten i,j folgende Gleichung für die Bestimmung des Laplaceoperators an dieser Stelle ergibt:

\begin{equation} \Delta u^{(1)}_{i,j} = \frac{u^{(1)}_{i-1,j}+u^{(1)}_{i+1,j}+u^{(1)}_{i,j-1}+u^{(1)}_{i,j+1} - 4 u^{(1)}_{i,j}}{h^2} \label{eq:laplace} \end{equation}

Wenn wir jetzt die Gleichungen \(\ref{eq:timediff}\) und \(\ref{eq:laplace}\) in Gleichung \(\ref{eq:wave2d}\) einsetzen, dann erhalten wir eine diskrete Form der Wellengleichung:

\begin{equation} \frac{u^{(0)}_{i,j} - 2 u^{(1)}_{i,j} + u^{(2)}_{i,j}}{t^2} - \alpha^2 * \left( \frac{u^{(1)}_{i-1,j}+u^{(1)}_{i+1,j}+u^{(1)}_{i,j-1}+u^{(1)}_{i,j+1} - 4 u^{(1)}_{i,j}}{h^2} \right)= 0 \end{equation}

Das ganze stellen wir jetzt nach \(u^{(0)}_{(i,j)}\) um, denn uns interessiert ja der Feldwert in der Zukunft:

\begin{equation} u^{(0)}_{i,j} = \left( \frac{t \alpha}{h}\right)^2 \left( u^{(1)}_{i-1,j}+u^{(1)}_{i+1,j}+u^{(1)}_{i,j-1}+u^{(1)}_{i,j+1} - 4 u^{(1)}_{i,j} \right) + 2 u^{(1)}_{i,j} - u^{(2)}_{i,j} \label{eq:findiff} \end{equation}

Mit dieser Gleichung können wir jeden Gitterknoten in \(u^{(0)}\) aus den beiden vorangegangenen Zeitschritten \(u^{(1)}\) und \(u^{(2)}\) berechnen. Dafür iterieren wir über alle Knoten des Feldes. Die Werte in den Feldern \(u^{(1)}\) und \(u^{(2)}\) werden von uns als Anfangsbedingungen vorgegeben. Bei diesem Verfahren handelt es sich um ein explizites Euler-Verfahren zweiter Ordnung.

Python Quellcode

Eine einfache Lösung der Wellengleichung mittels der Methode der Finiten Differenzen lässt sich in nur wenigen Zeilen Python-Quellcode implementieren. Der Quellcode für ein Beispielimplementierung unter Verwendung von Orts- und Zeitableitungen zweiter Fehlerordnung und mit statischen Randbedingungen befindet sich in der Datei waves_2d.py des Github Archivs dieser Artikelserie.

waves_2d.py 

import pygame
import numpy as np
import random

h = 1        # spatial step width
k = 1        # time step width
dimx = 300   # width of the simulation domain
dimy = 300   # height of the simulation domain
cellsize = 2 # display size of a cell in pixel

def init_simulation():
    u = np.zeros((3, dimx, dimy))         # The three dimensional simulation grid 
    c = 0.5                               # The "original" wave propagation speed
    alpha = np.zeros((dimx, dimy))        # wave propagation velocities of the entire simulation domain
    alpha[0:dimx,0:dimy] = ((c*k) / h)**2 # will be set to a constant value of tau
    return u, alpha

def update(u, alpha):
    u[2] = u[1]
    u[1] = u[0]

    # This switch is for educational purposes. The fist implementation is approx 50 times slower in python!
    use_terribly_slow_implementation = False
    if use_terribly_slow_implementation:
        # Version 1: Easy to understand but terribly slow!
        for c in range(1, dimx-1):
            for r in range(1, dimy-1):
                u[0, c, r]  = alpha[c,r] * (u[1, c-1, r] + u[1, c+1, r] + u[1, c, r-1] + u[1, c, r+1] - 4*u[1, c, r])
                u[0, c, r] += 2 * u[1, c, r] - u[2, c, r]
    else:
        # Version 2: Much faster by eliminating loops
        u[0, 1:dimx-1, 1:dimy-1]  = alpha[1:dimx-1, 1:dimy-1] * (u[1, 0:dimx-2, 1:dimy-1] + 
                                            u[1, 2:dimx,   1:dimy-1] + 
                                            u[1, 1:dimx-1, 0:dimy-2] + 
                                            u[1, 1:dimx-1, 2:dimy] - 4*u[1, 1:dimx-1, 1:dimy-1]) \
                                        + 2 * u[1, 1:dimx-1, 1:dimy-1] - u[2, 1:dimx-1, 1:dimy-1]

    # Not part of the wave equation but I need to remove energy from the system. 
    # The boundary conditions are closed. Energy cannot leave and the simulation keeps adding energy.
    u[0, 1:dimx-1, 1:dimy-1] *= 0.995

def place_raindrops(u):
    if (random.random()<0.02):
        x = random.randrange(5, dimx-5)
        y = random.randrange(5, dimy-5)
        u[0, x-2:x+2, y-2:y+2] = 120

def main():
    pygame.init()
    display = pygame.display.set_mode((dimx*cellsize, dimy*cellsize))
    pygame.display.set_caption("Solving the 2d Wave Equation")

    u, alpha = init_simulation()
    pixeldata = np.zeros((dimx, dimy, 3), dtype=np.uint8 )

    while True:
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                pygame.quit()
                return

        place_raindrops(u)
        update(u, alpha)

        pixeldata[1:dimx, 1:dimy, 0] = np.clip(u[0, 1:dimx, 1:dimy] + 128, 0, 255)
        pixeldata[1:dimx, 1:dimy, 1] = np.clip(u[1, 1:dimx, 1:dimy] + 128, 0, 255)
        pixeldata[1:dimx, 1:dimy, 2] = np.clip(u[2, 1:dimx, 1:dimy] + 128, 0, 255)

        surf = pygame.surfarray.make_surface(pixeldata)
        display.blit(pygame.transform.scale(surf, (dimx * cellsize, dimy * cellsize)), (0, 0))
        pygame.display.update()

if __name__ == "__main__":
    main()

Dieses Python-Script verwendet die Bibliotheken random, numpy und pygame. Die Bibliotheken numpy und pygame müssen über Pythons Paketmanager PIP installiert worden sein, bevor das Programm ausgeführt werden kann. Nach den import-Anweisungen folgt die Definition einiger globaler Variablen. Es werden Zeit- sowie Ortsschrittweiten definiert und in der Funktion init_simulation werden zwei Felder angelegt: Ein dreidimensionales Feld für die Feldgröße u (2 Ortsdimensionen und eine Zeitdimension) und ein zweidimensionales Feld, in dem die Wellenausbreitungsgeschwindigkeit innerhalb der Simulationsdomäne gespeichert wird.

In dieser Version der Simulation ist die Wellenausbreitungsgeschwindigkeit im gesamten Simulationsbereich konstant. Die Anfangsbedingungen des Feldes u sind im gesamten SImulationsgebiet 0. Wenn das Programm startet, wird also erst einmal nichts passieren, denn das initiale Wellenfeld ist überall konstant Null. Erst nach einiger Zeit werden in der Funktion place_raindrops zufällig "Regentropfen" platziert. 

import pygame
import numpy as np
import random

h = 1        # spatial step width
k = 1        # time step width
dimx = 300   # width of the simulation domain
dimy = 300   # height of the simulation domain
cellsize = 2 # display size of a cell in pixel

def init_simulation():
    u = np.zeros((3, dimx, dimy))         # The three dimensional simulation grid 
    c = 0.5                               # The "original" wave propagation speed
    alpha = np.zeros((dimx, dimy))        # wave propagation velocities of the entire simulation domain
    alpha[0:dimx,0:dimy] = ((c*k) / h)**2 # will be set to a constant value of tau
    return u, alpha

Das Herzstück der Simulation befindet sich in der Funktion update. Dort wird die Wellengleichung gelöst. Am Beginn der Funktion werden die Felder der letzten beiden Zeitschritte umkopiert bzw. "weitergerückt". Das Feld mit dem Index 0 wird die neuen Ergebnisse aufnehmen, dessen Zustand wird in das Feld mit dem Zeitindex 1 kopiert und die Werte von dort werden zum Zeitindex 2 verschoben. Ältere Werte werden verworfen. 

def update(u, alpha):
    u[2] = u[1]
    u[1] = u[0]

    # This switch is for educational purposes. The fist implementation is approx 50 times slower in python!
    use_terribly_slow_implementation = False
    if use_terribly_slow_implementation:
        # Version 1: Easy to understand but terribly slow!
        for c in range(1, dimx-1):
            for r in range(1, dimy-1):
                u[0, c, r]  = alpha[c,r] * (u[1, c-1, r] + u[1, c+1, r] + u[1, c, r-1] + u[1, c, r+1] - 4*u[1, c, r])
                u[0, c, r] += 2 * u[1, c, r] - u[2, c, r]
    else:
        # Version 2: Much faster by eliminating loops
        u[0, 1:dimx-1, 1:dimy-1]  = alpha[1:dimx-1, 1:dimy-1] * (u[1, 0:dimx-2, 1:dimy-1] + 
                                            u[1, 2:dimx,   1:dimy-1] + 
                                            u[1, 1:dimx-1, 0:dimy-2] + 
                                            u[1, 1:dimx-1, 2:dimy] - 4*u[1, 1:dimx-1, 1:dimy-1]) \
                                        + 2 * u[1, 1:dimx-1, 1:dimy-1] - u[2, 1:dimx-1, 1:dimy-1]

    # Not part of the wave equation but I need to remove energy from the system. 
    # The boundary conditions are closed. Energy cannot leave and the simulation keeps adding energy.
    u[0, 1:dimx-1, 1:dimy-1] *= 0.995

Die Berechnung findet in nur 4 Zeilen Quellcode statt. Diese Zeilen sind die direkte Umsetzung der Gleichung \(\ref{eq:findiff}\). Sie funktionieren, machen die Simulation aber unerträglich langsam, denn in Python ist das abarbeiten von for-Schleifen sehr langsam. Wer dieses Problem in einer anderen Programmiersprache nachprogrammieren möchte, kann gerne mit Schleifen programmieren. Für Python hat dieser Code nur akademischen Wert. 

def update(u, alpha):
    for c in range(1, dimx-1):
        for r in range(1, dimy-1):
            u[0, c, r]  = alpha[c,r] * (u[1, c-1, r] + u[1, c+1, r] + u[1, c, r-1] + u[1, c, r+1] - 4*u[1, c, r])
            u[0, c, r] += 2 * u[1, c, r] - u[2, c, r]

Für eine benutzbare Python-Version müssen wir den Code umschreiben. Wir verwenden dafür ein Python-Feature, die sogenannte Vektorisierung. Der folgende Code ist äquivalent zum obigen Quellcode, kommt aber gänzlich ohne for-Schleifen aus. 

	u[0, 1:dimx-1, 1:dimy-1] = alpha[1:dimx-1, 1:dimy-1] * ( \
		u[1, 0:dimx-2, 1:dimy-1] + \
		u[1, 2:dimx,   1:dimy-1] + \
		u[1, 1:dimx-1, 0:dimy-2] + \
		u[1, 1:dimx-1, 2:dimy] - 4*u[1, 1:dimx-1, 1:dimy-1]) \
		+ 2 * u[1, 1:dimx-1, 1:dimy-1] - u[2, 1:dimx-1, 1:dimy-1]

Vektorisierte Operationen in Python arbeiten mit kompletten Bereichen von Eingabevektoren. Somit ist es möglich die explizite Iteration über den Simulationsbereich zu vermeiden. Man kann mit einer einzigen Gleichung tausende Werte berechnen. Dort wo in der Gleichung oben einzelne Zellen adressiert wurde, werden nun ganze Indexbereiche auf einmal adressiert.

Mit u[0, 1:dimx-1, 1:dimy-1] addressieren wir das komplette Feld u, außer seinem Randbereich. Das ist äquivalent zur Verwendung von u[0, c, r] in der Doppelschleife oben. Wenn wir u[0,c-1,r] adressieren möchten, müssen wir einfach die "-1" vom Indexbereich der Spalte 1:dimx-1 abziehen und schreiben in der vektorisierten Form u[0, 0:dimx-2, 1:dimy-1]. Das sieht wesentlich komplizierte aus, ist aber um Größenordnungen (geschätzt 50x) schneller in der Ausführung.

Im Rest des Programmes befindet sich die bereits erwähnte Funktion zum zufälligen platzieren von Regentropfen und die Hauptschleife mit dem Code für die Visualisierung. Für die Darstellung der Ergebnisse wird der Wert des Feldes u in den drei Zeitschritten einfach in Graustufenwerte umgerechnet und mittels pygame direkt in das Ausgabefenster geschrieben.

Verwendung genauere Ortsableitungen

Wer nur einen Wassereffekt programmieren möchte, den wird es nicht stören, dass das explizite Euler-Verfahren nicht besonders genau ist. Es ist das mit Abstand am einfachsten zu implementierende Verfahren und daher auch das am häufigsten beschriebene. Der Effekt sieht auf den ersten Blick realistisch aus. Die Methode hat jedoch eine geringe Fehlerordnung O(h2). Das führt insbesondere bei Unstetigkeiten in den Anfangsbedingungen dazu, dass hinter diesen Wellen entstehen, die es eigentlich nicht geben dürfte. In der Realität würde ein, in der Mitte des Gitters platzierter Wellenberg ringförmig auseinander laufen. Man würde einen Ring sehen, sonst nichts. In unserer Simulation folgt dem Ring aber ein ganzer Wellenzug.

Eine einfache Möglichkeit die Rechengenauigkeit zu erhöhen, ist die Verbesserung der, vom Laplace-Operator berechneten, Ortsableitungen. Dies kann durch die Verwendung eines größeren Laplace-Operators erreicht werden. Anstelle einer Version mit 3x3 Einträgen (siehe Gleichung \(\ref{eq:laplace}\)), kann man eine mit 5x5, 7x7 oder auch 9x9 Verwenden. Die hier gezeigten Versionen von diskreten Laplace-Operatoren repräsentieren finite Differenzenverfahren höherer Ordnung für die Ortsableitung [3]:

\(O(h4)\): \begin{equation} \nonumber \mathbf{D}^2_{(i,j)}=\frac{1}{12}\begin{bmatrix}0 & 0 & -1 & 0 & 0\\0 & 0 & 16 & 0 & 0\\-1 & 16 & -60 & 16 & -1\\0 & 0 & 16 & 0 & 0\\0 & 0 & -1 & 0 & 0\end{bmatrix} \end{equation} \(O(h6)\): \begin{equation} \nonumber \mathbf{D}^2_{(i,j)}=\frac{1}{180}\begin{bmatrix}0 & 0 & 0 & 2 & 0 & 0 & 0\\0 & 0 & 0 & -27 & 0 & 0 & 0\\0 & 0 & 0 & 270 & 0 & 0 & 0\\2 & -27 & 270 & -980 & 270 & -27 & 2\\ 0 & 0 & 0 & 270 & 0 & 0 & 0\\ 0 & 0 & 0 & -27 & 0 & 0 & 0\\0 & 0 & 0 & 2 & 0 & 0 & 0\end{bmatrix} \end{equation}

Eine Version der Differenzengleichung \(\ref{eq:findiff}\) mit einem Ortsableitungsterm der Fehlerordnung 4 sieht so aus:

\begin{equation} u^{(0)}_{i,j} = \left( \frac{t \alpha}{h}\right)^2 \left(-u^{(1)}_{i,j-2}+16u^{(1)}_{i,j-1}-u^{(1)}_{i-2,j}+16u^{(1)}_{i-1,j}-60u^{(1)}_{i,j}+16u^{(1)}_{i+1,j}-u^{(1)}_{i+2,j}+16u^{(1)}_{i,j+1}-u^{(1)}_{i,j+2} \right) + 2 u^{(1)}_{i,j} - u^{(2)}_{i,j} \label{eq:findiff6} \end{equation}

Die Lösung der zweidimensionalen Wellengleichung unter Verwendung von Ortsableitungen bis zur Fehlerordnung 8 wird in der Datei waves_2d_improved.py des begleitenden GitHub Archivs demonstriert.

waves_2d_improved.py

Die Verwendung von Ortsableitungen höherer Fehlerordnung führt zu einer deutlichen Reduktion der numerischen Artefakte, kann diese aber nicht komplett unterdrücken. Das Video unten zeigt den Vergleich eines expliziten Euler Verfahrens mit Ortsableitung der Fehlerordnung 2 mit einem Verfahren dessen Ortsableitung die Fehlerordnung 8 hat (9x9 Laplace Operator). Man erkauft sich Genauigkeit mit höherem Rechenaufwand, reduziert aber auch die Stabilität des Verfahrens, weshalb die maximal mögliche Zeitschrittweite ein wenig geringer ist.

Numerische Lösung der zweidimensionalen Wellengeleichung unter Verwendung eines expliziten Euler-Verfahrens mit Ortsableitungen verschiedener Fehlerordnungen.

Absorbierende Randbedingungen

Die Simulation, wie bisher beschrieben, hat ein entscheidendes Problem. Alle Wellen werden am Rand reflektiert. Die Energie kann den Simulationsbereich nicht verlassen. Das wird, wenn immer neue Wellenberge platziert werden, nach kurzer Zeit dazu führen, dass man nichts mehr erkennt. Die Wellenenergie kann das Simulationsgebiet nicht verlassen.

Ein 3x3 Laplace-Operator benötigt für die Berechnung der Orstableitungen immer einen Wert in allen Richtungen neben dem zu berechnenden Gitterpunkt. Aus diesem Grund kann er auf den Gitterpunkten am Rand nicht angewendet werden. Wir haben deshalb bislang die Randwerte fest mit 0 vorgegeben. Das sind sogenannte Dirichlet-Randbedingungen. Deshalb kann keine Energie nach außen transferiert werden, weshalb auftreffenden Wellen zurückgeworfen werden.

Um einen offenen Rand zu simulieren, müssen wir auf diesem spezielle Differentialgleichungen lösen. Solche Bedingungen werden auch "abstrahlende" (engl. radiating), "absorbierende" (engl. absorbing oder ABC) oder "Mur Randbedingungen" genannt. [4], [5] Bei abstrahlenden Randbedingungen gelten an den Enden des Simulationsbereiches in X-Richtung folgende Differentialgleichungen:

\begin{align} \begin{alignedat}{3} \frac{\partial u}{\partial x} &&\biggr\rvert_{x=0} &= -&&\frac{1}{c} \frac{\partial u}{\partial t} \label{eq:mur1} \nonumber\\ \frac{\partial u}{\partial x} &&\biggr\rvert_{x=N} &= &&\frac{1}{c} \frac{\partial u}{\partial t} \label{eq:mur2} \end{alignedat} \end{align}

In Y-Richtung gelten analoge Gleichungen. Daraus ergeben sich für die Ränder x=0, x=N, y=0 und y=N folgende Differenzengleichungen für die Randgitterpunkte:

\begin{align} \begin{alignedat}{3} u^{(0)}_{0,j} &= u^{(1)}_{1,j} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{1,j} &&- u^{(1)}_{0,j}) \nonumber\\ u^{(0)}_{N,j} &= u^{(1)}_{N-1,j} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{N-1,j} &&- u^{(1)}_{N,j}) \nonumber\\ u^{(0)}_{i,0} &= u^{(1)}_{i,1} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,1} &&- u^{(1)}_{i,0}) \nonumber\\ u^{(0)}_{i,N} &= u^{(1)}_{i,N-1} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,N-1} &&- u^{(1)}_{i,N}) \end{alignedat} \end{align}

wobei N für den Index der letzten Gitterknotens steht und \(\kappa\) sich aus wie folgt aus Wellenausbreitungsgeschwingigkeit, Ortsschrittweite und Zeitschrittweite berechnet:

\begin{equation} \kappa = \alpha \frac{t}{h} \end{equation}

Diese Differenzengleichungen werden auf alle Gitterpunkte am Rand angewendet. Bei Verwendung eines größeren diskreten Laplace-Operators (5x5, 7x7 oder 9x9) muss der Randbereich breiter gewählt werden, damit alle Gitterknotenpunkte die der Laplace-Operator nicht abdecken kann, mit einbezogen werden. Die Implementierung von abstrahlenden Randbedingungen nach Mur wird in der Datei waves_2d_improved.py des begleitenden GitHub Archivs demonstriert.

waves_2d_improved.py

Literaturnachweis

  1. J. Douglas Faires, Richard L. Burden: "Numerische Methoden: Näherungsverfahren und ihre praktische Anwendung." Spektrum Akademischer Verlag; ISBN 3-86025-332-8; 1995
  2. Hans Stefani, Gerhard Kluge: "Theoretische Mechanik" Spektrum Akademischer Verlag; ISBN 3-86025-284-4; 1995
  3. Bengt Fornberg: "Generation of finite difference formulas on arbitrarily spaced grids." In: Mathematics of Computation. Band 51, Nr. 184, 1988, Seite 702
  4. Gerrit Mur: "Absorbing Boundary Conditions for the Finite-Difference Approximation of the Time-Domain Electromagnetic-Field Equations" IEEE Transactions on Electromagnetic Compatibility: Vol EMC-23, November 1981; Page 377-382
  5. unbekannter Autor: "Finite Difference Methods" online; unbekannter Herausgeber