Processingでシミュレーション～波動方程式

はじめに

Processingを使って波動方程式をシミュレーションします。

波動方程式

波動方程式は波の現象を表す方程式です。
次の微分方程式を満たす$u(t,x,y,z)$が波動を表します。
$\Delta_l$はラプラシアンです。$\Delta_l = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$です。

\frac{1}{s^2} \frac{\partial^2u}{\partial t^2} = \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} + \frac{\partial^2u}{\partial z^2} = \Delta_l u

差分化・離散化

2次元の波動方程式を扱います。
差分化

\frac{1}{s^2Δt^2}(u(t,x,y)-2u(t-Δt,x,y)+u(t-2Δt,x,y)) \\ 
= \frac{1}{Δx^2}(u(t-Δt,x-Δx,y)-2u(t-Δt,x,y)+u(t-Δt,x+Δx,y)) + \\
  \frac{1}{Δy^2}(u(t-Δt,x,y-Δy)-2u(t-Δt,x,y)+u(t-Δt,x,y+Δy))

離散化
Δt=1, t=0,1,2,...
Δx=1, x=0,1,2,...,len
Δy=1, y=0,1,2,...,len
とします。

\frac{1}{s^2}(u(t,x,y)-2u(t-1,x,y)+u(t-2,x,y)) \\ 
= (u(t-1,x-1,y)-2u(t-1,x,y)+u(t-1,x+1,y)) + \\
  (u(t-1,x,y-1)-2u(t-1,x,y)+u(t-1,x,y+1))

t=0のとき

u(0,x,y) = 2u(-1,x,y)-u(-2,x,y)) + \\ 
 s^2(u(-1,x-1,y)-2u(-1,x,y)+u(-1,x+1,y) + \\ 
     u(-1,x,y-1)-2u(-1,x,y)+u(-1,x,y+1))

u(-1,x,y),u(-2,x,y)は初期条件です。

u(t,-1,y),u(t,len+1,y),u(t,x,-1),u(t,x,len+1)
u(t,-2,y),u(t,len+2,y),u(t,x,-2),u(t,x,len+2)
は境界条件です。

Processingで可視化

t2は現在の分布、t1は1時刻前の分布、t0は2時刻前の分布です。
境界条件はすべて値0に固定しています。
時刻0でt2の中央に値を持たせています。波紋が広がっていく様子が分かります。

float s2 = 0.2;
float[][] u0;
float[][] u1;
float[][] u2;
PImage img;
int step = 1000;

void setup() {
  size(500, 500);
  init();
  //noSmooth();
  frameRate(30);
}

void draw() {
  for (int i=0; i<step; i++)
    for (int j=0; j<step; j++)
      img.pixels[i*step+j] = color(255-u2[i][j], 100, 100);
  img.updatePixels();
  image(img, 0, 0, width, height);
  for (int i=0; i<4; i++) // speed up
    updateTime();
}

void init()
{ 
  u0 = new float[step][step];
  u1 = new float[step][step];
  u2 = new float[step][step];
  u2[step/2][step/2] = 10000;
  img = createImage(step, step, RGB);
  colorMode(HSB, 255*100/68, 100, 100);
}

void updateTime()
{
  for (int i=0; i<step; i++) {
    for (int j=0; j<step; j++) {
      u0[i][j] = u1[i][j];
      u1[i][j] = u2[i][j];      
    }
  }

  for (int i=0; i<step; i++) {
    u2[0][i] = 0;
    u2[step-1][i] = 0;
    u2[i][0] = 0;
    u2[i][step-1] = 0;
  }
  
  for (int i=1; i<step-1; i++) {
    for (int j=1; j<step-1; j++) {
          u2[i][j] = 2*u1[i][j] - u0[i][j] +
      s2*(u1[i-1][j]-2*u1[i][j]+u1[i+1][j] +
          u1[i][j-1]-2*u1[i][j]+u1[i][j+1]);
    }
  }
}