フーリエ・スペクトル法による常微分方程式の周期解の数値計算

$\newcommand{\im}{\mathrm{i}}$ 常微分方程式の一つであるvan der Pol方程式の周期解の精度保証付き数値計算を紹介するため、フーリエ・スペクトル法による周期解の数値計算を行う。

van der Pol 方程式

van der Pol方程式とは、以下のような方程式である。

$$ \frac{d^2 x}{dt^2} - \mu (1-x^2)\frac{dx}{dt} + x = 0. $$

未知関数は $x(t)$ で、$\mu>0$ は非線形の減衰の強さを表すパラメータである。フーリエ・スペクトル法の計算に必要な参照軌道をまず得るために、van der Pol方程式をJulia言語のDifferentialEquations.jlというパッケージを使って解の挙動を数値計算する。

まず、van der Pol 方程式を次の連立常微分方程式系にしてDifferentialEquations.jlのODEソルバーで数値計算する。

$$ \begin{cases} \dot{x} = y\\ \dot{y} = \mu (1-x^2)y - x \end{cases} $$

初期値 $x(0)=0$, $y(0)=2$とし, $\mu=1$ のときの数値計算は以下のように実行される。

数値計算で得た結果、ある程度の時間が過ぎると下図のような周期軌道（に近い解の軌道）が現れる。

Newton法の初期値の設定

周期解をフーリエ級数で現わし、その係数と周期を求めるため、van der Pol方程式の周期解の周期を大まかに求める。

得た近似周期軌道と近似周期を使って、軌道のフーリエ補間を計算する。詳細はフーリエ級数を参照。

Newton法を用いた周期解の計算

van der Pol方程式は、$\dot{x} = \frac{dx}{dt}$ とおくと、以下のように表すことができる。

$$ \ddot{x} - \mu(1-x^2)\dot{x} + x = 0. $$

あとの計算のために、式を整理すると、

$$ \ddot{x} - \mu\dot{x} + \frac{\mu}{3} \dot{(x^3)} + x = 0. $$

周期解 $x(t)$ を周期 $L$ の周期関数とし、$\omega = \frac{2\pi}{L}$ とおくと、$x(t)$ とその微分やべき乗はフーリエ級数を使って、

$$ \begin{align*} x(t) &= \sum_{k \in \mathbb{Z}} a_k e^{\im k\omega t}\\ \frac{dx(t)}{dt} &= \sum_{k \in \mathbb{Z}}(\im k \omega) a_k e^{\im k \omega t} \\ \frac{d^2 x(t)}{dt^2} &= \sum_{k \in \mathbb{Z}} (- k^2 \omega^2 )a_k e^{\im k\omega t} \\ x(t)^3 &= \sum_{k \in \mathbb{Z}} (a * a * a)_k e^{\im k \omega t} \end{align*} $$

と書くことができる。ここで

$$ (a * a * a)_k := \sum_{\substack{k_1+k_2+k_3 = k\\k_i\in\mathbb{Z}}} a_{k_1}a_{k_2} a_{k_3},\quad k\in\mathbb{Z} $$

は3次の離散たたみこみである。

$\newcommand{\bx}{\mathrm x}$

そしてフーリエ係数に関する式を立てる。係数 $a = (a_k)_{k \in \mathbb{Z}}$ に対して、van der Pol方程式に求めたフーリエ級数を代入すると、

$$ f_k(a) := -k^2\omega^2 a_k - \mu\im k \omega a_k + \frac{\mu }{3}(\im k \omega)(a*a*a)_k + a_k $$

となる点列 $\left(f_k(a)\right)_{k\in\mathbb{Z}}$がえられる。そして、各 $k\in\mathbb{Z}$ について

$$ f_k(a) = 0 $$

となる点列 $a$ が得られれば、van der Pol方程式の解のフーリエ係数になる。未知数は周波数 $\omega$ と点列 $a$ であり、これらを並べて $\bx = (\omega,a)$ と書くことにする。未知数 $\bx$ に対して、$f_k(a) = 0$ という方程式だけでは不定な方程式になるため、解の形を一つに定める事ができない（並進対称性がある）。そこで、位相条件

$$ \begin{align*} \eta (a) &:= \sum_{|k|<N} a_k - \eta_0=0, \quad \eta_0 \in \mathbb{R} \end{align*} $$

を加える。この条件は $x(t)$ の切片 $x(0) = \eta_0$ を表している。最終的に van der Pol 方程式の周期解の求解は次の代数方程式を解くことに帰着される。

