Processingでシミュレーション～境界付近で波を消す

はじめに

波が境界に達したとき反射しないで消え去ってほしいことがあります。
これを実現するPerfectly matched layerという手法について記載します。

Perfectly matched layer(PML)

次のPaperが分かりやすかったです。この内容を記載します。

Notes on Perfectly Matched Layers (PMLs)

波動方程式

$u(t,x,y,z)$は波を表します。
$\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z})$です。"$\cdot$"は内積です。
以下が波動方程式です。少し不自然な形式に見えますが2つの方程式に分ける便宜上このようにします。

\nabla \cdot (a \nabla u) = \frac{1}{b} \frac{\partial^2 u}{\partial t^2}

微分の階数を減らす(2階→1階)ために2つの方程式に分解します。
$\vec{v}(t,x,y,z)$です。

\frac{\partial u}{\partial t} = b \nabla \cdot \vec{v} \\
\frac{\partial \vec{v}}{\partial t} = a \nabla u

上記方程式をまとめます。単なる書き換えです。
$\vec{w}$は$u$と$\vec{v}$を並べたベクトルです。$\hat{D}$は4x4の演算子です。

\frac{\partial \vec{w}}{\partial t} = 
\frac{\partial}{\partial t}
\left(
    \begin{array}{r}
      u \\
      \vec{v} \\
    \end{array}
\right) = 
\left(
    \begin{array}{rr}
       & b\nabla \cdot \\
      a\nabla & \\
    \end{array}
\right)
\left(
    \begin{array}{r}
      u \\
      \vec{v} \\
    \end{array}
\right) =
\hat{D} \vec{w} \\

\hat{D} = 
\left(
    \begin{array}{rr}
       & b\nabla \cdot \\
      a\nabla & \\
    \end{array}
\right)

次の仮定を置きます。原文抜粋

• We will assume that the space far from the region of interest is homogeneous.
• We will assume that the space far from the region of interest is linear and time-invariant.

上記仮定のもとで$\vec{w}(\vec{x},t)$は次のような解になるとのこと。plane waveの重ね合わせ表現になっています。
$\omega$は角周波数,$\vec{k}$は波数ベクトルです。$\vec{W}$は振幅です。

\vec{w}(\vec{x},t) = \sum_{\vec{k},\omega}{\vec{W}_{\vec{k},\omega}exp(i(\vec{k} \cdot \vec{x} - \omega t))}

上記例として以下が示されています。式(A)とします。

\vec{W}(y,z) exp(i(kx - \omega t))

解析接続

式(A)はxについて解析関数です。したがって解析接続できます。
このことは微分方程式を満たしながら変数が動く経路を変えられることを意味しています。
また、exp()であることから振幅を指数関数的に減衰させることができます。
これが反射を抑えて波を急激に減衰させられる理由になります。

分かりやすい図が記載されていましたので抜粋します。
$exp(ikx)$のグラフです。振幅1, 時間項は省略されています。

上側はxが実数成分のみの場合です。減衰はみられません。
下側はxの虚数成分を少しずつ増やした場合です。波形が減衰していることが分かります。
このようにして虚数成分を増やすことで波形を減衰させることができます。

複素数への変換

$\tilde{x} = x + if(x)$とします。
$\partial \tilde{x} = (1 + i \frac{df}{dx})\partial x$になります。$\frac{df}{dx}=\frac{\sigma_x(x)}{\omega}$とおきます。
$\sigma_x(x)$は虚数部を増加させる関数です。

微分演算子を次のように置き換えることでPMLにすることができます。変換(B)とします。

\frac{\partial}{\partial x} \to \frac{1}{1+i\frac{\sigma_x(x)}{\omega}} \frac{\partial}{\partial x}

1次元PML

$u(t,x),v(t,x)$について次の式が1次元波動方程式です。式(C)とします。
$u$と$v$の時間依存の部分は$exp(-i\omega t)$とします。

\frac{\partial u}{\partial t} = b \frac{\partial v}{\partial x} = -i\omega u \\
\frac{\partial v}{\partial t} = a \frac{\partial u}{\partial x} = -i\omega v \\

