畳み込み演算を利用して波動方程式をシミュレーション

  • 16
    Like
  • 0
    Comment
More than 1 year has passed since last update.

波動方程式とは

波動方程式は一般に以下の偏微分方程式で記述されます.

\frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u

ここで, $u$は波の高さ, $c$は波の速さに対応します.
なお, ここでは2次元空間上の波動方程式を取り扱うものとします.

波動方程式の離散化

波動方程式の両辺は2階の微分となっているので, 中心差分に基づいて離散化を行います.

\frac{u^{n+1}_{i, j} - 2u^n_{i, j} + u^{n-1}_{i, j}}{(\Delta t)^2} = c^2 \left( \frac{u^{n}_{i+1, j} - 2u^n_{i, j} + u^{n}_{i-1, j}}{(\Delta x)^2} + \frac{u^{n}_{i, j+1} - 2u^n_{i, j} + u^{n}_{i, j-1}}{(\Delta y)^2} \right)

ここで, $n$は時間ステップ数, $i, j$はそれぞれメッシュの$x$方向, $y$方向の番号, $\Delta t, \Delta x, \Delta y$はそれぞれ離散化する際の時間, $x$方向, $y$方向の細かさです.
この式から, $n+1$ステップ目の$u$を求めるには, $n, n-1$ステップの$u$の値を用いて, $n$ステップの値については$i, j$の周りの4マスを参照して計算を行うことになります.
以下では, 簡単のため$\Delta x = \Delta y = h$とすることにします.

MATLABでどう実装する?

for文が遅いことで有名なMATLABですが, 上の波動方程式を高速にシミュレーションできないでしょうか?
MATLABで高速に計算するためには, あらゆる演算を行列形式で計算することが必要です.
上の計算を行列形式で考えるとどうなるでしょうか?

空間座標と行列の要素の対応

ここでは, ある$i, j$マス目の波の高さの値を, 行列の$j$行$i$列目に対応させるように行列を構成するとうまく行きます.
すなわち,

U^n = \begin{bmatrix}
u^n_{1,1} & u^n_{2,1} & \cdots & u^n_{N_x, 1}\\
u^n_{1,2} & u^n_{2,2} & \cdots & u^n_{N_x, 2}\\
\vdots & \vdots & \ddots & \vdots \\
u^n_{1,N_y} & u^n_{2,N_y} & \cdots & u^n_{N_x, N_y}\\
\end{bmatrix}

のように行列を定義してみましょう. すると, 先ほど離散化した波動方程式が以下のように記述できます.

U^{n+1} = 2U^n - U^{n-1} + c^2\left(U^n \ast K \right)(\Delta t)^2
K = \frac{1}{h^2}\begin{bmatrix}
0 & 1 & 0\\
1 & -4 & 1\\
0 & 1 & 0
\end{bmatrix}

ここで, $\ast$は畳み込み演算を表す記号であり, $K$は畳み込みカーネルとなります.

MATLABによる畳み込み処理

ここで行う畳み込み処理は, 実は画像処理の分野ではよく行われている畳み込みフィルタリングに対応しています.
MATLABで畳み込みフィルタリングを行う関数として, imfilterと呼ばれる関数があります.
この関数を利用すると, 簡単に畳み込み演算が行えるだけでなく, 境界の処理も指定できます.

実際にシミュレーション

waveeqn.m
%% 変数の初期化とfigureのclose
clear
close

%% 定数の指定
nx = 128; % x方向のメッシュ数
ny = nx;  % y方向のメッシュ数
l = 10;   % 実空間での全長
h = l/(nx-1); % 1メッシュあたりの長さ
dt = 0.005; % 時間積分の時間の長さ
c = 2; % 波の速さ

x = [0:h:l]; % x方向のメッシュの座標情報
y = [0:h:l]; % y方向のメッシュの座標情報

F = [0 1 0;
     1 -4 1;
     0 1 0]./(h^2); % 畳み込みカーネル

%% 初期条件
[X, Y] = meshgrid(x, y); % 座標情報テーブルを作る
U = zeros(ny, nx)+0.5; % とりあえず水面は0.5の高さ
% (5, 5)を中心に, 半径sqrt(2)以内の円を1.5の高さに設定
U(((X-5).^2+(Y-5).^2)<2)=1.5;


%% シミュレーション条件を設定
loops = 300; % シミュレーションする時間ステップ数

%% 初期状態を表示
figure(1);
mesh(X, Y, U);

% 以下は表示に関するオプション
axis([0 10 0 10 -2 3]);
caxis([0 1.5]);
cmap(1,:) = [0 0 1];
cmap(2,:) = [0 0.5 1];
colormap(cmap);
text(2, -2,-0.2, 't = 0.000[s]' ,'FontSize', 20);
U1 = U;
U2 = U;
drawnow

%% シミュレーション開始
for i = 1:loops
    % 畳み込みによりn+1ステップ目の波の高さを計算
    % オプションで値を指定すると固定値境界に, replicateを指定すると値のコピー, circularを指定すると周期境界になる.
    U = 2.*U1 - U2  + c^2.*imfilter(U1, F, 'replicate').*dt*dt;
    % U1(U^n), U2(U^{n-1})を更新
    U2 = U1;
    U1 = U;

    % 表示
    mesh(X, Y, U);
    axis([0 10 0 10 -2 3]);
    caxis([0 1.5]);
    colormap(cmap);
    text(2, -2,-0.2, ['t =' num2str(i*dt, '%4.3f') '[s]'] ,'FontSize', 20);
    drawnow

    disp(['Loop...' num2str(i)]);
end

以上のように, 簡単にシミュレーションが可能となります.
いくつか試してみたので, 動画を貼っておきます.