$$ F(\bx) := \begin{bmatrix} \eta (a) \\ (f_k(a)_{k\in\mathbb{Z}} \end{bmatrix} =0. $$

以下、この零点探索問題 $F(\bx)=0$ について Newton 法で解を得ることを考える。まず $N$ をフーリエ係数の打ち切り番号（最大波数:$N-1$）とし、周期解の近似を次のように構成する。

$$ \bar{x}(t) = \sum_{|k|<N} \bar{a}_k e^{\im k \bar\omega t}. $$

このとき、フーリエ係数と（近似）周期をならべた

$$ \bar{\bx} = (\bar\omega, \bar{a}_{-N+1},\dots,\bar{a}_{N-1})\in \mathbb{C}^{2N} $$

を近似解とよぶ。近似解 $\bar \bx$ の項数は $2N$ 個。そして $f_k(a)=0$ を $|k|<N$ の範囲で打ち切る方程式

$$ F^{(N)}(\bx^{(N)}) = \begin{bmatrix} \eta (a^{(N)}) \\ (f_k(a^{(N)}))_{k < |N|} \end{bmatrix} =0 $$

を考える。ここで $a^{(N)} = (a_k)_{|k|<N}$, $\bx^{(N)} = (\omega,a^{(N)})$ をそれぞれ表し、$F^{(N)}:\mathbb{C}^{2N}\to \mathbb{C}^{2N}$ である。したがって $F^{(N)}(\bx^{(N)}) = 0$ という有限次元の非線形方程式を解くことで、近似解 $\bar \bx$ が得られる。

Newton法による周期解の数値計算

実際にNewton法を用いて、周期解の数値計算を行っていく。Newton法では、ある適当な初期値 $\bx_0$ を最初に定め、以下の反復計算によって計算できる。

$$ \bx_{n+1} = \bx_n - DF^{(N)}(\bx_n)^{-1} F^{(N)}(\bx_n),\quad n=0,1,\dots $$

このことから $DF^{(N)}(\bx_n)^{-1}$ と $F^{(N)}(\bx_n)$ を計算することができれば、近似解を得ることができる。はじめに離散畳み込みの関数powerconvfourierを用意する。詳細は離散畳み込みの精度保証を参照。

次のF_fourierは $F^{(N)}(\bx_n)$ を計算している。

ヤコビ行列の計算式

$F^{(N)}(\bx^{(N)})$のヤコビ行列は、

$$ DF^{(N)}(\bx^{(N)}) = \left[\begin{array}{c|ccc} 0 & 1 & \dots & 1\\\hline \vdots & &\vdots&\\ \partial_{\omega}f_k& \dots & \partial_{a_j}f_k & \dots\\ \vdots & &\vdots& \end{array}\right]\in\mathbb{C}^{2N\times 2N}\quad (|k|,|j|<N). $$

ここで、

$$ \begin{cases} \partial_\omega f_k = (-2k^2 \omega - \mu \im k) a_k + \frac{\mu \im k}{3}(a*a*a)_k & (|k|<N)\\ \partial_{a_j} f_k = (-k^2 \omega^2 - \mu \im k \omega + 1) \delta_{kj} + \mu \im k \omega (a*a)_{k-j}&(|k|,|j|<N) \end{cases},\quad \delta_{kj} = \begin{cases} 1 & (k=j)\\ 0 & (k\neq j) \end{cases} $$

である。ヤコビ行列の各要素との対応は

$$ \left(DF^{(N)}(\bx^{(N)})\right)_{\ell,m} = \begin{cases} 0 \ &(\ell=m=1) \\ 1 \ &(\ell=1, m = 2 \cdots 2N) \\ \partial_\omega f_k &(\ell = 2 \cdots 2N, m = 1,~\mbox{i.e.},~\ell = k + N + 1~\mbox{for}~|k|<N)\\ \partial_{a_j} f_k &(\ell,m = 2 \cdots 2N,~\mbox{i.e.},~\ell = k + N + 1~\mbox{for}~|k|<N,~m = j + N + 1~\mbox{for}~|j|<N). \end{cases} $$

Newton法の反復

Newton法の反復を以下のコードで実行する。

Newton法の3回の反復でJuliaのDifferentialEquations.jlで計算した軌道を初期値からvan der Pol方程式の近似周期軌道を得ることに成功した。この軌道のプロファイル、フーリエ係数の情報は以下で確認できる。

本資料は以下のような文献・Web ページ等を参考に書いています。

参考文献

  1. DifferentialEquations.jl: Scientific Machine Learning (SciML) Enabled Simulation and Estimation (https://diffeq.sciml.ai/stable/).
    (Julia言語のDifferentialEquations.jlのマニュアル。今回取り上げた常微分方程式だけでなく、遅延微分方程式、偏微分方程式、微分代数方程式の解法なども実装されている。)
  2. Jan Bouwe van den Berg, Jean-Philippe Lessard, Proceedings of Symposia in Applied Mathematics 74, American Mathematical Society, 2018.
    (力学系におけるNewton-Kantorovich型定理を使った精度保証付き数値計算の教科書、フーリエ級数・チェビシェフ級数の畳み込みを使った実践例をいくつか紹介している)
  3. 高橋和暉, 高安亮紀, フーリエ級数, 2022.
    （FFTを利用したフーリエ補間の計算方法をJuliaで実装する例を紹介している）
  4. 高橋和暉, 高安亮紀, 離散畳み込みの精度保証, 2022.
    （FFTを利用した離散畳み込みの計算方法とその精度保証付き数値計算方法をJuliaで実装する例を紹介している）
高橋和暉, 高安亮紀，2022年3月7日