Einführung

Hallo. In letzter Zeit ist Cartelet süchtig nach dem unklaren Spiel, wie man eine diskrete Version einer Differentialgleichung in Form eines Matrizenprodukts in einem Kurzcode schreibt. Dieser Artikel stammt aus dem Artikel von @ simonritchie, Zeichnen Sie eine Musteranimation mit dem Gray-Scott-Modell Es geht nur darum, kurz mit einer Matrix zu schreiben. Bitte beachten Sie, dass diesmal der einzige Vorteil der Verwendung einer Matrix darin besteht, dass der Code verkürzt wird. Außerdem habe ich mich auf nichts bezogen und laufe auf meiner eigenen Route, daher kann es viele seltsame Punkte geben, aber bitte verzeihen Sie mir.

Für eine detaillierte Erklärung des Gray-Scott-Modells lesen wir den obigen Artikel und stellen sofort eine Warteschlange auf.

Diskriminierung der Differenzierung

Erstens lautet die ursprüngliche Formel

\frac{\partial u}{\partial t}
= D_u \Delta u - uv^2 + f(1 - u)

\frac{\partial v}{\partial t}
= D_v \Delta v + uv^2 - v(f + k)

Es zeigt, wie sich die Dichten der beiden durch u und v dargestellten Substanzen im Laufe der Zeit ändern, während sie miteinander interagieren. Du, Dv sind Diffusionskoeffizienten, f ist ein Akronym für feed und k ist ein Akronym für kill. Ich werde.

Machen wir es knusprig.

Zuerst ist die linke Seite. Die linke Seite ist ein partielles Differential in Bezug auf t, und da wir keine zukünftigen Informationen kennen, nehmen wir die Vorwärtsdifferenz.

\frac{\partial u}{\partial t}
\approx
\frac{u(t+\Delta t)-u(t)}{\Delta t}

\frac{\partial v}{\partial t}
\approx
\frac{v(t+\Delta t)-v(t)}{\Delta t}

Als nächstes auf der rechten Seite Für den ersten Term wird der Diffusionskoeffizient vorerst ignoriert und Laplace wird unterschieden. Zweidimensionaler Laplace

\Delta 
= \nabla^2 
= \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}

Die Differenz zweiter Ordnung ergibt sich aus der Summe der Taylor-Erweiterungen der Vorwärtsdifferenz und der Rückwärtsdifferenz.

\frac{\partial^2 u}{\partial x^2}
\approx
 \frac{u(x+\Delta x) - 2u(x) + u(x-\Delta x)}{\Delta x^2}

Weil es als angenähert werden kann

\Delta u
=
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}
\approx
 \frac{u(x+\Delta x) - 2u(x) + u(x-\Delta x)}{\Delta x^2}
+  \frac{u(y+\Delta y) - 2u(y) + u(y-\Delta y)}{\Delta y^2}

\Delta v
=
\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}
\approx
 \frac{v(x+\Delta x) - 2v(x) + v(x-\Delta x)}{\Delta x^2}
+  \frac{v(y+\Delta y) - 2v(y) + v(y-\Delta y)}{\Delta y^2}

Es wird sein. Außerdem werden nur die Variablen, an die ich denke, in Klammern angezeigt, aber $ u (x + \ Delta x) $ usw. steht für $ u (t, x + \ Delta x, y) $ usw.

Matrix

Nachdem wir die Differenzierungsbegriffe getrennt haben, lassen Sie uns sie in die Form einer Matrix bringen.

Seit dieser Zeit gehen wir von einer Implementierung mit einem zweidimensionalen Array aus. $ X \ pm \ Delta x $ und $ y \ pm \ Delta y $ repräsentieren benachbarte Positionen oben, unten, links und rechts. Wenn Sie die ursprüngliche Formel mit der bisherigen Formeltransformation aufschreiben (alle folgenden werden in gleicher Anzahl geschrieben),

u(t+\Delta t,x,y)
=
u(t,x,y)
+\Biggl[
D_u\Biggl(
 \frac{u(t,x+\Delta x,t) - 2u(t,x,y) + u(t,x-\Delta x,y)}{\Delta x^2}
+
\frac{u(t,x,y+\Delta y) - 2u(t,x,y) + u(t,x,y-\Delta y)}{\Delta y^2}
\Biggr)
-
u(t,x,y)v(t,x,y)^2
+
 f\bigr(1 - u(t,x,y)\bigl)
\Biggr]
\Delta t

v(t+\Delta t,x,y)
=
v(t,x,y)
+\Biggl[
D_v\Biggl(
 \frac{v(t,x+\Delta x,t) - 2v(t,x,y) + v(t,x-\Delta x,y)}{\Delta x^2}
+
\frac{v(t,x,y+\Delta y) - 2v(t,x,y) + v(t,x,y-\Delta y)}{\Delta y^2}
\Biggr)
+
u(t,x,y)v(t,x,y)^2
-
 v(t,x,y)\bigr(f + k \bigl)
\Biggr]
\Delta t

