こんにちは。
今回は「チューリングパターンをつくる」という題でやっていきます。
皆さんは「ゲームさんぽ」という企画をご存じでしょうか。
様々なもののプロ・専門家とともにゲームをプレイするというものなのですが、専門知識が飛び出したり、ものの見方の違いを知ることができて面白いです。
その中のひとつに生物学者と一緒に釣りスピをやるという回があったのですが、そこで「チューリングパターン」という単語が出てきました。
大昔に「所さんの目がテン！」でも桝太一アナがチューリングパターンを描く企画をやっていたような記憶があるので、そちらで知ってるという方もいるかもしれません。
チューリングパターンは以下のような数式で表され、これによって生き物の模様を描くことができるようです。
$$ \begin{align} \frac{ \partial u}{ \partial t} & = f(u, v) + D_u \Delta u \\ \frac{ \partial v}{ \partial t} & = g(u, v) + D_v \Delta v \end{align} $$
この式は、ある成分$u$と$v$が$f(u,v)$または$g(u,v)$で反応し、$D_u \Delta u$または$D_v \Delta v$で拡散するという意味で、「反応拡散方程式」と呼ばれるようです。
（式の詳しい中身は後に説明します。）
これを見てしまった自分はprocessingで書くしかないと思い、この記事を書いている次第です。
長編になってしまったので目次を置いておきます。
拡散方程式
チューリングパターンの前に、練習がてら「一次元での拡散方程式」について書いてみます。
拡散方程式は以下のように表される式で、$u$が多いところは減り、$u$が少ないところは増えるというように、ある量$u$が均されていく、というような描像になります。
$$ \frac{ \partial u}{ \partial t} = D \Delta u $$
具体的には、物質の拡散や、熱の拡散（熱伝導）などが当てはまります。
$\Delta$は「ラプラシアン」や「ラプラス演算子」と呼ばれるもので、三次元でいえば以下のようになります。
$$ \begin{align} \Delta = \nabla^{2} & = \frac{ \partial^{2} }{ \partial x^{2} } + \frac{ \partial^{2} }{ \partial y^{2} } + \frac{ \partial^{2} }{ \partial z^{2} } \\ \nabla & = \frac{ \partial }{ \partial x } + \frac{ \partial }{ \partial y } + \frac{ \partial }{ \partial z } \end{align} $$
ここで、$\displaystyle \frac{ \partial }{ \partial x }$は微分（正確には偏微分）のことなので、$\nabla$（ナブラ）は「傾き」を表し、$\Delta$は「傾きの傾き」を表すといえます。
なぜ「ある量$u$が均される」という描像になるのか
この部分は高校数学の復習的な話になるため、わかるよという方は飛ばしてください。
（ここを押すと終わりまで飛びます。）
また、おおまかにわかれば良いので、偏微分$\displaystyle \frac{ \partial }{ \partial x }$と微分$\displaystyle \frac{ \mathrm{d} }{ \mathrm{d} x }$は同じものとして説明します。
ある関数$y = x^{3} - x$について考えます。
この関数について、一階微分$\displaystyle \frac{ \mathrm{d} }{ \mathrm{d} x } y$と二階微分$\displaystyle \frac{ \mathrm{d}^{2} }{ \mathrm{d} x^{2} } y$はそれぞれ、
$$ \frac{ \mathrm{d} }{ \mathrm{d} x } y = 3 x^{2} - 1~ , ~ \frac{ \mathrm{d}^{2} }{ \mathrm{d} x^{2} } y = 6 x $$
確か二階微分は数学IIIで習うことかと思いますが、要は微分した関数をもう一度微分すればいい話なので、難しくはないと思います。
グラフを描くというときには「増減表」をつくります。
（$y' = \frac{ \mathrm{d} }{ \mathrm{d} x } y$、$y'' = \frac{ \mathrm{d}^{2} }{ \mathrm{d} x^{2} } y$です。）
$$ \begin{array}{c|ccccccc} x & \cdots & - \frac{\sqrt{3}}{3} & \cdots & 0 & \cdots & \frac{\sqrt{3}}{3} & \cdots \\ \hline y' & + & 0 & & – & & 0 & + \\ \hline y'' & & - & & 0 & & + \\ \hline y & \nearrow & \frac{2 \sqrt{3}}{9} & \searrow & 0 & \searrow & - \frac{2 \sqrt{3}}{9} & \nearrow \end{array} $$
これを元にグラフを描くと、
ここからわかるように、グラフが凹んでいる（下に凸）のところで二階微分の値が正 (+) になっており、グラフが突き出ている（上に凸）のところで二階微分の値が負 (-) になっていることがわかります。
つまり先に挙げた拡散方程式では、$u$が小さい（下に凸）ところでは時間経過に伴って$u$の値に$D \Delta u$だけ足し（二階微分が正）、$u$が大きい（上に凸）ところでは時間経過に伴って$D \Delta u$だけ引く（二階微分が負）という意味となり、この結果$u$の大きさが同じになっていく（均される）ということです。
$$ \frac{ \partial u}{ \partial t} = D \Delta u $$
余計わからなくなってしまったかもしれませんが、百聞は一見に如かずということで、以降でつくっていくプログラムで理解してもらえればと思います。
終わり
拡散方程式の描画
それでは、拡散方程式についてプログラムにしてみます。
一次元で
簡単のためにまずは一次元でつくっていき、説明しやすいので熱伝導ということにします。
一次元での拡散方程式は、
$$ \frac{ \partial u}{ \partial t} = D \frac{ \partial^{2} }{ \partial x^{2} } u $$
目標としては、何か棒のようなものを横から見て、端から加熱したときに、棒の温度分布がどのように変わっていくかを図示する、というような形を目指します。
まずは、棒の温度が適当な温度分布を持っているとしておきます。
int n = 11; //棒上にとる点の数 float l; //点の間隔 PVector[] point = new PVector[n]; //点の座標（y成分が温度にあたる） void setup(){ size(500, 500); l = width/(n-1); float xoff = 0; for(int i = 0; i < n; i ++){ //初状態での温度分布 point[i] = new PVector(l*i, height - 100*noise(xoff)); xoff += 0.1; } } void draw(){ background(255); for(int i = 0; i < n-1; i ++){ //分布の図示 line(point[i].x, point[i].y, point[i+1].x, point[i+1].y); } }
続いて、$\Delta u$の計算をします。
draw()内に、
float[] laplacian = new float[n]; for(int i = 1; i < n-1; i ++){ laplacian[i] = (point[i+1].y-2*point[i].y+point[i-1].y)/(l*l); }
微分の定義が以下のようになることからもわかる通り、微分は傾きを求める計算です。
$y = f(x)$という関数について、微分の定義は以下のようになります。
$$ \frac{ \mathrm{d} }{ \mathrm{d} x } f(x) = \lim_{ h \to 0 } \frac{ f(x+h) - f(x) }{ h } $$
（$h$を限りなく0に近づけることで、$x=x$のときの一点の傾きに近づいていく。）
$$ \begin{align} \Delta y & = \frac{ \mathrm{d}^{2} }{ \mathrm{d}^{2} x } y \\ & = \frac{ \mathrm{d} }{ \mathrm{d} x } \left( \frac{ \mathrm{d} }{ \mathrm{d} x } f(x) \right) \\ & = \frac{ \frac{ \mathrm{d} }{ \mathrm{d} x } f \left( x+\frac{ \Delta x }{2} \right) - \frac{ \mathrm{d} }{ \mathrm{d} x } f \left( x-\frac{ \Delta x }{2} \right) }{\Delta x} \end{align} $$
ここで、
$$ \begin{align} \frac{ \mathrm{d} }{ \mathrm{d} x } f \left( x+\frac{ \Delta x }{2} \right) & = \frac{ f(x+\Delta x) - f(x) }{\Delta x} \\ \frac{ \mathrm{d} }{ \mathrm{d} x } f \left( x-\frac{ \Delta x }{2} \right) & = \frac{ f(x) - f(x-\Delta x) }{\Delta x} \end{align} $$
よって、
$$ \Delta y = \frac{ f(x+\Delta x) - 2f(x) + f(x-\Delta x) }{ \Delta x^{2} } $$
これが、コード中の
laplacian[i] = (point[i+1].y - 2*point[i].y + point[i-1].y) / (l*l);
の部分に当たります。
最後に、拡散方程式に当てはめて温度の時間変化を計算し、温度を更新する処理を追加します。
draw()内に、
float D = 1; float[] du = new float[n]; for(int i = 0; i < n; i ++){ du[i] = D * laplacian[i]; point[i].y += du[i]; }
拡散係数$D$はとりあえず1としておきました。
拡散方程式の左辺は$\displaystyle \frac{ \partial u}{ \partial t}$となっていますが、processingではdraw()を実行するたびに画面が更新されていくので、それを時間経過とみなすことができます。
更新速度を2倍にしたければ、duの計算を
du[i] = D * laplacian[i] * 2;
とするか、setup()に
frameRate(120); //デフォルトでframeRate(60)
と追加することでできます。
以上で一応動きはしますが、ひとつ問題があります。
ラプラシアンの計算のところを見ると[i+1]や[i-1]の点が必要なので、処理の都合上、両端の点（point[0]とpoint[n-1]）では$\Delta u$が計算できず、したがって温度変化も無いことになります。
そのため、両端の点での温度はその1つ内側の点の温度と同じであるということにしておきます。
（いわゆる「境界条件」というやつで、ここをどうするかでも結果は変わってきます。）
draw()内に、
point[0].y = point[1].y; point[n-1].y = point[n-2].y;
と追加すればOKです。
これにて、一次元での熱伝導を図示することができるようになりました。
例えば、左端から加熱した（point[0]とpoint[1]での温度が高温で一定）とすると、
int n = 101; float l; PVector[] point = new PVector[n]; float D; void setup(){ size(500, 500); l = width/(n-1); float xoff = 0; for(int i = 0; i < n; i ++){ point[i] = new PVector(l*i, height - 100*noise(xoff)); xoff += 0.1; } D = 1; } void draw(){ background(255); point[0].y = 100; point[1].y = 100; for(int i = 0; i < n-1; i ++){ line(point[i].x, point[i].y, point[i+1].x, point[i+1].y); } float[] laplacian = new float[n]; for(int i = 1; i < n-1; i ++){ laplacian[i] = (point[i+1].y-2*point[i].y+point[i-1].y)/(l*l); } float[] du = new float[n]; for(int i = 0; i < n; i ++){ du[i] = D * laplacian[i]; point[i].y += du[i]; } point[0].y = point[1].y; point[n-1].y = point[n-2].y; }
これを実行すると、
noise()でつくった凹凸が均されていくとともに、右側へだんだん温度が伝わっていくのがわかります。
そして最終的には全ての点で加熱している点と同じ温度になるはずです。
二次元で
ではこれを二次元に拡張します。
二次元での拡散方程式は、
$$ \frac{ \partial u}{ \partial t} = D \left( \frac{ \partial^{2} }{ \partial x^{2} } + \frac{ \partial^{2} }{ \partial y^{2} } \right) u $$
3Dで描画して三次元のグラフっぽくしても良いですが、後で反応拡散方程式についてプログラムをつくる上では二次元の絵としてつくった方が見やすいと思うので、そうします。
ある点$(i, j)$での$u$の量を、その点の色の濃さで表現するというような形です。
（テレビの砂嵐みたいな見た目になります。）
流れは一次元のときとほとんど同じです。
まずは、何か板のようなものを考え、それが適当な温度分布を持っているとしておきます。
float[][] u; void setup(){ size(500, 500); u = new float[width][height]; for(int i = 0; i < width; i ++){ for(int j = 0; j < height; j ++){ u[i][j] = random(0, 1); //各点でのuの初期量を設定 } } } void draw(){ background(255); loadPixels(); //画面上の点をピクセル単位で読み込む for(int i = 0; i < width; i ++){ for(int j = 0; j < height; j ++){ int pos = i + j * width; //(i,j)が何番目のピクセルかを計算 pixels[pos] = color(u[i][j] * 255); //ピクセルの色を変える } } updatePixels(); //画面上の点を更新する }
pixels[]は画面上の点の色情報を保存しておく配列で、color型の値で操作できます。
i + j * width
である点$(i, j)$がpixels[]の中で何番目の値なのかを計算できます。
$u$の値は0~1の間にしてあるので、
color(u[i][j] * 255);
で$u$の量によってグレースケールで色をつけられます。
（$u$が多いところでは白く、少ないところでは黒く）
続いて、$\Delta u$の計算をします。
draw()内に、
float[][] laplacianX = new float[width][height]; float[][] laplacianY = new float[width][height]; for(int i = 1; i < width-1; i ++){ for(int j = 1; j < height-1; j ++){ laplacianX[i][j] = u[i+1][j] - 2*u[i][j] + u[i-1][j]; laplacianY[i][j] = u[i][j+1] - 2*u[i][j] + u[i][j-1]; } }
一次元のときにやったことを、x方向とy方向それぞれでやればいいだけです。
最後に、拡散方程式に当てはめて温度の時間変化を計算し、温度を更新する処理を追加します。
draw()内に、
float D = 0.1; float[][] du = new float[width][height]; for(int i = 0; i < width; i ++){ for(int j = 0; j < height; j ++){ du[i][j] = D * (laplacianX[i][j] + laplacianY[i][j]); u[i][j] += du[i][j]; } }
ここまででも一応動きはしますが、これも一次元のときと同様、端の点での温度変化がありませんので、その1つ内側の点と同じだとしておきます。
draw()内に、
for(int i = 0; i < width; i ++){ u[i][0] = u[i][1]; u[i][height-1] = u[i][height-2]; } for(int j = 0; j < height; j ++){ u[0][j] = u[1][j]; u[width-1][j] = u[width-2][j]; }
以上で、二次元に拡張することができました。
実行してみると、
（ファイルサイズが大きくなってしまってアップロードできなかったので静止画で。）
砂嵐がだんだんぼけていって、最終的には同じ色になりました。
反応拡散方程式
初めに挙げた通り、反応拡散方程式は以下のように表される式で、反応項と拡散項からなります。
$$ \begin{align} \frac{ \partial u}{ \partial t} & = f(u, v) + D_u \Delta u \\ \frac{ \partial v}{ \partial t} & = g(u, v) + D_v \Delta v \end{align} $$
ここで挙げた式で、反応項は$f(u, v)$や$g(u, v)$のように一般化した形で書いてありますが、実際に計算するときは具体的に数式を入れて計算します。
その数式の種類によって呼び名が異なり、例えば、
- 拡散方程式
$$ f(u,v) = g(u,v) = 0 $$
- FitzHugh-南雲モデル
$$ \begin{align} f(u,v) & = u - u^{3} - v \\ g(u,v) & = \gamma \left( u - \alpha v - \beta \right) \end{align} $$
- Gray-Scottモデル
$$ \begin{align} f(u,v) & = - u v^{2} + f \left( 1 - u \right) \\ g(u,v) & = u v^{2} - v \left( f + k \right) \end{align} $$
（$f \left( 1 - u \right)$の$f$は定数で、関数を表しているわけではありません。）
今回は、Gray-Scottモデルで描画してみようと思います。
細かいことは省略しますが、Gray-Scottモデルは以下のような仮定の下でつくられたようです。
- 成分$u$、$v$は拡散方程式に従って拡散する
- 化学反応$u+2v \to 3v$が起こる
- $u$はその濃度に応じて自然発生する
- $v$はその濃度に応じて自然消滅する
反応拡散方程式の描画
さて、やっと本題の「チューリングパターン」の描画に移ります。
二次元で
拡散方程式のプログラムをちょっといじれば良いだけなので、直で二次元からやります。
拡散方程式からの変更点は、
- 反応項があること
- 成分が2つであること
反応項は、先に挙げた式をduの計算のところに加えれば良いだけですし、成分が2つになって新たに$v$が加わるなら、$u$と同じことを$v$に対してもやれば良いだけです。
float[][] u, v; float Du, Dv, f, k; void setup(){ size(300, 300); //uのみを一様に広げる u = new float[width][height]; v = new float[width][height]; for(int i = 0; i < width; i ++){ for(int j = 0; j < height; j ++){ u[i][j] = 1; v[i][j] = 0; } } //真ん中に20*20の大きさでvをまく int x = ceil(width/2); int y = ceil(height/2); for(int i = x-10; i < x+10; i ++){ for(int j = y-10; j < y+10; j ++){ u[i][j] = 1; v[i][j] = 1; } } Du = 0.2; Dv = 0.1; f = 0.046; k = 0.063; } void draw(){ background(255); loadPixels(); for(int i = 0; i < width; i ++){ for(int j = 0; j < height; j ++){ int pos = i + j * width; pixels[pos] = color(u[i][j] * 255); } } updatePixels(); float[][] laplacianXu = new float[width][height]; float[][] laplacianYu = new float[width][height]; float[][] laplacianXv = new float[width][height]; float[][] laplacianYv = new float[width][height]; for(int i = 1; i < width-1; i ++){ for(int j = 1; j < height-1; j ++){ laplacianXu[i][j] = u[i+1][j] - 2*u[i][j] + u[i-1][j]; laplacianYu[i][j] = u[i][j+1] - 2*u[i][j] + u[i][j-1]; laplacianXv[i][j] = v[i+1][j] - 2*v[i][j] + v[i-1][j]; laplacianYv[i][j] = v[i][j+1] - 2*v[i][j] + v[i][j-1]; } } float[][] du = new float[width][height]; float[][] dv = new float[width][height]; for(int i = 0; i < width; i ++){ for(int j = 0; j < height; j ++){ //反応項を追加 du[i][j] = Du * (laplacianXu[i][j] + laplacianYu[i][j]) - u[i][j]*v[i][j]*v[i][j] + f*(1-u[i][j]); dv[i][j] = Dv * (laplacianXv[i][j] + laplacianYv[i][j]) + u[i][j]*v[i][j]*v[i][j] - (k+f)*v[i][j]; u[i][j] += du[i][j]; v[i][j] += dv[i][j]; } } for(int i = 0; i < width; i ++){ u[i][0] = u[i][1]; u[i][height-1] = u[i][height-2]; v[i][0] = v[i][1]; v[i][height-1] = v[i][height-2]; } for(int j = 0; j < height; j ++){ u[0][j] = u[1][j]; u[width-1][j] = u[width-2][j]; v[0][j] = v[1][j]; v[width-1][j] = v[width-2][j]; } }
これまでは、初状態では砂嵐状に分布させていましたが、チューリングパターンにおいては、成分$u$のみを全体に一様に広げておいて、ある区画に成分$v$を加えておくという形にします。
機械的に追加するだけなので簡単だと思います。
実行すると、
不思議な模様を描きながら広がっていく様子が見えます。
パラメータはDu、Dv、f、kなので、これらを変えることで様々な模様を描くことができます。
0.001変えるだけでも結果はかなり変わるので、無限通りの組み合わせがあるといっても良いくらいです。
また、最初に成分$v$をまいておく場所を複数にすることでも模様が描け、それが生き物の模様に似ているね、というのがチューリングパターンということになります。
（こんな模様の魚とかいた気がしませんか…？）
「畳み込み」について
ここはおまけ的なところになります。
チューリングパターンを描くにあたっていろいろとウェブページを見ていましたが、ラプラシアンの計算時に本ブログでやったような、微分の定義から直接計算する、という方法をとっているものが少なかったように感じました。
それらは代わりに「畳み込み」という方法を使っていました。
初めて聞いた単語なので解釈が間違っているかもしれませんが、畳み込み自体は、1ステップ後の値の計算をある範囲の値それぞれに重みをつけて足して行う、というような処理だと捉えています。
例えば、「参考」の項の1番目にあるサイトではラプラシアンの計算を以下のようにしています。
float laplaceA = 0; laplaceA += a*-1; laplaceA += prev[i+1][j].a*0.2; laplaceA += prev[i-1][j].a*0.2; laplaceA += prev[i][j+1].a*0.2; laplaceA += prev[i][j-1].a*0.2; laplaceA += prev[i-1][j-1].a*0.05; laplaceA += prev[i+1][j-1].a*0.05; laplaceA += prev[i-1][j+1].a*0.05; laplaceA += prev[i+1][j+1].a*0.05;
図的に表すと、
0.2とか0.05とかいう値がどこからきたのかわかりませんが、やっていることは似ていると思うので、単純化することで処理が早くなるとか、何か利点があるのだと思います。
自分は今のところ、微分の定義から計算する方法でうまくいっている（と感じている）ことと、理解もしやすいということから、畳み込みを無理に使うことはないと思います。
まとめ
かなり長くなってしまいましたが、以上でチューリングパターンの描画ができました。
Gray-Scottモデルであれば、パラメータを少し変えるだけでもかなり違う模様になるので、ぜひ自分だけの模様を作ってみてください。
FitzHugh-南雲モデルで描画してみても面白いと思います。
