Introducing the integral operator (c) 2009 Akio Wakabayashi

4.5. .予測子と修正子
. . 微分演算子は微分方程式の解を得、その演算は(4.35)で表される。その演算は初期値を除いて(4.31)を用いて実行されねばならない。被積分ベクトル Y₀' が y₀' から y_p−1' 迄の p 個の値で構成されているなら、(4.28)に示すように、解 ΔY₀ は y₀ から y_p 迄の p+1 個の値により構成される。従って、最後の点 t_p における y_p' の値は結果に影響を与えない。その結果、被積分関数が p−1 次関数なら結果は剰余項のない p 次関数になる。しかし、被積分関数が任意の連続関数なら、y_p' の影響のある解の方が被積分関数が高次関数であるから y_p' を使用しない解よりよい精度を持つであろう。
. . 著者は Δ を掛けた積分演算子、即ち、(4.28)の hD⁻¹ を予測子と名付ける。その理由は、被積分関数の t₀ から t_p−1 迄を用いて t₀ から t_p 迄の積分を求めるからである。従来の予測子は y_p を y₀', y₁', ………, y_p−1' で表せば逐次の p の値に対して以下に示すように得られる。
. . 行列 D⁻¹ が１次、即ち、p=1 のときは、行列演算式は Δy₀=h[1]y₀' に同じであるから Euler 法である。行列 D⁻¹ が２次、即ち、p=2 のときは次に示すように中点法である。

行列 D⁻¹ が３次のときは次のようになる。
Eq4_e2

行列 D⁻¹ が４次のときは次に示す Milne の予測子になる。

. . 予測子 hD⁻¹ により積分されたベクトル Y₀ の最後の要素 Δ^py₀ を持たないサブベクトルを考えると、それは積分範囲に等しい区間の被積分関数を用いた t₀ から t_p−1 迄の積分を表す。故に、著者は行列 D⁻¹ の最後の行の全要素の値を零にしたとき、この演算子 hD⁻¹ を修正子と名付ける。従って、p 次の予測子に対する修正子は p+1 次行列のサブ行列であり、(p×(p+1))次の矩形行列で表される。従来の修正子は y_p を y₀', y₁', ………, y_p−1' 及び y_p' で表すことにより、逐次の p の値に対して以下のように得られる。
. . 行列 D⁻¹ が(1×2)次のとき Euler 法の修正子であり、次に示すように台形法である。
Eq4_e4

行列 D⁻¹ が(2×3)次のとき中点法の修正子であり、次に示すように Simpson の 1/3 則である。
Eq4_e5

Runge-Kutta 法は第２章 1.4 節に示したように Euler 法と中点法の二つの予測子とこの修正子からなる。
. . 行列 D⁻¹ が(3×4)次のとき Simpson の 3/8 則である。即ち、

行列 D⁻¹ が(4×5)次のとき Newton-cotes の５点法である。即ち、

. . これが Milne の予測子に対する修正子であるべきである。従来の予測子の修正子として Simpson の 1/3 を用いるべきというのは奇妙な考えである。われわれが顕微鏡視的な現象をより精密に観測するとき、拡大率を上げなければならない。同様に、より正確な結果を得るためには予測子より高次の修正子を用いなければならない。
. . 微分方程式 y' =−ty の Milne 法による解を反復修正毎にプリントアウトすると、それらは t=0 の近傍では１回ないし２回の反復修正で正確な解に漸近するが、t=5 の近傍では修正子を何回も反復しても正確な解に漸近しないことが分かる。その原因は t=0 の近傍では Simpson の 1/3 を修正子として使用することは十分であるが、 t=5 の近傍では十分でないことである。第２章１節及び２節で述べたように、中点法を用いることで十分であるから Milne の予測子を用いることは効果がない。
. . (4.31)において、演算子は p 次の正方行列であるからベクトル ΔY₀ 及び Y₀' は同数個の要素を持つ。修正子は t₀ から t_p 迄の積分を被積分関数の t₀ から t_p 迄を用いて求めなければならないから、最後の要素 Δ^py₀' をベクトル Y₀' に追加しなければならない。その場合、演算子は p+1 次の行列になり、ベクトル ΔY₀ は最後の要素 Δ^p+1y₀ を持つ。それは ΔY₀ に Δ^T を左より掛けることにより取り除かれ、ベクトル Δ^TΔY₀ は初期値零を第１要素として持つ。(4.33)に示したように、この演算は y' (t) の t から t+h 迄の積分を y' (t) の t₀ から t+h 迄の積分から減算することを意味し、その行列演算式は(4.34)で表される。初期値 y₀ を(4.33)の両辺に加えると解 y(t) が得られる。それは(4.35)で表され、また、次の式でも表される。