上記式に変換(B)を行います。

b \frac{\partial v}{\partial x} = -i\omega u + \sigma_x u \\
a \frac{\partial u}{\partial x} = -i\omega v + \sigma_x v \\

上記式と式(C)から$-i\omega u$と$-i\omega v$を消去します。$-i \omega$が$\frac{\partial}{\partial t}$になると考えてもよいです。
これがtime-domainの1次元PMLです。

\frac{\partial u}{\partial t} = b \frac{\partial v}{\partial x} - \sigma_x u \\
\frac{\partial v}{\partial t} = a \frac{\partial u}{\partial x} - \sigma_x v \\

差分化・離散化

差分化

\frac{1}{Δt}(u(t+Δt,x)-u(t,x)) = b \frac{1}{Δx}(v(t, x)-v(t, x-Δx)) - σ_x u(t, x) \\
\frac{1}{Δt}(v(t+Δt,x)-v(t,x)) = a \frac{1}{Δx}(u(t+Δt, x+Δx)-u(t+Δt, x)) - σ_x v(t, x) \\

(*)
2番目の式の$u(t+Δt,x+Δx),u(t+Δt,x)$の項で時間・空間差分を1番目の式に対して1つずらしています。
シミュレーション結果が不安定になったためこのようにしました。

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

u(t+1,x) = u(t,x) + b(v(t,x)-v(t,x-1)) - σ_x u(t,x) \\
v(t+1,x) = v(t,x) + a(u(t+1,x+1)-u(t+1,x)) - σ_x v(t,x) \\

t=0のとき

u(1,x) = u(0,x) + b(v(0,x)-v(0,x-1)) - σ_x u(0,x) \\
v(1,x) = v(0,x) + a(u(1,x+1)-u(1,x)) - σ_x v(0,x) \\

初期条件
$u(0,x),v(0,x)$

境界条件
$u(t,0),u(t,len)$

減衰条件
$σ_x(x)$

Processingで可視化 PML版1次元波動方程式

u1,v1は現在の分布、u0,v0は1時刻前の分布です。
左端のu0をsin関数で振っています。
u1,v1の右端数点は0固定にしています。float範囲で収束しきれない値が反射するのを強制的に抑えています。
右側半分が減衰領域になります。

float a = 0.7;
float b = 0.7;
float[] u0;
float[] v0;
float[] u1;
float[] v1;
float[] s; // sigma
int step = 100;
float w = 0;
int time = 0;

void setup() {
  size(600, 300);
  init();
  //noSmooth();
  frameRate(20);
  background(0);
  fill(0, 0, 0, 0);
  stroke(255, 255, 255);
}

void draw() {
  background(0);
  for (int i=0; i<step; i++) {
    float x1,x2,y1,y2;
    x1 = w*i;
    y1 = height/2+u1[i];
    ellipse(x1, y1, 2, 2);
    if (i < step-1) {
      y2 = height/2+u1[i+1];
      x2 = w*(i+1);
      line(x1, y1, x2, y2);
    }
  }
  for (int i=0; i<1; i++) // speed up
    updateTime();
}

void init()
{
  w = width/(step-1);
  u0 = new float[step];
  v0 = new float[step];
  u1 = new float[step];
  v1 = new float[step];
  s = new float[step];
  int llen = step/2; // length of layer
  float inc_step = 0.5/llen;
  for (int i=step-llen; i<step; i++) {
    s[i] = inc_step*(i-(step-llen));
  }
}

void updateTime()
{
  int force_zero_len = 5;
  u1[0] = 100*sin(2*PI/20*time);

  for (int i=1; i<step-1; i++) {
    u1[i]  = u0[i] + b*(v0[i] - v0[i-1]) - s[i]*u0[i];
  }
  for (int i=0; i<step-1; i++) {
    v1[i]  = v0[i] + a*(u1[i+1] - u1[i]) - s[i]*v0[i];
  }
  for (int i=0; i<step-force_zero_len; i++) {
    u0[i] = u1[i];
    v0[i] = v1[i];
  }

  time++;
}

