反応拡散系

2020年11月1日 (日) 09:49時点における版

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

2変数の連立偏微分方程式を考える。

解説

• https://www.karlsims.com/rd.html
• https://mrob.com/pub/comp/xmorphia/ogl/index.html

シミュレーション

• http://pmneila.github.io/jsexp/grayscott/

ソース

• https://github.com/MStrandh/gray_scott_reaction_diffusion
• https://github.com/Nekodigi/Reaction-Diffusion-Algorithm

作品例

• 芸術学部フェスタ2020「反応と拡散」

int M = 640;
int N = 480;

//System parameters
double diffU;
double diffV;
double paramF;
double paramK;

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];

//int[][] offset = new int[N][2];

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();
   colorMode(HSB,1.0);
   
   //Set default parameters;
   diffU = 0.16;
   diffV = 0.08;
   paramF = 0.035;
   paramK = 0.06;
   
   rndInitCondition = true;
   
   //Populate U and V with initial data
   generateInitialState();
   
}

void timestep(double F, double K, double diffU, double diffV) {
      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 + F*(1 - u);
          dV[i][j] = diffV*lapV + uvv - (K+F)*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 k = 0; k < 10; k++) {
       timestep(paramF, paramK, diffU, diffV);
    }
   
    // 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 ));
//        set(i, j, color(0.74, 0.87, (float)(1-U[i][j])));
//        set(i, j, color(0.5, 0.6, (float)(V[i][j])));
      }
    }
}


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();
  }
}

参考

スケーラブルアート論