お話：数値解析 第 4 回
数値積分
長田直樹

1

今回は補間型公式としてニュートン・コーツ公式、はじめに変数変換型公式として二重指数関数型公式を取り上定積分や無限積分の値を数値として求めることがげる。ロンベルク積分法は次回に話す。

数値積分である。通常の定積分 (特異点を持たない有限区間の積分) は、被積分関数 f (x) の原始関数 F (x)が分かれば、∫ bf (x)dx = F (b) − F (a)2ニュートン・コーツ公式

2.1

補間型公式

a

関数 f (x) は区間 [a, b] で連続であるとする。[a, b]により解決する。原始関数が分からなくても、MAT-の n + 1 個の点 x0 < · · · < xn における関数値
LAB や Octave などの数値解析ソフトで値を求める
f (x0 ), . . . , f (xn ) が与えられているとき、pn (xj ) =こともできる。さらに、大学 1-2 年程度の微積分学にf (xj ) (0 ≤ j ≤ n) を満たす高々n 次 (n 次以下のこ登場する積分のほとんどは、Mathematica や Mapleとを高々n次という) の多項式を f (x) のn次補間多などの数式処理ソフトを利用すると解決できる。た項式という。
∫bとえば、∫ ∞f (x)dx の定積分を
2ae−x cos xdx∫ b0pn (x)dxの値は容易に求めることができる。
aしかしながら、自然科学や工学に現れる積分は、原始関数が初等関数で表せないことや、表せても手で近似するのが、補間型数値積分公式である。補間型公式は間がかかることが多い。数値解析ソフトや数式処理

n∑

ソフトを用いる場合でも、プログラミングを行わなWj f (xj )j=0ければ解決できないケースもある。またこれらのソの形をしている。x0 , . . . , xn を分点、W0 , . . . , Wn はf (x) に無関係な定数で重みという。

フトが利用できないこともあるだろう。このような場合に利用するのが数値積分である。
数値積分公式は次のように分類することができる。
代表的な公式も挙げておく。

2.2

1. 補間型 (複合公式を含む)

ニュートン・コーツ公式[a, b] を N 等分しニュートン・コーツ公式、ガウス型公式

xj = a + j

2. 変数変換型

b−a,
Nj = 0, . . . , Nを分点にとる公式を (閉じた) N + 1点ニュートン・

IMT 公式、二重指数関数型公式

コーツ公式という。

3. 加速法の利用
ロンベルク積分法、広義積分に対するリチャー[a, b] を N + 2 等分し、両端を除いた

ドソン補外

xj = a + (j + 1)
1b−a,
N +2

j = 0, . . . , N

を分点にとる公式を開いたN + 1点ニュートン・コー

開いた 1 点ニュートン・コーツ公式を求める。

ニュートン・コーツ公式を実際上用いるのは、台形則、中点則、シンプソン則がほとんどであるので、これらのみを説明する。

p0 ((a + b)/2) = f ((a + b)/2) を満たす 0 次多項式は()
a+bp0 (x) = f2
だから

台形則

2.3

中点則

2.4

ツ公式という。

