Die 3D-Wellengleichung

Simulation der dreidimensionalen Wellengleichung mit der Methode der Finiten Differenzen in Python.

Die Grundlagen der Methode der Finite Differenzen

Simulation stehender Wellen durch numerisches Lösen der dreidimensionalen Wellengleichung in Python.

Interferenz und Beugung einer Wellenfront an zwei runden Durchbrüchen.

Nachdem ich in einem früheren Artikel die zweidimensionale Wellengleichung besprochen habe, kommt nun die logische Fortsetzung. Wir bauen einen dreidimensionalen Wellentank! Das Prinzip ist das Gleiche, wie bei der zweidimensionalen Wellengleichung, nur jetzt in drei Raumdimensionen.

Das ist die dreidimensionale Wellengleichung:

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

Wer den Artikel zur 2-D Wellengleichung gelesen hat, weiß bereits, dass der eigentliche Zweck dieser Formulierung ist, Studenten zu erschrecken und das funktioniert prima. Der Ein oder Andere angehende Naturwissenschaftler wird vom Versuch die Gleichung analytisch zu lösen posttraumatisches Belastungsstörungen davongetragen haben, denn dafür benötigt man die diskrete Fouriertransformation. Wer möchte, kann in diesem Link eine Beschreibung des Lösungsweges für nur eine Raumdimension nachlesen.

Jean Baptiste Joseph Fourier war zweifelsohne ein Genie, denn rein mathematisch betrachtet ist sein Lösungsweg eine elegante und geniale Methode aber darum geht es hier nicht, denn dieser Artikel richtet sich an "Nicht-Genies" mit einem Computer.

Gleichung \(\ref{eq:pde1}\) 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 den Luftdruck stehen.

Wenn Physiker untereinander reden, verwenden sie häufig eine andere Formulierung der Gleichung:

\begin{equation} u_{tt} - \alpha^2 \Delta u = 0 \label{eq:wave2d} \end{equation}

Hier fällt zunächst auf, das dies exakt die selbe Gleichung, wie im Artikel über die zweidimensionale Wellengleichung ist. Das liegt daran, dass das "Dreieck" in der Gleichung für den Laplaceoperator steht, dessen Dreidimensionalität hier einfach "mitgedacht" wird. Es spielt keine große Rolle ob das Problem zweidimensional oder dreidimensional ist. Im dreidimensionalen Fall muss man einfach ein paar Terme mehr für die Z-Richtung auf das Problem werfen, damit es weggeht. Für die Lösung benötigen wir also die zweite Ableitung der Feldgröße nach der Zeit und das Ergebnis der Anwendung des, jetzt dreidimensionalen, Laplace-Operators.

Da wir numerisch vorgehen, müssen wir das Problem diskretisieren. Unsere Simulationsdomäne wird ein dreidimensionales Gitter sein. Wir benötigen drei 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,z)\) an der Stelle \(u(i*h, j*h, k*h)\).

Wir beginnen 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,k)}\) einen Wert aus der Zukunft \(u^{(0)}_{(i,j,k)}\) und einen aus der Vergangenheit \(u^{(2)}_{(i,j,k)}\).

\begin{equation} u^{(1)}_{{tt}_{i,j,k}} = \frac{u^{(0)}_{i,j,k} - 2 u^{(1)}_{i,j,k} + u^{(2)}_{i,j,k}}{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 gehen. Neben der zweiten Zeitableitung benötigen wir noch den Laplace-Operator für die 2. Ortsableitung in x-, y- und z-Richtung. Die hier verwendete Version des dreidimensionalen, diskreten Laplace-Operators \(\mathbf{D}^3_{i,j,k}\) verwendet für die Berechnung die Werte in den 6 direkten Nachbarzellen in x, y und z-Richtung mit einfacher Wichtung und subtrahiert den Wert des Feldes im Knotenpunktes i,j,k mit sechsfacher Wichtung:

\begin{equation} \quad \mathbf{D}^3_{i,j,k}= \left\{ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \begin{bmatrix} 0 & 1 & 0 \\ 1 & -6 & 1 \\ 0 & 1 & 0 \end{bmatrix}, \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} \right\} \end{equation}

Das sieht kompliziert aus, aber es ist einfach die Erweiterung des zweidimensionalen Laplaceoperators um eine Dimension, denn in die Berechnung der zweiten Ableitungen sollen nun auch die beiden Nachbarzellen in Z-Richtung eingehen. Was hier in Matrixform niedergeschrieben wurde, ist der Faltungskernel des diskreten dreidimensionalen Laplaceoperators. Was es bedeutet ist, dass sich für den Gitterknoten i,j,k folgende Gleichung für dessen Bestimmung an dieser Stelle ergibt:

