９：N次指数関数法と呼んでもいいっすか？

前回、必要な安定領域は実軸方向で\([-3,0]\)、虚軸方向で\([-\pi i,+\pi i]\)の範囲と勝手に決めつけたわけですが、この領域で安定でまともな計算精度がある数値解法はあるのでしょうか？
いずれじっくりやりたいと思いますが、こういう場合の定番の後退微分法（BDF）やAdams-Moulton法が候補に入らなことは事は図1^[1]いずれじっくりやりますので、図が甘いのはご容赦くださいで示しておきます。
後退微分法（Gear法とも言う）は3段^[2]Wikipediaの訳に従い”段”としましたが、英語ではStepです。Runge-Kutta法のときは、Stageを段と訳しているので、非常に紛らわしいの場合は（緑の線）と虚軸がわずかに安定領域から外れ、4段になるとオイシイ領域が不安定になり、6段を超えると不安定という、実はけっこう痛い奴です^[3]Stiff-Solverの代表格である後退微分法は3次以上はA安定じゃない、というのが衝撃の事実。

おや、台形積分も、後退Eulerも、Adams-Moultonも、Gearも、全部だめなら、候補は残るのか？

それがあります。線形式に限定すれば、ありえます。

最初の「常微分方程式と状態方程式」で書きましたが、線形の常微分方程式である状態方程式
\begin{equation}
\frac{dX(t)}{dt}=AX(t)+BU(t)
\label {eq:StateEq_9}
\end{equation}
の場合、\(X_i=X(t_0)、　X_{i+1}=X(h+t_0)\)とすれば、
\begin{equation}
X_{i+1}=e^{A h}X_i + (e^{A h} – I )A^{-1} BU_i
\label{eq:diffExp_9}
\end{equation}
となるので式\eqref{eq:diffExp_9}を計算すれば、\(X_iからX_{i+1}\)が即座に得られますよね、と。

ですが、行列指数関数\(e^{A}\)は、
\begin{equation}
e^{A}:=I+A+\frac{A^2}{2!}+…+\frac{A^N}{N!}= \sum_{k=0}^{N} \frac{A^k}{k!}
\label{eq:defMatExp_9}
\end{equation}
において\(N\rightarrow\infty\)なんですが、これはどうやって計算するのでしょうか？
式\eqref{eq:defMatExp_9}をそもまま頑張って計算する？　\(N=\infty\)まで？

たしかに、\(\infty\)まではできませんが、この手の式の高次項は収束するのがお約束なので適当なところまで頑張ればOK、のはずです。適当なところまで頑張ってみましょう。

しかし、「台形積分法に言いがかりをつける」でも書きましたが、\(h\lambda=3\)の場合\(( 3^k/k! )\)が１より小さくなるのは\(k=7\)の時です。従って、少なくとも\(N=7\)を超えないとマズそうなので、イロを付けて\(N=10\)まで頑張ってみたのが図2です。

図2　N次で打ち切った指数関数の安定領域

式\eqref{eq:defMatExp_9}の行列指数関数を\(N\)次で打ち切ったものを使用して、式\eqref{eq:diffExp_9}を計算した場合の安定領域は、実はRunge-Kutta法と同じ^[4]「Taylor展開とRunge-Kutta法」の時にやりましたなので、4次まではよく見かけると思いますが、5次から先は初めて見たのではないでしょうか？
まちの制御屋もこの図2が出てくる論文^[5]”On the absolute stability regions corresponding to partial sums of the exponential function” D.Ketchenson et al. IMA Journal of Numerical Analysis, Volume 35, Issue 3, July 2015, Pages … Continue readingを見つけたときは、胸が張り裂けそうでした。

式\eqref{eq:diffExp_9}は厳密解なので、\(N=\infty\)まで頑張ればA安定のはずですが、たかだか\(N=10\)程度では、安定領域はかなり小さいです。ですが\(N=8\)を超えれば、まちの制御屋が意味のある安定領域とした領域を含みます。さらに\(N=10\)であれば^[6] … Continue reading、計算精度は\(h^{10}\)まで適合していて、なんとワンダフル。

• 

• 

　図3a　\(h\lambda=\pm\pi i \)の場合　　　　図3b　\(h\lambda=-3\)（重根）の場合

さらに、式\eqref{eq:defMatExp_9}は順番に計算するだけなので、陽解放（Explicit-Method)です。反復計算は必要ありません。

をぉぉ、これこそがまちの制御屋が求めていた数値解法ではないですか？

この手法は、あまりに当たり前なのか、式\eqref{eq:diffExp_9}を計算するのがあまりに面倒なのか^[7] … Continue reading、線形に限定するのは数値計算屋の沽券にかかわるのか、教科書では見かけない手法です。
特に名称も決まっていないようです。
確かに式\eqref{eq:StateEq_9}をTaylor展開すれば得られるのでTaylor展開法と呼んでもよいのでしょうが、偏微分が出てこないTaylor展開法なんて。。。という気もします。

ですので、「N次指数関数法」と、呼ぶことにします。
あぁ、「N次指数関数法」に栄光あれ。

という、画期的な数値解法を掘り起こしたわけですが、
残念ながら、この手法は線形方程式限定です。

「線形方程式限定？　お客を選んだら、商売にならへんのや」
数値解法君の嘆く声が聞こえます^[8] … Continue reading。

やっぱ、任意関数でどうにかならんですかね？