수치해석에서의 스펙트럴 메서드
메서드
편미분방정식의 수치적인 풀이를 위해, 편미분방정식 자체를 푸는 대신 해의 후보를 푸리에 급수로 두고 푸리에 계수의 상미분방정식으로써 우회하는 기법을 스펙트럴 메서드spectral method라 한다.
설명
그게 어떤 편미분방정식이든 관계 없이, 우리는 해의 푸리에 전개가 가능하다는 가정 하에서 다음과 같이 충분히 큰 $N$ 에 대해 $u$ 의 푸리에 급수로써 해를 근사시킨다는 아이디어에서 시작하려고 한다. $$ u ( t , x ) = \mathcal{F} (u) \approx \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) $$ 여기서 우리가 알 수 있는 것과 알 수 없는 것이 무엇인지 살펴보면, 우선 $\exp \left( i k x \right)$ 는 $k$번째 푸리에 계수에 대응되어 계산할 수 있고 $x$ 는 격자점의 위치에 해당하는 값으로써 언제나 알 수 있는 것에 반해 $k$번째 $c_{k}$ 가 시간에 종속된 $c_{k} (t)$ 는 알 수 없다. 그나마 만약 이들이 다음과 같이 상미분방정식으로 표현되는 지배 방정식을 따른다면, 수치해석적인 솔버를 통해 비교적 간단하게 풀 수 있을 것이다. $$ \begin{align*} {\frac{ d c_{-N} }{ dt }} =& f_{-N} \left( c_{-N} , \cdots , c_{N} \right) \\ \vdots & \\ {\frac{ d c_{N} }{ dt }} =& f_{N} \left( c_{-N} , \cdots , c_{N} \right) \\ \implies \dot{\mathbf{c}} =& \mathbf{f} \left( \mathbf{c} \right) \end{align*} $$
이제 우리가 모르는 모든 것은 시스템을 표현하는 $\mathbf{f}$ 와 초기값 $\mathbf{c} (0)$ 으로 요약되었다고 볼 수 있다. 예로써 쿠라모토-시바신스키 방정식을 스펙트럴 메서드로 푼다고 생각해보자.
방정식은 다음과 같이 시작한다. 개인적으로 이 방정식은 선형적인 항이 두 개 이상, 거기에 비선형적인 항까지 포함되어 있어 스펙트럴 메서드를 공부하기에 탁월한 예시라고 생각한다. $$ {\frac{ \partial u }{ \partial t }} + {\frac{ \partial }{ \partial x }} \left( {\frac{ 1 }{ 2 }} u^{2} \right) + {\frac{ \partial^{2} u }{ \partial x^{2} }} + {\frac{ \partial^{4} u }{ \partial x^{4} }} = 0 $$
선형항 1
일단은 비선형항인 $u^{2}$ 은 없다고 치자. $u = \sum c \exp$ 이라 가정했으므로 그대로 대입하면 다음과 같다.
$$ \begin{align*} & {\frac{ d }{ d t }} \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) \\ =& - {\frac{ \partial^{2} }{ \partial x^{2} }} \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) \\ & - {\frac{ \partial^{4} }{ \partial x^{4} }} \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) \end{align*} $$
편미분을 모조리 처리하자.
$$ \begin{align*} & \sum_{k = -N}^{N} \dot{c}_{k} (t) \exp \left( i k x \right) \\ =& - \sum_{k = -N}^{N} \left( ik \right)^{2} c_{k} (t) \exp \left( i k x \right) \\ & - \sum_{k = -N}^{N} \left( ik \right)^{4} c_{k} (t) \exp \left( i k x \right) \end{align*} $$
푸리에 기저의 직교성: 두 집합 $\left\{ e^{inx} \right\}_{n=-\infty}^\infty$와 $\left\{ \cos nx\ \right\}_{n=0}^\infty \cup \left\{ \sin nx \right\}_{n=1}^\infty$는 $L^{2}(-\pi,\ \pi)$의 정규직교기저 이다. 또한 $\left\{ \cos nx \right\}_{n=0}^{\infty}$와 $\left\{ \sin nx \right\}_{n=1}^{\infty}$는 $L^{2}(0,\ \pi)$의 정규직교기저이다.
덧셈으로 표현된 식이지만, 만약 KS 방정식에 비선형 항이 없었다면 푸리에 기저의 직교성에 따라 $k$ 에 대해 다음과 같이 따로따로 풀 수 있다. $$ \dot{c}_{k} (t) = k^{2} c_{k} (t) - k^{4} c_{k} (t) $$
수학적으로야 이것이 몇 차원이 되든 어려운 미분방정식은 아니지만, 수치해석적인 관점에서는 $k^{4}$ 라는 항은 꽤나 부담이 될 수 있다. 가령 $N = 32$ 정도만 되어도 $k^{4} = 1,048,576 \approx 10^{6}$ 으로써 스티프한 문제가 될 것이다. 물론 실제로는 공간의 크기에 따라 달라지겠지만, 푸리에 전개의 편미분에서 네제곱 항이 떨어지는 것은 어떤 방식으로든 응어리로 남는 게 당연하다.
상식적으로 생각해봐도 아무리 놀라운 트릭이 있다고 한들 그 어려운 편미분방정식이 그냥 상미분방정식으로 떨어질 리는 없다. 본디 유한차분법과 같은 방식에서 미분계수를 근사하면서 가해졌어야 할 부하가, 로젠브록 메서드처럼 한 스텝마다의 큰 비용으로 전가되는 것이라 보면 된다.
비선형항
푸리에 변환의 성질: (c) $f$가 연속이면서 조각마다 매끄럽다고 하자. 그러면 $$ \mathcal{F} \left[ f^{\prime} \right] (\xi) = i \xi \mathcal{F} f (\xi) $$
푸리에 변환에서 $k$번째 항에 대해 $\left( u^{2} / 2 \right)_{x}$ 은 다음과 같이 우변에 반영되어야 한다. $$ \mathcal{F} \left[ {\frac{ \partial }{ \partial x }} \left( {\frac{ 1 }{ 2 }} u^{2} \right) \right] (k) = i k \mathcal{F} \left[ {\frac{ 1 }{ 2 }} u^{2} \right] (k) $$ 따라서 비선형항까지 포함하는 미분방정식은 $k$번째 항에 대해 다음과 같이 나타난다. 우변에 $u$ 가 남은 것처럼 보이지만, 실제로는 푸리에 변환 $\mathcal{F}$ 을 통해 주파수 도메인으로 넘어가서 계산된다. $$ \dot{c}_{k} (t) = - ik \mathcal{F} \left[ {\frac{ 1 }{ 2 }} u^{2} \right] (k) + \left[ k^{2} - k^{4} \right] c_{k} (t) $$
코드
위의 영상은 $L = 22$ 이고 그리드의 수가 $64$ 이면서 시간 간격 $h = 0.01$ 일 때 $t \in [0, 200]$ 만큼의 쿠라모토-시바신스키 방정식을 스펙트럴 메서드로 푸는 것을 줄리아로 구현한 것을 보여주고 있다.
using DifferentialEquations, FFTW, SparseArrays, LinearAlgebra, Random, Plots, Arpack
function ks_spectral_rhs!(du_hat, u_hat, p, t)
k, Lk = p.k, p.Lk
u = real.(ifft(u_hat))
u22 = (u .^ 2) / 2
nonlinear_term = -im .* (k .* fft(u22))
linear_term = Lk .* u_hat
du_hat .= linear_term .+ nonlinear_term
return nothing
end
상미분방정식의 우변을 보면 본질적으로 수식으로 살핀 것과 정확히 같음을 쉽게 확인할 수 있다. 계속해서 $u$ 에 대한 정보를 주파수 도메인에 업데이트 해주어야 하기에 매 스텝마다 푸리에 변환이 사용된다.
L = 22.0; Q = 64; dt = 0.01; T = 200
tspan = 0:dt:T
x = (0:Q-1) .* (L/Q)
k = vcat(collect(0:Q÷2-1), 0, collect(-Q÷2+1:-1)) .* (2pi / L)
Lk = k .^ 2 .- k .^ 4
Random.seed!(0)
prob = ODEProblem(ks_spectral_rhs!, fft(randn(Q)), extrema(tspan), (; k, Lk))
@time sol = solve(prob, ROS2(autodiff=false), saveat=collect(tspan), reltol=1e-4, abstol=1e-6);
u = stack(real.(ifft.(sol.u)))
실제로 스티프한 상미분방정식의 풀이에는 2단계 로젠브록 메서드가 사용되었다. 푸리에 역변환은 모든 풀이가 끝난 뒤에 단 한 번만 쓰면 된다.
_p1 = heatmap(u, xlabel = "t", ylabel = "grid", xticks = [])
@time anime = @animate for tk in 1:100:size(u, 2)
p1 = vline(_p1, [tk], label = :none, color = :white)
p2 = plot(u[:, tk], label = :none, ylims = [-3, 3], color = :black, xlabel = "grid", ylabel = "u")
plot(p1, p2, layout = (2, 1), size = [600, 600])
end
mp4(anime, "ks_pde.mp4")
이는 단지 애니메이션을 만들기 위한 시각화 코드다.