📂수치해석

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

메서드

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

$$\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$일 때 수렴한다.

2차원

2차원인 경우에 시공간을 다음과 같이 쪼갰다고 하자.

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

그러면 근사식은 다음과 같다.

$$\dfrac{\partial ^{2} u}{\partial x^{2}} \approx \dfrac{u_{i+1,j}^{n} - 2u_{i,j}^{n} + u_{i-1,j}^{n}}{h^{2}}$$

$$\dfrac{\partial ^{2} u}{\partial y^{2}} \approx \dfrac{u_{i,j+1}^{n} - 2u_{i,j}^{n} + u_{i,j-1}^{n}}{h^{2}}$$

$$\dfrac{\partial ^{2} u}{\partial t^{2}} \approx \dfrac{u_{i,j}^{n+1} - 2u_{i,j}^{n} + u_{i,j}^{n-1}}{\Delta t^{2}}$$

$$\begin{align*} \implies & \dfrac{\partial ^{2} u}{\partial t^{2}} - \left( \dfrac{\partial ^{2} u}{\partial x^{2}} + \dfrac{\partial ^{2} u}{\partial y^{2}} \right) \\ &\approx \dfrac{u_{i,j}^{n+1} - 2u_{i,j}^{n} + u_{i,j}^{n-1}}{\Delta t^{2}} - \dfrac{u_{i+1,j}^{n} - 2u_{i,j}^{n} + u_{i-1,j}^{n}}{h^{2}} - \dfrac{u_{i,j+1}^{n} - 2u_{i,j}^{n} + u_{i,j-1}^{n}}{h^{2}} = 0 \end{align*}$$

$$\implies u_{i,j}^{n+1} = 2(1 - 2\mu^{2})u_{i,j}^{n} - u_{i,j}^{n-1} + \mu^{2} \left( u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} \right)$$

여기서 $\mu = \dfrac{\Delta t}{h}$이다.

코드

경계 조건을 $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)