2． 数値積分の誤差
2.1　Simpson の積分公式の精度
. . Simpson の 1/3 則積分公式の精度は O(h5)であるとされているが、これは積分結果が４次関数近似の精度であり、従って、被積分関数を３次関数で近似して積分したことに相当するということである。任意の３次関数の任意の区間[x0, x2]を求積法で積分した値と２等分点 x1 を用いて Simpson の積分公式で求めた値とは正確に一致するが、３次関数上のこの３点を通る２次関数の区間[x0, x2]を求積法により積分した値とも正確に一致する。従って、これら３点が任意の関数上の３点であれば、その関数を３次関数で近似しても２次関数で近似しても積分精度は同じであり、３次関数で近似する効果はない。更に、３次関数で近似しようとしても３点だけでは不定である。それなのに積分公式が被積分関数を３次関数近似した積分精度を持つというのは疑問である。
. . x0，x1=x0+h，x2=x0+2h において y0，y1，y2 である３点を通る２次関数を次式で表して２区間を積分する。

Eq2_1___________(2.1)
Eq2_2
一方、この３点を通る３次関数を f(x)=ax3+bx2+cx+k とすると、
Eq2_3
上記の３点を通る２次関数は一意に定まるが３次関数は無数に存在し、この区間の積分値は３点の関数値のみにより決まるから全て同じ値となる。従って、これら３次関数は区間[x0, x1]で２次関数の上側にあれば、区間[x1, x2]では下側にあり、その隔たりがどんなに大きくても２次関数の積分値に対する出入りは完全に相殺され、被積分関数を３次関数で近似する効果は全くない。この様な意味のない３次の項を加えても同じ積分値を与えるからといって Simpson の公式は３次関数近似の積分精度があるというのは誤りである。３次関数近似の積分精度を論ずるためには３次関数を一意に確定しなければならないのであり、そのためにはもう一点の関数値を確定しなければならない。その４点を用いた積分公式の精度が３次関数近似の積分公式の精度であり、Simpson の公式より精度は良いのである。この３次関数近似の積分公式の精度が O(h5)である。Simpsonの公式の精度は２次関数近似積分の精度でしかないのであり、O(h4)である。
. . 例えば、被積分関数が指数関数(2.4)であるとき、３点 (0, 1), (1/2, √2), (1, 2)を通るから区間[0, 1]は２次関数(2.5)で近似でき、この近似関数の積分値は求積法でも Simpson の公式でも1.442809であり、真の積分値は1.442695……である。
Eq2_4___________(2.4)
Eq2_5___________(2.5)
この３点を通る３次関数の１つである
Eq2_6___________(2.6)
Fig2_1
図 2.1
の求積法による積分値も1.442809であるから積分精度は同じである。しかし、３次関数の場合は前半[0, 1/2]の積分値は 0.9089045、後半[1/2, 1]の積分値は 0.5339045 であり、真値はそれぞれ 0.59758385…，0.84511118… であるから著しく異なり、それぞれの差は略相殺されている。２次関数の積分値との差は完全に相殺される。図2.1 に示す様に被積分関数自体も全く近似とは言えない。この様な３次関数は無数に存在し、これを一意に確定するには上記３点以外に任意の１点を通ることを指定しなければならない。この１点を(2.4)上の何処に取るかにより区間[0, 1]の被積分関数近似の精度は変わるが、中点を使う限り、その１点を等間隔区分点としてこの区間の内部にとることは出来ない。中点 x=0.5 の使用を止めて区間[0, 1]を３等分した４点を使用することは出来るが Simpson の公式はこの区間内の２点を使用出来ない。従って、この１点は区間[0, 1]の外にとらなければならない。
. . この第４点を区間[0, 1]の外、x=−0.5 における(2.4)上に取った３次関数は一意に確定し(2.7)で表される。
Eq2_7___________(2.7)
積分範囲を Simpson の公式に合わせるために、この関数の区間[0, 1]を積分すると下記に示すように Simpson の公式と同じ、即ち、第４点を使用しない(2.6)の積分結果と同じになる。これは(2.7)が(2.6)と同じ３点を通る無数に存在する３次関数の一つに過ぎないからである。第４点を確定した効果を得るにはこの点を含む区間[−0.5, 1]を積分しなければならない。
Eq2_7a
第４点を含む区間[−0.5, 1]を積分すると下記に示す値になり、(2.4)のこの区間を積分した値、即ち、真値は 1.8652506… である。Simpson の1/3則ではこの４点を使用した積分は出来ないので両端と中点を用いて積分すると下記の S の値となり、３次関数の積分値の方が精度が良い。区間[0, 1]の内部に３等分点をとった３次関数と Simpson の1/3則の比較の場合も全く同じである。Simpson の1/3則は３次関数近似の精度はなく、２次関数近似の精度 O(h4)である。
Eq2_7b
Eq2_7c
下記に示す４点の関数値を用いた Simpson の3/8則が上記３次関数近似の積分値に一致する結果を与えるのである。Simpson の1/3則が３次関数近似の積分精度を持つという従来の議論は間違いである。
Eq2_7d
. . 図 2.1 の区分点を負側に３点増やし積分区間を[−2,1]とすると、同じ区分幅で 1/3則は３回で積分できて 2.524916、3/8則は２回で積分できて 2.525159となり、真値は 2.5247163…であるから 1/3則の方が精度がよいが、この比較は 1/3則に有利な比較で正しくない。１回の積分範囲が狭ければ被積分関数値の変化が小さいから関数の近似精度が向上し、積分精度がよくなるのは当然である。この関数値の変化を同じにする、即ち、１回の積分範囲を同じにしなければ対等の比較ではない。故に、両者共区間[−2,1]を２回又は３回で積分すれば 3/8則の方の精度がよくなる。これは微分方程式の数値解法の精度比較の際にも注意すべき重要なことである。等間隔区分点の解を求めるのに、Runge K.法は１区分幅をその中点を用いて積分するが、これと比較するとき、Milne 法の修正子は Simpson の 1/3則に基づいて２区分幅を積分するから積分範囲が２倍であり、従来の精度比較は Runge K.法に有利な比較で正しくない。
. . Newton の前進補間公式を２区間積分すると Δ3f0の項は消えるから剰余項はf(4)(x)で決まるということを前提として Simpson の 1/3則の剰余項は証明される。しかし、３点のみを与えられたとき、Δ3f0 は定義できないのであり、定義できないものを有ると仮定して剰余項を高次に設定するのは誤りである。３点のみ与えられた場合は Δ2f0 迄しか定義できないから剰余項は f(3)(x)を用いなければならない。中央差分を用いるものも定義できない奇数次項が消えることで剰余項を高次に設定する誤りを犯している。又、求積法で証明するものも、３次関数は３点のみでは一意に確定できないにも関わらず、できると仮定して剰余項を f(4)(x)に設定している上に、被積分関数に平均値の定理を適用し、それを積分したものが正しい積分になるという間違った計算をしている。これらは皆同じ結論に達するが皆同じ誤りを犯しているからであって正しいからではない。Simpson の 1/3則の正しい剰余項は第３章にて示す。

