1.1　演算子法の概略
. . 等間隔変数t_k=t₀+kh, k=0,1,2,3,.......における関数f(t)の値がf₀, f₁, f₂, f₃, ........であるとき、関数f(t)はニュートンの補間公式で表され、等間隔変数t_kに対する積分値y_kはスターリングの積分公式で表される。従って、y_kの階差はf_kの階差を用いた公式で表される。その公式は行列を用いた数式で(1.1)のように記述できる。公式の行列は上部三角行列であり、第i行は第i−1行を1列右に移して第1列に0をセットした構造を持つ。その第1行は次のようになる。

1.2　常微分方程式の数値解法
. . 微分方程式y' =f(y, t)及び初期値y₀が与えられたとき、(1.1)式の１行１列の部分行列を用いてΔy₀を予測する。次に、２行２列の部分行列とf(y, t)により求められたΔf₀を用いてΔy₀を修正し、Δ²y₀を予測する。最後に、２行３列の行列とf(y, t)により修正されたΔf₀とΔ²f₀を用いてΔy₀とΔ²y₀を数回修正する。次に、初期値点をt₂に移して同じ操作を繰り返す。

(1.1)

. . 等間隔区分点をn+1個使用する解法では積分範囲[t₀, t_n]の積分をn行n+1列の行列を用いた同様の方法で求め、初期値点をt_nへ移す。その新しい範囲の等間隔区分点数は、積分が収束する限り、新しい任意の区間長hを用いて決めてよい。その理由は、この解法は初期値のみから出発でき、その他の区分点の値は全て予測と修正により決定されるからである。
. . 全積分範囲をn区間に等分割し、その１区間を更に２等分するとき、(1.1)の解法をn回繰り返せば全積分範囲の解が得られる。この方法は、Runge-Kuttaより精度がよく、以下に示すような簡単にして新しい解法を導く。ここより、変数tは実数を、xはそれをある精度に丸めた値を取るものとする。

(1.2)

. . 手順[1]はEuler法によりy₁を予測する。手順[2]は台形法によりy₁を修正し、同じf₁を用いた中点法でy₂を予測する。手順[3]は負の区分幅を持つ３次のAdams-Moulton法によりy₁を修正し、同じf₁及びf₂を用いたSimpson法によりy₂を修正する。手順[3]は更に２回繰り返される。その理由は、f₁及びf₂は手順[3]で得たy₁及びy₂を用いた値でないから、これら及びy₁とy₂は修正されるべきであり、更に、もう一度修正したとき解は収束しているべきであるからである。もしそれが収束していなければ区分幅hが不適当である。積分区間h₂がHではなく、h₂=x₂−x₀である理由は後で述べる。
. . 関数z(t)により駆動されるn階の常微分方程式を解く方法は二通りある。その一つは微分行列Dの多項式を計算して得た行列の逆行列を計算しなければならないのでやや複雑である。もう一つは微分方程式が(1.3)で表せるので(1.2)の解法をn回繰り返せばよく、簡単である。

(1.3)

1.3　プログラムと実行結果
. . 著者の解法はBASICによりProgram 1 のように書くことが出来る。これは配列のサイズを６に設定しているので１階から５階迄の微分方程式を解くことが出来る。プログラムは行番号56と58に副プログラムを与えなければならない。前者は(1.3)の第１式の右辺を与え、後者は、全積分範囲、等間隔区分数、微分方程式の階数、及び、初期値を与える。行番号60は計算結果を比較するための真値を与える。しかし、行番号54の"X#=X2:GOSUB 60"及び変数Y#と共に削除することが出来る。解の値は末尾に単精度長の０を付加して２倍長の２進数に変換した後８桁の浮動少数に10進化している。

(1.4)

1.4　数値解の精度
. . １積分区間の２分の１をhで表すとき微分方程式(1.4)の解は手順[1][2][3]及び[3]により級数に展開され、x=2hにおいて(1.5)で表される。Runge-Kutta法の解も同様に級数に展開されるが最後の項が存在しない。最後の項の分母はTaylor展開では30である。Table 1 はこれら三者の関係を示している。

(1.5)

著者の解法でも手順[3]の繰返しがないと最後の項は現れない。手順[3]を２回より多く繰り返すと2hの６乗より高次の項が現れるがそれらは正しくない。これらが無視し得る値でないなら数値解は真の解とは少し異なった値に収束する。これは従来のP-C法と同じ現象であり、区分幅hが適当でない。これらが(1.4)のy"=f(y' )の解y' について手順[3]を３回繰り返す理由である。
. . Runge-Kutta法は１積分区間の中点に対してx₁, f₁, y₁及び区間幅hを用い、積分区間の端に対してx₂, f₂, y₂及び区間幅2hを用いると次のように記述できる。

