search
LoginSignup
32

More than 5 years have passed since last update.

posted at

updated at

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

IMAGE ALT TEXT HERE

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
32