introduction

Bonjour. Récemment, Cartelet est accro au jeu incompréhensible de l'écriture d'une version discrète d'une équation différentielle sous la forme d'un produit de matrices dans un code court. Cet article est extrait de l'article de @ simonritchie, Animation de motif de dessin avec le modèle Gray-Scott C'est juste une question d'écrire brièvement à l'aide d'une matrice. Veuillez noter que cette fois, le seul avantage de l'utilisation d'une matrice est qu'elle raccourcit le code. De plus, je n'ai fait référence à rien, et je cours sur ma propre route, donc il peut y avoir beaucoup de points étranges, mais pardonnez-moi s'il vous plaît.

Pour une explication détaillée du modèle Gray-Scott, lisons l'article ci-dessus et faisons une file d'attente immédiatement.

Discrimination de la différenciation

Premièrement, la formule originale est

\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)

Il montre comment les densités des deux substances représentées par u et v évoluent dans le temps tout en interagissant l'une avec l'autre. Du, Dv sont des coefficients de diffusion, f est un acronyme pour feed et k est un acronyme pour kill. Comme vous pouvez le voir dans la formule, plus f est grand, plus il est facile pour u d'augmenter en densité et v, et plus k est grand, plus il est probable que v diminue en densité. Je vais.

Rendons-le croustillant.

Le premier est le côté gauche. Le côté gauche est un différentiel partiel par rapport à t, et puisque nous ne connaissons pas les informations futures, nous prenons la différence avant.

\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}

Ensuite, sur le côté droit Pour le premier terme, le coefficient de diffusion est ignoré pour le moment et le laplacien est discriminé. Laplacien bidimensionnel

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

Le différentiel de second ordre provient de la somme des expansions de Taylor de la différence avant et de la différence arrière.

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

Parce qu'il peut être approximé comme

\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}

Ce sera. De plus, seules les variables auxquelles je pense sont indiquées entre parenthèses, mais $ u (x + \ Delta x) $ etc. représente $ u (t, x + \ Delta x, y) $ etc.

Matrice

Maintenant que nous avons séparé les termes de différenciation, mettons-les sous la forme d'une matrice.

Depuis ce temps, nous supposons une implémentation utilisant un tableau à deux dimensions, $ x \ pm \ Delta x $ et $ y \ pm \ Delta y $ représentent les positions adjacentes en haut, en bas, à gauche et à droite. Premièrement, si vous écrivez la formule originale en utilisant la transformation de formule jusqu'à présent (tout ci-dessous sera écrit en nombres égaux),

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

Ce sera. Si cela se résume à $ \ Delta x = \ Delta y \ equiv \ Delta r $ dans le même ordre et dans la même section de position,

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

(Les termes $ u (t, x, y) v (t, x, y) ^ 2 $ interagissent et sont séparés car ils doivent être calculés séparément à chaque fois). 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}

L'équation ci-dessus utilise une matrice

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

Peut être exprimé comme. C'est la fin de la transformation. Notez que la partie coefficient à deux items (non-u ou v) du produit de la matrice est transposée puis multipliée s'il ne s'agit pas d'une matrice symétrique. Puisqu'il est ajouté deux fois dans les directions x et y, la composante diagonale est ajustée de 1/2.

** Je n'ai rien fait de difficile ** C'est difficile à regarder, mais ne vous inquiétez pas, le programme le rend extrêmement simple à écrire.

la mise en oeuvre

Certaines expressions sont censées être exécutées dans Jupyter Notebook. La partie expliquée dans la transformation de formule ci-dessus ne comprend que les quatre lignes marquées d'une étoile (*).

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

#Paramètres
dt = 1
dr = 1e-2
Du = 2e-5 #Coefficient de diffusion(u)
Dv = 1e-5 #Coefficient de diffusion(v)
f = 0.022 #feed
k = 0.051 #kill

Du /= dr**2;Dv /= dr**2 #Pause à l'avance

n = 256 #Un côté de la taille de la scène
u = np.ones((n, n));v = np.zeros_like(u) #u,initialisation de v
eye = np.eye(n) #Matrice unitaire
diffu = Du * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Du + f) * eye * .5 #Matrix pour vous*
diffv = Dv * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Dv + f + k) * eye * .5 #Matrice pour v*
sqsize = 10 #La moitié d'un côté du carré initial
u[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .5 #Définition des conditions initiales(u)
v[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .25 #Définition des conditions initiales(v)

def update(i):
    global u,v
    for t in range(10): #La croissance est un peu lente, alors retournez-la 10 fois par mise à jour
        uvv = u * v * v #La partie interaction de l'exemple*
        u, v = u + (diffu@u + u@diffu - uvv + f)*dt,v + (diffv@v + v@diffv + uvv) * dt #u,v mise à jour*
    im.set_array(u) #Pour le graphique ci-dessous
    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()

Pour les tableaux numpy, utilisez l'opérateur numpy.dot () ou @ pour calculer le produit de la matrice, et l'opérateur numpy.multiply () ou * pour calculer le produit élément par élément (produit adamal). Je peux.

résultat

Résumé

C'est incroyable qu'un motif organique naisse sans permission ... Regarder la formule est encombrée, alors regardons le code.