2次元PML

$u(t,x,y),v_x(t,x,y),v_y(t,x,y)$について次の式が2次元波動方程式です。
$u$と$v_x,v_y$の時間依存の部分は$exp(−iωt)$とします。

\frac{\partial u}{\partial t} = b \nabla \cdot v = b \frac{\partial v_x}{\partial  x} + b \frac{\partial v_y}{\partial  y} = -i \omega u \\
\frac{\partial v_x}{\partial t} = a \frac{\partial u}{\partial x} = -i \omega v_x \\
\frac{\partial v_y}{\partial t} = a \frac{\partial u}{\partial y} = -i \omega v_y \\

上記式に変換(B)を行います。y方向にも変換します。
$\partial x \to (1+i \frac{\sigma_x}{\omega})\partial x$
$\partial y \to (1+i \frac{\sigma_y}{\omega})\partial y$

b \frac{\partial v_x}{\partial  x}(1+i\frac{\sigma_y}{\omega}) + b \frac{\partial v_y}{\partial  y}(1+i\frac{\sigma_x}{\omega}) =
 -i \omega u + (\sigma_x + \sigma_y) u + i\frac{\sigma_x \sigma_y}{\omega} u\\
a \frac{\partial u}{\partial x} = -i \omega v_x + \sigma_x v_x\\
a \frac{\partial u}{\partial y} = -i \omega v_y + \sigma_y v_y\\

2,3番目の式を$-i\omega$を$\frac{\partial}{\partial t}$に置き換えます。

\frac{\partial v_x}{\partial t} = a \frac{\partial u}{\partial x} - \sigma_x u \\
\frac{\partial v_y}{\partial t} = a \frac{\partial u}{\partial y} - \sigma_y u \\

1番目の式をtime-domainに書き換えるため、$\frac{i}{\omega}$を時間積分と考えます。
任意関数$\psi_x(t,x,y)$を導入します。1番目の式の$i\frac{\sigma_x}{\omega}$の項を次のようにします。

\psi_x = \frac{ib\sigma_x}{\omega} \frac{\partial v_y}{\partial y}　\\
-i \omega \psi_x = b\sigma_x \frac{\partial v_y}{\partial y}　\\
\frac{\partial \psi_x}{\partial t} = b\sigma_x \frac{\partial v_y}{\partial y}　\\

1番目の式の$i\frac{\sigma_y}{\omega}$の項についても関数$\psi_y(t,x,y)$を導入します。

\psi_y = \frac{ib\sigma_y}{\omega} \frac{\partial v_x}{\partial x}　\\
-i \omega \psi_y = b\sigma_y \frac{\partial v_x}{\partial x}　\\
\frac{\partial \psi_y}{\partial t} = b\sigma_y \frac{\partial v_x}{\partial x}　\\

1番目の式の$i\frac{\sigma_x \sigma_y}{\omega}$の項についても関数$\psi_{xy}(t,x,y)$を導入します。

\psi_{xy} = \frac{i\sigma_x \sigma_y}{\omega} u　\\
-i \omega \psi_{xy} = \sigma_x \sigma_y u　\\
\frac{\partial \psi_{xy}}{\partial t} = \sigma_x \sigma_y u　\\

以上をまとめると2次元のtime-domain PMLは次の式になります。

\frac{\partial u}{\partial t} = b \nabla \cdot v - (\sigma_x + \sigma_y)u + \psi_x + \psi_y - \psi_{xy} \\
\frac{\partial v_x}{\partial t} = a \frac{\partial u}{\partial x} - \sigma_x v_x \\
\frac{\partial v_y}{\partial t} = a \frac{\partial u}{\partial y} - \sigma_y v_y \\
\frac{\partial \psi_x}{\partial t} = b\sigma_x \frac{\partial v_y}{\partial y}　\\
\frac{\partial \psi_y}{\partial t} = b\sigma_y \frac{\partial v_x}{\partial x}　\\
\frac{\partial \psi_{xy}}{\partial t} = \sigma_x \sigma_y u　\\

