파동 방정식의 수치적 풀이: k-space method

メソッド

以下のような波動方程式が与えられたとする。

$$ \partial_{t}^{2} u(\mathbf{x}, t) = \Delta_{\mathbf{x}} u (\mathbf{x}, t), \qquad (\mathbf{x}, t) \in \mathbb{R}^{n} \times \mathbb{R} $$

両辺に変数$\mathbf{x}$に対するフーリエ変換を取ると次を得る。

$$ \partial_{t}^{2} \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) = \mathcal{F}_{\mathbf{x}}[\Delta u (\cdot, t)]u(\mathbf{k}, t) = -|\mathbf{k}|^{2} \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) $$

また、上記の式の左辺は有限差分法によって次のように近似される。

$$ \partial_{t}^{2} \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) \approx \dfrac{\mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t + \Delta t) - 2\mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) + \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t - \Delta t)}{(\Delta t)^{2}} $$

上記の二つの式から次を得る。

$$ \dfrac{\mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t + \Delta t) - 2\mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) + \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t - \Delta t)}{(\Delta t)^{2}} = -|\mathbf{k}|^{2} \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) $$

$$ \implies \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t + \Delta t) = 2\mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) - \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t - \Delta t) - (\Delta t)^{2}|\mathbf{k}|^{2} \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) $$

両辺にフーリエ逆変換を取ると、

$$ u(\mathbf{x}, t + \Delta t) = 2u(\mathbf{x}, t) - u(\mathbf{x}, t - \Delta t) - \mathcal{F}_{\mathbf{k}}^{-1}\Big[ (\Delta t)^{2}|\mathbf{k}|^{2} \mathcal{F}_{\mathbf{x}}u(\mathbf{k}, t) \Big](\mathbf{x}, t) \tag{1} $$

説明

$k$-spaceメソッドは、上記のようにフーリエ変換を使用して波動の伝播^propagationを計算する方法を指す。フーリエ変換を取ると空間変数$\mathbf{x}$が周波数変数$\mathbf{k}$に変わるため、この名前が付けられている。有限要素法で波動方程式を解く際には空間格子と時間格子の比率が重要だが、$k$-spaceメソッドではそのような考慮が必要ない。

また、$(1)$の最後の項で、$(\Delta t |\mathbf{k}|)^{2}$を$4 \sin^{2} \left( \dfrac{|\mathbf{k}|\Delta t}{2} \right)$に変える方が良い結果をもたらす^{1 2}。

コード

以下はジュリアで作成された$k$-spaceメソッドのコードと実行結果をアニメーションで示したものである。

using Plots
using FFTW

# Initialization
U = zeros(Float64, 256, 256, 256)

# x ∈ [-1, 1]
# t ∈ [0, 1]
Δt = 1/256

for i ∈ 1:256, j ∈ 1:256
    if 10^2 < (i-80-64)^2 + (j-80-64)^2 < 12^2
        U[i, j, 1] = 1
        U[i, j, 2] = 1
    end
    if (i-45-64)^2 + (j-45-64)^2 < 4^2
        U[i, j, 1] = 1
        U[i, j, 2] = 1
    end
end

# make k-space grid
k = LinRange(-128., 128, 256) * π
k1 = k * fill!(similar(k'), 1)
k2 = fill!(similar(k), 1) * k'

K = sqrt.(k1.^2 + k2.^2)
K = 4 * (sin.(Δt * K/2)).^2
K = ifftshift(K)

for i ∈ 3:256
    U[:, :, i] = 2U[:, :, i-1] - U[:, :, i-2] - real(ifft(K .* fft(U[:, :, i-1])))
end

w2b = cgrad([:white, :black])
anim = @animate for i ∈ 1:256
    heatmap(U[:,:,i], c=w2b, clim=(0,1), title="$i")
end

gif(anim)