📂数値解析

파동 방정식의 수치적 풀이: 유한차분법(FDM)

メソッド

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

$$ \dfrac{\partial^{2} u}{\partial t^{2}} = c^{2} \dfrac{\partial^{2} u}{\partial x^{2}}, \qquad 0 \le x \le 1, \quad t \ge 0 \tag{1} $$

我々の目的は有限な点で上記のソリューションを近似することである。時空間ドメインを次のように分割しよう。

$$ \left\{ (\ell \Delta x, n\Delta t) : \ell=0,1,\dots,d+1,\ n\ge 0 \right\}\quad \text{ where } \Delta x = \frac{1}{d+1} $$

$u(\ell \Delta x, n\Delta t)$の近似を$u_{\ell}^{n}$と表記しよう。上付き文字${}^{n}$が累乗ではなく時間のインデックスであることに注意すること。$(1)$の右辺と左辺はそれぞれ次のように近似される。

$$ \dfrac{\partial^{2} u(x,t)}{\partial x^{2}} = \dfrac{1}{(\Delta x)^{2}} \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] + \mathcal{O}\left( (\Delta x)^{2} \right) $$

$$ \dfrac{\partial^{2} u(x,t)}{\partial t^{2}} = \dfrac{1}{(\Delta t)^{2}} \big[ u(x, t-\Delta t) - 2u(x,t) + u(x, t+\Delta t) \big] + \mathcal{O}\left( (\Delta t)^{2} \right) $$

ここで$\mu = \dfrac{\Delta t}{\Delta x}$とする。二つの式を$(1)$に代入して整理すると次を得る。

$$ u(x, t-\Delta t) - 2u(x,t) + u(x, t+\Delta t) \approx \mu^{2} \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] $$

$$ \implies u(x, t+\Delta t) \approx 2u(x,t) - u(x, t-\Delta t) + \mu^{2} \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] $$

$$ \implies u(x, t+\Delta t) \approx 2(1 - \mu^{2})u(x,t) - u(x, t-\Delta t) + \mu^{2} \big[ u(x-\Delta x, t) + u(x+\Delta x, t) \big] $$

$$ \implies u_{\ell}^{n+1} = \mu^{2} u_{\ell+1}^{n} + 2(1 - \mu^{2})u_{\ell}^{n} + \mu^{2} u_{\ell-1}^{n} - u_{\ell}^{n-1},\quad \ell=1,\dots,d,\ n=0,1,\dots \tag{2} $$

ここで空間間隔と時間間隔の比率$\mu = \dfrac{\Delta t}{\Delta x}$はクーラン数^{Courant number}と呼ばれる重要な定数であり、$(2)$の収束可否を決定する。

説明

上記のメソッドは$\mu \le 1$のとき収束する。

コード

境界条件を$0$に設定して計算すると、次のような結果が得られる。

using Plots

x = LinRange(-1, 1, 300)
Δx = x[2] - x[1]

t = LinRange(0, 4, 600)
Δt = t[2] - t[1]

μ = Δt / Δx

U = zeros(length(x), length(t))

U[:, 1] = exp.(-30 * x.^2)
U[:, 2] = exp.(-30 * x.^2)

for i ∈ 3:length(t)
    U[2:end-1, i] = (μ^2)*U[3:end, i-1] + 2(1-μ^2)*U[2:end-1, i-1] + (μ^2)*U[1:end-2, i-1] - U[2:end-1, i-2]
end

heatmap(U)

anim = @animate for i ∈ 1:length(t)
    plot(x, U[:, i], xlims=(-1, 1), ylims=(-2,2), legend=false)
end

gif(anim, fps=50)