差分化・離散化

差分化

\frac{1}{Δt}(u(t+Δt,x,y)-u(t,x,y)) = \\ 
b(\frac{1}{Δx}(v_x(t,x,y)-v_x(t,x-Δx,y) + \frac{1}{Δy}(v_y(t,x,y)-v_y(t,x,y-Δy))) - (\sigma_x + \sigma_y) u(t,x,y) + \psi_x(t,x,y) + \psi_y(t,x,y) - \psi_{xy}(t,x,y) \\
\frac{1}{Δt}(v_x(t+Δt,x,y)-v_x(t,x,y)) = a(\frac{1}{Δx}(u(t+Δt,x+Δx,y)-u(t+Δt,x,y))) - \sigma_x v_x(t,x,y) \\
\frac{1}{Δt}(v_y(t+Δt,x,y)-v_y(t,x,y)) = a(\frac{1}{Δy}(u(t+Δt,x,y+Δy)-u(t+Δt,x,y))) - \sigma_y v_x(t,x,y) \\
\frac{1}{Δt}(\psi_x(t+Δt,x,y)-\psi_x(t,x,y)) = b \sigma_x(\frac{1}{Δy}(v_y(t+Δt,x,y+Δy)-v_y(t+Δt,x,y))) \\
\frac{1}{Δt}(\psi_y(t+Δt,x,y)-\psi_y(t,x,y)) = b \sigma_y(\frac{1}{Δx}(v_x(t+Δt,x+Δx,y)-v_x(t+Δt,x,y))) \\
\frac{1}{Δt}(\psi_{xy}(t+Δt,x,y)-\psi_{xy}(t,x,y)) = \sigma_x \sigma_y u(t,x,y) \\

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

u(t+1,x,y)-u(t,x,y) = b(v_x(t,x,y)-v_x(t,x-1,y) + v_y(t,x,y)-v_y(t,x,y-1)) - (\sigma_x + \sigma_y) u(t,x,y) + \psi_x(t,x,y) + \psi_y(t,x,y) - \psi_{xy}(t,x,y) \\
v_x(t+1,x,y)-v_x(t,x,y) = a(u(t+1,x+1,y)-u(t+1,x,y)) - \sigma_x v_x(t,x,y) \\
v_y(t+1,x,y)-v_y(t,x,y) = a(u(t+1,x,y+1)-u(t+1,x,y)) - \sigma_y v_y(t,x,y) \\
\psi_x(t+1,x,y)-\psi_x(t,x,y) = b \sigma_x(v_y(t+1,x,y+1)-v_y(t+1,x,y)) \\
\psi_y(t+1,x,y)-\psi_y(t,x,y) = b \sigma_y(v_x(t+1,x+1,y)-v_x(t+1,x,y)) \\
\psi_{xy}(t+1,x,y)-\psi_{xy}(t,x,y) = \sigma_x \sigma_y u(t,x,y) \\

t=0のとき

u(1,x,y) = u(0,x,y) + b(v_x(0,x,y)-v_x(0,x-1,y) + v_y(0,x,y)-v_y(0,x,y-1)) - (\sigma_x + \sigma_y) u(0,x,y) + \psi_x(0,x,y) + \psi_y(0,x,y) - \psi_{xy}(0,x,y) \\
v_x(1,x,y) = v_x(0,x,y) + a(u(1,x+1,y)-u(1,x,y)) - \sigma_x v_x(0,x,y) \\
v_y(1,x,y) = v_y(0,x,y) + a(u(1,x,y+1)-u(1,x,y)) - \sigma_y v_y(0,x,y) \\
\psi_x(1,x,y) = \psi_x(0,x,y) + b \sigma_x(v_y(1,x,y+1)-v_y(1,x,y)) \\
\psi_y(1,x,y) = \psi_y(0,x,y) + b \sigma_y(v_x(1,x+1,y)-v_x(1,x,y)) \\
\psi_{xy}(1,x,y) = \psi_{xy}(0,x,y) + \sigma_x \sigma_y u(0,x,y) \\

初期条件
$u(0,x,y),v_x(0,x,y),v_y(0,x,y),\psi_x(0,x,y),\psi_y(0,x,y),\psi_{xy}(0,x,y)$

境界条件
$u(t,x,0),u(t,x,len),u(t,0,y),u(t,len,y)$

減衰条件
$\sigma_x,\sigma_y$

Processingで可視化 PML版2次元波動方程式

float a = 0.5;
float b = 0.5;
float[][][] u;
float[][][] vx;
float[][][] vy;
float[][][] px; // psi
float[][][] py;
float[][][] pxy;
float[] sx; // sigma
float[] sy;
PImage img;
int step = 101;
int time = 0;

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[j*step+i] = color(127-u[i][j][1], 100, 100);
  img.updatePixels();
  image(img, 0, 0, width, height);
  for (int i=0; i<1; i++) // speed up
    updateTime();
}

