Julia で波動方程式シミュレーション

NOTE

本記事では, DifferentialEquations.jl の Tist5 で数値計算をしています. また, 計算結果を GIF で可視化していますが, ファイルサイズの都合上, 読み込みに時間がかかる可能性があります.

波動方程式

波動方程式とは, 次式で表される定数係数二階線形偏微分方程式である.

\begin{equation*} \dfrac{1}{c^{2}}\dfrac{\partial^{2}u}{\partial t^{2}}=\Delta u \end{equation*}

ここで, $c$ は波動の位相速度を表す係数である.

1次元波動方程式

1次元の場合, 次式で表される.

\begin{equation*} \dfrac{1}{c^{2}}\dfrac{\partial^{2}u}{\partial t^{2}}=\dfrac{\partial^{2}u}{\partial x^{2}} \end{equation*}

ここで, $0\leq x\leq L$ , $0\leq t$ .
また, 本節では, 境界条件はディリクレ境界条件とし, 初期変位はガウス分布とする.

wave_1d.jl

using DifferentialEquations
using Plots
using LaTeXStrings
plot_font = "Computer Modern"
default(fontfamily=plot_font, linewidth=1, label=nothing, grid=true)

L = 1.0
c = 1.0
nx = 100
dx = L / (nx - 1)
x = range(0, stop=L, length=nx)

u0 = exp.(-100 .* (x .- L/2).^2)
v0 = zeros(nx)

function wave_eq!(du, u, p, t)
    du[1:nx] .= u[nx+1:2*nx]
    for i in 2:nx-1
        du[nx+i] = c^2 * (u[i+1] - 2u[i] + u[i-1]) / dx^2
    end
    du[nx+1] = 0.0
    du[2*nx] = 0.0
end

u0_combined = [u0; v0]
tspan = (0.0, 2.0)

prob = ODEProblem(wave_eq!, u0_combined, tspan)
sol = solve(prob, Tsit5(), dt=0.01)

@gif for (i, u) in enumerate(sol.u)
    t = sol.t[i]
    plot(x, u[1:nx], ylim=(-1, 1), xlabel=L"x", ylabel=L"u(x,t)", title=L"t=" * string(round(t, digits=2)), lw=2, dpi=300)
end

2次元波動方程式

2次元の場合, 次式で表される.

\begin{equation*} \dfrac{1}{c^{2}}\dfrac{\partial^{2}u}{\partial t^{2}}=\dfrac{\partial^{2}u}{\partial x^{2}}+\dfrac{\partial^{2}u}{\partial y^{2}} \end{equation*}

ここで, $0\leq x\leq L$ , $0\leq y\leq L$ , $0\leq t$ .
また, 本節では, 境界条件はディリクレ境界条件とし, 初期変位はガウス分布とする.

wave_2d.jl

using DifferentialEquations
using Plots
using LaTeXStrings
plot_font = "Computer Modern"
default(fontfamily=plot_font, linewidth=1, label=nothing, grid=true)

Lx = 1.0
Ly = 1.0
c = 1.0
nx = 50
ny = 50
dx = Lx / (nx - 1)
dy = Ly / (ny - 1)
x = range(0, stop=Lx, length=nx)
y = range(0, stop=Ly, length=ny)

u0 = [exp(-100 * ((xi - Lx/2)^2 + (yi - Ly/2)^2)) for xi in x, yi in y]
v0 = zeros(nx, ny)

function wave_eq!(du, u, p, t)
    du[1:nx, 1:ny] .= u[nx+1:2*nx, 1:ny]
    for i in 2:nx-1
        for j in 2:ny-1
            du[nx+i, j] = c^2 * (
                (u[i+1, j] - 2u[i, j] + u[i-1, j]) / dx^2 +
                (u[i, j+1] - 2u[i, j] + u[i, j-1]) / dy^2
            )
        end
    end
    du[nx+1, :] .= 0.0
    du[2*nx, :] .= 0.0
    du[:, 1] .= 0.0
    du[:, ny] .= 0.0
end

u0_combined = [u0; v0]
tspan = (0.0, 2.0)

