Erlangで数値解析(ただし、1次元)

数値解析のプログラムを書いてて、Erlangで書けるもんなんだろうか、と思ってやってみた。題材は、電磁波の伝搬問題を解くFDTD(Finite Difference Time Domain)法。めんどくさいので、1次元で。

FDTD法とは

Maxwellの方程式を差分法で解く手法。


(上記は真空場合の式)

これを何も考えずに差分化して、更新式の形にすると(1次元のTEM波。波はX軸方向に進みます)、

となります。要は、ただの1次風上差分方程式です。

Erlangで書く

さて、上記の更新式をErlangで書くと、

calc_ey({Eps, Dt, Dx}, Ey, Hz) ->
    Diff_Hz = create_diff_list(Hz),   % Hzの差分リストを作成
    Cut_Ey = cut_both_end_element(Ey),   % Eyの両端をカット
    [E - Dt / (Eps * Dx) * DH || {E, DH} <- lists:zip(Cut_Ey, Diff_Hz)].   % Ey更新式

calc_hz({Mu, Dt, Dx}, Hz, Ey) ->
    Diff_Ey = create_diff_list(Ey),   % Eyの差分リストを作成
    [H - Dt / (Mu * Dx) * DE || {H, DE} <- lists:zip(Hz, Diff_Ey)].  % Hz更新式

と言う形になります。書き方は色々あると思いますが、いまいちErlangに慣れていないので、リスト内包を使いました。

ちなみにCで書くと、

void calc_ey(double* ey, double* hz, double* dx, double dt)
{
  int i;
  for (i = 1; i < NX; i++) {
    ey[i] = ey[i] - dt / (PERMITTIVITY * dx[i]) * (hz[i] - hz[i-1]);
  }
}

void calc_hz(double* hz, double* ey, double* dx, double dt)
{
  int i;
    for (i = 0; i < NX; i++) {
      hz[i] = hz[i] - dt / (PERMEABILITY * dx[i]) * (ey[i+1] - ey[i]);
    }
}

こうなります。

複雑な式ではないので、どちらで書いても同じような形にはなるのですが、Erlangの場合、

• forループが無い
• forループがないぶん、式を直接プログラム化したような感じになる。

と言うような特徴があると思います。

Erlangには全く慣れていないので、このような書き方しかできませんが、もっとエレガントな書き方があると思います。

今後

値の更新の箇所を並列化したいね〜、などと考えています。また時間ができたときにでもw

最後にソース

• Erlang版
• C版