_______(4.39). . この演算式は(4.38)で表された要素を a_n で表すと詳細な行列演算子式で次のように表される。
Eq4_39a

両辺に Δ を左より掛けると、
Eq4_40

____(4.40). . この行列と h の積が p 次の予測子に対する修正子である。この行列は最後の対角要素が零であることを除き、 p+1 次の行列 D に同じである。最後の行の全要素が零であるから、この行列は最後の行の表記を省略して(p×(p+1))次の矩形行列で表してよい。(4.40)は(4.39)に演算子 Δ を左より掛けたものと等価である。故に、

________(4.41) 従って、厳密に言えば ΔY₀ を与える修正子は ΔΔ^ThD⁻¹ と表すべきであるが、積 ΔΔ^T は最後の対角要素が零である以外単位行列と同じであり、省略しても混乱は生じないから通常は省略してよい。
. . (4.35)の積分演算子 h(ΔD)⁻¹ は、(4.39)を詳細に表した演算式に示したように、t₀ から t_p 迄の積分を被積分関数の t₀ から t_p 迄を用いて求めるから厳密に言えば修正子である。しかし、最後の要素 Δ^py₀' が未知のときそれを零であると仮定すれば、その演算子は t₀ から t_p 迄の積分を被積分関数の t₀ から t_p−1 迄を用いて求め、この仮定は予測を意味する。故に、著者は通常予測子と修正子の区別なしに(4.31)及び(4.35)の両方を使用する。
. . 与えられた関数の積分を得たいとき、通常は p 次の予測子より精度のよい結果を与える(p×(p+1))次の修正子を用いるのがよい。与えられた微分方程式の解を得たいときは、以下に示すように予測子と修正子の対を用いることができる。
. . 本解法においては必要とされる出発値は初期値一つだけであり、先ず最初に与えられた初期値 y₀ により y₀' を求める。それから逐次に、１次の予測子、(1×2)次の修正子、２次の予測子、(2×3)次の修正子、………、p 次の予測子、(p×(p+1))次の修正子と用いて必要な精度の p 個の解を求める。最後の修正子以外の修正子は全て省略することができる。何故なら、k 次の予測子は((k−1)×k)次の修正子を含んでいるからである。
. . 最後の修正子は次数 p を固定したまま何回も反復繰り返し使用してはならない。その理由は、２ないし３回の繰り返しを除き、より良い精度の結果を得るうえで効果がなく、従来の予測子・修正子法のように不安定現象となる原因であるからである。

4.6　２重積分演算子
. . 導関数 y' (t) は y(t)−y₀ の微分によって得られ、２次導関数 y"(t) は y' (t)−y₀' の微分によって得られる。故に、

____(4.42) 両辺を t₀ から t 迄積分すると、
Eq4_43

___(4.43) これが(4.42)の逆演算であり一意である。初期値 y₀ 及び y₀' が与えられるなら、関数 y(t) は２階の微分方程式の解であり、次のように表される。
Eq4_44

_____(4.44) この２重積分演算子を導くために、先ず最初に y(t) 及び y"(t) の逐次階差を用いて(4.43)を行列演算式で表す。

. . ［定理 4.9］. .(4.43)から y₀ 及び y₀' を消去すると、

_____(4.45) ［証明］. .積分変数を x 及び z で表すと、
Eq4_45a

(4.43)から y₀ を消去すると、
Eq4_45b

次に y₀' を消去すると、
Eq4_45c

[Q.E.D.]

［補足］

. . ［定理 4.10］. .m≥0 のとき、Δ^my"(t) の t から t+h 迄の２重積分は、
Eq4_46

_____(4.46) ［証明］. .定理 4.9 により m=0 のとき成立する。m=k のとき成立すると仮定すれば m=k+1 のときも次のように成立する。
Eq4_46a

[Q.E.D.]

補足

