Erstellen einer Reaktions-Diffusions-Simulation (German)

Alexander Schoch,Fri Aug 16 2024•math python german

In diesem Artikel schauen wir uns an, wie wir diese Animation (mit einer etwas tieferen Auflösung) erstellen können:

Diese Animation zeigt ein Reaktions-Diffusions-System, ein System mit mehreren Komponenten (z.B. Chemikalien oder Bakterien), welche miteinander Reagieren können und über den Raum diffundieren. Für die Simulation in diesem Artikel schauen wir uns dabei spezifisch das Gray-Scott-Modell an.

Das Gray-Scott-Modell beschreibt ein System mit den folgenden Regeln:

Das System kennt zwei Komponenten, $A$ und $B$ , mit deren Konzentrationen $a(x, y, t)$ und $b(x, y, t)$ an jedem Pixel,
Die Komponente $A$ wird dem System ständig zugeführt,
Die Komponente $B$ wird ständig aus dem System entfernt,
Ein Teilchen $A$ und zwei Teilchen $B$ reagieren zu drei Teilchen $B$ ,
Beide Komponenten diffundieren mit unterschiedlicher Rate.

Dieses System kann mittels der folgenden partiellen Differenzialgleichungen augedrückt werden (Quelle):

\frac{\partial a}{\partial t} = D_\text{A} \nabla^2 a - ab^2 + F(1 - a)

\frac{\partial b}{\partial t} = D_\text{B} \nabla^2 b + ab^2 - (F + k) b

Brechen wir diese Gleichungen auf und betrachten wir ihre einzelnen Teile:

$\partial a / \partial t$ ist die Änderung der Konzentration von $A$ über die Zeit.
$D_\text{A} \partial^2 a / \partial x^2$ ist die Diffusion von $A$ und entspricht dem Zweiten Fick’schen Gesetz (analog zur Wärmegleichung).
$ab^2$ ist die Reaktion, gegeben als Rate Equation (wir multiplizieren einfach die Konzentrationen mit dem stöchiometrischen Faktor als Exponent). Hier wird $A$ verbraucht (daher das $-$ ) und $B$ erzeugt (daher das $+$ ).
$F(1-a)$ ist die Zufuhr von $A$ . Der Faktor $(1-a)$ sorgt dafür, dass die Konzentration von $A$ den Wert 1 nicht überschreitet, indem er diesen Term für grössere Werte von $a$ kleiner macht.
$(F+k)b$ ist das Entfernen von $B$ . Der Term $(F+k)$ sorgt dafür, dass die Entfernungsrate von $B$ immer grösser ist als die Zufuhr von $A$ , und der Skalierungsfaktor $b$ erhöht die Entfernung für höhere Werte von $b$ .

Diese Gleichung kann für unsere Zwecke diskretisiert werden (Quelle):

\Delta a = (D_\text{A} \nabla^2 a - ab^2 + F(1-a)) \Delta t

\Delta b = (D_\text{B} \nabla^2 b + ab^2 - (F+k)b) \Delta t

wobei $\Delta a$ die Änderung in den Konzentration von $A$ in einem Zeitschritt ist. $\nabla^2$ (Laplace-Operator, zweite Ableitung im Ort) kann über einen 3x3 Gaussischen Kernel (auch als Faltung oder Blur bezeichnet) berechnet werden. Wir schauen uns das nachher noch an.

Schauen wir uns an, wie wir das in Python implementieren.

Der Code

Zunächst erstellen wir unsere virtuelle Umgebung und installieren alle erforderlichen Pakete:

python -m venv venv
source venv/bin/activate
 
pip install numpy matplotlib alive_progress

Als nächstes importieren wir diese Pakete in unsere neu erstellte Python-Datei:

import numpy as np
import matplotlib.pyplot as plt
import os
 
from alive_progress import alive_bar

Definiere die Parameter, die das Ergebnis der Simulation steuern. Ich habe hier die Standardwerte von Karlsims übernommen. R_init beschreibt den Radius des Kreises in Pixeln, den wir als anfänglichen Seed verwenden werden. num_frames definiert, wie viele Frames wir rendern wollen. Hier ergeben 1200 Bilder mit einer Framerate von 60 FPS eine Animation von 20 Sekunden.

N = 400           # Anzahl Pixel in beiden Dimensionen
T = 5000          # Anzehl Zeitschritte
DA = 1            # Diffusivität von A
DB = 0.5          # Diffusivität von B
f = 0.055         # Zufuhrrate von A
k = 0.062         # Entfernunsrate von B
dt = 1            # Zeitschritt
R_init = 5        # Radius des Seed-Kreises
num_frames = 1200 # Anzahl der zu generierenen Frames