void init()
{
  u   = new float[step][step][2];
  vx  = new float[step][step][2];
  vy  = new float[step][step][2];
  px  = new float[step][step][2];
  py  = new float[step][step][2];
  pxy = new float[step][step][2];
  sx  = new float[step];
  sy  = new float[step];
  int llen = 10; // length of layer
  float gain = 0.5;
  float inc_step = gain/llen;

  for (int i=step-llen; i<step; i++) {
      sx[i] = sy[i] = inc_step*(i-(step-llen));
  }
  for (int i=0; i<llen; i++) {
      sx[i] = sy[i] = gain - inc_step*i;
  }

  img = createImage(step, step, RGB);
  colorMode(HSB, 255*100/68, 100, 100);
}

void updateTime()
{
  u[step/2][step/2][0] = 3000*sin(2.0*PI/25.0*(float)time);

  /* du/dt */
  for (int i=1; i<step; i++) {
    for (int j=1; j<step; j++) {
      u[i][j][1] = u[i][j][0] + b*(vx[i][j][0]-vx[i-1][j][0]+vy[i][j][0]-vy[i][j-1][0])
        -(sx[i]+sy[j])*u[i][j][0]+px[i][j][0]+py[i][j][0]-pxy[i][j][0];
    }
  }
  /* dvx/dt */
  for (int i=0; i<step-1; i++) {
    for (int j=0; j<step; j++) {
      vx[i][j][1] = vx[i][j][0] + a*(u[i+1][j][1]-u[i][j][1])-sx[i]*vx[i][j][0];
    }
  }
  /* dvy/dt */
  for (int i=0; i<step; i++) {
    for (int j=0; j<step-1; j++) {
      vy[i][j][1] = vy[i][j][0] + a*(u[i][j+1][1]-u[i][j][1])-sy[j]*vy[i][j][0];
    }
  }
  /* dpx/dt */
  for (int i=0; i<step; i++) {
    for (int j=1; j<step; j++) {
      px[i][j][1] = px[i][j][0] + b*sx[i]*(vy[i][j][0]-vy[i][j-1][0]);
    }
  }
  /* dpy/dt */
  for (int i=1; i<step; i++) {
    for (int j=0; j<step; j++) {
     py[i][j][1] = py[i][j][0] + b*sy[j]*(vx[i][j][0]-vx[i-1][j][0]);
    }
  }
  /* dpxy/dt */
  for (int i=0; i<step; i++) {
    for (int j=0; j<step; j++) {
      pxy[i][j][1] = pxy[i][j][0] + sx[i]*sy[j]*u[i][j][0];
    }
  }

  for (int i=0; i<step; i++) {
    for (int j=0; j<step; j++) {
         u[i][j][0] =   u[i][j][1];
        vx[i][j][0] =  vx[i][j][1];
        vy[i][j][0] =  vy[i][j][1];
        px[i][j][0] =  px[i][j][1];
        py[i][j][0] =  py[i][j][1];
       pxy[i][j][0] = pxy[i][j][1];
    }
  }

  time++;
}