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