📂수치해석

열 방정식의 수치적 풀이: 유한차분법

열 방정식의 수치적 풀이¹

아래와 같은 1차원 열 방정식이 주어졌다고 하자.

$$\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}$일 때 수렴한다. 이는 공간 간격에 따라 한 번에 전진할 수 있는 타임 스텝의 upper bound가 존재함을 말해준다. 즉 공간 도메인을 촘촘하게 나눌수록, 한 번에 전진할 수 있는 타임스텝도 그만큼 줄어든다. 간단한 예제에서라면 문제 없겠지만, 시간 도메인의 길이가 길다면 계산 시간이 오래걸리는 원인이 될 수 있다. 아래의 예제에서만 봐도 공간 도메인을 크기가 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)$를 결합하면 증명 끝.

■

예제

줄리아로 이를 다음과 같이 구현할 수 있다.

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