FDTD法

簡単な原理

初心者向け講義資料（FDTD法による音波伝搬シミュレーション）に詳しく書いてありますので、そちらを参考にしてください。

matlabでの実装

以下は、一次元のFDTDシミュレーションプログラムです。

見ていただいて分かるとおり、ほとんどが初期条件設定及びアニメーション部分となっています。 シミュレーション条件はプログラムの最初の部分に書いてありますが、両端固定の弦振動となっています。 単純に両端固定の場合はプログラムのようにFDTDの計算部分で境界の計算部分を入れる必要はありません。 そのほかの境界の場合は、次の時刻の音圧（振動）分布を求めた後に境界条件を入れる必要があります。 matlabでは行列演算をするのが大切ですので、このようなプログラムになっています。 一次元ではfor文で回すのとあまり変わりませんが、二次元になると実行時間に大きな差ができます。

なお、二次元のFDTD法によるシミュレーションもほとんど変わりません。 fdtd_2D.mを見ていただければ分かると思います。

clear
%条件設定
nx=0.1;%弦の長さ(m)
nt=0.1;%計算時間(s)
c=100;%音速(m/s)
dx=1e-3;%距離刻み(m)
dt=1e-5;%時間刻み(s)

%定数等決定
A=(c*dt/dx).^2;%差分計算の時に用いる係数
p=zeros(nx/dx+1,nt/dt+1);%音圧p(x,t)
w=nx/20;%ガウス分布のパルス幅（分散）                           
shift=nx/2;%ガウス分布の中心(平均値)
p(2:nx/dx,1)=exp(-((dx:dx:nx-dx) - shift).^2/(2*w^2));%ガウス分布型初期値
p(:,2)=p(:,1);%初期条件から静かに離す

%FDTD計算
for t=2:nt/dt
    p(2:nx/dx,t+1)=A*(p(1:nx/dx-1,t)-2*p(2:nx/dx,t)+p(3:nx/dx+1,t))-p(2:nx/dx,t-1)+2*p(2:nx/dx,t);
    %両端固定である場合は境界条件を入れる必要はない
    %固定境界以外はここで境界条件を入れる。
end

% ムービーの作成および表示
for t = 1:nt/dt+1
    plot(p(:,t));
    axis([1 nx/dx+1 -1 1])
    M(t) = getframe;
end
movie(M)