反応拡散系(3)
前回、Gray-Scott モデルにおける反応拡散方程式を見ました。以下のようなものでした。
vt = dvΔv - u^2v + F(1-v)
ut = duΔu + u^2v - (F+k)u
dv : V の拡散率
du : U の拡散率
F : 原材料 V の外部からの供給率＆中間生成物 U の外部への流出率
k : 中間生成物 U の最終生成物 P への転換率（U の除去率）
これを実際の計算式に落とし込むとどうなるのか、というのが今回の話。
計算を実行する前準備。
この化学反応は2次元平面上で発生し、その2次元平面の各ピクセルごとに V の濃度と U の濃度の初期値を与え、各ピクセルごとに方程式を計算し続けることで、それぞれの物質の濃度に変化が生じる、ということをします。
1ピクセルごとの計算ですが、第2項、第3項はそのまま乗算、加算、減算をすればいいとして、第1項はどうするのか。
これはラプラシアンというもので。近傍4ピクセル（上、下、右、左）の濃度の合計値から、カレントピクセルの濃度の4倍の値を引いた値になるという。
ConvolutionFilter のマトリックス風に表現すると以下のような感じですね。
[
0, 1, 0,
1 ,-4, 1,
0, 1, 0
];
なおラプラシアンは △u や △v の値であり、第1項の計算値は、そのラプラシアンに dv なり du を掛けた値になります。
U と V の初期値、そして vt と ut の値は 0.0 ～ 1.0 の間に収まるようにしました。
実際の Gray-Scott モデルではその範囲を超える値を許容しているのか否か分かりませんが（化学反応のモデルなのでマイナス値は許容していないように思えますが）、パラメータ値によっては結果のビジュアライズがおかしくなるのは、0.0 ～ 1.0 に収めていないからっぽいので、そのようにしました。これはコード中の calc 関数で定義しています。
var len:int = NUM_OF_PIXELS;
for (var i:int = 0; i < len; i++) {
// カレントピクセルの座標
var posX:int = i % IMAGE_WIDTH;
var posY:int = i / IMAGE_WIDTH >> 0;
// 近傍4ピクセルのリスト上の位置を計算
var west:int = posX == 0 ? i : i - 1; // 左
var east:int = posX == IMAGE_WIDTH_FOR_LIST ? i : i + 1; // 右
var north:int = posY == 0 ? i : i - IMAGE_WIDTH; // 上
var south:int = posY == IMAGE_HEIGHT_FOR_LIST ? i : i + IMAGE_WIDTH; // 下
// カレントの濃度値
var currentV:Number = vList_[i];
var currentU:Number = uList_[i];
// 拡散項の計算
var diffusionV:Number = dv_ * (vList_[west] + vList_[east] + vList_[north] + vList_[south] - 4 * currentV);
var diffusionU:Number = du_ * (uList_[west] + uList_[east] + uList_[north] + uList_[south] - 4 * currentU);
// 反応項の計算(1)
var reaction1:Number = currentV * currentU * currentU; // 2U + V -> 3U
// 反応項の計算(2)
var reaction2V:Number = f_ * (1 - currentV); // 原材料 V の外部からの供給
var reaction2U:Number = (f_ + k_) * currentU; // U -> P （U の除去）
// 反応拡散の計算
currentV += (diffusionV - reaction1 + reaction2V); // reaction1 は除去、reaction2V は供給：原材料
currentU += (diffusionU + reaction1 - reaction2U); // reaction1 は供給、reaction2U は除去：中間生成物
if (currentV < 0.0) currentV = 0.0;
if (currentV > 1.0) currentV = 1.0;
if (currentU < 0.0) currentU = 0.0;
if (currentU > 1.0) currentU = 1.0;
vList_[i] = currentV;
uList_[i] = currentU;
}
そして、投稿したコードでは、原材料 V の濃度をビジュアライズしています。V が U と反応して、濃度が薄くなるほど、画面上では黒くなるということです。draw 関数で定義しています。
// 原材料の濃度を色に変換し、リストに格納
var len:int = NUM_OF_PIXELS;
for (var i:int = 0; i < len; i++) {
var gray:int = vList_[i] * 255 >> 0;
viewList_[i] = gray << 16 | gray << 8 | gray;
}
// BitmapData に反映
viewBmd_.lock();
viewBmd_.setVector(RECT, viewList_);
viewBmd_.unlock();
ビューア上でドラッグすると、当該ピクセルと近傍4ピクセルの V の濃度を 1.0、U の濃度を 0.5 + Math.random() * 0.5 にします。それによって、黒い点を描いたようになったり、すでに黒くなったところは白く消えたりします。interference 関数で定義しています。
なお、初期値は V の濃度は 1.0、U の濃度は基本的には 0.0 ですが 0.1％の確率で 0.5 + Math.random() * 0.5 にしています。これは init 関数で定義。
なお、今までの流れから分かるように、原材料 V は中間物質 U と反応して減少していく一方なので、外部から常に供給され続けないと、この反応は継続しません。
その外部からの供給を示すのが方程式の中の F(1-v) という項ですが、ここにロジスティック写像が存在しているのが肝ですね。
Comments
Tell me what you're thinking...
and oh, if you want a pic to show with your comment, go get a gravatar!