파동 방정식의 수치적 풀이: 유한차분법(FDM)
📂数値解析파동 방정식의 수치적 풀이: 유한차분법(FDM)
メソッド
下記のような一次元波動方程式が与えられたとする。
∂t2∂2u=c2∂x2∂2u,0≤x≤1,t≥0(1)
我々の目的は有限な点で上記のソリューションを近似することである。時空間ドメインを次のように分割しよう。
{(ℓΔx,nΔt):ℓ=0,1,…,d+1, n≥0} where Δx=d+11
u(ℓΔx,nΔt)の近似をuℓnと表記しよう。上付き文字nが累乗ではなく時間のインデックスであることに注意すること。(1)の右辺と左辺はそれぞれ次のように近似される。
∂x2∂2u(x,t)=(Δx)21[u(x−Δx,t)−2u(x,t)+u(x+Δx,t)]+O((Δx)2)
∂t2∂2u(x,t)=(Δt)21[u(x,t−Δt)−2u(x,t)+u(x,t+Δt)]+O((Δt)2)
ここでμ=ΔxΔtとする。二つの式を(1)に代入して整理すると次を得る。
u(x,t−Δt)−2u(x,t)+u(x,t+Δt)≈μ2[u(x−Δx,t)−2u(x,t)+u(x+Δx,t)]
⟹u(x,t+Δt)≈2u(x,t)−u(x,t−Δt)+μ2[u(x−Δx,t)−2u(x,t)+u(x+Δx,t)]
⟹u(x,t+Δt)≈2(1−μ2)u(x,t)−u(x,t−Δt)+μ2[u(x−Δx,t)+u(x+Δx,t)]
⟹uℓn+1=μ2uℓ+1n+2(1−μ2)uℓn+μ2uℓ−1n−uℓn−1,ℓ=1,…,d, n=0,1,…(2)
ここで空間間隔と時間間隔の比率μ=ΔxΔtはクーラン数Courant numberと呼ばれる重要な定数であり、(2)の収束可否を決定する。
説明
上記のメソッドはμ≤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)