数値解析におけるスペクトル法
メソッド
偏微分方程式の数値的な解法として、偏微分方程式自体を解く代わりに解の候補をフーリエ級数と置き、フーリエ係数の常微分方程式として迂回する手法を スペクトル法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)$ に要約されたと見なせる。例として Kuramoto–Sivashinsky 方程式をスペクトル法で解くことを考えてみる。
方程式は次のように始まる。個人的にはこの方程式は線形項が二つ以上あり、さらに非線形項も含むためスペクトル法を学ぶには格好の例だと思う。 $$ {\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]$ だけの Kuramoto–Sivashinsky 方程式をスペクトル法で解く Julia による実装を示している。
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")
これは単にアニメーションを作成するための可視化コードである。