. . ［定理 4.11］. . y(t)∈C^p+1[a, b] 及び 1≤m≤p+1 と仮定すると次のような ξ が存在する。
Eq4_47

______(4.47) ［証明］. . Δ^my(t) は y(t+kh) の減算で構成されるからΔ^my(t)∈C^p+1[a, b] が成立する。平均値の定理により、

1). .	m=1 のとき、. .
2)	m=2 のとき、. . 但し、. .
3)	m=n のとき成立すると仮定すると m=n+1 のときは、 . .[Q.E.D.]

. . ［定理 4.12］. . y(t)∈C^p+1[a, b] と仮定すると、Δ^my(t) は t=t₀+rh∈[a, b] を満たす実数 r に対して(4.48)で表される。
Eq4_48

___(4.48) ［証明］

1)	. .m=0 のとき、(4.48)は定理 3.1 に示した y(t) の Newton の補間公式である。
2)	. . m=1 のとき、定理 3.3 に示したようにその補間公式の階差である。Δy(t)∈C^p+1[a, b] であるが、剰余項が零になる等間隔区分点の数は r=0 から p−1 迄の p 点に減少する。
3)	. .m=n のとき成立すると仮定すれば、定理 4.11 により次のような値 ξ が存在する。これは関数 Δⁿy(t) についての Newton の p−n 次の補間公式の剰余項であり、r=0 から p−n 迄の p−n+1 点で零になる。故に、その補間公式及びそれをシフトした公式は次のように表せる。これらは定理 3.3 の証明の y(t) 及び p にそれぞれ Δⁿy(t) 及び p−n を代入した公式と同じものである。 . . 関数 Δⁿy(t₀+uh) と定数 R_q を持つ近似関数との残差として関数 F(u) を考えると、関数 G(u) 即ち F(u+1) は関数 Δⁿy{t₀+(u+1)h} と R_h の代わりに定数 R_q を持つ近似公式との残差である。G(u)−F(u) は u=0 から p−n−1 及び r の p−n+1 点で零になることを用いると定理 3.3 の証明と同様にして次の結果が得られる。定理 4.11 により、を満たす ξ_r が存在する。. .[Q.E.D.]

. . 上記証明の 3) は直接(4.48)の m に n を代入した式を用いて以下のように示してもよい。
Eq4_48e

但し、剰余 R は、平均値の定理により ξ_r<ξ<ξ_r+h を満たす ξ が存在するから次のようになる。 Eq4_48f

. . ２次微分演算の逆演算は(4.43)で表される。それを行列演算式で表すとき左辺の逐次階差を求めなければならない。t=t₀ なら左辺は y(t)−y₀−y₀' (t−t₀)=0 になるが、これは２重積分の結果ではない。t=t₀+h なら左辺は２重積分の結果であるが、左辺の２つの値の差は Δy₀−y₀' h で表されるから２重積分に関係のない値を含んでおり、ベクトル Y₀ の要素ではない。従って、これらは本演算子法においては逆演算の結果ではない。k≥2 なら左辺の第 k 階差は Δ^ky₀ になる。従って、逆演算の結果は Δ²Y₀ と表され、逆演算は(4.45)に等価である。

. .［定理 4.13］. . y"(t)∈C^p−1[a, b] のとき、それを Newton の補間公式で表すと、２重積分の結果の逐次階差は、

___(4.49)
［証明］. .y(t) と y"(t) の関係は(4.46)で表される。それに t=t₀ を代入すると、Δ^m+2y₀ が得られる。
Eq4_49a

積分変数を x=t₀+rh 及び z=t₀+qh で表すと、
Eq4_50

____(4.50) 条件 y"(t)∈C^p−1[a, b] が存在するから、定理 4.12 の y 及び p の表記をそれぞれ y" 及び p−2 に変えると、被積分関数 Δ^my"(t₀+rh) は次のようになる。
Eq4_50a

剰余項は導関数 y"(t) の p−1 次導関数、即ち、y^(p+1)(t) により決まり、変数 r に関して左辺が連続であるから連続関数である。従って、この式を(4.50)に代入すれば(4.49)が得られる。

[Q.E.D.]

. . ［定理 4.14］. .(4.49)の Δ^k+my₀" の係数にある２重積分を a_k²で表すと、
Eq4_51