\begin{equation} \Delta u^{(1)}_{i,j,k} = \frac{u^{(1)}_{i-1,j,k}+u^{(1)}_{i+1,j,k}+u^{(1)}_{i,j-1,k}+u^{(1)}_{i,j+1,k} +u^{(1)}_{i,j,k-1}+u^{(1)}_{i,j,k+1} - 6 u^{(1)}_{i,j,k}}{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 dreidimensionalen Wellengleichung:

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

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

\begin{equation} u^{(0)}_{i,j, k} = \left( \frac{t \alpha}{h}\right)^2 \left( u^{(1)}_{i-1,j,k}+u^{(1)}_{i+1,j,k}+u^{(1)}_{i,j-1,k}+u^{(1)}_{i,j+1,k} +u^{(1)}_{i,j,k-1}+u^{(1)}_{i,j,k+1} - 6 u^{(1)}_{i,j,k} \right) + 2 u^{(1)}_{i,j,k} - u^{(2)}_{i,j,k} \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. Zum Beginn der Berechnung werden die Werte in den Feldern \(u^{(1)}\) und \(u^{(2)}\) als Anfangsbedingungen vorgegeben. Bei diesem Verfahren handelt es sich um ein explizites Euler-Verfahren zweiter Ordnung.

Verwendung genauere Ortsableitungen

Das explizite Euler-Verfahren 2. Ordnung ist nicht besonders gut und wird in der Regel nicht für ernsthafte Berechnungen verwendet. Eine einfache Möglichkeit die Rechengenauigkeit zu erhöhen, ist die Verbesserung der, vom Laplace-Operator berechneten, Ortsableitungen. Dies kann durch die Verwendung eines Laplace-Operators mit größerem Faltungskernel erreicht werden. Anstelle einer Version mit 3x3x3 Einträgen (siehe Gleichung \(\ref{eq:laplace}\)), kann man eine mit 5x5x5, 7x7x7 oder auch 9x9x9 Verwenden. Dieses Versionen der Faltungskernel des diskreten Laplace-Operators können leicht aus den bereits gezeigten zweidimensionalen Versionen höherer Ordnung abgeleitet werden. Exemplarisch sei das hier am Laplaceoperator 4. Ordnung demonstriert:

\(O(h4)\): \begin{equation} \quad \mathbf{D}^3_{i,j,k}= \frac{1}{12} \left\{ \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & 16 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & -1 & 0 & 0\\0 & 0 & 16 & 0 & 0\\-1 & 16 & -90 & 16 & -1\\0 & 0 & 16 & 0 & 0\\0 & 0 & -1 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & 16 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix}, \begin{bmatrix}0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0\\0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0\end{bmatrix} \right\} \end{equation}

Die Erweiterung auf drei Dimensionen besteht darin, eine Z-Dimension mit 4 weiteren Matrizen hinzuzufügen und im Zentrum der Matrizen die gleichen Werte einzutragen, wie in x- und y-Richtung der mittleren Matrix. Der Wert im Zentrum der 5x5x5 Matrix muss von -60 auf -90 gesenkt werden, weil die Summe der Elemente in den, für die Z-Richtung angefügten, Matrizen gleich 30 ist. Eine Version der Differenzengleichung \(\ref{eq:findiff}\) mit einem Ortsableitungsterm der Fehlerordnung 4 sieht so aus:

\begin{equation} \begin{aligned} u^{(0)}_{i,j,k} = \left( \frac{t \alpha}{h}\right)^2 \Bigl( & -u^{(1)}_{i-2,j,k}+16u^{(1)}_{i-1,j,k}+16u^{(1)}_{i+1,j,k}-u^{(1)}_{i+2,j,k} \\ & -u^{(1)}_{i,j-2,k}+16u^{(1)}_{i,j-1,k}+16u^{(1)}_{i,j+1,k}-u^{(1)}_{i,j+2,k} \\ & -u^{(1)}_{i,j,k-2}+16u^{(1)}_{i,j,k-1}+16u^{(1)}_{i,j,k+1}-u^{(1)}_{i,j,k+2} \\ & -90u^{(1)}_{i,j,k} \Bigr) + 2 u^{(1)}_{i,j,k} - u^{(2)}_{i,j,k} \label{eq:findiff16} \end{aligned} \end{equation}

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

waves_3d_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. 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. Ein generelles Problem ist, dass nur die Genauigkeit der partiellen Differenzen der Ortsableitungen verbessert wurde. Die Zeitableitung wird weiterhin nur mit einem Verfahren zweiter Ordnung berechnet.