2.2　Taylor 展 開
. . 常微分方程式の数値解の精度検討に於いて Taylor 展開及びその Lagrange の剰余がよく用いられるが、理論的には(2.8)の x−a を必要なだけ小さく仮定し、必要なだけの次数を仮定できるので問題はなくても、実際の数値計算ではできるだけ広い区間を低次の関数で近似することが計算の簡単化と累積誤差を避けるために必要であり、このような数値計算の精度を Taylor 展開で論ずることは適当でない。剰余項は更に高次の微係数と積分の平均値を与える点 x=ξ を必要とするが、実際にこれらを計算するのは殆ど不可能であり、剰余項は零として近似しなければならないからである。

Eq2_8___________(2.8)
y(x)のｎ次以下を yn(x)と置けば平均値の定理により区間 b−a=nhの積分誤差は、
Eq2_9___________(2.9)
(2.8)を見れば明らかなように、Taylor 展開の問題点は Lagrange の剰余が n+1 点で零になることが一般的には保証されないことにある。yn(x)は真値に対して常に剰余項に相当する差違があり、ｎ及び x−a が適当な条件を満たし、剰余項が必要な精度以内であれば真値と一致すると見做すことができることを保証するだけである。そしてｎが小さいときそのような条件を満たす x−a の範囲は極めて狭く、y(x)及び展開点 ａ により著しく異なる。一方、数値積分の公式や常微分方程式の数値解法では y(x)上の n+1 点の値を用いているから、被積分ｎ次近似関数はこれらの点で y(x)と一致することを前提としているのであって、ｎ及び x−a に関する条件に言及せずに、これら前提の成立しない Taylor 展開で論ずることは種々の問題を生ずる。
. . (2.4)を n=2, a=0 で Taylor 展開すると、剰余項を R(n+1)で表せば、
Eq2_9a
Eq2_10___________(2.10)
３点を通る２次関数(2.5)は、___ Eq2_10a
であるから係数は１桁しか一致しないし、g(x)は係数有効桁数と同じ７桁が真値と一致する２点(0.5, 1.414213), (1, 2)を通るが、y(x)の２次関数は真値と上位３桁しか一致しない２点(0.5, 1.406630), (1, 1.933374)を通り、近似精度が悪い。
. . (2.4)を n=3, a=−1/2 で Taylor 展開すると、
Eq2_11___________(2.11)
３点を通る３次関数(2.6)は、Eq2_11a
であるから高次の項の係数程不一致が著しく、f(x)は(0, 1)の他 g(x)と同じ２点を通るが、y(x)の３次関数はそれぞれ(0, 0.9995436), (0.5, 1.406349), (1, 1.956958)の３点を通り、やはり近似精度は悪い。
. . これら近似関数の 0.1 間隔の関数値を表 2.1 に、求積法による積分値を表 2.2 に示す。但し、上記 Taylor 展開の３次関数を区間[0, 1]で用いるのは定数項を見れば明らかなように、精度が悪すぎるので a=0 で展開した下式を用いている。
Eq2_12___________(2.12)

