Introduction

Hello. Recently, Cartelet is addicted to the incomprehensible play of how to write a discretized differential equation in the form of a matrix product in a short code. This article is from @ simonritchie's article, Draw pattern animation with Gray-Scott model It's just a matter of writing shortly using a matrix. Please note that this time, the only advantage of using a matrix is that it shortens the code. Also, I didn't refer to anything, and I'm running on my own route, so there may be many strange points, but please forgive me.

For a detailed explanation of the Gray-Scott model, let's read the above article and let's make a line immediately.

Discretization of differentiation

First, the original formula is

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

It shows how the densities of the two substances represented by u and v change with time while interacting with each other. Du, Dv are diffusion coefficients, f is an acronym for feed and k is an acronym for kill. As you can see from the formula, the larger f is, the more likely the density of u is to increase and the density of v is likely to decrease, and the larger k is, the more likely the density of v is to decrease. I will.

Let's break it up crisply.

First is the left side. The left side is the partial derivative with respect to t, and since we do not know the future information, we take the forward difference.

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

Next, on the right side, For the first term, the diffusion coefficient is ignored for the time being, and the Laplacian is discretized. Two-dimensional Laplacian

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

The second derivative is from the sum of the Taylor expansions of the forward difference and the backward difference.

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

Because it can be approximated to

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

It will be. Also, only the variables I'm thinking of are shown in parentheses, but $ u (x + \ Delta x) $ etc. represents $ u (t, x + \ Delta x, y) $ etc.

Matrix

Now that the differential term has been discretized, let's put it in the form of a matrix.

Since we are assuming implementation using a two-dimensional array this time, $ x \ pm \ Delta x $ and $ y \ pm \ Delta y $ represent adjacent positions in the vertical and horizontal directions. First, if you write down the original formula using the formula transformation so far (all below will be written with equal signs),

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

It will be. If this is summarized as $ \ Delta x = \ Delta y \ equiv \ Delta r $ in the same order and same position section,

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

(The terms $ u (t, x, y) v (t, x, y) ^ 2 $ are interacting and are separated because they need to be calculated separately each time). 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}

Expressed as, the above equation uses a 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

Can be expressed as. This is the end of the transformation. Note that the two-item coefficient (non-u or v) part of the matrix product is transposed and then multiplied if it is not a symmetric matrix. Since it is added twice in the x direction and the y direction, the diagonal component is adjusted by 1/2.

** I haven't done anything difficult ** It's awkward to even look at, but don't worry, the program makes it extremely simple to write.

Implementation

There are some expressions that are supposed to be executed in Jupyter Notebook. The part explained in the above formula transformation is only 4 lines marked with a star (*).

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

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

Du /= dr**2;Dv /= dr**2 #Divide in advance

n = 256 #One side of the stage size
u = np.ones((n, n));v = np.zeros_like(u) #u,initialization of v
eye = np.eye(n) #Identity matrix
diffu = Du * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Du + f) * eye * .5 #Matrix for u*
diffv = Dv * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Dv + f + k) * eye * .5 #Matrix for v*
sqsize = 10 #Half of one side of the initial square
u[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .5 #Setting initial conditions(u)
v[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .25 #Setting initial conditions(v)

def update(i):
    global u,v
    for t in range(10): #Growth is a little slow, so turn it 10 times per update
        uvv = u * v * v #Interaction part of the example*
        u, v = u + (diffu@u + u@diffu - uvv + f)*dt,v + (diffv@v + v@diffv + uvv) * dt #u,v update*
    im.set_array(u) #For plot below
    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()

For numpy arrays, use the numpy.dot () or @ operator to calculate the matrix product, and the numpy.multiply () or * operator to calculate the element-by-element product (Hadamard product). I can.

result

Summary

It's amazing that organic patterns are created without permission ... Looking at the formula is cluttered, so let's look at the code.