Es wird sein. Wenn dies als $ \ Delta x = \ Delta y \ equiv \ Delta r $ in derselben Reihenfolge und demselben Positionsabschnitt zusammengefasst ist,

u(t+\Delta t,x,y)
=
u(t,x,y)
+
\Biggl[
-
u(t,x,y)\bigr(
4\frac{D_u}{\Delta r^2}+f
\bigl)
+
\frac{D_u}{\Delta r^2}\Biggl(
 u(t,x+\Delta x,t) + u(t,x-\Delta x,y)+u(t,x,y+\Delta y) + u(t,x,y-\Delta y)
\Biggr)
+
f
-
u(t,x,y)v(t,x,y)^2
\Biggr]
\Delta t

v(t+\Delta t,x,y)
=
v(t,x,y)
+
\Biggl[
-
v(t,x,y)\bigr(
4\frac{D_v}{\Delta r^2}+f+k
\bigl)
+
\frac{D_v}{\Delta r^2}\Biggl(
 v(t,x+\Delta x,t) + v(t,x-\Delta x,y)+u(t,x,y+\Delta y) + v(t,x,y-\Delta y)
\Biggr)
+
u(t,x,y)v(t,x,y)^2
\Biggr]
\Delta t

(Die Terme $ u (t, x, y) v (t, x, y) ^ 2 $ interagieren und werden getrennt, da sie jedes Mal separat berechnet werden müssen.) u, v

\begin{bmatrix}
u_{00} & u_{01} & \ldots & u_{ 0n }\\
u_{10} & u_{11} & \ldots & u_{ 1n }\\
\vdots & \vdots & \ddots & \vdots \\
u_{n0} & u_{n1} & \ldots & u_{ nn }\\
\end{bmatrix}

\begin{bmatrix}
v_{00} & v_{01} & \ldots & v_{ 0n }\\
v_{10} & v_{11} & \ldots & v_{ 1n }\\
\vdots & \vdots & \ddots & \vdots \\
v_{n0} & v_{n1} & \ldots & v_{ nn }\\
\end{bmatrix}

Die obige Gleichung verwendet eine Matrix

u(t+\Delta t,x,y)
=
u(t,x,y)
+
\Biggl[
\begin{bmatrix}
-\frac{(4\frac{D_u}{\Delta r^2}+f)}{2} & \frac{D_u}{\Delta r^2} & 0 & \ldots & \frac{D_u}{\Delta r^2}\\
\frac{D_u}{\Delta r^2} & -\frac{(4\frac{D_u}{\Delta r^2}+f)}{2} & \frac{D_u}{\Delta r^2} & \ldots & 0\\
0 & \frac{D_u}{\Delta r^2} & -\frac{(4\frac{D_u}{\Delta r^2}+f)}{2} & \ldots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\frac{D_u}{\Delta r^2} & 0 & 0 & \ldots & -\frac{(4\frac{D_u}{\Delta r^2}+f)}{2}\\
\end{bmatrix}
\begin{bmatrix}
u_{00} & u_{01} & \ldots & u_{ 0n }\\
u_{10} & u_{11} & \ldots & u_{ 1n }\\
\vdots & \vdots & \ddots & \vdots \\
u_{n0} & u_{n1} & \ldots & u_{ nn }\\
\end{bmatrix}
+
\begin{bmatrix}
u_{00} & u_{01} & \ldots & u_{ 0n }\\
u_{10} & u_{11} & \ldots & u_{ 1n }\\
\vdots & \vdots & \ddots & \vdots \\
u_{n0} & u_{n1} & \ldots & u_{ nn }\\
\end{bmatrix}
\begin{bmatrix}
-\frac{(4\frac{D_u}{\Delta r^2}+f)}{2} & \frac{D_u}{\Delta r^2} & 0 & \ldots & \frac{D_u}{\Delta r^2}\\
\frac{D_u}{\Delta r^2} & -\frac{(4\frac{D_u}{\Delta r^2}+f)}{2} & \frac{D_u}{\Delta r^2} & \ldots & 0\\
0 & \frac{D_u}{\Delta r^2} & -\frac{(4\frac{D_u}{\Delta r^2}+f)}{2} & \ldots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\frac{D_u}{\Delta r^2} & 0 & 0 & \ldots & -\frac{(4\frac{D_u}{\Delta r^2}+f)}{2}\\
\end{bmatrix}
+
f
-
u(t,x,y)v(t,x,y)^2
\Biggr]
\Delta t

