logo

파동 방정식의 수치적 풀이: k-space method 📂数値解析

파동 방정식의 수치적 풀이: 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)

  1. T Douglas Mast, Laurent P Souriau, D-LD Liu, Makoto Tabei, Adrian I Nachman, and Robert C Waag. A k-space method for large-scale models of wave propagation in tissue. IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 48(2):341–354, 2001. ↩︎

  2. T Douglas Mast, Laurent P Souriau, D-LD Liu, Makoto Tabei, Adrian I Nachman, and Robert C Waag. A k-space method for large-scale models of wave propagation in tissue. IEEE transactions on ltrasonics, ferroelectrics, and frequency control, 48(2):341–354, 2001. ↩︎