数理解析研究所講究録 1295 巻 2002 年 51-55

51

遅延微分方程式の級数による解法
平山弘 (Hiroshi Hirayama)’

神奈川工科大学1

佐藤創大郎 (Soutarou Satou)\dagger
神奈川工科大学

はじめに関数を Taylor 級数展開することは、すでに多くの人によって研究されている自動微分と呼ばれる方法 [12]と数学的には同じものである。 この自動微分を利用した常微分方程式の数値解法については、 伊理、小野、
戸田 $[7][10]$ によって、 詳しく論じられている。 これら論文では、 自動微分を利用して、 微分係数が得られるので、 これを利用すると Runge-Kutta 法を改良することができることが示されている。 しカルながら、

この方法は、 まだ一般化していないプリ・プロセッサの存在を仮定している。 また、高次の公式を得るには、 通常の Runge-Kutta の公式 [14][15] と同様にかなり難しいにように思われる。

常微分方程式の解を、 任意の次数まで Taylor 級数展開できることは、 数学的には、 常微分方程式の級数展開による解法と同じであり、 その可能性については、 久保田 [8] などで述べられており、 Corliss andChang[l] によって、 任意の次数の Taylor 級数展開式を利用した常微分方程式を解くためのパツケージ・ソフトウエアが発表されている。 Corliss and Chang のソフトウエアは、 上記の自動微分と同じようにプリ.
プロセッサを利用しているパッケージ・ソフトウエアである。
著者等 [15] は、 プログラムでよく使われる演算子 ( ,-,*,/など) を、 被演算の型が異なる場合、 別の意$+$味を与えることができる機能、 Fortran 90 や C+十言語 [2] の機能 (オペレータ・オーバロード、 operatoroverload) を使い、 有限項で打ち切った Taylor 級数問の四則演算、 Taylor 級数の関数演算を定義する。 この機能を使うと、 プログラムの形で与えられた任意の関数を Taylor 級数展開することができる。 常微分方程式の初期値問題の解を任意次数の Taylor 級数展開の形で、容易に得られる。得られた Taylor 級数解を数値計算に利用すれば、 任意の次数の数値解法が得られる。 この Taylor 級数展開をある次数の Pade’ 展開式に変形すると、 任意次数の A 安定な数値計算方法を与えることを示した。
法) と比較する
この計算方法は、微分方程式の数値解法 [9] に、 よく使われる陽的 Runge-Kutta 法 (と、 計算次数が任意とれる点や A 安定な計算法であるため、 計算のステップサイズを大きくとっても不安定性が生じない等の利点がある。
陰的 Runge-Kutta 法 (IRK 法) と比較すると、 A 安定な計算法となるが、 IRK 法は計算式が一般に非線形連立方程式となるため、各計算ステップ毎に非線形連立方程式を解く必要がある。 Newton 法等の反復法$\mathrm{R}\mathrm{K}$が必要になるので計算が複雑になり、 計算時間も多くかかる。
A 安定な計算法は、 陽的計算法ではあり得ないとの証明があるが、 この計算法は、 Pade’ 展開するために連立一次方程式を解くことが陽的計算となっている。 この連立一次方程式は、 一般的な連立一次方程式ではなく、特別な形をした Toeplitz 型となる。 この連立一次方程式は、通常の連立一次方程式は、未知数の数をとすると $O(n^{3})$ の計算量が必要であるが、 Toeplitz 型は Levinson-Durbin の方法 [11] を使えば、 $O(n^{2})$$n$’hirayama@sd kanagawa-it.ac jP$\uparrow \mathrm{s}\mathrm{o}\mathrm{u}@\mathrm{p}\mathrm{r}\mathrm{o}.\mathrm{s}\mathrm{d}$.kanagawa-it.ac$.\mathrm{j}\mathrm{p}$VII-I52の計算量で計算できる。 さらに、 $O(n\log n)$ の計算量で計算できることが知られている。 この事実を考慮すれば、 ほぼ陽的な計算法と変わらない計算で任意次数の A 安定な計算ができることがわかる。

本論文では、 この計算法を遅延微分方程式 [3] に適用することを提案する。 ここで提案する計算法は、現在最もよく使われるプログラム言語の一つであり、入手が容易な C+十言語を使う。 このような計算方法は、数値計算でよく使われる Fortran 90 でも可能である。 オペレータ・オーバロードの機能がない語やFortran 77 では、 このような計算は不可能ではないが、 現在までこのような計算が行われていない事実を$\mathrm{C}-\overline{\equiv}$見ると、 実際上不可能と思われる。

2Taylor 展開
この節では、 関数を Taylor 展開するための基本的な考え方を説明し、 その計算方法について簡単に論じる。詳しくは、 RallHenrici[4] や平山 [5] などに述べられている。
$[12]_{\text{、}}$

2.1

Taylor 級数の四則演算

Taylor 級数の四則計算のプログラムは、 以下のように簡単に作ることができる。 平行移動によって、 展開位置を原点へ移すことができるので一般性を失うことなしに、原点で展開した式だけを扱うことができる。 この級数を次のように定義する。

2.1.1$f(x)=f\mathrm{o}+f_{1}x+f_{2}x^{2}+f_{3}x^{3}+f_{4}x^{4}+\cdots$(1)$g(x)=g0+g_{1}x+g_{2}x^{2}+g_{3}x^{3}+g_{4}x^{4}+\cdots$(2)$h(x)=h_{0}+h_{1}x+h_{2}x^{2}+h_{3}x^{3}+h_{4}x^{4}+\cdots$(3)加減算$h(x)$ が $f(x)$と $g(x)$ の和差のとき、 $f,$ $g$ および $h$ の係数は、 次のような関係になる。

$h(x)=f(x)\pm g(x)$212:(4)

乗算$h(x)$ が $f(x)$と $g(x)$ の積のとき、 $f,$ $g$ および $h$ の係数は、 次のような関係 (こなる。

$h(x)=f(x)g(x)$213$h_{:}=f_{i}\pm g$$h_{n}= \sum_{k=0}^{n}f_{k}g_{n-k}$(5)除算$h(x)$ が $f(x)$と $g(x)$ の商のとき、 すなわち

$h(x)= \frac{f(x)}{g(x)}$VII-2(6)53であるとき、 $f,$$g$および$h$の係数には、 以下のような関係が成り立つ。

(2)、 および (3) を代入して、 展開する。両辺の同じ次数(7) の式は、 (6) において、 両辺にを掛け、 (1)の係数が等しいことを利用して得られる。

2.2(7)$h_{n}= \frac{1}{g_{0}}(f_{n}-\sum_{k=0}^{n-1}h_{k}g_{n-k})$$h_{0}= \frac{f_{0}}{\mathit{9}0}$

Taylor 級数の関数

他の多くの初等関数でも同様に係数間の漸化式を導くことができる。 ここでは、 Taylor 級数の指数関数の計算方法を示す。 この方法は後に述べる常微分方程式の解の Taylor 展開の単純な例にもなっている。
$h(x)=e^{f(x)}$

とおくと、

(8)

$\frac{dh}{dx}=h\frac{df}{dx}$である。 (8)$[]_{\llcorner}-(1)_{\text{、}}$(3) を代入して、 両辺の係数を比較することによって、 次のような関係が得られる。
$h_{0}=e^{f\mathrm{o}}$$h_{n}= \frac{1}{n}\sum_{k=1}^{n}kh_{n-k}f_{k}$(9)$(n\geq 1)$対数関数や三角関数の場合も同様な公式を得ることができる。式 (9) のように、 Taylor 展開式の指数関数の計算には、 指数関数の計算は、 1 回しか入らないので、 Taylor 展開の指数関数の計算量は、通常の四則演算とそれほど大きな違いにはならない。 このことは、 対数関数や三角関数などでも同様である。

3遅延微分方程式の解の Taylor 展開
この節では、 Taylor 展開の方法を利用して、 遅延微分方程式の解を任意の次数まで Taylor 展開する方法を示す。 Taylor 展開された解は、 任意次数の一段公式として、 Runge-Kutta 公式の代わりとして利用できる。

次のように$\frac{d}{d}\mathrm{X}x$について解かれた形になっている遅延微分方程式について考える。
$\frac{d\mathrm{y}}{dx}=\mathrm{f}(x, \mathrm{y}(x),$ここで、$\mathrm{y}$および$\mathrm{f}$(10)$\mathrm{y}(x-\tau))$

は、 一般にベクトルである。 初期条件は$\mathrm{y}(x)=\mathrm{g}(x)$(11)$x_{0}\geq x\geq x_{0}-\tau$で与えられる。 この微分方程式の解の Taylor 展開は、 以下に示す Picard の逐次計算法 [13] を使うことによって計算する。

(12)$\mathrm{y}_{n}=\mathrm{y}0+\int_{0}^{x}\mathrm{f}(x,\mathrm{y}_{n-1}(t),\mathrm{y}(t-\tau))dt$

ただし、 (12) の計算はは、 過去の次以上の高次の項は省略して計算する。 この計算において、解であるから、既知関数である。補間計算等を行わないようにするため、 計算ステップサイズとして、 を$\mathrm{y}(t-\tau)$$n$$\tau$整数個に分割した大きさをとる。 このようにとれば、補間計算を行う必要がないため、誤差が入りにくくなるため、 計算精度の向上につながる。 この積分を何回も繰り返すことによって、解の Taylor 展開式を得ることが出来る。得られた Taylor 展開を利用して、 次の計算位置の定数項を計算する。 この定数項を初期値VII-354として、 (12) の計算を繰り返すことによって、 次の計算位置での Taylor 展開を得る。 この操作を繰り返すことによって、 微分方程式を解くことができる。

Taylor 展開を利用して、 次の位置における定数項の計算を行う楊合、 この Taylor 展開式を Pade’ 展開式に変換し、 それを利用して、 次の位置における定数項を計算すると安定な計算法になる。 以下の手順を$\mathrm{A}$まとめると次のようになる。

(1) 計算ステップサイズを(2) 初期値関数$\mathrm{g}(x)$$\tau$の整数分の 1 の大きさにする。

を利用して、 微分方程式の初期値 (yo) を計算する。

(3) (12) を使って、 Taylor 展開式を求める。
(4) 得られた Taylor 展開式を可能ならば、 Pade’ 展開式に変換する。
(5)

この Pade’ 展開式または Taylor 展開式を利用して、 次の計算 Taylor 展開位置の定数項を計算する。

(6) さらに計算する場合、 (3) に戻る。

4

数値例
厳密解が知られている、 次の方程式 [16] を解く。 $\varphi(x)=exp(2+cos(x)^{2})$ のとき,$\Delta_{=-y(x-1)}ddx(1+y(x)^{2})+\varphi(x-1)(1+\varphi(x)^{2})+\varphi’(x)$$y(x)=\varphi(x)$,$2\geq x\geq 0$$0\geq x\geq-1$(13)を解くことを考える。 この時の厳密解は、 $y=\varphi(x)$ である。 この微分方程式を $x=0$ から、 計算ステップを $h=1/1024$ としてとして解く。 $x=0$ における 6 次の Taylor 展開式を求めると、$y=20.0855-20$ $.0855*x^{2}+16.7379*x^{4}-5.49829e+13*x^{6}$
この式は、 分子分母 3 次の Pade’ 展開できないので、 この Taylor 展開式を利用して、 $x=$ 旧こおける$y=20.0854$ を計算する。 この値を利用して、 $x=$ 眉こおける 6 次の Taylor 展開式を求めると、$y=20.0854-0.0655854*x$ -22.6421’ $x^{2}+337.936*x^{3}$ -33376.9’ $x^{4}+2.63471e+06*x^{5}+6.01748e+07*x^{6}$
これを y=num/&n と Pade’ 展開すると$num=20.0855-2814.47*x-105800*x^{2}+6.84746e+06*x^{3}den=1-140$$.124*x-5266.46*x^{2}+340720*x^{3}$となる。 この Pade’ 展開式を利用して、 $x=2h$ における$y$の値を計算すると $y=20.0854$ となる。 $x=2h$における 6 次の Taylor 展開式を求めると、$20.0854-0.0655854*x-22.6421*x^{2}+337.936*x^{3}-33376.9*x^{4}+2.63471e+06*x^{5}+6.01748e+07*x^{6}$となる。 6 次の係数を見ると非常に大きな係数になっているが、計算ステップ十分に小さいので Pade’ 展開を利用しないでも十分高精度で計算できる。 この計算を 1 万ステップ計算を進めても、 6 桁以上の精度を保ち非常に安定した計算ができる。

VII-4

555

まとめ
本論文では、 Taylor 展開を利用して、遅延方程式を解く方法を提案した。 この方法は通常の微分方程式を解く場合と同じように、高精度で安定した計算が可能でる。
遅延微分方程式では、 硬い (Stiff) 状態が発生し、 その場合本方法が有用であることを示すことができると期待したが、残念ながら、 そのような例を挙げることができなかった。
この方法と他の方法を比較し、 どのような場合、 どちらの計算方法が良いか判断すべきであるが、厳密解がわかっているような問題では、 その差が現れるような現象が見受けられなかった。

いろいろな問題を調べ、 本方法を優劣を調べることがこれからの課題である。

参考文献[1] Corliss G. and Chang Y. F., Solving Ordinary Differential Equations Using Taylor Series, ACM
bans. Math. Soft.,. 114-144 (1982)$\mathrm{V}\mathrm{o}\mathrm{l}.8,$$\mathrm{p}\mathrm{p}$[2] Ellis M. A. and Stroustrup B., The Annotated$\mathrm{C}++\mathrm{R}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}$Manual, Addison-Wesley,(1990)[3] Hairer E., Wanner G., Solving Ordinary Differential Equations , Springer-Verlag, (1987)$\mathrm{I}$[4] Henrici, P., Applied Computational Complex Analysis, Vol. 1, John Wiely&Sons, New
1, (1974)$\mathrm{Y}\mathrm{o}\mathrm{r}\mathrm{k},\mathrm{C}\mathrm{h}\mathrm{a}\mathrm{p}$[5] 平山弘, C+十言語によるべき型特異点をもつ関数の数値積分, 日本応用数理学会論文誌,
257-266 (1995)[6] 平山, 小宮, 佐藤, Taylor 級数法による常微分方程式の解法, 日本応用数理学会論文誌,$\mathrm{V}\mathrm{o}\mathrm{l}.5,$$\mathrm{V}\mathrm{o}\mathrm{l}.12$$\mathrm{p}\mathrm{p}$.
.

, (2002)

出版予定[7] 伊理, 小野, 戸田, 合成関数の高速微分法とその導関数を含む Runge-Kutta 系の常微分方程式数値解法
公式への応用, 情報処理学会論文誌, 沙 26, pp. 389-396(1986)[8] 久保田光一, コンピューテングの玉手箱続. 自動微分 Taylor 展開, bit,$\mathrm{V}\mathrm{o}\mathrm{l}.22$, pp. 860-861 (1991)[9] 三井斌友、 数値解析入門、 朝倉書店、 (1985)[10] 小野, 戸田, 伊理, 微分係数を用いた埋め込み型 Runge-Kutta 系 2 段公式について, 情報処理学会論文誌, 沙 28, PP. 807-814(1987), Flannery, Teukolsky S.A., Vetterling W.T., Numerical Recipes in , Cambridge
University Press, Cambridge, (1988) (邦訳 : 丹慶、 奥村、 佐藤、 小林訳、 技術評論社、 (1993))[11] Press[12]$\mathrm{W}.\mathrm{H}.$$\mathrm{R}\mathrm{a}11,\mathrm{L}.$ence,$\mathrm{B}.\mathrm{P}.$$\mathrm{C}$B. , Automatic Differentiation-Technique and Applications, Lecture Notes in Computer Sci, Springer-Verlag, Berlin-Heidelberg-New York(1981)$\mathrm{V}\mathrm{o}\mathrm{l}.120$[13] 佐野理, キーポイント微分方程式, 岩波書店, 東京, (1993)[14] 田中, 高山, 安定性のよい 9 段数 7 次陽的 Runge-Kutta 法について, 情報処理学会論文誌,
52-61 (1993)[15] 田中, 高山,2 段数陰的 Runge-Kutta 法について, 情報処理学会論文誌,$\mathrm{V}\mathrm{o}\mathrm{l}.36,$[16] 小藤, 梁, 遅延微分方程式の半群的解法、 常微分方程式の数値解法とその周辺,
ベンチャービジネスラボラトリ, pp. 145-158(2000)

VII-5$\mathrm{p}\mathrm{p}$沙 34,$\mathrm{p}\mathrm{p}$.

. 226-234 (1995)
名古屋大学