v(t+\Delta t,x,y)
=
v(t,x,y)
+
\Biggl[
\begin{bmatrix}
-\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2} & \frac{D_v}{\Delta r^2} & 0 & \ldots & \frac{D_v}{\Delta r^2}\\
\frac{D_v}{\Delta r^2} & -\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2} & \frac{D_v}{\Delta r^2} & \ldots & 0\\
0 & \frac{D_v}{\Delta r^2} & -\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2} & \ldots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\frac{D_v}{\Delta r^2} & 0 & 0 & \ldots & -\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2}\\
\end{bmatrix}
\begin{bmatrix}
v_{00} & v_{01} & \ldots & v_{ 0n }\\
v_{10} & v_{11} & \ldots & v_{ 1n }\\
\vdots & \vdots & \ddots & \vdots \\
v_{n0} & v_{n1} & \ldots & v_{ nn }\\
\end{bmatrix}
+
\begin{bmatrix}
v_{00} & v_{01} & \ldots & v_{ 0n }\\
v_{10} & v_{11} & \ldots & v_{ 1n }\\
\vdots & \vdots & \ddots & \vdots \\
v_{n0} & v_{n1} & \ldots & v_{ nn }\\
\end{bmatrix}
\begin{bmatrix}
-\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2} & \frac{D_v}{\Delta r^2} & 0 & \ldots & \frac{D_v}{\Delta r^2}\\
\frac{D_v}{\Delta r^2} & -\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2} & \frac{D_v}{\Delta r^2} & \ldots & 0\\
0 & \frac{D_v}{\Delta r^2} & -\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2} & \ldots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\frac{D_v}{\Delta r^2} & 0 & 0 & \ldots & -\frac{(4\frac{D_v}{\Delta r^2}+f+k)}{2}\\
\end{bmatrix}
+
u(t,x,y)v(t,x,y)^2
\Biggr]
\Delta t

Kann ausgedrückt werden als. Dies ist das Ende der Transformation. Es ist zu beachten, dass der zweiteilige Koeffiziententeil (nicht u oder v) des Produkts der Matrix transponiert und dann multipliziert wird, wenn es sich nicht um eine symmetrische Matrix handelt. Da es zweimal in x- und y-Richtung hinzugefügt wird, wird die diagonale Komponente um 1/2 eingestellt.

** Ich habe nichts Schwieriges getan ** Selbst das Anschauen ist zu einem Problem geworden, aber keine Sorge, das Programm macht das Schreiben extrem einfach.

Implementierung

Es gibt einige Ausdrücke, die in Jupyter Notebook ausgeführt werden sollen. Der in der obigen Formeltransformation erläuterte Teil besteht nur aus den vier mit einem Stern (*) gekennzeichneten Zeilen.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
%matplotlib nbagg

#Parameter
dt = 1
dr = 1e-2
Du = 2e-5 #Diffusionskoeffizient(u)
Dv = 1e-5 #Diffusionskoeffizient(v)
f = 0.022 #feed
k = 0.051 #kill

Du /= dr**2;Dv /= dr**2 #Pause im Voraus

n = 256 #Eine Seite der Bühnengröße
u = np.ones((n, n));v = np.zeros_like(u) #u,Initialisierung von v
eye = np.eye(n) #Einheitsmatrix
diffu = Du * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Du + f) * eye * .5 #Matrix für u*
diffv = Dv * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Dv + f + k) * eye * .5 #Matrix für v*
sqsize = 10 #Die Hälfte einer Seite des ursprünglichen Quadrats
u[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .5 #Anfangsbedingungen einstellen(u)
v[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .25 #Anfangsbedingungen einstellen(v)

def update(i):
    global u,v
    for t in range(10): #Das Wachstum ist etwas langsam, also drehen Sie es 10 Mal pro Update
        uvv = u * v * v #Der Interaktionsteil des Beispiels*
        u, v = u + (diffu@u + u@diffu - uvv + f)*dt,v + (diffv@v + v@diffv + uvv) * dt #u,v Update*
    im.set_array(u) #Für Handlung unten
    return im,

fig = plt.figure()
im = plt.imshow(u,vmin = 0, vmax = 1, cmap = "twilight")
anim = FuncAnimation(fig, update, interval = 1, blit = True)
plt.show()

Verwenden Sie für numpy-Arrays den Operator "numpy.dot ()" oder "@", um das Matrixprodukt zu berechnen, und den Operator "numpy.multiply ()" oder "*", um das Element-für-Element-Produkt (Adamal-Produkt) zu berechnen. Ich kann.

Ergebnis

Zusammenfassung

Es ist erstaunlich, dass ein organisches Muster ohne Erlaubnis geboren wird ... Das Betrachten der Formel ist unübersichtlich. Schauen wir uns also den Code an.