#はじめに
こんにちは。最近、微分方程式を離散化したものを行列の積の形にして如何に短いコードで書くか、という意味不明な遊びにハマっているCarteletです。
本記事は、こちらの@simonritchieさんの記事、
Gray-Scottモデルで模様のアニメーションを描く
を行列を使って短く書いてやろうというだけのモノです。
今回、行列を使うことの利点はコードが短くなるくらいしかないのでご注意ください。
また、何かを参考にしたわけでもなく、独自路線をひた走っていますのでおかしな点が多々あるかもしれませんがご容赦ください。

Gray-Scottモデルについての細かい説明は上記の記事を読んでもらうことにして、早速行列にしていきましょう。

#微分の離散化
まず、もととなる式は

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

で、それぞれ u, v で表される2つの物質の密度が相互作用しながら時間変化する様子を表しています。
Du, Dvは拡散係数、
f は feed、k は kill の頭文字で、式から分かる通り、f が大きいほど u の密度が増えやすく、v の密度が減少しやすくなり、k が大きいほど v の密度が減少しやすくなっています。

ではサクサク離散化していきましょう。

まずは左辺です。
左辺は t に関する偏微分で、未来の情報は分からないので前進差分を取ります。

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

次に右辺、
第一項はとりあえず拡散係数は無視して、ラプラシアンの離散化をします。2次元のラプラシアンは、

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

であり、2階微分は前進差分と後進差分のテイラー展開の和から

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

と近似できるので、

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

となります。
また、考えている変数のみを括弧内に示しましたが　$u(x+\Delta x)$　などは　$u(t,x+\Delta x,y)$　などを表しています。

#行列化
微分の項の離散化ができたので次は行列の形に落とし込んでいきましょう。

今回は２次元配列を用いての実装を想定しているため、　$x\pm\Delta x$　や　$y\pm\Delta y$　は上下左右で隣り合った位置を表しています。
まず、ここまでの式変形を用いてもとの式を書き下すと（以下全部等号で書きます）、

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

となります。これを　$\Delta x=\Delta y \equiv\Delta r$　として、同次、同位置の項でまとめると、

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

となり（　$u(t,x,y)v(t,x,y)^2$　の項は相互作用しており、毎回別個計算する必要があるため分離しています）、
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}

と表すと、上式は行列を用いて

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

と表すことができます。
変形はここまでです。行列の積の二項目の係数（uやvでない方）部分は対称行列でない場合は転置してから掛ける点に注意してください。
x方向、y方向の2回分足し合わせているため対角成分を1/2倍して調節しています。

何も難しいことをしていないのに見ることすら面倒くさい形になってしまいましたが、プログラムでは至極シンプルに書けますのでご安心ください。

#実装
一部Jupyter Notebookでの実行を想定した表現があります。
上の式変形で説明した部分は星(*)印のついた4行だけです。

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

#パラメータ
dt = 1
dr = 1e-2
Du = 2e-5 #拡散係数(u)
Dv = 1e-5 #拡散係数(v)
f = 0.022 #feed
k = 0.051 #kill

Du /= dr**2;Dv /= dr**2 #予め割っておく

n = 256 #ステージサイズの一辺
u = np.ones((n, n));v = np.zeros_like(u) #u, vの初期化
eye = np.eye(n) #単位行列
diffu = Du * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Du + f) * eye * .5 #u用の行列*
diffv = Dv * (np.roll(eye, -1, axis=1) + np.roll(eye, 1, axis=1)) - (4 * Dv + f + k) * eye * .5 #v用の行列*
sqsize = 10 #初期条件の正方形の一辺の半分
u[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .5 #初期条件の設定(u)
v[n//2-sqsize:n//2+sqsize, n//2-sqsize:n//2+sqsize] = .25 #初期条件の設定(v)

def update(i):
    global u,v
    for t in range(10): #成長が若干遅いので一更新に付き10回回す
        uvv = u * v * v #例の相互作用の部分*
        u, v = u + (diffu@u + u@diffu - uvv + f)*dt,v + (diffv@v + v@diffv + uvv) * dt #u,vの更新*
    im.set_array(u) #以下プロット用
    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()

numpy配列ではnumpy.dot()または@演算子を用いて行列の積を、numpy.multiply()または*演算子で要素ごとの積（アダマール積）を計算することができます。

#結果

#まとめ
勝手に有機的な模様生まれるのしゅごい…
式を見ると雑然としているのでコードを見ましょう。