📂数値解析

熱方程式の数値解：有限差分法

数値的な熱方程式の解¹

以下のような一次元の熱方程式が与えられたとしよう。

$$ \dfrac{\partial u}{\partial t} = \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 u(x,t)}{\partial t} = \dfrac{1}{\Delta t} \big[ u(x,t+\Delta t) - u(x,t) \big] + \mathcal{O}(\Delta t) $$

これで、$\mu = \dfrac{\Delta t}{(\Delta x)^{2}}$としよう。 この二つの式を $(1)$に代入して、$\Delta t = \mu (\Delta x)^{2}$を掛けると次を得る。

$$ \dfrac{u(x,t+\Delta t) - u(x,t)}{\Delta t} = \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) $$ $$ \implies u(x,t+\Delta t) = u(x,t) + \dfrac{\mu (\Delta x)^{2}}{(\Delta x)^{2}} \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] + \mathcal{O}\left( (\Delta x)^{4} \right) $$ $$ \implies u(x,t+\Delta t) = u(x,t) + \mu \big[ u(x-\Delta x, t) - 2u(x,t) + u(x+\Delta x, t) \big] + \mathcal{O}\left( (\Delta x)^{4} \right) $$ $$ \implies u_{\ell}^{n+1} = u_{\ell}^{n} + \mu \left( u_{\ell-1}^{n} - 2u_{\ell}^{n} + u_{\ell+1}^{n} \right),\quad \ell=1,2,\dots,d,\ n=0,1,\dots \tag{2} $$

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

説明

以下の定理から、$(2)$は $\mu \le \frac{1}{2}$の時に収束する。これは、空間間隔に応じて、一度に進められるタイムステップの上限が存在することを意味する。つまり、空間ドメインを細かく分割するほど、一度に進められるタイムステップもそれだけ小さくなる。単純な例では問題ないかもしれないが、時間ドメインの長さが長い場合、計算時間が長くなる原因になり得る。例えば、空間ドメインを大きさ100のベクトルで分割した場合、時間ドメインはなんと1733個に分割される。

収束性

$\mu \le \dfrac{1}{2}$であれば、数値解$(2)$は収束する。

証明

$t^{\ast} \gt 0$を任意の定数としよう。各点での近似値と実際のソリューションのエラーを次のように定義しよう。

$$ e_{\ell}^{n} := u_{\ell}^{n} - u(\ell\Delta x, n\Delta t),\qquad \ell = 0,1,\dots,d+1,\quad n=0,1,\dots,n_{\Delta t} $$

ここで、$n_{\Delta t} = \lfloor t^{\ast}/{\Delta t} \rfloor = \lfloor t^{\ast} / (\mu (\Delta x)^{2}) \rfloor$は$n$の終点である。 方法が収束するとは、エラー$e_{\ell}^{n}$に対して次の式が成立することとして表現できる。

$$ \lim\limits_{\Delta x \to 0} \left[ \max\limits_{n=0,1,\dots,n_{\Delta t}} \left( \max\limits_{\ell=0,1,\dots,d+1} \left| e_{\ell}^{n} \right| \right) \right] = 0 $$

括弧内の値を次のようにしよう。

$$ \eta^{n} := \max\limits_{\ell=0,1,\dots,d+1} \left| e_{\ell}^{n} \right|,\quad n=0,1,\dots,n_{\Delta t} \tag{3} $$

その場合、収束に対する表現は、

$$ \lim\limits_{\Delta x \to 0} \left( \max\limits_{n=0,1,\dots,n_{\Delta t}} \eta^{n} \right) = 0 \tag{4} $$

今、$\tilde{u}_{\ell}^{n} = u(\ell\Delta x, n\Delta t)$を$u_{\ell}^{n}$の真値としよう。

$$ \begin{align*} u_{\ell}^{n+1} &= u_{\ell}^{n} + \mu \left( u_{\ell-1}^{n} - 2u_{\ell}^{n} + u_{\ell+1}^{n} \right) \\ \tilde{u}_{\ell}^{n+1} &= \tilde{u}_{\ell}^{n} + \mu \left( \tilde{u}_{\ell-1}^{n} - 2\tilde{u}_{\ell}^{n} + \tilde{u}_{\ell+1}^{n} \right) + \mathcal{O}\left( (\Delta x)^{4} \right) \end{align*} \qquad \begin{array}{l} \ell = 0,1,\dots,d+1 \\ n=0,1,\dots,n_{\Delta t} \end{array} $$

