反応拡散系

Gray-Scott Reaction and Diffusion（反応拡散系）モデル

参照資料

github

概要

物質は空間内で濃度の濃いところから薄いところへ拡散し、最終的に均一になります。２つの物質が反応しながら拡散するとき、物質の濃淡の波ができ、生物に見られる多種多様な模様が生成されます。これはチューリングパターンとして知られていて、シマウマの縞模様やヒョウの斑模様などがよい例です。

反応と拡散のアルゴリズム

物質A：一定の割合で発生する(Feed)
物質B：一定の割合で消滅する(Kill)
1つのAと2つのBとが反応するとAがBに変化する。

反応拡散方程式

2つの物質AとBが反応しながら、空間に拡散していく様子を反応拡散方程式で表わす。
Laplacian(ラプラシアン) = (上下左右の4方向から流入する量) - (4方向に流出する量)

描画

物質A,Bの濃淡を模様として描画する。
物質Aの発生率と物質Bの消滅率を微妙に変化させるだけで、現れる模様は大きく変化する。

ソース例

int M = 400;
int N = 400;

//System parameters
double diffU; // U(A) diffuse rate
double diffV; // V(B) diffuse rate
double paramF; // Feed rate
double paramK; // Kill rate

boolean rndInitCondition;

double[][] U = new double[M][N];
double[][] V = new double[M][N];

double[][] dU = new double[M][N];
double[][] dV = new double[M][N];

void settings() {
  size(M,N);
}

void generateInitialState() {
    for (int i = 0; i < M; i++) {
      for (int j = 0; j < N; j++) { 
        U[i][j] = 1.0;
        V[i][j] = 0.0;
      }
    }
    
    if (rndInitCondition) {
        for (int i = M/3; i < 2*M/3; i++) {
            for (int j = N/3; j < 2*N/3; j++) {     
               U[i][j] = 0.5*(1 + random(-1, 1));
               V[i][j] = 0.25*( 1 + random(-1, 1));
          }
        }
    } else {
      for (int i = M/3; i < 2*M/3; i++) {
            for (int j = N/3; j < 2*N/3; j++) {     
               U[i][j] = 0.5;
               V[i][j] = 0.25;
          }
        }
    }
}

void setup() {
   frameRate(48);
//   smooth();  Processing 4以降ではエラーになる
   colorMode(HSB,1.0);
   
   //Set default parameters;
   diffU = 0.16;
   diffV = 0.08;
   paramF = 0.0545;
   paramK = 0.062;
   
   rndInitCondition = true;
   
   //Populate U and V with initial data
   generateInitialState();
}

void timestep() {
      for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
//          int p = i + j*N;
          
          double u = U[i][j];
          double v = V[i][j];
          double uvv = u*v*v;
          
          int left = (i-1+M) % M;
          int right = (i+1) % M;
          int up = (j-1+N) % N;
          int down = (j+1) % N;
       
          double lapU = (U[left][j] + U[right][j] + U[i][up] + U[i][down] - 4*u);
          double lapV = (V[left][j] + V[right][j] + V[i][up] + V[i][down] - 4*v);
          
          dU[i][j] = diffU*lapU  - uvv + paramF*(1 - u);
          dV[i][j] = diffV*lapV + uvv - (paramK+paramF)*v;
        }
      }
    for (int i= 0; i < M; i++) {
      for (int j = 0; j < N; j++){
          U[i][j] += dU[i][j];
          V[i][j] += dV[i][j];
      }
    }
}

void draw(){ 
    for (int i = 0; i < 20; i++) {
       timestep();
    }
   
    // Draw points
    for (int i = 0; i < M; i++) {
      for (int j = 0; j < N; j++) {
        set(i, j, color((float)(1-U[i][j]),0.9, 0.5 ));
      }
    }
}

キーボードでパラメータを変更

void keyPressed() {
  switch (key) {
    case '1':
          diffU = 0.16;
          diffV = 0.08;
          paramF = 0.035;
          paramK = 0.06;
          generateInitialState();
          break;
    case '2':
          diffU = 0.16;
          diffV = 0.08;
          paramF = 0.042;
          paramK = 0.065;
          generateInitialState();
          break;
    case '3':
          diffU = 0.18;
          diffV = 0.13;
          paramF = 0.025;
          paramK = 0.056;
          generateInitialState();
          break;
    case '4':
          diffU = 0.18;
          diffV = 0.09;
          paramF = 0.02;
          paramK = 0.056;
          generateInitialState();
          break;
    case '5':
          diffU = 0.14;
          diffV = 0.06;
          paramF = 0.035;
          paramK = 0.065;
          generateInitialState();
          break;
    case '6':
          diffU = 0.19;
          diffV = 0.09;
          paramF = 0.062;
          paramK = 0.062;
          generateInitialState();
          break;
     case '7':
          diffU = 0.16;
          diffV = 0.08;
          paramF = 0.05;
          paramK = 0.065;
          generateInitialState();
          break;
    case 'r':
          rndInitCondition = true;
          generateInitialState();
          break;
    case 'n':
          rndInitCondition = false;
          generateInitialState();
  }
}