prob = ODEProblem(wave_eq!, u0_combined, tspan)
sol = solve(prob, Tsit5(), dt=0.01)

@gif for (i, u) in enumerate(sol.u)
    t = sol.t[i]
    heatmap(x, y, u[1:nx, 1:ny], clim=(-1, 1), xlabel=L"x", ylabel=L"y", title=L"t=" * string(round(t, digits=2)), dpi=300)
end

3次元波動方程式

3次元の場合, 次式で表される.

\begin{equation*} \dfrac{1}{c^{2}}\dfrac{\partial^{2}u}{\partial t^{2}}=\dfrac{\partial^{2}u}{\partial x^{2}}+\dfrac{\partial^{2}u}{\partial y^{2}}+\dfrac{\partial^{2}u}{\partial z^{2}} \end{equation*}

ここで, $0\leq x\leq L$ , $0\leq y\leq L$ , $0\leq z\leq L$ , $0\leq t$ .
また, 本節では, 境界条件はディリクレ境界条件とし, 初期変位はガウス分布とする.

wave_3d.jl

using DifferentialEquations
using Plots
using LaTeXStrings

plot_font = "Computer Modern"
default(fontfamily=plot_font, linewidth=1, label=nothing, grid=true)

Lx = 1.0
Ly = 1.0
Lz = 1.0
c = 1.0
nx = 30
ny = 30
nz = 30
dx = Lx / (nx - 1)
dy = Ly / (ny - 1)
dz = Lz / (nz - 1)
x = range(0, stop=Lx, length=nx)
y = range(0, stop=Ly, length=ny)
z = range(0, stop=Lz, length=nz)

u0 = [exp(-100 * ((xi - Lx/2)^2 + (yi - Ly/2)^2 + (zi - Lz/2)^2)) for xi in x, yi in y, zi in z]
v0 = zeros(nx, ny, nz)

u0_flat = reshape(u0, :)
v0_flat = reshape(v0, :)
u0_combined = [u0_flat; v0_flat]

function wave_eq!(du, u, p, t)
    u_field = reshape(u[1:nx*ny*nz], nx, ny, nz)
    v_field = reshape(u[nx*ny*nz+1:end], nx, ny, nz)
    
    du_field = zeros(nx, ny, nz)
    dv_field = zeros(nx, ny, nz)

    for i in 2:nx-1
        for j in 2:ny-1
            for k in 2:nz-1
                du_field[i, j, k] = v_field[i, j, k]
                dv_field[i, j, k] = c^2 * (
                    (u_field[i+1, j, k] - 2u_field[i, j, k] + u_field[i-1, j, k]) / dx^2 +
                    (u_field[i, j+1, k] - 2u_field[i, j, k] + u_field[i, j-1, k]) / dy^2 +
                    (u_field[i, j, k+1] - 2u_field[i, j, k] + u_field[i, j, k-1]) / dz^2
                )
            end
        end
    end

    du_field[:, 1, :] .= 0.0
    du_field[:, ny, :] .= 0.0
    du_field[:, :, 1] .= 0.0
    du_field[:, :, nz] .= 0.0

    du_combined = [reshape(du_field, :); reshape(dv_field, :)]
    du .= du_combined
end

tspan = (0.0, 2.0)

prob = ODEProblem(wave_eq!, u0_combined, tspan)
sol = solve(prob, Tsit5(), dt=0.01)

@gif for (i, u) in enumerate(sol.u)
    t = sol.t[i]
    u_field = reshape(u[1:nx*ny*nz], nx, ny, nz)
    
    surface(x, y, u_field[:, :, Int(floor(nz/2))], clim=(-1, 1), zlims=(-1, 1), xlabel=L"x", ylabel=L"y", zlabel=L"u", title=L"t=" * string(round(t, digits=2)), dpi=300)
end

まとめ

今回は, 1～3次元の波動方程式をシミュレーションし, 可視化した. 数値解法を自分で実装することも数値解法への理解を深めるためには必要であるが, 小規模なシミュレーションならば既存のライブラリ (ODEソルバ) で十分¹.

Footnotes

実装時間は短縮できる. ↩