____(4.51)
［証明］. .Stirling 数を用いると、
Eq4_51a

j=0 のとき、k=0 に対する Stirling 数は１であり、その他の k の値に対しては零である。

[Q.E.D.]

. . Table 3.3 の Stirling 数を用いると、(4.51)の値は k=0 から 10 迄に対して次のように求められる。
Eq4_51b

故に、(4.45)の Δ²y(t) と y"(t) の関係は(4.49)の m の値が零から p−2 迄の式を並べることにより次の行列演算式で表される。
Eq4_52

__(4.52) 但し、 Eq4_52a

. . m の値が m=p−1 のとき、(4.49)は剰余項のみとなる。
Eq4_53

_____(4.53) 平均値の定理を２度用いると、
Eq4_53a

故に、________ Eq4_53b

従って、それは行列演算式(4.52)に表されない。

. .［定理 4.15］. .h→0 なら(4.52)の剰余項は h^p+1R→0 である。
［証明］. .ベクトル h^p+1R を展開した行列の対角要素 h^p+1R_m は(4.45)の被積分関数を Newton の補間公式で表した式の t=t₀+mh における剰余項に同じである。(4.45)の積分変数 x, z を x=t₀+rh 及び z=t₀+qh となる r, q に変換すると、

u に m を代入すると、Δ²y_m は次のようになる。
Eq4_54

___(4.54) m≤p−2 且つ r≤q+1≤p であるから y^(p+1)(ξ_r) の係数は定理 4.6 により p 以下である。|y^(p+1)(t)|≤M となる値 M が存在するから、
Eq4_54a

[Q.E.D.]

. .従って、(4.52)の行列を D⁻² で表すと、２重積分は次のように表される。

_______(4.55) Δ²Y₀ に第 2.6 節の定理 2.12 に述べた Δ⁻² を左から掛けると、
Eq4_56

___(4.56) これが(4.44)で表される２重積分に対応する演算式である。故に、著者は演算子 h²Δ⁻²D⁻² を２重積分演算子と名づける。この演算子は二つの初期値 y₀, y₀' の代わりに y₀ 及び Δy₀ の二つの値が既知であることを必要とする。従って、この演算子により２階の微分方程式を解くためには初期値 y₀' は階差 Δy₀ に変換されなければならない。
. . 初期値 y₀' は微分 Y₀' =h⁻¹DΔY₀ の行列演算式の第１行で表される。即ち、
Eq4_57

_____(4.57) 従って、初期値 y₀' 及び(4.55)により得られるベクトル Δ²Y₀ の要素を用いてこの式から必要な値 Δy₀ は得られる。
. . 行列演算式(4.55)は(4.52)に示すように被積分関数の t₀ から t_p−2 迄を積分して t₀ から t_p 迄の結果を得る。故に、演算子 h²D⁻² は p−1 次の行列であり、t_p−2 から t_p 迄の２区間を予測する予測子である。一方、行列演算式(4.56)は被積分関数が y₀" から Δ^py₀" 迄の p+1 要素を持つベクトルで表され得るならその t₀ から t_p 迄を積分して t₀ から t_p 迄の結果を得る。行列 (Δ^T)²D⁻² は、第１行及び２行は全て零の要素を持ち、第３行から p+1 行迄は p+1 次の行列 D⁻² の第１行から p−1 行迄と同じである。(4.56)に Δ² を左から掛けると、
Eq4_58

_______(4.58) 演算子 (ΔΔ^T)²h²D⁻² は２重積分の修正子であり、演算子 (ΔΔ^T)² は最後の二つの対角要素が零であることを除いて単位行列と同じであるから修正子の最後の２行は全て零である。故に、通常は演算子 (ΔΔ^T)² は省略して、２重積分の修正子は((p−1)×(p+1))次の矩形行列で表される。

. .［定理 4.16］. .次の関係は剰余項に無関係に常に成立する。

D⁻²=D⁻¹×D⁻¹_________(4.59) ［証明］. .y"(t)∈C^p−1[a, b] のとき、0≤m≤p−2 に対する p−1 等間隔区分点 t=t₀+mh で関数 y"(t) に一致する p−2 次の関数 z"(t) を考えると、
Eq4_59a_j

