9. 数値積分 (1) 補間型積分公式

目次

数値積分 (1): 概要, 補間型積分公式, ロンバーグ積分, 数値実験, 積分公式の係数表
数値積分 (2): 台形則の誤差, 変数変換型積分公式, 誤差の推定, 自動積分, 積分公式の係数表
数値積分 (3): 実用ルーチンのベンチマーク (実用ルーチンの選択, 性能比較)

概要補間型積分公式ロンバーグ積分数値実験積分公式の係数表

概要

有限区間あるいは無限区間における定積分の近似計算方法を数値積分法という.

次の積分を数値計算で求める例を見てみる.

  
  I = ∫ excos(x)dx [0.2, 1] = 1.1582

 
最も簡単な方法は, 下図の網掛けのように積分範囲の下端と上端を結ぶ直線で囲まれる領域で近似する方法である. これを S で表すと, S = (f(a) + f(b))(b – a)/2 となり, これは台形の面積でこの方法は台形則とよばれる. 実際に f(x) = excos(x), a = 0.2, b = 1 として計算すると S = 1.0663 (誤差 -7.9%) となる.

次に, もっと精度をよくするために a と b に加え中点 x1 = (a + b)/2 の3点を使って, 下図の網掛けのように積分範囲の下端と上端を結び中点を通る二次曲線で囲まれる領域で近似する方法で計算してみる. S = (f(a) + 4f(x1) + f(b))(b – a)/6 となる. これはシンプソン則とよばれる. 計算値は S = 1.1575 (誤差 -0.06%) となった.

一般的に, 数値積分においては積分区間 [a, b] を n分割する適当な分点列 a = x0 < x1 < ⋯ < xn = b をとり, 分点(標本点ともいう)における関数値f(xi)の重み付きの和を積分の近似値とする.

  
  I = ∫ f(x)dx [a, b] ≒ In = Σwif(xi) (i = 0~n)

 
ここで wi は xi における重みである.

分点xiにおいて関数値f(xi)をとる補完多項式を求めその積分という形で上の公式を作るものを補間型積分公式とよぶ. 上の台形則とシンプソン則は補間型積分公式の代表的な例である. 他にこれに属するより高次な公式も使われている.

積分区間を複数の小区間に分割し, その小区間ごとに公式を適用するものを複合公式という. 例えば, 積分区間を下図のように4分割しそれぞれに台形則を適用すれば単体の台形則より高精度が期待できそうである.

これを複合台形則というが, 台形則の場合には実用上複合公式を使うのが普通であるためこれを単に台形則ということが多い(シンプソン則も同様に複合シンプソン則を指すことがある). m分割したときの複合台形則の計算式は次のようになる.

  
  Im = h(f(x0)/2 + f(x1) + f(x2) + ・・・ + f(xm)/2), h = xi+1- xi (i = 0~m)

 
h は小区間の幅である. ここでは h を一定 (= (b – a)/m) としているが可変にする方法も考えられる. 実際に4分割した場合の計算をしてみると S = 1.1523 (誤差 -0.5%) となり精度は大幅に改善されたが5点を使っているにもかかわらず3点しか使っていないシンプソン則におよばない.

それではシンプソン則よりもさらに3次, 4次, … と高次の公式を使えばよいかというと, 場合によっては桁落ちによりむしろ精度が悪くなる恐れがある. そのため適当な次数に抑えた公式を適用した複合公式を使用するのが一般的である.

数値積分を効率よく行うために種々の技法とこれらの公式を組み合わせて使うことも多い. 例えば, 加速法と台形則を組み合わせるロンバーグ積分, 変数変換と台形則を組み合わせる変数変換型積分公式, 複合公式(区間分割)と高次公式を組み合わせる適応自動積分などがある. 変数変換型積分公式および適応自動積分については「数値積分(2)」で説明する.

実用の観点からは, 関数値がとびとびの点(等間隔, 不等間隔)におけるデータとしてしか与えられない場合と, 任意の点における関数値が計算できるサブルーチンとして与えられる場合とで使用できる公式やプログラムが異なる. ここでは後者を想定している. 前者については「6. 補間法」で説明されるスプライン補間が参考になるかもしれない.

補間型積分公式

積分区間 [a, b] を n分割する分点列を a = x0 < x1 < ⋯ < xn = b, 各分点における重みを wi, 分点における関数値を f(xi) として積分公式を以下に再掲する.

  
  I = ∫ f(x)dx [a, b] ≒ In = Σwif(xi) (i = 0~n)

 
分点が積分区間の端点に一致しない公式の場合には, N個の分点 a < x1 < x2 < ⋯ < xN < b により次のように表すことにする.

  
  I = ∫ f(x)dx [a, b] ≒ IN = Σwif(xi) (i = 1~N)

 
補間型積分公式では重みと分点のとり方により次のように種々の公式が導かれる.