Das hier verwendete Modell hat eine, in Z-Richtung, linear zunehmender Wellenausbreitungsgeschwindigkeit und eine Zone konstant niedriger Geschwindigkeit in der Mitte (z.B. Hohlraum bei der Seismik). Der Einfluss der Genauigkeit der Ortsableitungen zeigt sich im direkten Vergleich von verschiedenen Implementierungen mit O(h^2) (links), O(h^4) (mitte) und O(h^6) (rechts). Insbesondere in der Simulation links zeigt sich deutlich mehr numerische Dispersion im Ergebnis.

Absorbierende Randbedingungen

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

Ein 3x3x3 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- und Z-Richtung gelten analoge Gleichungen. Daraus ergeben sich für die Ränder x=0, x=N, y=0, y=N, z=0 und z=N folgende Differenzengleichungen für die Randgitterpunkte:

\begin{align} \begin{alignedat}{3} u^{(0)}_{0,j,k} &= u^{(1)}_{1,j,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{1,j,k} &&- u^{(1)}_{0,j,k}) \\ u^{(0)}_{N,j,k} &= u^{(1)}_{N-1,j,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{N-1,j,k} &&- u^{(1)}_{N,j,k}) \\ u^{(0)}_{i,0,k} &= u^{(1)}_{i,1,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,1,k} &&- u^{(1)}_{i,0,k}) \\ u^{(0)}_{i,N,k} &= u^{(1)}_{i,N-1,k} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,N-1,k} &&- u^{(1)}_{i,N,k}) \\ u^{(0)}_{i,j,0} &= u^{(1)}_{i,j,1} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,j,1} &&- u^{(1)}_{i,k,0}) \\ u^{(0)}_{i,j,N} &= u^{(1)}_{i,j,N-1} &&+ \frac{\kappa-1}{\kappa+1} (u^{(0)}_{i,j,N-1} &&- u^{(1)}_{i,j,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 (5x5x5, 7x7x7 oder 9x9x9) 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_3d_improved.py des begleitenden GitHub Archivs demonstriert.

waves_3d_improved.py

Python Quellcode

Die hier gezeigte kurze Version des Quellcodes verwendet Ortsableitungen der Fehlerordnung 2 und absorbierende Randbedingungen 1. Ordnung nach Mur. Er basiert auf dem der zweidimensionalen Verison. Für die Darstellung wird die mayavi Bibliothek verwendet.

waves_3d.py

import numpy as np
import math

from mayavi import mlab
from tvtk.util.ctf import *

hs = ts = 1        # time and space step width
dimx = dimy = 100  # dimensions of the simulation domain 
dimz = 60           

u = np.zeros((3, dimx, dimy, dimz))   # three state arrays of the field variable

# Initial conditions
vel = np.full((dimx, dimy, dimz), 0.3)
vel[0:dimx, 0:dimy, 4:6] = 0.0

for i in range(dimx):
    for j in range(dimy):
        for k in range(4, 6):
            distance1 = np.linalg.norm(np.array([i, j, k]) - np.array([40, 50, 5]))
            distance2 = np.linalg.norm(np.array([i, j, k]) - np.array([60, 50, 5]))
            if distance1 <= 5 or distance2 <= 5:
                vel[i, j, k] = 0.3

tau = ( (vel*ts) / hs )**2
k = ts * vel / hs  

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

    # the 3D Wave equation
    u[0, 1:dimx-1, 1:dimy-1, 1:dimz-1] = tau[1:dimx-1, 1:dimy-1, 1:dimz-1] \
        * (   u[1, 0:dimx-2, 1:dimy-1, 1:dimz-1] +     u[1, 1:dimx-1, 0:dimy-2, 1:dimz-1]
            + u[1, 1:dimx-1, 1:dimy-1, 0:dimz-2] - 6 * u[1, 1:dimx-1, 1:dimy-1, 1:dimz-1]
            + u[1, 2:dimx  , 1:dimy-1, 1:dimz-1] +     u[1, 1:dimx-1, 2:dimy,   1:dimz-1]
            + u[1, 1:dimx-1, 1:dimy-1, 2:dimz] ) + 2 * u[1, 1:dimx-1, 1:dimy-1, 1:dimz-1] - u[2, 1:dimx-1, 1:dimy-1, 1:dimz-1]

    # Absorbing boundary conditions for all 6 sides of the simulation domain
    c = dimx - 1
    u[0, dimx-2:c, 1:dimy-1, 1:dimz-1] = u[1, dimx-3:c-1, 1:dimy-1, 1:dimz-1] \
        + (k[dimx-2:c, 1:dimy-1, 1:dimz-1]-1)/(k[dimx-2:c, 1:dimy-1, 1:dimz-1]+1) \
        * (u[0, dimx-3:c-1, 1:dimy-1, 1:dimz-1] - u[1, dimx-2:c,1:dimy-1, 1:dimz-1])
    c = 0
    u[0, c:1, 1:dimy-1, 1:dimz-1] = u[1, 1:2, 1:dimy-1, 1:dimz-1] \
        + (k[c:1, 1:dimy-1, 1:dimz-1]-1)/(k[c:1, 1:dimy-1, 1:dimz-1]+1) \
        * (u[0, 1:2,1:dimy-1, 1:dimz-1] - u[1, c:1,1:dimy-1, 1:dimz-1])
    r = dimy-1
    u[0, 1:dimx-1, dimy-2:r, 1:dimz-1] = u[1, 1:dimx-1, dimy-3:r-1, 1:dimz-1] \
        + (k[1:dimx-1, dimy-2:r, 1:dimz-1]-1)/(k[1:dimx-1, dimy-2:r, 1:dimz-1]+1) \
        * (u[0, 1:dimx-1, dimy-3:r-1, 1:dimz-1] - u[1, 1:dimx-1, dimy-2:r, 1:dimz-1])
    r = 0
    u[0, 1:dimx-1, r:1, 1:dimz-1] = u[1, 1:dimx-1, 1:2, 1:dimz-1] \
        + (k[1:dimx-1, r:1, 1:dimz-1]-1)/(k[1:dimx-1, r:1, 1:dimz-1]+1) \
        * (u[0, 1:dimx-1, 1:2, 1:dimz-1] - u[1, 1:dimx-1, r:1, 1:dimz-1])
    d = dimz-1
    u[0, 1:dimx-1, 1:dimy-1, dimz-2:d] = u[1, 1:dimx-1, 1:dimy-1, dimz-3:d-1] \
        + (k[1:dimx-1, 1:dimy-1, dimz-2: d]-1)/(k[1:dimx-1, 1:dimy-1, dimz-2:d]+1) \
        * (u[0, 1:dimx-1, 1:dimy-1, dimz-3:d-1] - u[1, 1:dimx-1, 1:dimy-1, dimz-2:d])
    d = 0
    u[0, 1:dimx-1, 1:dimy-1, d:1] = u[1, 1:dimx-1, 1:dimy-1, d+1:2] \
        + (k[1:dimx-1, 1:dimy-1, d:1]-1)/(k[1:dimx-1, 1:dimy-1, d:1]+1) \
        * (u[0, 1:dimx-1, 1:dimy-1, d+1:2] - u[1, 1:dimx-1, 1:dimy-1, d:1])

@mlab.animate(delay=20)
def update_loop():
    fig = mlab.figure(bgcolor=(0, 0, 0))
    src = mlab.pipeline.scalar_field(u[0])
    volume = mlab.pipeline.volume(src)
    tick = 0
    while True:
        tick += 0.1
        mlab.view(azimuth=4*tick, elevation=80, distance=200, focalpoint=(int(dimx/2), int(dimy/2), 20))

        u[0:2, 1:dimx-1, 1:dimy-1, 1:3] += 2 * math.sin(tick*1.5)
        update(u)

        src.mlab_source.scalars = u[0]
        absmax = 1 

        # define opacity function
        otf = PiecewiseFunction()
        for val, opacity in [(absmax, 1), (absmax * 0.2, 0), (-absmax, 1), (-absmax * 0.2, 0)]:
            otf.add_point(val, opacity)
        volume._volume_property.set_scalar_opacity(otf) 

        # define color function (red-black-green)
        ctf = ColorTransferFunction()
        for p in [(absmax, 0, 1, 0), (0, 0, 0, 0), (-absmax, 1, 0, 0)]:
            ctf.add_rgb_point(*p)
        volume._volume_property.set_color(ctf)

        yield

if __name__ == "__main__":
    animate = update_loop()
    mlab.show()

Literaturnachweis

J. Douglas Faires, Richard L. Burden: "Numerische Methoden: Näherungsverfahren und ihre praktische Anwendung." Spektrum Akademischer Verlag; ISBN 3-86025-332-8; 1995
Hans Stefani, Gerhard Kluge: "Theoretische Mechanik" Spektrum Akademischer Verlag; ISBN 3-86025-284-4; 1995
Bengt Fornberg: "Generation of finite difference formulas on arbitrarily spaced grids." In: Mathematics of Computation. Band 51, Nr. 184, 1988, Seite 702
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
unbekannter Autor: "Finite Difference Methods" online; unbekannter Herausgeber