トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS

シミュレーションプログラム/FDTD法 の変更点


*FDTD法 [#s5a18688]

**簡単な原理 [#a7aa9b4d]
初心者向け講義資料([[FDTD法による音波伝搬シミュレーション>シミュレーション/FDTD法による音波伝搬シミュレーション]])にもう少し詳しく書いてありますが、簡単に原理を説明します。
初心者向け講義資料([[FDTD法による音波伝搬シミュレーション>初心者向け講義資料/FDTD法による音波伝搬シミュレーション]])にもう少し詳しく書いてありますが、簡単に原理を説明します。

FDTD法(Finite-difference time-domain method;時間領域差分法)は、名前に差分法とついていることから分かるように、差分近似によって計算をします。
具体的にはどういうことかと言いますと、微分を差分で置き換えるということです。~
&mimetex( \frac{d}{dx}f(x) = \lim_{h \to 0} \frac{f(x+h)-f(x-h)}{h}. );(1)~
式(1)が微分の定義ですが、右辺のリミットを外して単純な差としてしまうのが差分法です。

FDTD法による音波伝搬シミュレーションでは、波動方程式を差分法により近似します。
音速&mimetex(c);で音圧&mimetex(p);の波動方程式は~
&mimetex( \frac{\partial^2 p}{\partial t^2}p = c^2 \nabla^2 p,);(2)~
で与えられますが、簡単のため一次元の波動方程式~
&mimetex( \frac{\partial^2 p}{\partial t^2} = c^2 \frac{\partial^2 p}{\partial x^2},);(3)~
を考え、差分法による近似をします。~
&mimetex( p(x, t+\Delta t) = \frac{c^2 \Delta t^2}{\Delta x^2} \{ p(x+\Delta x,t)+p(x-\Delta x,t)-2p(x,t) \} -p(x,t-\Delta t)+2p(x,t) .);(4)~
&mimetex(x);方向の刻みが&mimetex(\Delta x);で、&mimetex(t);方向の刻みが&mimetex(\Delta t);となっています。
式(4)を見ていただくと、現時点の音波の情報と過去の音波の情報から未来の情報が得られるということが分かります。
従って、初期値として時刻0の音波情報と時刻&mimetex(\Delta t);での音波情報を与えれば音波伝搬シミュレーションができます。


**matlabで実装したFDTD法の音波伝搬シミュレーションプログラム [#f61f179a]
以下は、一次元のFDTDシミュレーションプログラムです。

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

なお、二次元のFDTD法によるシミュレーションもほとんど変わりません。
[[fdtd_2D.m:http://www.tohlab.net/cont/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)