表 2.1_____近似関数値 (x の太字は区分点)
x Taylor 展開 y(x)の３点を通る 真__値
２次式３次式２次式３次式
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0 		1.0
1.071717
1.148238
1.229565
1.315695
1.406630
1.502370
1.602914
1.708263
1.818416
1.933374		 1.0
1.071772
1.148683
1.231063
1.319247
1.413568
1.514359
1.621952
1.736681
1.858878
1.988878		 1.0
1.069117
1.145097
1.227939
1.317645
1.414214
1.517645
1.627939
1.745097
1.869117
2.0		 1.0
1.071529
1.148312
1.230753
1.319253
1.414214
1.516036
1.625125
1.741880
1.866704
2.0		 1.0
1.0717734
1.1486983
1.2311444
1.3195079
1.4142135
1.5157165
1.6245047
1.7411011
1.8660659
2.0
. . Taylor 展開の２次関数は x=0.1 では y(x)上の３点を通る３次関数より精度がよいけれど x=0.4 では３点を通る２次関数より劣り急激に精度は悪化していく。更に、Taylor 展開の３次関数でも x=0.8 以上では３点を通る２次関数に及ばない。Taylor 展開は y(x)との一致点を展開点たゞ一点に集中することにより局所性を強くして展開点近傍の精度を向上させているのに対し、y(x)上の n+1 点を通るｎ次関数は Newton の前進補間公式に等しく、一致点を n+1 点に分散することで精度の平均化を図り、適用範囲を広く、かつ、明確にしたものである。これを示すのが両者の剰余項の違いである。
. . Taylor 展開の強い局所性は、展開点における関数の形を求めればよい数式計算では優れている。又、ただ一点の関数値を計算すればよい場合も展開点を近づけることができれば次数固定でも必要なだけの精度の解を得ることができるという利点もある。例えば、２次の Taylor 展開(2.10)で x=0.01 とすれば単精度で７桁正しい関数値を求めることができるし、更に x を小さくし２倍精度で計算すれば２倍精度の関数値を得ることもできる。常微分方程式の数値解法では Euler 法の精度がよい理由であるが、長区間の計算では区分点の数が膨大になること、一定区分幅では場所により局所性が著しく異なることにより大きな誤差を生ずる。被積分関数の高次微分が可能であれば広い区分幅でも精度のよい解を得ることはできるが、一般的には高次微分は簡単ではない。
表 2.2_____近似関数の積分値 (x の太字は区分点)
x Taylor 展開 y(x)の３点を通る 真__値
２次式３次式２次式３次式
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0 		0.0
0.103546
0.214504
0.333354
0.460577
0.596653
0.742063
0.897287
1.062806
1.239100
1.426649		 0.0
0.103547
0.214526
0.333466
0.460932
0.597520
0.743861
0.900619
1.068489
1.248204
1.440525		 0.0
0.103399
0.214052
0.332647
0.459869
0.596405
0.742940
0.900162
1.068757
1.249410
1.442809		 0.0
0.103534
0.214481
0.333385
0.460834
0.597451
0.743905
0.900901
1.069186
1.249546
1.442809		 0.0
0.1035472
0.2145263
0.3334708
0.4609524
0.5975838
0.7440217
0.9009699
1.0691829
1.2494690
1.4426950
. . Taylor 展開の欠点は高次近似式が簡単に使えないことと、単に O(hp+1)といっても実際の精度は ｈ 及び展開点 ａ により著しく異なり必要な精度を満たす範囲が明確でないことである。これを積分した場合にも同様のことが言える。表 2.2 は Taylor 展開の２次関数を積分すると y(x)の３点を通る２次関数の積分より精度のよい範囲が関数値の場合に比べ拡大することを示している。積分は平均化作用を持つので局所性が緩和されるからであり、(2.8)の剰余項と(2.9)を比較すると x−a>η−a であることがこれを示している。しかし、O(hp+1)を正確に計算するのは一般的には困難で、必要とする精度の保証される範囲は明確にできない。
. . これに対し、p+1 点を通る ｐ 次関数は Newton の ｐ 次補間公式に等しくその精度が保証される区間は[a, a+ph]であることが明確である。Newton の補間公式は ph=x−a を固定して h→0 にすれば Taylor 展開に収束するから、h→0 を仮定した理論的な精度の検討では両者は同じであるが、ｈ が有限の大きさである実際の数値計算の精度は大きな相違を生ずる。