これを引くと、次を得る。

$$ \begin{align*} e_{\ell}^{n+1} &= e_{\ell}^{n} + \mu \left( e_{\ell-1}^{n} - 2e_{\ell}^{n} + e_{\ell+1}^{n} \right) + \mathcal{O}\left( (\Delta x)^{4} \right) \\[1em] &\le e_{\ell}^{n} + \mu \left( e_{\ell-1}^{n} - 2e_{\ell}^{n} + e_{\ell+1}^{n} \right) + c (\Delta x)^{4},\qquad c \gt 0 \end{align*} $$

不等式は漸近記法の定義によって成立する。三角不等式により、次を得る。

$$ \begin{align*} \left| e_{\ell}^{n+1} \right| &\le \left| e_{\ell}^{n} + \mu \left( e_{\ell-1}^{n} - 2e_{\ell}^{n} + e_{\ell+1}^{n} \right) + c (\Delta x)^{4} \right| \\ &\le \left| e_{\ell}^{n} + \mu \left( e_{\ell-1}^{n} - 2e_{\ell}^{n} + e_{\ell+1}^{n} \right) \right| + c (\Delta x)^{4} \\ &\le \mu \left| e_{\ell-1}^{n} \right| + \left| 1 - 2\mu \right| \left| e_{\ell}^{n} \right| + \mu \left| e_{\ell+1}^{n} \right| + c (\Delta x)^{4} \\ &\le (2\mu + \left| 1 - 2\mu \right|)\eta^{n} + c (\Delta x)^{4},\qquad n=0,1,\dots,n_{\Delta t}-1 \end{align*} $$

最後の不等式は$(3)$によって成立する。$\mu \le 1/2$であるので、$(2\mu + \left| 1 - 2\mu \right|)=1$であり、再び$(3)$から次を得る。

$$ \eta^{n+1} = \max\limits_{\ell=0,1,\dots,d+1} \left| e_{\ell}^{n+1} \right| \le \eta^{n} + c (\Delta x)^{4},\qquad n=0,1,\dots,n_{\Delta t}-1 $$

これを帰納的に繰り返すと、

$$ \eta^{n+1} \le \eta^{n} + c (\Delta x)^{4} \le \eta^{n-1} + 2c (\Delta x)^{4} \le \eta^{n-2} + 3 c (\Delta x)^{4} \le \cdots $$ $$ \implies \eta^{n} \le \eta^{0} + n c (\Delta x)^{4},\qquad n=0,1,\dots,n_{\Delta t} $$

初期値が$\eta^{0} = 0$であり、$n(\Delta x)^{2} = n\Delta t / \mu = t^{\ast} / \mu$であるので、

$$ \eta^{n} \le \dfrac{ct^{\ast}}{\mu}(\Delta x)^{2},\qquad n=0,1,\dots,n_{\Delta t} $$

したがって、すべての$n$に対して、

$$ \eta^{n} \to 0\quad \text{ as } \Delta x \to 0 $$

これと$(4)$を組み合わせると証明終わり。

■

例

これをJuliaで次のように実装できる。

1次元

今、次のような1次元の熱方程式を考えよう。

$$ \begin{array}{lc} \text{Domain:} & \\ & \Omega = [-1,1]\times [0,0.35] \\ \text{Heat equation:} & \\ & u_{t} - u_{xx} = 0 \\[0.5em] & u(x,0) = \sin (\pi x) \\[0.5em] & u(-1, t) = 0,\quad u(1, t) = 0 \\[0.5em] \text{Exact solution:}& \\ & u(x,t) = \sin(\pi x) e^{-\pi^2 t} \end{array} $$

ここで、$\Delta x = \dfrac{1-(-1)}{99}$そして$\Delta t = \mu (\Delta)^{2} = 0.5(\Delta x)^{2}$として設定し、$(2)$で解くと、次のような結果が得られる。

一方、$\Delta t = \mu (\Delta)^{2} = 0.51(\Delta x)^{2}$として解くと、解が正しく得られないことが分かる。

using Plots
using LinearAlgebra