関数 z"(t) の積分も２重積分も z^(p+1)(t)≡0 であるから剰余項を持たない。従って、(4.32)及び(4.45)により、
Eq4_59b

Δ²y(t₀+uh) の剰余項は u=m に対しては(4.54)のそれと同じであるから２番目の式は行列演算式(4.52)に等価である。故に、
Eq4_59d

_____________[Q.E.D.]
. . 従って、h²D⁻² による２重積分は hD⁻¹ による２回の積分と等価である。後者の場合、単積分のプログラムを２回使うことができ、初期値 y₀' を(4.57)により階差 Δy₀ に変換する必要はない。何故なら、階差 Δy₀ は ΔY₀=hD⁻¹Y₀' で表される２回目の積分で初期値 y₀' を用いて得られるからである。しかし、それは前者の場合の２倍の計算時間を必要とし、定理 4.16 が D⁻¹の要素の丸め誤差のため正確に成立しないことにより結果の精度が多少失われる場合もあり得る。

4.7　n 重責分演算子
. . 導関数 y' (t) は y(t)−y₀ の微分により得られ、２次導関数 y"(t) は y' (t)−y₀' の微分により得られる。同様に、n 次導関数は y⁽ⁿ⁻¹⁾(t)−y₀⁽ⁿ⁻¹⁾ の微分により得られる。従って、

___(4.60) 両辺を t₀ から t 迄逐次に積分すると、
Eq4_60a

____(4.61) これが(4.60)の逆演算であり一意である。n 個の初期値 y₀ から y₀⁽ⁿ⁻¹⁾ が与えられていれば、関数 y(t) は n 階の微分方程式の解であり、次のように表される。
Eq4_62

___(4.62) これは Taylor 展開と呼ばれている。剰余項は通常は以下に示す n 重積分の平均値の定理を用いて表される。この解はまた初期値を上記の各積分の右辺へ移項して次のように積分で表すこともできる。
Eq4_63

___(4.63)

. .［定理 4.17］. .関数が y(t)∈C^p[a, b] なら以下を満たす ξ が存在する。

_____(4.64) ［証明］. .p が１のときは良く知られているように成立する。p=n のとき成立すると仮定すれば、(4.63)を単積分の平均値の定理を用いて積分すると、
Eq4_64a_j

仮定及び(4.62)により、
Eq4_64b

p=n+1 の場合は、平均値 y⁽ⁿ⁾(ξ) は次のように表される。
Eq4_64c

これを(4.62)の n 重積分に代入すると、y⁽ⁿ⁺¹⁾(t) の n+1 重積分で表した解 y(t) に同じである。従って、y⁽ⁿ⁺¹⁾(t) の n+1 重積分の場合は平均値 y⁽ⁿ⁺¹⁾(ξ₁) が存在し、係数 φ は n/(n+1) である。何故なら、
Eq4_64d

故に、_________ Eq4_64e

. . n 重積分演算子を導くためには、(4.61)は先ず y(t) 及び y⁽ⁿ⁾(t) の逐次階差を用いた行列演算式で表されなければならない。従って、初期値は y(t) 及び y⁽ⁿ⁾(t) の両ベクトルの要素ではないから全て消去されねばならない。

. .［定理 4.18］. .(4.61)から全ての初期値を消去すると、

________(4.65) ［証明］. .n=2 のときは定理 4.9 の場合である。n=k のとき成立すると仮定すると、n=k+1 のときは次のようになる。
Eq4_65a

[Q.E.D.]

. . ［定理 4.19］. .整数 m が m≥0 なら Δ^m+ny(t) は次のようになる。
Eq4_66

________(4.66) [Proof]. .m=0 のときは定理 4.18 の場合である。m=k のとき成立すると仮定すると m=k+1 のときは次のようになる。
Eq4_66a

[Q.E.D.]

. . ［定理 4.20］. . y⁽ⁿ⁾(t)∈C^p−n+1[a, b] のとき、それを Newton の剰余項付き補間公式で表すと、t₀ から t 迄の n 重積分の逐次階差は次のようになる。但し、0≤m≤p−n である。
Eq4_67

____(4.67) ［証明］. .(4.66)の Δ^m+ny(t) に t₀ を代入すると、Δ^m+ny₀ であり、次のように表される。
Eq4_67a