(1.6)

これは著者の解法とは次の点で異なっている。手順[2']は後退Euler法によりy₁を修正し、更に、f₁を修正している。手順[3']は２つのf₁の平均値を用いたSimpson法でy₂を修正しているがy₁は修正しない。手順[3']を繰り返すことが意味があるか否かは不明である。
. . 例えば、dy/dt =−y, y(0)=1の解は４次のTaylor展開に展開されるが、もう一度手順[3']を繰り返しても何の違いも得られない。その理由は、k₄が中点における関数fの値にのみ依存するy₀+k₃により決定されるからである。しかし、手順[3']のようにy₂でy₀+k₃を置き換えるなら、それは(1.6)に示すようにx=2hに対して予測した値でもよいから繰返しにより2hの５乗の項が生ずるが正しくない。更に、4乗の項が正しくないものに変わってしまう。そして、数値計算結果はTable 2 に2h=0.2の場合を示すが精度が悪くなる。

Table 2. The results repeated [3] and [3']
Author's method		Runge-Kutta		True
1×[3] 2×[3] 3×[3]	0.8186667 0.8187333 0.8187311	1×[3'] 2×[3'] 3×[3']	0.8187333 0.8187089 0.8187097	0.81873075

. . 著者の解法もhの４乗の項はTaylor展開の項に必ずしも一致しない。著者の解法はTaylor展開を求めるのではなく、Newtonの補間公式の積分を求めている。区分幅hが大きい場合には前者は後者より精度のよい近似法ではない。３点が等間隔であるとき、関数f(t)は剰余項つきの２次のTaylor展開で表され、その積分は剰余項の積分を剰余項とする３次のTaylor展開で表され、その剰余項はh⁴以上の項からなる。著者の解法の剰余項はTaylor展開とは少し異なったh⁴の項を持っているが、剰余項全体はRunge-Kuttaより精度良く近似されている。従って、著者の解法の数値解はTable 1, Table 2 及びTable 3 に示すように精度がよいのである。
. . 解が４次以上のTaylor展開で近似されるべきであるとき、例えばdy/dt=y⁻¹, y(0)=0.5, 2h=0.125の場合、著者の解法の解はTable 3 に示すように７次のTaylor展開より約３桁精度がよい。Simpsonの1/3則によるy⁻¹の積分はTaylor展開より約２桁精度がよい。即ち、４次の項がTaylor展開と一致することは重要ではない。３次の項より高次の項の全体を正確に近似することが重要である。従って、数値解法はStirlingの積分公式により導くべきである。

Table 3. dy/dt=y−1, y(0)=0.5, 2h=0.125
t	Author's	Simpson	Runge-K	Taylor(7)	True
0.125	0.7071319	0.7072122	0.7071896	0.7106933	0.70710678

. . Runge-Kuttaの解法は４次の項までTaylor展開と一致し、それより高次の項は無視し得るものとして導かれているが、その解は７次のTaylor展開より約２桁精度がよい。その理由は、その手順が(1.6)に示すように著者の解法に似ているからである。しかし、著者の演算子法でその解法の正当性は証明されない。

1.5　数値解の存在条件
. . (1.2)に示した手順[1][2][3][3][3]により得られるy₁及びy₂の値はそれぞれが微分方程式の解に収束する数列である。それは(1.3)の最初の微分方程式のy⁽ⁿ⁻¹⁾についても同様である。何故ならその微分方程式はy⁽ⁿ⁻¹⁾に関する１階の微分方程式と考えてよいからである。故に、その微分方程式の右辺はLipschitzの条件を満たさねばならないがそれだけでは十分ではない。

. . ［定理 1.1］. .積分演算子のノルムは区分幅hに等しい。
［証明］. .t₀<t_a<tであるt_aに対して、f(t)の積分がf(t_a)(t−t₀)に等しくなるようなt_aが存在する。このとき(1.7)を満たすようなある数Mが存在する。

. . 関数f(t)がy(t)の導関数のとき、上に述べた積分は定積分、即ち、y(t₀+h)−y(t₀)又はy(t₀)=0の場合はy(t₀+h)である。f(t)がy(t)のn次導関数でk=0からn−1に対してy^(k)(t₀)=0のとき、区間hの間でf(t)をn重積分する演算子のノルムは次のように導かれる。

. . ［定理 1.2］. .n重積分演算子のノルムはhⁿ/n!である。
［証明］. .f(t)のn重積分は(1.9)の全ての初期値を０にし、y⁽ⁿ⁾(t)をf(t)に置き換えた式で表される。

. .(1.9)

(1.7)及び(1.8)と同様に、. .

. .［証明終］

. . 微分の定義は下記の第１式で表され、導関数の積分は第２式で表される。何故なら、積分演算はh→0の逆演算を含むべきだからである。故に、積分演算子のノルムは積分区間hである。

. . ［定理 1.3］. .差分法の計算においてはn重積分演算子のノルムはhⁿである。
［証明］. .Eq.(1.9)のy⁽ⁿ⁾(t)を置き換えたf(t)のn重積分をNewtonの補間公式で表すためには全ての初期値は消去されねばならない。Eq.(1.9)は剰余項つきn−1次のTaylor展開に同じであるからn階の階差を求めれば全ての初期値は消去される。故に、Eq.(1.9)のn階の階差は剰余項の階差になり、以下のように求められる。

(1.10)

		. .(1.11)
∴. .	∴. . Minf=hn	[ Q.E.D.]

. . 差分法においてはEq.(1.10)により、n重積分は n-1次より高次の階差Δ^my(t₀)を求める。従って、y(t)を求めるためにはn次より低次の階差が初期値の代わりに与えられねばならない。Eq.(1.10)はまた積分毎に積分範囲がhづつ増加することも示している。
. . 被積分関数がf(t)=t³であるとき、これを階差で表すには４点の等間隔点を必要とする。f(t)の積分はもう一点を必要とし、その4hの範囲で積分演算のノルムはhである。従って、０から4h迄の区間のノルムはh×(4h)³=(4h)⁴/4である。区間が０から3hならノルムは(3h/4)×(3h)³=(3h)⁴/4である。故に、ノルムの計算は結果が正になることを除いて、関数の計算及び被積分関数に対する積分演算と同様である。
. . 著者の演算子法はこれらを行列の数式表現に明示する。例えば、(1.1)のベクトルFは区間2hを必要とし、２階以上の階差は無視し得ることを必要とする。ベクトルΔYは行列Dの第３行を省略しない表現ではΔ³y₀を持つから区間3hを必要とし、行列Dのノルムは１であるから積分のノルムはhである。

. . ［定理 1.4］. .微分方程式(1.3)の解は(1.12)が満たされるとき収束する。

(1.12)

(1.13)

(1.14)

　［証明終］

. . 条件式(1.12)にあるtの値は出発点t₀であってはならない。その理由は、(1.14)及び、y_i(t₀)=y(t₀)は固定点であることにある。tの値はt₀+hでなければならない。関数列y_i(t)を構成するため、及び、t₀とt₀+hの間に極大値が存在しても階差はその値を使用しないからである。
. . この収束条件はどのような数値解法でも満たされなければならない。例えば、Euler法は最初の近似解をy(t)=y₀と仮定してy₂(t)を得る。真の解y₁(t)とこれら２つの近似解の関係は(1.14)で表され、y₂(t)がより良い近似解であるためにはこの収束条件が満たされなければならない。
. . 積分演算子のノルムは使用されている最小の区間幅に等しい。従って、Runge-Kutta法のノルムは刻み幅の半分に等しく、従来行われているP-C法の結果との精度の比較の方法は公正でない。
. . 微分方程式(1.4)を解くプログラムで使用している区分幅h=0.05は(1.12)の条件h(2+2h)<1、即ち、h<0.366を十分に満たしている。微分方程式が(1.15)であるとき、収束条件はh(20−2h)<1、即ち、h<0.05となる。

1.6　硬い微分方程式と従来の問題点
. . 微分方程式(1.16)は硬い微分方程式といわれる。これは収束条件100h<1 を持ち、著者の方法ではh=0.01でも精度は極めて悪いが解を得る。しかし、h=0.015では解は得られない。

(1.20) (1.21)

(1.22) (1.23)

1.7　数値解の安定性
. . 数値解の不安定は幾つかの原因による。数値解の収束条件が十分に満たされないのはその原因の一つである。微分方程式(1.24)の数値解の収束条件は(1.25)である。従って、区分幅hは出発点t₀が０から遠くなるに従って小さくしなければならない。t₀が０に近いなら解の収束条件は(1.26)によりh<tである。これは臨海条件であり、かつ、y' (0)=0であるけれども解は容易に正確に求まる。それはy(0)が零でなく、y(h)もy' (h)も零でないからである。例えば、hが0.1なら１ステップの積分はt=10に近づくに従ってかなりの精度を失うので、全積分範囲が０から10までの場合には終端近傍では近似値が得られないこととなる。

(1.27)
(1.28)

(1.29)

y" +ty' +y=0,. . y(0)=10, y' (0)=0,. . h(t+h)<1. . (1.32)
(1.33)

Table 8. Valid h and True solutions.
	t=12.31400…	t=12.31415926535…	t=12.31500…
	4.48846403…×10−6	4.48774923…×10−6	4.48397781…×10−6
	1.18294852…×10−33	1.18063078…×10−33	1.16847043…×10−33

目　次　へ