function Solver_Heat_1d(Δt, T, x, f, g)
    # Δt: time step,    T: last time,   x: space vector
    # f: initial value,  g: boundary value
    
    t = Vector(0:Δt:T)
    n = length(t)
    ℓ = length(x)
    Δx = x[2]-x[1]
    μ = Δt/((Δx)^2)
    μ = Δt/((Δx)^2)

    U = zeros(ℓ,n)
    U[1,:] .= g[1]
    U[end,:] .= g[end]

    U[:,1] = f

    for i ∈ 2:n
        U[2:end-1, i] = (1-2μ)*U[2:end-1, i-1] + μ*(U[1:end-2, i-1] + U[3:end, i-1])                     
    end
    return U
end

function ref_sol(x, t, exact_sol)
    xt =[kron(ones(length(t)),x);; kron(t,ones(length(x)))]
    ref_U = reshape(exact_sol(xt),length(x),length(t))
    return ref_U
end

function results(t, x, U, ref_U, μ)
    h1 = heatmap(t, x, U, colormap=:viridis, clim=(-1,1))
    # xlabel!("time [0, $(round(t[end], digits=2)) ]")
    # ylabel!("space [$(x[1]),  $(x[end])]")
    title!("FDM solution")

    h2 = heatmap(t, x, ref_U, colormap=:viridis, clim=(-1,1))
    title!("exact solution")

    h3 = heatmap(t, x, U - ref_U, colormap=:viridis)
    title!("error")

    h4 = heatmap(framestyle=:none)
    annotate!([0.5],[0.5], text("μ = $(round(μ, digits=3))", 30))

    plot(h1, h2, h3, h4, size=(1500,750))
end

for μ ∈ [0.5, 0.51]
    x = Vector(LinRange(-1., 1, 100))
    Δx = x[2]-x[1]
    initial_f = sin.(π*x)
    Δt = μ*((Δx)^2)
    T = 0.35
    bdry_g = [0.0, 0.0]
    t = Vector(0:Δt:T)

    exact_sol(xt) = sin.(π*xt[:,1]) .* exp.(- (π)^2 * xt[:,2])
    ref_u = ref_sol(x, t, exact_sol)
    U = Solver_Heat_1d(Δt, T, x, initial_f, bdry_g)
    results(t, x, U, ref_u, μ)
    savefig("mu=" * string(μ) * ".png") 
end

2次元

$$ \begin{array}{lc} \text{Domain:} & \\ & \Omega = [-2,2] \times [-2,2] \times [0,1] \\ \text{Heat equation:} & \\ & u_{t} - u_{xx} = 0 \\[0.5em] & u(x,y,0) = \frac{1}{4} e^{-x^{2} -y^{2}} \\[0.5em] & u(x, y, t) = 0\qquad \text{ on } \partial \left( [-2,2] \times [-2,2] \right) \end{array} $$

同じ方法で上の微分方程式を解くと、

using Plots

cd = @__DIR__

function Solver_Heat_2d(Δt, T, x, f)
    # Δt: time step,    T: last time,   x: space vector
    # f: initial value
    
    t = Vector(0:Δt:T)
    n = length(t)
    ℓ = length(x)
    Δx = x[2]-x[1]
    μ = Δt/((Δx)^2)

    U = zeros(ℓ,ℓ,n)
    U[:,:,1] = f

    for i ∈ 2:n
        U[2:end-1, 2:end-1, i] = (1-4μ)*U[2:end-1, 2:end-1, i-1] + μ*(U[1:end-2, 2:end-1, i-1] + U[3:end, 2:end-1, i-1] + U[2:end-1, 1:end-2, i-1] + U[2:end-1, 3:end, i-1])                     
    end

    return U
end

x = range(-2., stop=2, length=100)'
y = range(-2., stop=2, length=100)

initial_f = @. (1/4)exp(-x^2) * exp(-y^2) 

Δx = x[2]-x[1]
Δt = 0.49*((Δx)^2)/2
T = 1.0

U = Solver_Heat_2d(Δt, T, x, initial_f)

anim = @animate for i ∈ range(1, stop=size(U)[3], step=100)
    surface(U[:,:,i], zlims=(0,0.5), camera=(30,30), colorbar=:none, showaxis = false, grid=false, size=(500,500))
end

gif(anim, cd*"/2d.gif", fps=10)

環境

• OS: Windows11
• Version: Julia v1.8.3, Plots v1.38.6