公式 nまたはN 重みwi 分点xi 参考
台形則 n = 1 w0 = w1 = (b – a)/2 x0 = a, x1 = b ニュートン・コーツ公式で n = 1 とした場合に一致.
シンプソン則 n = 2 w0 = w2 = h/3
w1 = 4h/3
x0 = a, x1 = a + h, x2 = a + 2h = b
(等間隔: h = (b – a)/2)
ニュートン・コーツ公式で n = 2 とした場合に一致.
ニュートン・コーツ公式 n = 1 ~ 7, 9 有理数(分数)で表される.
n = 8 および n ≧ 10 では負の値が現れ桁落ちの原因になるので使わないほうがよいとされる.
x0 = a, x1 = a + h, … , xn = b
(等間隔: h = xi+1– xi)
nが奇数のときはn次, nが偶数のときはn+1次までの多項式について正確な積分値を与える.
n = 3 の場合をシンプソンの3/8則ともよぶ.
チェビシェフ公式 N = 2 ~ 7, 9 一定. a < x1 < … < xN < b (不等間隔) N = 8 および N ≧ 10 では分点の値が複素数になるため通常は使われない.
ガウス型公式 N = 2 ~ できるだけ精度が高くなるように分点と重みの両方を選ぶため, 重みは一定ではなく分点は不等間隔 (a < x1 < … < xN < b)
ガウス・ルジャンドル公式(ガウス公式), ガウス・ラゲール公式, ガウス・エルミート公式などの種類がある.
2N – 1 次までの多項式について正確な積分値を与える.
クレンショウ・カーチス公式 n = 2 ~ チェビシェフ級数展開より定める. xi = (((b – a)/2)cos(πi/n) + (a + b)/2 (i = 0~n) 特異性の強い関数への適用可能性が特長.

ニュートン・コーツ公式

区間[a, b]をn等分して分点を等間隔にとる方法である. n = 1 の場合は台形則, n = 2 の場合はシンプソン則とよばれる.

分点は a = x0, x1, ・・・, xn = b で間隔 h = (b – a)/n である. すなわち, xk = a + kh となる.

被積分関数f(x)を近似するラグランジュ補間多項式(*) pn+1(x) は次のように与えられる.

  
  pn+1(x) = ΣLk(n)(x)f(xk) (k = 0~n)

 
(*) ラグランジュ補間多項式については「6. 補間法」を参照せよ.

これを積分すると,

  
  In = ∫ pn+1(x)dx [a, b]
  = Σ(∫ Lk(n)(x)dx [a, b])f(xk) (k = 0~n)
  = (b - a)/n Σwkf(xk) (k = 0~n)
  ただし,
    wk = n/(b - a) ∫ Lk(n)(x)dx [a, b]

 
と積分公式の形になる. これをn次のニュートン・コーツ公式という. この公式はn次までの多項式に対して正確な値を与える.

重み wk は多項式の計算により有理数(分数)で求めることができる.

n = 1~10 のときの重みの値を「積分公式の係数表」に掲載した. 表のように n = 8 および n = 10 (実際には n ≧ 10) では負の値が現れる. これは桁落ちの原因になるので n = 8 および n ≧ 10 の公式は使わないほうがよいとされる.

一般的には, nを大きくして精度を上げるのではなく, 積分区間を小区間に分割してそれぞれに公式を適用するニュートン・コーツ複合公式が使用される. 例えば, n = 2 (シンプソン則)の場合, [a,b]をm個の小区間に分割し, 各小区間でシンプソン則を適用すると次の計算式が得られる.

  
  I2m = h/3(f(x0) + 2Σf(x2j) + 4Σf(x2j-1) + f(x2m)) (最初のΣは j = 1~m-1, 2番目のΣは j = 1~m)
  ただし h = (b - a)/2m, xj = a + jh (j = 0, 1, ..., 2m)

 
これは複合シンプソン公式とよばれる.

チェビシェフ公式

ニュートン・コーツ公式とは逆に重みが一定になるように分点を選ぶ方法である.

チェビシェフ公式は次のように表される.

  
  IN = (b - a)/N Σf((a + b)/2 + (b - a)xk/2) (k = 1~N)

 
xk は分点で, 多項式 pN(x) = Π(x – xk) (k = 1~N) = ΣaixN-i (i = 0~N) をゼロとおいた方程式のN個の解として求められる. 展開した多項式の係数 ai は次のように与えられる(導出の詳細は文献[2]を参照).

  
  a0 = 1
  ai = -N/i Σai-j/(j + 1) (Σは j = 2, 4, ..., i) (i = 2, 4, ..., N2)
  ai = 0  (i = 1, 3, ..., N1)

  ただし, N1 = N - 1, N2 = N (Nが偶数のとき), N1 = N, N2 = N - 1 (Nが奇数のとき)

 
N = 2~7, 9 のときの分点の値を「積分公式の係数表」に掲載した. N = 8 と N ≧ 10 は方程式の解に複素数が含まれ, 計算が複雑になるうえ精度の点でも不利なため通常は使われない.

ガウス型公式

できるだけ精度が高くなるように分点と重みの両方を選ぶ方法である. 被積分関数を直交多項式補間(*)することにより求められ, 2N – 1 次までの多項式に対して正確な値を与える.

(*) 直交多項式および直交多項式補間については「6. 補間法」を参照せよ.

区間[a, b]における f(x)w(x) の積分を考える. w(x) は密度関数である.

  
  I = ∫ f(x)w(x)dx [a, b]

 
f(x) を直交多項式補間 fN(x) で置き換えると次のようになる.

  
  IN = ∫ fN(x)w(x)dx [a, b]
    = Σ(1/λj){Σwkpj(xk)f(xk) (k = 1~N)} ∫ pj(x)w(x)dx [a, b] (j = 0~N-1)

 
pj(x) の直交性より上式の積分の部分は次のように簡単になる.

  
  ∫ pj(x)w(x)dx [a, b] =
    λ00 (j = 0)
    0       (j ≠ 0)

 
すなわち, IN の式において外側のΣは j = 0 以外の項は0になり, また p0(x) = μ0 であるから, 次のガウス型積分公式が得られる.

  
  IN = (1/λ0)(λ00)Σwkp0(xk)f(xk) (k = 1~N)
     = Σ wkf(xk) (k = 1~N)

 
ここで, 分点 xk は直交多項式のゼロ点である. 重み wk は直交多項式補間の wk に等しく, 次式で与えられる.

  
  wk = μNλN-1/(μN-1pN-1(xk)p'N(xk))

 
直交多項式の種類に応じていくつかの積分公式があり, 代表的なものは次のとおりである.

直交多項式 積分公式 区間 [a,b] 密度関数 w(x)
ルジャンドル多項式 Pn(x) ガウス・ルジャンドル公式 (ガウス公式) [-1, 1] 1
ラゲール多項式 Ln(x) ガウス・ラゲール公式 [0, +∞] exp(-x)
エルミート多項式 Hn(x) ガウス・エルミート公式 [-∞, +∞] exp(-x2)

ガウス・ルジャンドル公式 (ガウス公式)

ルジャンドル多項式を使用して, pn(x) = Pn(x), w(x) = 1, [a, b] = [-1, 1] とするとガウス・ルジャンドル公式が得られる. ガウス・ルジャンドル公式は単にガウス公式とよばれることが多い.

  
  I = ∫f(x)dx [-1, 1] ≒ IN = Σwkf(xk) (k = 1~N)

 
ただし, xk は Pn(x) = 0 のN個のゼロ点である. wk は, λN = 2/(2N + 1), μN = Π(2k – 1)/k (k = 1~N) より, 次のように求めることができる.

  
  wk = 2/(nPN-1(xk)P'N(xk)) = 2(1 - xk2)/(nPN-1(xk))2

 
N = 2~10 のときの分点および重みの値を「積分公式の係数表」に掲載した.

任意の区間[a, b]の積分を求めるときは次のように変数変換して計算するとよい.

  
  I = ∫f(x)dx [a, b] ≒ IN = ((b - a)/2)Σwkf((a + b)/2 + ((b - a)/2)xk) (k = 1~N)

 

ガウス・ラゲール公式

ラゲール多項式を使用して, pn(x) = Ln(x), w(x) = exp(-x), [a, b] = [0, +∞] とするとガウス・ラゲール公式が得られる.

  
  I = ∫f(x)exp(-x)dx [0, +∞] ≒ IN = Σwkf(xk) (k = 1~N)

 
ただし, xk は Ln(x) = 0 のN個のゼロ点である. wk は, λN = 1, μN = (-1)n/n! より, 次のように求めることができる.

  
  wk = -1/(NLN-1(xk)L'N(xk)) = xk/(n2 LN-1(xk)2)

 
N = 2~10 のときの分点および重みの値を「積分公式の係数表」に掲載した.

ガウス・エルミート公式

エルミート多項式を使用して, pn(x) = Hn(x), w(x) = exp(-x2), [a, b] = [-∞, +∞] とするとガウス・エルミート公式が得られる.

  
  I = ∫f(x)exp(-x2)dx [-∞, +∞] ≒ IN = Σwkf(xk) (k = 1~N)

 
ただし, xk は Hn(x) = 0 のN個のゼロ点である. wk は, λN = π1/22nn!, μN = 2n より, 次のように求めることができる.

  
  wk = π1/22n(n - 1)!/(HN-1(xk)H'N(xk)) = π1/22n-1(n - 1)!/(nHN-1(xk)2)

 
N = 2~10 のときの分点および重みの値を「積分公式の係数表」に掲載した.

クレンショウ・カーチス公式

滑らかな関数の積分

被積分関数 f(x) をチェビシェフ級数展開し, 項別に積分する方法である. チェビシェフ級数の収束が速い関数であれば少ない関数計算回数で積分を求めることができる.

区間[-1, 1]で滑らかな関数 f(x) をチェビシェフ級数展開する.

  
  f(x) = (1/2)a0 + ΣakTk(x) (k = 1~∞)
  Tk(x) = cos(kθ), x = cos(θ)

  ak = 2/π∫f(cos(θ))cos(kθ)dθ [0, π]

 
係数 ak の積分を台形則で計算すると次のようになる.

  
  ak(n) = (1/n)f(1) + (2/n)Σf(cos(jπ/n))cos(kjπ/n) + (-1)k(1/n)f(-1) (Σ は j = 1~n-1)

 
右上添字の(n)は, n小区間の複合台形則を使ったことを示す. この係数はFFT(コサイン変換)ルーチンを使って効率的に計算できる(XLPack であれば Cost1f が使える).

f(x) のチェビシェフ級数展開を第n項までで打ち切って積分を次のように近似する.

  
  I = ∫f(x)dx [a, b]
  ≒ In = (1/2)a0(n)∫T0(x)dx[-1, 1] + Σak(n)∫Tk(x)dx[-1, 1] + (1/2)an(n)∫Tn(x)dx[-1, 1] (Σ は k = 1~n-1)

 
このように項別の積分の形になる. チェビシェフ多項式の項の積分は次の関係を使って計算することができる.

  
  ∫Tkdx =
    T1(x)  (k = 0)
    (1/4)(T2(x) + 1)  (k = 1)
    (1/2)(Tk+1(x)/(k + 1) - Tk-1(x)/(k - 1))  (k ≧ 2)
  (積分定数は省略)

 
以上の関係式を使って f(x) の[-1, 1]における積分を計算することができる(この場合, kが奇数の ∫Tkdx [-1, 1] は 0 になる).

任意の区間[a, b]の積分を求めるには次のように変数変換して計算するとよい.

  
  x = ((b - a)/2)cos(θ) + (a + b)/2

 

特異性が強い関数の積分

クレンショウ・カーチス公式で f(x)w(x) の形の関数を積分することを考える. w(x) は特異性が強いような関数を想定して, f(x) だけをチェビシェフ級数展開し w(x) を含む項別積分は別途行うことにする.

  
  I = ∫f(x)w(x)dx [-1, 1]
  ≒ In = (1/2)a0(n)∫T0(x)w(x)dx[-1, 1] + Σak(n)∫Tk(x)W(x)dx[-1, 1] + (1/2)an(n)∫Tn(x)w(x)dx[-1, 1] (Σ は k = 1~n-1)
  = (1/2)a0(n)μ0 + Σak(n)μk + (1/2)an(n)μn (Σ は k = 1~n-1)

  ただし μk = ∫Tk(x)w(x)dx [-1, 1]

 
μk はモーメントとよばれる. モーメントが安定的に計算できる(例えば, 解析的に)ような w(x) を含む積分については特異性が強くてもクレンショウ・カーチス公式で計算できることになる.

次の積分を考える.

  
  I = ∫f(x)/(x - c)dx [-1, 1]

 
これはコーシーの主値積分である.

w(x) = 1/(x – c) についてモーメントは次のように計算できる(文献[3]参照).

  
  μk = ∫Tk(x)/(x - c)dx [-1, 1]

  μ0 = ln|(1 - c)/(1 + c)|
  μ1 = 2 + cμ0
  μk+1 - 2cμk + μk-1 =
    0 (k が奇数)
    4/(1 - k2) (k が偶数)

 
c の値が [-1-ε, 1+ε] (ε はおおむね 0.1) の中にあれば上式は安定に計算できる. そうでなければ, c は積分範囲から十分に遠くにあるとみなして通常のガウス公式などで計算できる.

文献[3]では他に以下の関数についての計算法も示されており, これらの計算法はQUADPACKのサブルーチンで採用されている.

・ sin(ωx) および cos(ωx)
・ (x – a)α(b – x)βlogμ(x – a)logν(b – x) (μ, ν = 0 または 1)

(3) ロンバーグ積分

数列 ak が定数aに収束するものとして, 収束率λを

  
  λ = lim (ak+1 - a)/(ak - a)  [k→∞]

 
と定義する. λの値が事前にわかっているとして, これをaについて解いて求めた値

  
  bk = (ak+1 - λak)/(1 - λ)

 
による数列 {bk} は元の数列 {ak} よりも速く収束すると予想される. これをリチャードソンの加速法という.

複合台形則において分点数を倍々にした系列 I1, I2, ・・・ を作れば, h は半々になってゆき, 誤差は1段進むごとに1/4になる(すなわち収束率が1/4: 台形則の誤差の項を参照せよ). これにリチャードソンの加速法を適用した補外を行い,

  
  Ik(1) = (Ik+1 - λIk)/(1 - λ) = (4Ik+1 - Ik)/(4 - 1)

 
とした系列 { I1(1), I2(1), ・・・ } を得ることができる. これはシンプソン則に一致していることがわかる. 収束率λ(1)は1/42になる. 更に収束を加速するために新しい系列に第2回補外を行うと

  
  Ik(2) = (Ik+1(1) - λ(1)Ik(1))/(1 - λ(1)) = (42Ik+1(1) - Ik(1))/(42 - 1)

 
とした系列 { I1(2), I2(2), ・・・ } を得る. このように補外を繰り返して収束を加速する方法をロンバーグ積分法という.

複合台形則において分点数を倍々にしていく際には, 毎回半分は倍にする前の分点と一致するので半分だけ計算すればよいことになる.

数値実験

本文で説明した補間型積分公式についていくつかの積分関数を使って計算精度を比較する. なお, 計算はVBAのDouble型を使って倍精度(15~16桁)で行った.

有限区間の積分

(1) 性質の良い関数

次の積分を考える.

  
  I = ∫ ex cos x dx [0, 1] = 1.37802 46135 47364

 
これは性質の良い(解析的な)関数で, グラフに表すと次のような形をしている.

ex1f

これを種々の積分公式で分点数N(横軸)を変えて計算し, その時の相対誤差(縦軸)を対数グラフにすると次のようになった.

台形則は小区間数を N – 1 とした複合公式の結果を示した. シンプソン則は小区間数を (N – 1)/2 とした複合公式の結果を示した. ニュートン・コーツ公式およびクレンショウ・カーチス公式では N = n + 1 としてプロットした.

ニュートン・コーツ公式およびチェビシェフ公式では分点数が偶数よりも奇数のほうが成績がよいので, これらの公式を使う際には分点数が奇数の公式を使った方がよさそうである.

ガウス公式は他に比べて精度がよく, 可能であればこの公式を使うのがよさそうである. N = 7 以上で頭打ちなのは倍精度計算の丸め誤差のためである.

分点を倍々に増やしていったときの台形則の計算結果に着目すると次のようになっている.

分点数 積分値(計算値) 積分値(真値) 誤差 収束率
2 1.34061800327106 1.37802461354736 0.037406610
4 1.36858238253106 1.37802461354736 0.009442231 1/3.961
8 1.37565843490021 1.37802461354736 0.002366179 1/3.990
16 1.37743271822098 1.37802461354736 0.000591895 1/3.997
32 1.37787661780930 1.37802461354736 0.000147996 1/3.999
収束率はほぼ1/4になっており, ロンバーグ積分の項で示した収束率が成り立っていることがわかる.

ロンバーグ積分は簡単な仕組みで, 台形則やシンプソン則と比べると手間はあまり変わらないが精度はよいので, こちらを使ったほうがよいかもしれない.

(2) 積分区間外に特異点を持つ関数

積分区間外に特異点を持つ次の関数の積分を考える.

  
  I = ∫ 1/(1+x)dx [0, 4] = ln(1+4) = ln(5)

 
次のように積分区間 [0, 4] では滑らかな関数に見えるが, 近くの x = -1 に特異点を持つ.

ex2f

これを上と同じように種々の積分公式で計算し誤差を対数グラフにする.

図のようにどの公式でも精度が落ちたのがわかる. ガウス公式でも N = 18 以上でないと最高精度は得られない.

次に以下の積分を考える.

  
  I = ∫ 1/(1+x2)dx [0, 4] = arctan(4)

 
これは実数軸上には特異点がないが, 複素平面では近く(x = i) に特異点がある.

ex3f

計算結果は次のとおりである.

これも一つ前の例と同じようにどの公式でも精度が落ちたのがわかる. ロンバーグ積分が失速しているが, 特異性のために加速が成り立たなくなっているものと思われる.

(3) 端点に特異点を持つ関数

次に端点に特異点があるとどうなるか見てみる.

  
  I = ∫ √(1-x2)dx [0, 1] = π/4

 
x = 1 が特異点である.

ex4f

計算結果は次のとおりである.

このように端点とはいえ区間内に特異点がある場合にはどの公式でもまともな積分値は得られなかった.

この例における台形則の収束率は次のようになった.

分点数 積分値(計算値) 積分値(真値) 誤差 収束率
2 0.683012701892219 0.785398163397448 0.10238546
4 0.748927267025610 0.785398163397448 0.036470896 1/2.807
8 0.772454786089293 0.785398163397448 0.012943377 1/2.818
16 0.780813259456935 0.785398163397448 0.0045849039 1/2.823
32 0.783775605719283 0.785398163397448 0.0016225577 1/2.826
収束率は1/4ではなく1/2.8にしかなっておらず, 正常な動きをしていないことを示唆している. これに伴ってロンバーグ積分の加速の効果が小さくなっていることがわかる.

無限区間, 半無限区間の積分

(1) 半無限区間の積分

次の積分を考える.

  
  I = ∫ e-x/x [a, +∞] = E1(a)

 
E1(a) はよく使われる特殊関数である指数積分である. グラフに表すと次のような形をしている.

ガウス・ラゲール公式による計算結果は次のとおりである.

(2) 無限区間の積分

次の積分を考える.

  
  I = ∫ e-x2/(1 + x2) [-∞, +∞] = 1.3432934216

 
関数のグラフは次のとおりである.

ガウス・エルミート公式による計算結果は次のとおりである.

積分公式の係数表

ニュートン・コーツ公式

  
  ∫ f(x)dx [a, b] = h Σwkf(xk) (k = 0~n)
  ただし, h = (b - a)/n

 
ニュートン・コーツ公式における係数 wi を以下に示す. 本文で説明した式を数式処理ソフトを使って有理数のまま計算した.

  
         1  1
nc(1) = [-, -]
         2  2

         1  4  1
nc(2) = [-, -, -]
         3  3  3

         3  9  9  3
nc(3) = [-, -, -, -]
         8  8  8  8

         14  64  8   64  14
nc(4) = [--, --, --, --, --]
         45  45  15  45  45

         95   125  125  125  125  95
nc(5) = [---, ---, ---, ---, ---, ---]
         288  96   144  144  96   288

         41   54  27   68  27   54  41
nc(6) = [---, --, ---, --, ---, --, ---]
         140  35  140  35  140  35  140

         5257   25039  343  20923  20923  343  25039  5257
nc(7) = [-----, -----, ---, -----, -----, ---, -----, -----]
         17280  17280  640  17280  17280  640  17280  17280

         3956   23552    3712   41984    3632  41984    3712   23552  3956
nc(8) = [-----, -----, - -----, -----, - ----, -----, - -----, -----, -----]
         14175  14175    14175  14175    2835  14175    14175  14175  14175

         25713  141669  243   10881  26001  26001  10881  243   141669  25713
nc(9) = [-----, ------, ----, -----, -----, -----, -----, ----, ------, -----]
         89600  89600   2240  5600   44800  44800  5600   2240  89600   89600

          80335   132875    80875  28375    24125  89035    24125  28375    80875  132875  80335
nc(10) = [------, ------, - -----, -----, - -----, -----, - -----, -----, - -----, ------, ------]
          299376  74844     99792  6237     5544   12474    5544   6237     99792  74844   299376

 

チェビシェフ公式

  
  ∫ f(x)dx [a, b] = (b - a)/N Σf((a + b)/2 + (b - a)xi/2) (i = 1~N)

 
チェビシェフ公式における分点 xi を以下に示す. 拡張倍精度(80ビット浮動小数)の計算により本文で説明した多項式の方程式の解を求め, 結果の16桁を表示した.

  
N = 2
  xi: -5.773502691896258E-0001  5.773502691896258E-0001
N = 3
  xi: -7.071067811865475E-0001  0.000000000000000E+0000  7.071067811865475E-0001
N = 4
  xi: -7.946544722917661E-0001 -1.875924740850799E-0001  1.875924740850799E-0001
       7.946544722917661E-0001
N = 5
  xi: -8.324974870009819E-0001 -3.745414095535811E-0001  0.000000000000000E+0000
       3.745414095535811E-0001  8.324974870009819E-0001
N = 6
  xi: -8.662468181078206E-0001 -4.225186537611115E-0001 -2.666354015167047E-0001
       2.666354015167047E-0001  4.225186537611115E-0001  8.662468181078206E-0001
N = 7
  xi: -8.838617007580490E-0001 -5.296567752851568E-0001 -3.239118105199076E-0001
       0.000000000000000E+0000  3.239118105199076E-0001  5.296567752851568E-0001
       8.838617007580490E-0001
N = 9
  xi: -9.115893077284345E-0001 -6.010186553802381E-0001 -5.287617830578800E-0001
      -1.679061842148039E-0001  0.000000000000000E+0000  1.679061842148039E-0001
       5.287617830578800E-0001  6.010186553802381E-0001  9.115893077284345E-0001

 

ガウス・ルジャンドル公式 (ガウス公式)

  
  ∫ f(x)dx [a, b] = ((b - a)/2)Σwkf((a + b)/2 + ((b - a)/2)xk) (k = 1~N)

 
ガウス・ルジャンドル公式における分点 xi および重み wi を以下に示す. 分点としてルジャンドル多項式のゼロ点を求め, それから重みを計算した. 計算は拡張倍精度(80ビット浮動小数)により行い結果の16桁を表示した.

  
N = 2
  xi: -5.773502691896258E-0001  5.773502691896258E-0001
  wi:  1.000000000000000E+0000  1.000000000000000E+0000
N = 3
  xi: -7.745966692414834E-0001  0.000000000000000E+0000  7.745966692414834E-0001
  wi:  5.555555555555556E-0001  8.888888888888889E-0001  5.555555555555556E-0001
N = 4
  xi: -8.611363115940526E-0001 -3.399810435848563E-0001  3.399810435848563E-0001
       8.611363115940526E-0001
  wi:  3.478548451374539E-0001  6.521451548625461E-0001  6.521451548625461E-0001
       3.478548451374539E-0001
N = 5
  xi: -9.061798459386640E-0001 -5.384693101056831E-0001  0.000000000000000E+0000
       5.384693101056831E-0001  9.061798459386640E-0001
  wi:  2.369268850561891E-0001  4.786286704993665E-0001  5.688888888888889E-0001
       4.786286704993665E-0001  2.369268850561891E-0001
N = 6
  xi: -9.324695142031520E-0001 -6.612093864662645E-0001 -2.386191860831969E-0001
       2.386191860831969E-0001  6.612093864662645E-0001  9.324695142031520E-0001
  wi:  1.713244923791703E-0001  3.607615730481386E-0001  4.679139345726910E-0001
       4.679139345726910E-0001  3.607615730481386E-0001  1.713244923791703E-0001
N = 7
  xi: -9.491079123427585E-0001 -7.415311855993944E-0001 -4.058451513773972E-0001
       0.000000000000000E+0000  4.058451513773972E-0001  7.415311855993944E-0001
       9.491079123427585E-0001
  wi:  1.294849661688697E-0001  2.797053914892767E-0001  3.818300505051189E-0001
       4.179591836734694E-0001  3.818300505051189E-0001  2.797053914892767E-0001
       1.294849661688697E-0001
N = 8
  xi: -9.602898564975362E-0001 -7.966664774136267E-0001 -5.255324099163290E-0001
      -1.834346424956498E-0001  1.834346424956498E-0001  5.255324099163290E-0001
       7.966664774136267E-0001  9.602898564975362E-0001
  wi:  1.012285362903763E-0001  2.223810344533745E-0001  3.137066458778873E-0001
       3.626837833783620E-0001  3.626837833783620E-0001  3.137066458778873E-0001
       2.223810344533745E-0001  1.012285362903763E-0001
N = 9
  xi: -9.681602395076261E-0001 -8.360311073266358E-0001 -6.133714327005904E-0001
      -3.242534234038089E-0001  0.000000000000000E+0000  3.242534234038089E-0001
       6.133714327005904E-0001  8.360311073266358E-0001  9.681602395076261E-0001
  wi:  8.127438836157441E-0002  1.806481606948574E-0001  2.606106964029355E-0001
       3.123470770400028E-0001  3.302393550012598E-0001  3.123470770400028E-0001
       2.606106964029355E-0001  1.806481606948574E-0001  8.127438836157441E-0002
N = 10
  xi: -9.739065285171717E-0001 -8.650633666889845E-0001 -6.794095682990244E-0001
      -4.333953941292472E-0001 -1.488743389816312E-0001  1.488743389816312E-0001
       4.333953941292472E-0001  6.794095682990244E-0001  8.650633666889845E-0001
       9.739065285171717E-0001
  wi:  6.667134430868814E-0002  1.494513491505806E-0001  2.190863625159820E-0001
       2.692667193099964E-0001  2.955242247147529E-0001  2.955242247147529E-0001
       2.692667193099964E-0001  2.190863625159820E-0001  1.494513491505806E-0001
       6.667134430868814E-0002

 

ガウス・ラゲール公式

  
  ∫ f(x)exp(-x)dx [0, +∞] = Σwkf(xk) (k = 1~N)

 
ガウス・ラゲール公式における分点 xi および重み wi を以下に示す. 分点としてラゲール多項式のゼロ点を求め, それから重みを計算した. 計算は拡張倍精度(80ビット浮動小数)により行い結果の16桁を表示した.

  
N = 2
  xi: 5.857864376269050E-0001 3.414213562373095E+0000
  wi: 8.535533905932738E-0001 1.464466094067262E-0001
N = 3
  xi: 4.157745567834791E-0001 2.294280360279042E+0000 6.289945082937479E+0000
  wi: 7.110930099291730E-0001 2.785177335692408E-0001 1.038925650158614E-0002
N = 4
  xi: 3.225476896193923E-0001 1.745761101158347E+0000 4.536620296921128E+0000
      9.395070912301133E+0000
  wi: 6.031541043416336E-0001 3.574186924377997E-0001 3.888790851500538E-0002
      5.392947055613275E-0004
N = 5
  xi: 2.635603197181409E-0001 1.413403059106517E+0000 3.596425771040722E+0000
      7.085810005858838E+0000 1.264080084427578E+0001
  wi: 5.217556105828087E-0001 3.986668110831759E-0001 7.594244968170760E-0002
      3.611758679922048E-0003 2.336997238577623E-0005
N = 6
  xi: 2.228466041792607E-0001 1.188932101672623E+0000 2.992736326059314E+0000
      5.775143569104511E+0000 9.837467418382590E+0000 1.598287398060170E+0001
  wi: 4.589646739499636E-0001 4.170008307721210E-0001 1.133733820740450E-0001
      1.039919745314907E-0002 2.610172028149321E-0004 8.985479064296212E-0007
N = 7
  xi: 1.930436765603624E-0001 1.026664895339192E+0000 2.567876744950746E+0000
      4.900353084526485E+0000 8.182153444562861E+0000 1.273418029179781E+0001
      1.939572786226254E+0001
  wi: 4.093189517012739E-0001 4.218312778617198E-0001 1.471263486575053E-0001
      2.063351446871694E-0002 1.074010143280746E-0003 1.586546434856420E-0005
      3.170315478995581E-0008
N = 8
  xi: 1.702796323051010E-0001 9.037017767993799E-0001 2.251086629866131E+0000
      4.266700170287659E+0000 7.045905402393466E+0000 1.075851601018100E+0001
      1.574067864127800E+0001 2.286313173688926E+0001
  wi: 3.691885893416375E-0001 4.187867808143430E-0001 1.757949866371718E-0001
      3.334349226121565E-0002 2.794536235225673E-0003 9.076508773358213E-0005
      8.485746716272532E-0007 1.048001174871510E-0009
N = 9
  xi: 1.523222277318082E-0001 8.072200227422558E-0001 2.005135155619347E+0000
      3.783473973331233E+0000 6.204956777876613E+0000 9.372985251687576E+0000
      1.346623691109209E+0001 1.883359778899170E+0001 2.637407189092738E+0001
  wi: 3.361264217979625E-0001 4.112139804239844E-0001 1.992875253708856E-0001
      4.746056276565160E-0002 5.599626610794583E-0003 3.052497670932106E-0004
      6.592123026075352E-0006 4.110769330349548E-0008 3.290874030350708E-0011
N = 10
  xi: 1.377934705404924E-0001 7.294545495031705E-0001 1.808342901740316E+0000
      3.401433697854900E+0000 5.552496140063804E+0000 8.330152746764497E+0000
      1.184378583790007E+0001 1.627925783137810E+0001 2.199658581198076E+0001
      2.992069701227389E+0001
  wi: 3.084411157650201E-0001 4.011199291552736E-0001 2.180682876118094E-0001
      6.208745609867775E-0002 9.501516975181101E-0003 7.530083885875388E-0004
      2.825923349599566E-0005 4.249313984962686E-0007 1.839564823979631E-0009
      9.911827219609009E-0013

 

ガウス・エルミート公式

  
  ∫ f(x)exp(-x2)dx [-∞, +∞] Σwkf(xk) (k = 1~N)

 
ガウス・エルミート公式における分点 xi および重み wi を以下に示す. 分点としてエルミート多項式のゼロ点を求め, それから重みを計算した. 計算は拡張倍精度(80ビット浮動小数)により行い結果の16桁を表示した.

  
N = 2
  xi: -7.071067811865475E-0001  7.071067811865475E-0001
  wi:  8.862269254527580E-0001  8.862269254527580E-0001
N = 3
  xi: -1.224744871391589E+0000  0.000000000000000E+0000  1.224744871391589E+0000
  wi:  2.954089751509193E-0001  1.181635900603677E+0000  2.954089751509193E-0001
N = 4
  xi: -1.650680123885785E+0000 -5.246476232752903E-0001  5.246476232752903E-0001
       1.650680123885785E+0000
  wi:  8.131283544724518E-0002  8.049140900055128E-0001  8.049140900055128E-0001
       8.131283544724518E-0002
N = 5
  xi: -2.020182870456086E+0000 -9.585724646138185E-0001  0.000000000000000E+0000
       9.585724646138185E-0001  2.020182870456086E+0000
  wi:  1.995324205904591E-0002  3.936193231522412E-0001  9.453087204829419E-0001
       3.936193231522412E-0001  1.995324205904591E-0002
N = 6
  xi: -2.350604973674492E+0000 -1.335849074013697E+0000 -4.360774119276165E-0001
       4.360774119276165E-0001  1.335849074013697E+0000  2.350604973674492E+0000
  wi:  4.530009905508846E-0003  1.570673203228566E-0001  7.246295952243925E-0001
       7.246295952243925E-0001  1.570673203228566E-0001  4.530009905508846E-0003
N = 7
  xi: -2.651961356835233E+0000 -1.673551628767471E+0000 -8.162878828589647E-0001
       0.000000000000000E+0000  8.162878828589647E-0001  1.673551628767471E+0000
       2.651961356835233E+0000
  wi:  9.717812450995192E-0004  5.451558281912703E-0002  4.256072526101278E-0001
       8.102646175568073E-0001  4.256072526101278E-0001  5.451558281912703E-0002
       9.717812450995192E-0004
N = 8
  xi: -2.930637420257244E+0000 -1.981656756695843E+0000 -1.157193712446780E+0000
      -3.811869902073221E-0001  3.811869902073221E-0001  1.157193712446780E+0000
       1.981656756695843E+0000  2.930637420257244E+0000
  wi:  1.996040722113676E-0004  1.707798300741348E-0002  2.078023258148919E-0001
       6.611470125582413E-0001  6.611470125582413E-0001  2.078023258148919E-0001
       1.707798300741348E-0002  1.996040722113676E-0004
N = 9
  xi: -3.190993201781528E+0000 -2.266580584531843E+0000 -1.468553289216668E+0000
      -7.235510187528376E-0001  0.000000000000000E+0000  7.235510187528376E-0001
       1.468553289216668E+0000  2.266580584531843E+0000  3.190993201781528E+0000
  wi:  3.960697726326438E-0005  4.943624275536947E-0003  8.847452739437657E-0002
       4.326515590025558E-0001  7.202352156060510E-0001  4.326515590025558E-0001
       8.847452739437657E-0002  4.943624275536947E-0003  3.960697726326438E-0005
N = 10
  xi: -3.436159118837738E+0000 -2.532731674232790E+0000 -1.756683649299882E+0000
      -1.036610829789514E+0000 -3.429013272237046E-0001  3.429013272237046E-0001
       1.036610829789514E+0000  1.756683649299882E+0000  2.532731674232790E+0000
       3.436159118837738E+0000
  wi:  7.640432855232621E-0006  1.343645746781233E-0003  3.387439445548106E-0002
       2.401386110823147E-0001  6.108626337353258E-0001  6.108626337353258E-0001
       2.401386110823147E-0001  3.387439445548106E-0002  1.343645746781233E-0003
       7.640432855232621E-0006

 


参考文献

[1] 森正武「数値解析 (第2版)」(2002) 共立出版
[2] 森口繁一「数値計算工学」(1989) 岩波書店
[3] Piessens, de Doncker-Kapenga, Überhuber, Kahaner: “QUADPACK: A subroutine package for automatic integration”, Springer-Verlag. (1983)