」<前回

      次回>「

5:数値解法の次数とTaylor展開

このBlogではStiffな常微分方程式の数値解法を扱っているのですが、多くの数値解法は、計算式の段階で近似式です。
それは、離散化した差分方程式だから、ではなく、差分方程式であっても厳密解を得る計算式はあるが、そうなっていない、という事です。

では差分方程式であっても厳密解を求める式とは、つまり\(x_{i+1}=x(t_0 + h)  x_i=x(t_0)\)とした場合、\(x_iからx_{i+1}\)を厳密に求める式とはなんでしょうか?
それはTaylor展開です。

厳密な定義はWikipediaを参照してもらうとしてTaylorの定理-Wikipedia
\begin{align}
x(h + t_0)=x(t_0)&+\left.h\frac{dx}{dt} \right|_{t=t_0} +\left. \frac{h^2}{2} \frac{d^2x} {dt^2} \right|_{t=t_0}\nonumber \\
&+…+ \left. \frac{h^\infty}{\infty!}\frac{d^\infty x}{dt^\infty} \right|_{t=t_0}= \sum_{k=0}^{\infty} \frac{h^k}{k!}\left.\frac{d^kx}{dt^k} \right|_{t=t_0}\\
\label{eq:TaylorExpantion}
\end{align}
が成立します[1]厳密には、\(x(t)\)が無限に微分可能(=\(C^\infty\)級)とかの条件が付くはずです。記号の使い方、あってますかね?

つまり、\(x_i=x(t_0)\)の値と、そこでの高次の導関数(\(\frac{d^k}{dt^k}x(t)\)等)があれば、\(x_{i+1}=x(t_0+h)\)が得られる、と、


対象とする常微分方程式を
\begin{equation}
\frac{d}{dt}x(t) = f(x(t))
\end{equation}
とすれば、
\begin{align}
x(t_0 + h)=x(t_0)&+ h f(t_0) + \left. \frac{h^2}{2} \frac{d f} {dt} \right|_{t=t_0}\nonumber\\
&+…+ \left. \frac{h^{k+1}}{(k+1)!}\frac{d^kf}{dt^k} \right|_{t=t_0}+…+ \left. \frac{h^\infty}{\infty!}\frac{d^\infty f}{dt^\infty} \right|_{t=t_0}
\label{eq:TylorExpantODE}
\end{align}
ですね。

しかしですね、「\(\frac{d^k}{dt^k}f(t)\)ってどこの子ですか?」と言いたいです。

そんな知らない子は無視すると、
\begin{eqnarray}
x(t_0 + h)&=&x(t_0)+h f(x(t_0)) \nonumber\\
x_{i+1}&=&x_i+h f(x_i)
\end{eqnarray}
ですね。これは前進Euler法です。

前進Euler法は1次の精度を持つ、あるいは1次で適合している[2]他の文献の言い方は多少違いますが、参考文献にあげた「常微分方程式と微分方程式の数値解法」ではこの言い方ですとは、厳密解を与えるTaylor級数の\(h\)の1次の項まで一致している(そしてその後は無視している)という意味であることが解りました。

一方、式\eqref{eq:TylorExpantODE}は\(h\)がマイナスであっても通用するので、
\begin{eqnarray}
x(t_0 – h)&=&x(t_0) – h f(x(t_0)) \nonumber\\
x_{i-1}+h f(x_{i}) &=&x_{i} \nonumber\\
x_{i+1} &=&x_i + h f(x_{i+1})
\end{eqnarray}
後退Euler法と一致します。後退Euler法もやはり、1次でしか適合しない、ことが解りました。

では足して2で割る台形積分法はどうでしょうか?
展開の中心を\(x_{i+1/2}=x(t_0 + \frac{h}{2})\)として\(\pm{\frac{h}{2}}\)で式\eqref{eq:TylorExpantODE}を使い展開すると、
\begin{eqnarray}
x_{i+1}&=&x_{i+1/2}+\frac{h}{2} f(x_{i+1/2})+\frac{h^2}{2\cdot2^2} \frac{df}{dt}+ \frac{h^3}{6\cdot2^3} \frac{d^2f}{dt^2}+…\\
\label{eq:TrapezoidP}
x_{i}&=&x_{i+1/2}-\frac{h}{2} f(x_{i+1/2})+\frac{h^2}{2\cdot2^2} \frac{df}{dt} – \frac{h^3}{6\cdot2^3} \frac{d^2f}{dt^2}+…
\end{eqnarray}
式(6)ー式(7)を計算すると、\(h^2\)の項が消去されるので、
\begin{equation}
x_{i+1}-x_{i}= h f(x_{i+1/2}) +2\frac{h^3}{6\cdot2^3} \frac{d^3f}{dt^2}
\label{eq:TrapIMD1}
\end{equation}
ここで、さらに\( f(x_{i+1/2})\)を\(\pm{\frac{h}{2}}\)で展開し、
\begin{eqnarray}
f(x_{i+1})&=&f(x_{i+1/2}) + \left.\frac{h}{2} \frac{df}{dt} \right|_{x=x_{i+h/2}} +\frac{h^2}{2\cdot2^2}\frac{d^2f}{dt^3}+… \\
 f(x_{i})&=&f(x_{i+1/2}) -\left.\frac{h}{2} \frac{df}{dt} \right|_{x=x_{i+h/2}} + \frac{h^2}{2\cdot2^2}\frac{d^2f}{dt^3}+…
\label{eq:trapIMD2}
\end{eqnarray}
その和をとると
\begin{eqnarray}
 f(x_{i})+f(x_{i+1}) =2f(x_{i+1/2}) +\left.\frac{h^2}{2^2} \frac{d^2f}{dt^2}\right|_{x=x_{i+h/2}} +… \nonumber\\
f(x_{i+1/2}) =\frac{1}{2}\left(f(x_{i})+f(x_{i+1}) \right) -\left.\frac{h^2}{2^3} \frac{d^2f}{dt^2}\right|_{x=x_{i+h/2}} +…
\label{eq:TrapIMD3}
\end{eqnarray}
式\eqref{eq:TrapIMD3}を式\eqref{eq:TrapIMD1}に代入し、知らない子が関わっている\(h^3\)に続く項を無視すれば、
\begin{equation}
 x_{i+1}=x_i+\frac{h}{2} \left( f(x_{i}) +f(x_{i+1}) \right)
\end{equation}
という台形積分法の式が得られます。

\(h^3\)以上の項を切り捨てているので、2次の式ということになります。

以上、前進Euler法、後退Euler法、台形積分法の計算式がTaylor級数の一部と一致しているが、同時にTaylor級数のかなり多くの項を無視していることがわかりました。ですがが、無視されるのは寂しいので、もっと多くの項まで面倒を見ることはできないでしょうか?

そのためには、知らない子\(\frac{d^k}{dt^k}f(t)\)をなんとかしなくてはいけません。

参考文献
「常微分方程式と微分方程式の数値解法」 U.M. アッシャー,L.R.ペツォルド 共著,中村 眞理雄 監訳 嘉村友作・三ツ間 均 共訳 2006年 培風館
「微分方程式による計算科学入門」    三井 斌友 小藤 俊幸 斉藤 善弘 2004年共立出版

脚注
脚注
1 厳密には、\(x(t)\)が無限に微分可能(=\(C^\infty\)級)とかの条件が付くはずです
2 他の文献の言い方は多少違いますが、参考文献にあげた「常微分方程式と微分方程式の数値解法」ではこの言い方です

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です