2.3　局所相対誤差
. . 数値積分や常微分方程式の数値解法で用いられる積分公式は固定次数近似、固定区分幅であるのが常識的となっている。最近では自動区分幅調節の解法も行われてはいるが台形法や Simpson 等の低次の固定次数近似である。積分公式の精度は O(hp)と一般的に表されるが x の関数でもあり、被積分関数によっては x の影響の方が大きく１ステップの積分で有効桁がなくなるほどの誤差を生ずる場合がある。このような x の近傍では低次の固定次数近似の場合 h を極端に小さくする必要があり、一定区間の積分に要するステップ数が膨大となることにより大きな累積誤差を生ずる。
. . 連続関数 f(x)を任意の点 x0 から x1=x0+h 迄、台形法で積分する場合、f(x)の直線近似は、

Eq2_13___________(2.13)

この積分を z(x)で表し、真の積分を y(x)で表すと、
Eq2_14___________(2.14)

真値は、___________Eq2_15___________(2.15)
差分の相対誤差は、______Eq2_16___________(2.16)
関数値の相対誤差は、______Eq2_17___________(2.17)
. . 被積分関数が f(x)=−ae−x のとき、積分は微分方程式 y' =−y の解 y(x)=ae−x であり、
Eq2_18___________(2.18)

表 2.3___(2.18)の場合の z1 の局所相対誤差
y(真値)fz1(2.17)の値
x0=0
x1=0.1		 10
9.04837		−10
−9.04837
9.04758
−0.876278×10−4
x0=10
x1=10.1		 4.53999×10−4
4.10795×10−4		 −4.53999×10−4
−4.10795×10−4
4.10759×10−4
−0.876278×10−4