積分変数 t を t=t₀+rh の変数 r に変えると、積分範囲 t₀ から t₀+h は r の 0 から 1 に変わり、積分範囲 t から t+h は r から r+1 に変わる。
Eq4_68

________(4.68) 条件 y⁽ⁿ⁾(t)∈C^p−n+1[a, b] があるから定理 4.12 の y 及び p の表記をそれぞれ y⁽ⁿ⁾ 及び p−n に替えることにより、被積分関数は次のように表される。
Eq4_68a

剰余項は導関数 y⁽ⁿ⁾(t) の p−n+1 次導関数、即ち、y^(p+1)(t) により決まり、上式の左辺は連続であるから剰余項は r の連続関数である。従って、上式を(4.68)に代入すれば(4.67)を得る。

[Q.E.D.]

. . ［定理 4.21］. .整数 m 及び n が 0≤m<n を満たすとき、
Eq4_69

_____(4.69) ［証明］. .1). .m=0 のとき式[1]は成立する。何故なら、
Eq4_69a

2). .m=1 のときも成立する。何故なら、
Eq4_69b

3). .m=k のとき成立すると仮定すると、m=k+1 の場合は、
Eq4_69c

故に、式[1]は常に成立する。式[2]の左辺を s_n と置くと、
Eq4_69d

. . ［定理 4.22］. .(4.67)の Δ^k+my₀⁽ⁿ⁾ の係数を a_kⁿ と表すと、
Eq4_70

_____(4.70) ［証明］. .二項係数を Stirling 数 α_j^(k) を用いて表すと、
Eq4_71

______(4.71) 任意の整数 j≥0 に対して、
Eq4_71a

故に、r^j の r から r+1 迄の n 重積分が
Eq4_72

_____(4.72) と表されると仮定すると、n+1 重積分は次のようになる。
Eq4_72a

故に、(4.72)は整数 n≥1 に対して常に成立する。従って、(4.71)を n 重積分して最外側の積分の積分範囲の r に零を代入すると証明は完了する。但し、整数 j は Stirling 数により k=0 の場合にのみ零に成ることができ、(4.72)の総和は r=0 の場合には定理 4.21 により n! になるから a₀ⁿ は１になる。

[Q.E.D.]

. . 定理 4.20 及び 4.22 により、導関数 y⁽ⁿ⁾(t) のベクトルと Δⁿy(t) のベクトルの関係は(4.67)の 0≤m≤p−n の各式を並べて次の行列演算式で表すことができる。

__(4.73)
但し、. . Eq4_73a

. . 左辺のベクトルの要素 Δ^p+1y₀ は(4.67)の総和が m=p−n+1 のときは零になるから次のように剰余項のみで表される。
Eq4_74

______(4.74) これはベクトル Y₀⁽ⁿ⁾ の要素に無関係であるから行列演算式には現れない。定理 4.11 の t=t₀ の場合により、右辺の n 重積分が y^(p+1)(ξ) となるような値 ξ が存在する。それは n 重積分の平均値の定理によっても次のように示される。
Eq4_74a

従って、関数 y(t) 及び y⁽ⁿ⁾(t) がベクトルで表されるとき Δ^p+1y₀ は消える。

. .［定理 4.23］. .両関数 y(t) 及び g(t) が連続で、t₀ から t の区間で符号を変えないなら、次の関係を満たす値 ξ が存在する。

____(4.75) ［証明］. .g(t)>0 及び L≤y(t)≤G と仮定すると、
Eq4_75a

L=G の場合は(4.75)は明らかに成立する。L≠G の場合は、p=1 なら、
Eq4_75b

故に、

中間値の定理により次のような ξ が存在する。
Eq4_75d

. . (4.75)が p=n のとき成立すると仮定すると、p=n+1 のときは、
Eq4_75e

故に、

中間値の定理により次を満たす ξ₁ が存在する。
Eq4_75g

g(t)<0 が負の場合も同様にして証明される。______[Q.E.D.]
. . 剰余項 Δ^mR₀ の一番内側の積分を f(r) で表すと、定理 4.17 により次式を満たす u が存在する。
Eq4_75h

関数 f(u) は u<j<u+1 を満たすある整数 j において被積分関数の二項係数の符号を変える。定理 4.23 により、二つの値 ξ₁, ξ₂ が存在して次式が成立するする。
Eq4_75i