Nun berechnen wir einige Parameter, die wir für die Bilderzeugung verwenden werden. frame_interval ist eine ganze Zahl, die das Intervall für die Bilderzeugung speichert (z.B. “wir wollen jedes 12. Frame ein Bild ausgeben”). // ist die Python Ganzzahldivision, die sicherstellt, dass frame_interval eine Ganzzahl ist. num_digits speichert die maximale Anzahl der Ziffern, die die Bildnummer enthalten kann.

frame_interval = T // num_frames
frame_nr = 0
num_digits = int(np.log10(num_frames)) + 1 # Für die Generation von Filenames

Für unsere Anfangsbedingung definieren wir zwei Matrizen A und B, wobei A mit 1 und B mit 0 gefüllt wird. Dann erstellen wir eine Maske, die immer dann True ist, wenn das entsprechende Pixel weniger als R_init vom Zentrum des Bildes entfernt ist, und ansonsten False. Dann setzen wir die Konzentration von $B$ in diesem Kreis auf 1.

A = np.ones((N, N))
B = np.zeros((N, N))
X, Y = np.meshgrid(range(N), range(N))
mask = np.sqrt((X - N//2)**2 + (Y - N//2)**2) < R_init
B[mask] = 1

Eine kurze Demonstration, wie das genau funktioniert:

In [1]: import numpy as np
 
In [2]: X, Y = np.meshgrid(range(5), range(5))
 
In [3]: X
Out[3]: 
array([[0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4]])
 
In [4]: Y
Out[4]: 
array([[0, 0, 0, 0, 0],
       [1, 1, 1, 1, 1],
       [2, 2, 2, 2, 2],
       [3, 3, 3, 3, 3],
       [4, 4, 4, 4, 4]])
 
In [5]: N = 5
 
In [6]: np.sqrt((X - N//2)**2 + (Y - N//2)**2) < 2
Out[6]: 
array([[False, False, False, False, False],
       [False,  True,  True,  True, False],
       [False,  True,  True,  True, False],
       [False,  True,  True,  True, False],
       [False, False, False, False, False]])

Um sicherzugehen, dass dies auch wirklich so funktioniert, wie wir es uns wünschen, plotten wir das kurz.

plt.imshow(B)
plt.colorbar()
plt.show()

Anfangsbedingung

Nun müssen wir den oben beschriebenen mathematischen Ausdruck umsetzen. Beachte, dass es wichtig ist, zuerst dA und dB zu berechnen und sie anschliessend zu A und B zu addieren. Andernfalls würden wir A bereits überschreiben, während wir die alten Werte noch für die Berechnung von B benötigen.

t % frame_interval == 0 prüft, ob die aktuelle Bildnummer ein Vielfaches von frame_interval ist. Wenn ja, wird ein Bild der aktuellen Konzentration gespeichert, berechnet als $B / (A + B)$ . Dies ergibt eine Zahl zwischen 0 und 1 für jedes Pixel.

Die Funktionen convolution und save_image werden wir weiter unten schreiben.

with alive_bar(T) as bar: # Fortschritts-Balken
    for t in range(T):
        dA = DA * convolution(A) - A * B * B + f * (1 - A)
        dB = DB * convolution(B) + A * B * B - (k + f) * B
        A += dA * dt
        B += dB * dt
 
        if(t % frame_interval == 0):
            save_image(B / (A + B), frame_nr)
            frame_nr += 1
 
        bar()

Faltung

Der oben erwähnte 3x3-Gauss-Kernel beschreibt im Grunde den gewichteten Durchschnitt aller Zellen innerhalb eines 3x3-Quadrats um ein beliebiges Pixel. Hier sind die spezifischen Gewichte

-1 für das Pixel selbst,
0.2 für die vier direkten Nachbarn, und
0.05 für die vier diagonalen Pixel.

Beispielsweise hat das mittlere Pixel des Arrays

[[ 1, 2, 3],
 [ 4, 5, 6],
 [ 7, 8, 9]]

eine Faltung von

\text{Faltung}(\text{Mitte}) = -5 + 0.2(2 + 4 + 6 + 8) + 0.05(1 + 3 + 7 + 9) = 0

Bevor wir dies umsetzen, gibt es noch eine weitere Sache zu bedenken: Was machen wir am Rande des Gebiets? Am äussersten linken Rand gibt es keinen linken Nachbarn, also müssen wir uns überlegen, was in diesem Fall zu tun ist. Zwei einfache Möglichkeiten sind

Ignorieren der Randzellen und Anwenden des Algorithmus auf alle nicht randständigen Pixel
die Domäne herumloopen, so dass der linke Nachbar einer linken Randzelle eine rechte Randzelle ist

In diesem Artikel werden wir Letzteres verwenden, da es eine sehr praktische Funktion für genau diesen Zweck gibt: np.roll. Diese Funktion nimmt ein Array, eine Verschiebung und eine Achse als Input und rollt die Werte um diese Verschiebung in der Dimension axis.

In [1]: X
Out[1]: 
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
 
In [2]: np.roll(X, 1, axis=0) # Rolle vertikal
Out[2]: 
array([[7, 8, 9],
       [1, 2, 3],
       [4, 5, 6]])
 
In [3]: np.roll(X, 1, axis=1) # Rolle horizontal
Out[3]: 
array([[3, 1, 2],
       [6, 4, 5],
       [9, 7, 8]])

Damit können wir nun unsere Faltugnsfunktion schreiben. Beachte, dass diese Funktion oberhalb des Main Loop definiert werden muss.

def convolution(X):
    return -X
    + 0.2 * (
        np.roll(X, 1, axis=0) 
        + np.roll(X, 1, axis=1) 
        + np.roll(X, -1, axis=0) 
        + np.roll(X, -1, axis=1)) 
    + 0.05 * (
        np.roll(X, (1, 1), axis=(0, 1)) 
        + np.roll(X, (-1, 1), axis=(0, 1)) 
        + np.roll(X, (1, -1), axis=(0, 1)) 
        + np.roll(X, (-1, -1), axis=(0, 1)))

Erstellen des Videos

Zu guter Letzt plottet die Funktion save_image die Konzentration von $B$ unter Verwendung von matplotlib.pyplot.imshow mit einer schönen Color Map und speichert sie z.B. in frames/0001.png. Die oberen vier Zeilen der Funktion entfernen den weissen Rand um den Plot und stellen ein Seitenverhältnis von 1:1 ein. Diese Funktion muss ebenfalls oberhalb des Main Loop definiert werden. Stelle ausserdem sicher, dass der Ordner frames existiert!

def save_image(X, i):
    fig = plt.figure(figsize=(1, 1))
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)
 
    ax.imshow(X, cmap='inferno')
    filename = "%s.png" % str(i).zfill(num_digits)
    plt.savefig(os.path.join('frames', filename), dpi=N)
    plt.close('all')

Nachdem wir das Skript ausgeführt und alle Bilder erzeugt haben, können wir sie mit ffmpeg zu einem einzigen Video zusammenfügen.

ffmpeg -r 60 -pattern_type glob -i 'frames/????.png' -c:v libx264 video.mp4

Das war’s!

Die Variation der Parameter $D_\text{A}$ , $D_\text{B}$ , $F$ und $k$ ermöglicht eine sehr grosse Bandbreite an Mustern. Es ist auch möglich, die Parameter über die Zeit oder den Raum zu verändern, um sehr hübsche Animationen zu erhalten!

Teile gerne deine fertige Animation mit mir! Du kannst sie mir per E-Mail schicken (schochal[at]ethz.ch) oder sie auf Mastodon posten und mich taggen (@alexander_schoch@linuxrocks.online). Ich würde mich auch über Feedback zu diesem Beitrag freuen!

Kompletter Code

import numpy as np
import matplotlib.pyplot as plt
import os
 
from alive_progress import alive_bar
 
N = 400
T = 20000
DA = 1
DB = 0.5
f = 0.055
k = 0.062
dt = 1
R_init = 5
num_frames = 1200
 
frame_interval = T // num_frames
frame_nr = 0
num_digits = int(np.log10(num_frames)) + 1 # For filename generation
 
A = np.ones((N, N))
B = np.zeros((N, N))
X, Y = np.meshgrid(range(N), range(N))
mask = np.sqrt((X - N//2)**2 + (Y - N//2)**2) < R_init
B[mask] = 1
 
def convolution(X):
    return -X
    + 0.2 * (
        np.roll(X, 1, axis=0) 
        + np.roll(X, 1, axis=1) 
        + np.roll(X, -1, axis=0) 
        + np.roll(X, -1, axis=1)) 
    + 0.05 * (
        np.roll(X, (1, 1), axis=(0, 1)) 
        + np.roll(X, (-1, 1), axis=(0, 1)) 
        + np.roll(X, (1, -1), axis=(0, 1)) 
        + np.roll(X, (-1, -1), axis=(0, 1))) 
 
def save_image(X, i):
    fig = plt.figure(figsize=(1, 1))
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)
 
    ax.imshow(X, cmap='inferno')
    filename = "%s.png" % str(i).zfill(num_digits)
    plt.savefig(os.path.join('frames', filename), dpi=N)
    plt.close('all')
 
 
with alive_bar(T) as bar:
    for t in range(T):
        dA = DA * convolution(A) - A * B * B + f * (1 - A)
        dB = DB * convolution(B) + A * B * B - (k + f) * B
        A += dA * dt
        B += dB * dt
 
        if(t % frame_interval == 0):
            save_image(B / (A + B), frame_nr)
            frame_nr += 1
 
        bar()

Comments

Subscribe via RSS.