. . 従って、任意の x0 における近似積分 Δz0 の相対誤差 E は全区間に渡り一定で h のみにより決まり、h=0.1 ならば 0.8331945×10−3 であるから、Δz0 は全区間に渡って有効桁３桁となる。関数値 z1 の相対誤差は表 2.3 に示すように(2.17)により約 −0.105 倍になるから１桁精度は良くなる。固定次数固定区分幅解法でこの微分方程式が安定に解ける理由である。しかし、有効桁７桁の初期値で出発しても、次の区間は有効桁４桁の初期値で解かなければならないから台形法では精度はどんどん悪化していく。h=0.01 で E=0.8333319×10−5、関数値 z1 の相対誤差はこの約 −0.0105 倍になるから有効桁は７桁になる。従って、これより小さな区分幅にしなければならない。
. . 被積分関数が Eq2_18aのとき、積分は微分方程式 y' =−xy の解 Eq2_18bであるから、任意の x0 における真の積分 Δy0 及び近似積分 Δz0 の相対誤差 E は下記のようになり、h=0.1 のとき x0=0 における Δz0 の有効桁は３桁である。
Eq2_19___________(2.19)
関数値 z1 の相対誤差は下記に示すように E の約 −0.005 倍になり有効桁は５桁である。
Eq2_20___________(2.20)

表 2.4___(2.20)の場合の z1 の局所相対誤差
y(真値)fz1(2.20)の値
x0=0
x1=0.1		 10
9.95013		0
−9.95013×10−1
9.95025
1.25200×10−5
x0=10
x1=10.1		 1.92875×10−21
7.06007×10−22		 −1.92875×10−20
−7.13067×10−21
6.07840×10−22
−1.39046×10−1

しかし、x0=10 では E=0.08028511、z1 の相対誤差も表 2.4 の値であるから Δz0 も z1 も真値から出発して１ステップで有効桁１桁しか得られない。１区間の積分精度を有効桁７桁にするためには h=0.001 とする必要があり、このとき E=8.084153×10−6、z1 の相対誤差は −0.01005067×E=−0.8125117×10−7 である。僅か 0.1 区間の積分に 100 ステップを要し、その積分値は 7.05997×10−22 で有効桁は５桁に減少する。
. . Simpson の 1/3 則を用いる場合には２区間を１ステップで積分するので任意の x0 に対する積分値の相対誤差は次式で表される。
Eq2_21___________(2.21)

. . Eq2_21aの場合は、
Eq2_22___________(2.22)
従って、この場合は台形法と同様に全域に渡って相対誤差は同じで、h のみにより定まり、h=0.1 では有効桁７桁(−1.22855×10−7)である。
Eq2_22aの場合は、
Eq2_23___________(2.23)
従って、相対誤差は x0 の値により異なり、h=0.1 では x0=0 においては有効桁７桁である(−1.67501×10−7)が、x0=10 では１ステップで有効桁２桁(−3.03927×10−2)しかない。h=0.01 とすれば x0=10 において有効桁７桁に改善され、５ステップで 0.1 の区間の積分が出来、その積分値は 7.06007×10−22 で表 2.4 の y に一致する。即ち、台形法に比べ区分幅は 10 倍、ステップ数は 1/20 で台形法より精度の良い積分値が得られる。
. . 従来常識的となっている固定区分幅、固定次数の数値積分方式では、これ等の例に見られるように、常に正しい結果を得られることは保証されない。可変区分幅では著しいステップ数の増加を伴うので累積誤差を考慮しなければならない。それよりも可変次数とする方が精度の改善効果は顕著である。微分方程式の数値解法において、Runge K.法は安定であるが Milne 法は不安定であると言われるが無意味である。安定、不安定の問題ではなく、十分な精度の解を得るために必要な値に局所相対誤差を抑えることの出来る区分幅を用いているか否かの問題と、有効桁を失ったデータを生かしておく計算法に起因するものである。

目　次　へ