u<u₁<j, j<u₂<u+1 を満たす二つの値 u₁, u₂ が存在して次の関係が成立する。
Eq4_75j

従って、h^p+1|Δ^mR₀|→0 が成立し、定理 1.7 により剰余項は、
Eq4_75k

. . 故に、(4.73)の行列を D⁻ⁿ で表すと、この行列演算式は次のように表される。
Eq4_76

______(4.76) これは(4.65)で表される n 重積分の演算子法表記である。ΔⁿY₀ に左より第 2.6 の定理 2.12 で述べた Δ⁻ⁿ を掛けると、
Eq4_77

_____(4.77) これは(4.62)で表された n 重積分に対応する演算子法表記である。故に、著者は hⁿΔ⁻ⁿD⁻ⁿ を n 重積分演算子と呼ぶ。この演算子は y₀ から y₀⁽ⁿ⁻¹⁾ 迄の n 個の初期値の代りに y₀ から Δⁿ⁻¹y₀ 迄の n 個の値が既知であることを必要とする。従って、n 階の微分方程式をこの演算子で解くためには、初期値 y₀' から y₀⁽ⁿ⁻¹⁾ は Δy₀ から Δⁿ⁻¹y₀ 迄の逐次階差に変換されなければならない。
. . この変換は著者がこの演算子法を発見した切っ掛けであった。この変換はこの演算子法に類似の行列演算式で表されるが、その行列はこの演算子法の演算子ではない。何故なら、その被演算子は要素が逐次階差ではなく、初期値点における逐次微分値であるベクトルだからである。この行列は通常の数値計算をこの演算子法に結び付ける演算子である。その詳細は第四章で述べる。
. . (4.76)の演算子 hⁿD⁻ⁿ は行列 D⁻ⁿ が p−n+1 次であるから t₀ から t_p−n 迄の被演算子を用いて t₀ から t_p 迄の n 重積分を得る。故に、これは t_p−n から t_p 迄の n 重積分を予測する予測子である。(4.77)の行列 D⁻ⁿ が p+1 次のとき、(4.77)に左から Δⁿ を掛けると、
Eq4_78

________(4.78) この演算子は t₀ から t_p 迄の被演算子を用いて t₀ から t_p 迄の n 重積分を得るから n 重積分の修正子である。演算子 (ΔΔ^T)ⁿ は最後の n 個の対角要素が零であることを除いて単位行列に同じである。故に、通常はこの演算子を省略し、n 重積分の修正子は ((p−n+1)×(p+1)) 次の矩形行列で表す。

. .［定理 4.24］. .次の関係は剰余項の有無にかかわらず常に成立する。

_______(4.79) ［証明］. .y⁽ⁿ⁾(t)∈C^p−n+1[a, b] のとき、関数 y⁽ⁿ⁾(t) と 0≤m≤p−n に対する t=t₀+mh の p−n+1 点で等しい p−n 次関数 z⁽ⁿ⁾(t) を考えると、
Eq4_79a

z^(p+1)(t)≡0 であるから z⁽ⁿ⁾(t) の逐次積分は剰余項を持たない。故に、定理 4.18 により z' (t) の n−1 次導関数のn−1 重積分及び y⁽ⁿ⁾(t) の n 重積分は次のようになる。
Eq4_79b

逐次階差を計算すると、定理 4.19 及び 4.20 により次の行列演算式で表される。
Eq4_79d

２番目の式は剰余項も含めて(4.73)に等しい。故に、
Eq4_79e

同様に、Δz⁽ⁿ⁻¹⁾(t₀+rh) 及び Δⁿy(t₀+rh) を使えば、
___________ Eq4_79f

________[Q.E.D.]
. . 従って、n 重積分は単積分プログラムを n 回使って行っても良い。その場合は y₀' から y₀⁽ⁿ⁻¹⁾ 迄の初期値を Δy₀ から Δⁿ⁻¹y₀ 迄の逐次階差に変換する必要はない。何故なら、逐次積分は次のようになるからである。
Eq4_80

_____(4.80) 従って、その方法は簡単であり、便利である。しかし、計算時間は n 重積分演算子を使用する場合のおよそ n 倍であり、D⁻¹ の n 乗迄の数値計算の丸め誤差により結果の精度が悪化する場合もあり得る。

目次へ