(∫ bp0 (x)dx = (b − a)f

閉じた 2 点ニュートン・コーツ公式を求めよう。
ap1 (a) = f (a), p1 (b) = f (b) を満たす高々1 次多項式はとなる。

(M = (b − a)ff (b) − f (a)(x − a) + f (a)b−a∫bである。定積分 a f (x)dx を
p1 (x) =∫ b
p1 (x)dx =a1(b − a) (f (a) + f (b))
2対し∫ b(1)aa<ξ<b∫ (b−a)/2 ()2b−a−t20))][ ((a+ba+b0000−t +f+ t dt× f22

命題 1 f (x) は [a, b] で C 2 級とする。台形則 T に対しaf 00 (ξ)(b − a)3 ,24証明

を台形則あるいは台形公式と呼ぶ。

f 00 (ξ)(b − a)3 ,12)となる ξ が存在する。

1(b − a) (f (a) + f (b))
2f (x)dx = T −(4)命題 2 f (x) は [a, b] で C 2 級とする。中点則 M にで近似する公式∫ ba+b2)を中点則あるいは中点公式と呼ぶ。

f (x)dx = M +T =a+b2a<ξ<bを 2 通りに変形する。[証明を完成せよ。] ¤を満たす ξ が存在する。

2.5

証明 部分積分を 2 回用いることにより∫ b(b − x)(x − a)f 00 (x)dx = 2T − 2a∫ bシンプソン則

閉じた 3 点ニュートン・コーツ公式を求める。

f (x)dx (2)p2 (a) = f (a), p2 ((a + b)/2) = f ((a + b)/2), p2 (b) =af (b) を満たす高々2 次多項式 p2 (x) の [a, b] での積分となる。一方、(b − x)(x − a) は [a, b] で非負だから、 値は積分の平均値の定理により∫ bb−aS=6(b − x)(x − a)f 00 (x)dx

a= f 00 (ξ)∫ b(b − x)(x − a)dx =a()a+bf (a) + 4f () + f (b)2である。[確かめよ。](b − a)3 00f (ξ)6(3)

S をシンプソン則あるいはシンプソン公式という。
台形則、中点則、シンプソン則の間にはS=となる ξ ∈ (a, b) が存在する。(2)(3) より求める結果が得られる。¤

2M + T3の関係がある。台形則と中点則の剰余項が消えるように重みをつけた平均になっている。

2

命題 3 f (x) は [a, b] で C 4 級とする。シンプソン則

S に対し∫ bf (4) (ξ)
f (x)dx = S −(b − a)5 ,2880aとなれば反復を終了し T2n を近似値に採用する。

Tn の分点は T2n の分点でもあるので、計算の無駄を省くため T2n の計算には Tn の値を利用する。
h = (b − a)/n とおくとn−1∑11
Tn = h  f (a) +f (a + jh) + f (b)
22j=1


2n−1
∑h 11hT2n =f (a) +f (a + j ) + f (b)2 222j=1

a<ξ<bを満たす ξ が存在する。
証明略

複合ニュートン・コーツ公式

3より、積分区間 [a, b] を n 等分し、各小区間にニュートン・コーツ公式を適用する公式を複合ニュートン・T2n =コーツ公式という。ニュートン・コーツ公式はほとんどの場合、複合公式として使われる。

nh∑
11
Tn +f (a + (j − )h)
22 j=12(6)と書ける。複合中点公式

3.1

複合台形公式

Mn = h[a, b] を n 等分し、h = (b−a)/n, xj = a+jh, (j =0, . . . , n) とおく。[xj , xj+1 ](j = 0, . . . , n − 1) に台形則を適用し、それらの和をとると複合台形公式

n−1
∑
11
Tn = h  f (x0 ) +f (xj ) + f (xn )
2
2
j=1n
∑

1f (a + (j − )h)2j=1を用いると、(6) は

T2n =1(Tn + Mn )2(7)と書ける。

が得られる。

3.2

複合台形公式の漸近公式が前回 (2008 年 7 月号) 話したオイラー・マクローリンの公式である。

複合シンプソン公式[a, b] を 2n 等分し、h = (b − a)/2n, xj = a +jh, (j = 0, . . . , 2n) とおく。[x2j , x2(j+1) ](j =2m+2関数 f (x) は区間 [a, b] で C級であるとする。 0, . . . , n − 1) にシンプソン則を適用し、それらの和∫ bをとると複合シンプソン公式
Tn −f (x)dxn−1na∑∑hm]∑Sn = f (x0 ) + 4f (x2j ) + f (x2n )
f (x2j−1 ) + 2B2k 2k [ (2k−1)
3=hf(b) − f (2k−1) (a) (5)j=1
j=1(2k)!
定理 1 (オイラー・マクローリンの公式)

k=1+O(h2m+2 ),(h → +0)が得られる。

Tn , Mn , Sn の間に証明 前号を見よ。

Sn =複合台形公式は、

2Mn + Tn
3(8)が成り立つ。さらに、(7)(8) より

T1 , T 2 , T 4 , T 8 , . . .と分割数を 2 倍づつ増やしてゆき、適当な終了条件

Sn =を満たすまで反復を繰り返す。たとえば、あらかじ

4T2n − Tn
3(9)が導ける。(5) と (9) より Sn の漸近公式が得られる。

め与えた ² > 0 に対し|Tn − T2n | < ²3

定理 2 (複合シンプソン公式の漸近公式)

に対する刻み幅 h の複合台形公式は

関数 f (x) は区間 [a, b] で C級であるとする。
n∑1
ck = B2k (f (2k−1) (b) − f (2k−1) (a))/(2k)! とおくと(13)f (a + jh)T (h) = h  f (a) +2j=1
h → +0 のとき∫ bm∑
4 − 22kである。n は (12) により定める。
f (x)dx =Sn −ck h2k + O(h2m+2 )3ak = 1, 2, . . . (あるいは最初のいくつかの k) に対し
k=2¯¯¯ B2k (2k−1)¯
証明 容易である。[導いてみよ。] ¤(2k−1)¯f(a) = 0, ¯f(a + nh)¯¯ < ² (14)(2k)!
複合台形公式と複合シンプソン公式の間の関係 (9)2m+2は次回ロンベルク積分にも登場する。

が成り立つとき、(13) は高精度の結果が期待できる。
例 1 f (x) = e−x cos x の [0, ∞) での積分に (13) を24急減衰する無限積分

4.1

全無限積分適用する。ライプニッツの公式より

f (2k−1) (0) = 0,k = 1, 2, . . .

がいえる。[確かめよ。] h = 1, ² = 10−16 とすると全無限積分∫ ∞∫ b|f (6)| + |f (7)| > ²,f (x)dx

f (x)dx = limb→∞a→−∞−∞

より n = 7 である。(13) は表 1 のようになる。h = 0.5に対する刻み幅 h の複合台形公式は

T (h) = hn∑

|f (7)| + |f (8)| < ²af (jh)で小数第 14 位まで正確な値が得られる。計算に使った分点数は 16(x = 0, 0.5, 1, 1.5, . . . , 6.5, 7 と 8) で(10)ある。

j=−mである。積分区間の両端の値は無視できるので、係表 1: I =数の 1/2 はつかない。m, n は与えた ² > 0 に対し、|f (−mh)| + |f (−(m + 1)h)| < ²(11)|f (nh)| + |f ((n + 1)h)| < ²(12)

h
1.000
0.500
0.250により定める。f (x) が振動しながら減衰する (0 に∫ ∞ −x2ecos xdx
0T (h)

T (h) − I

0.691021866829514
0.690194223521574
0.690194223521571

8.28 × 10−4
2.66 × 10−150

収束する) 関数でたまたま f (nh) = 0 などとなる場²∫∞合を除外するため、連続する 2 つの関数値を用いている。

0f (3) で近似しただけだが、誤差はなんと 8.3 × 10−4

k = 1, 2, . . . (あるいは最初のいくつかの k) に対し¯)¯¯¯ B2k ( (2k−1)(2k−1)¯¯<²f(nh)−f(−mh)¯ (2k)!
¯である。
図 1 に曲線と近似の折れ線を示す。区間 [0, 1] で過小に評価したものと [1, 2] で過大に評価したものがが成り立つときは、(10) はオイラー・マクローリン相殺し合ったということになるが、これほどの精度の公式により高精度の結果が期待できる。

4.2は図 1 からだけでは納得できないであろう。
なお、無限積分の値は√∫ ∞π −1/4−x2ecos xdx =e20半無限積分

半無限積分
∫ ∞

∫ bb→∞となることが留数定理を用いて証明できる。(一松信、f (x)dxf (x)dx = lim
a= 10−3 と す る と n = 3 で あ る 。
exp(−x2 ) cos xdx の値を 12 f (0) + f (1) + f (2) +講座複素解析、本誌 2007 年 12 月号 p.17 をみよ。)

a4

exp(−x x) cos(x)5.2

二重指数関数型公式

1

ある意味で最適の変数変換型公式が、1974 年
0.8に高橋秀俊と森正武が発表した二重指数関数型公式 (double exponential formula の頭文字を取り0.6DE 公式と略される) である。変換後の被積分関数
g(t) = f (φ(t))φ0 (t) が |t| → ∞ のとき二重指数関数 exp(− exp(|t|)) のオーダーで減衰するためその0.4ような名前がつけられている。((16) において α =−∞, β = ∞ としている。)

0.2

広義積分 (15) は変数変換 x = ((b−a)u+(b+a))/2000.51

図 1: ² = 101.5x−322.5により積分区間が (−1, 1) に移るので、u を x に置き3換えて∫ 1f (x)dx, h = 1 での複合台形公式

5二重指数関数型公式

5.1

変数変換型公式を考える。(17) に対し変換(π)
x = φ(t) = tanhsinh t ,2る公式を二重指数関数型公式という。(18) は DE 変

∫ b²→+0−∞ < t < ∞ (18)を行い、変換後の積分 (16) に複合台形公式を適用す関数 f (x) は a の近傍で有界ではないが極限値

I = lim(17)−1換という。

f (x)dx(18) の右辺に現れる双曲線関数はa+²が存在するとき、I を (a, ∞] における広義積分といet − e−t
et + e−t, cosh t =22sinh t
et − e−ttanh t == tcosh te + e−tsinh t =い、通常の定積分と同じ記号∫ bf (x)dxaにより定義される。(微積分学の教科書を見よ。)
で表す。a を特異点という。[a, b) についても同様で

φ0 (t) =ある。(a, c] と [c, b) における広義積分が存在するとき、それらの和を (a, b) における広義積分という。
広義積分

∫ bI=

なので、(17) に対する DE 公式は

f (x)dx(15)n(π))πh ∑ (cosh kh(π)
f tanhsinh kh22
cosh2
sinh khk=−m2(19)aは複合中点公式など開いた補間型公式で計算することはできるが、一般に収束は遅い。そこで、区間 (α, β)において連続微分可能で

lim φ(t) = a,t→αとなる。
被積分関数が 1 + x や 1 − x を因数に持つときは、

lim φ(t) = bt→β

1 ± tanh((π/2) sinh t)となる関数 φ(t) を用いて変数変換 x = φ(t) を行うと

∫ βI=cosh tπ)(2 cosh2 π sinh t2f (φ(t))φ0 (t)dtの計算において桁落ち (近い数同士の引き算で有効(16)桁が少なくなる現象) が生じやすい。そのため、(19)αのまま計算するのではなく、桁落ちを避ける変形をとなる。(16) に精度のよい数値積分公式を適用する

行った後にプログラミングする。

方法を変数変換型数値積分公式という。

5

例 2 ベータ関数∫ 1
B(p, q) =xp−1 (1 − x)q−1 dx,0

参考文献
(p, q > 0) (20) [1] 森正武、数値解析第 2 版、共立出版、2002

は 0 < p < 1 のとき x = 0 が特異点、0 < q < 1 の[2] 長田直樹、数値微分積分法、現代数学社、1987

とき x = 1 が特異点である。

(おさだなおき/東京女子大学)

B(1/4, 3/4) = 4.44288293815836610179の値を DE 公式で計算する。

x = 12 (u + 1) と変数変換すると() ∫ 113
1 3,
B=(1 + u)− 4 (1 − u)− 4 du
4 4−1である。f (x) = (1 + x)−3/4 (1 − x)−1/4 として (19)を適用する。

1 − tanh2 x =1,
cosh2 x

1 + tanh x =exp x
cosh xに注意すると

cosh((π/2) sinh x)
f (tanh((π/2) sinh x)) = √exp((π/2) sinh x)となるので、複合台形公式 T (h) はnπh ∑cosh kh√
2exp((π/2) sinh kh) cosh((π/2) sinh kh)k=−m(21)である。

² = 10−16 としたとき、(11)(12) で決まる m, n は
m = 5, n = 4 である。計算結果を表 2 に示す。
表 2: I = B(1/4, 3/4)

T (h)T (h) − I

1.000

4.445844600516824

2.96 × 10−30.500
0.250

4.442883163952324
4.442882938158366

2.26 × 10−70h

桁落ち防止 (21) を行わず、f (x) = (1+x)−3/4 (1−x)−1/4 に (19) を適用すると t = ±4 で f (φ(t))φ0 (t)が inf になり、計算できない。
無限積分に対する複合台形公式や、DE 公式は複素関数論を用いた誤差解析がなされている。開発者
自身の [1] を見よ。同書には無限区間の積分に対する DE 公式についても解説がある。

DE 公式の C 言語によるプログラム (例 2 で用いたもの) を

http://www.cis.twcu.ac.jp/~osada/rikei/
においておく。

6