数値演算法

(9) 数値積分法 -1-

積分法の歴史は古く、古代ギリシャの有名な数学者「アルキメデス ( Archimedes )」は、「取り尽くし法 ( Method of Exhaustion )」と呼ばれる方法を使って円や放物線などの面積を計算したことで知られています。その後、微分積分学は、「アイザック・ニュートン ( Isaac Newton )」と「ゴットフリート・ライプニッツ ( Gottfried Wilhelm Leibniz )」によって 17 世紀のヨーロッパで確立されました。今では、高校で多項式や三角関数・指数関数の原始関数を習い、それを利用すれば放物線の面積なども正確に求めることができます。
しかし、原始関数を求めることができる関数は限られていて、常に式の形が得られるわけではないので、定積分の計算には数値積分法がよく利用されます。今回は、数値積分法の手法について紹介したいと思います。

1) 区分求積法 ( Quadrature by Parts )

高校での微分・積分の授業では、最初に微分、次に積分を教えるのが通常のようで、自分も先に微分法を習いました。積分法では、「微分すると f(x) となる関数を f(x) の原始関数 ( Primitive Functions )という」という内容から始まって不定積分を覚え、次に定積分によって面積や体積を求める手法の説明がありました。今でもこの流れに変化はないと思いますが、このような教え方はあまり好評ではないそうで、実際には積分法の方が歴史的に古く、その考え方も最初は素朴なものでした。

図 1-1. 区分求積法

上図のような曲線 y = f(x) に対して、区間 [ a, b ] における x 軸に挟まれた部分の面積 S を求めたい時、正確にその値を求めるには、

のように定積分から得ることができます。しかし、そのためには f(x) の原始関数 F(x) を求める必要があり、それは常に得られるとは限りません。古代の数学者は、原始関数を求めるのではなく、面積を求める範囲を矩形によって細かく分割し、その面積の和を求めることによって近似値として得ることを考えました。左側から k 番目の矩形の高さを、その区間内から x 座標の代表値 x_k を適当に選んで f(x_k) としたとき、S の近似値 S' は

と表すことができます。分割する数を多くすればそれだけ得られる値の精度は向上し、n → ∞ ( Δx_k → 0 ) とすれば、S' → S となるので、

これが定積分の定義になります ( 補足 1 )。

矩形の高さは「適当に選ぶ」としましたが、実際にはランダムに選ぶことはなく、何らかの規則に従って決まった位置を選択することになります。例えば、S を n 等分する場合、S は Y 軸に平行な n - 1 本の直線によって分割されるので、それらの x 座標を { x_k } ( k = 1, 2, ... n - 1 ) として a = x₀, b = x_n とすれば、f(x₀) 〜 f(x_n-1) を選んだとき各矩形の左側、f(x₁) 〜 f(x_n) を選べば各矩形の右側を選択したことになります。また、矩形の幅 ( 分割時の区間の長さ ) も任意ではなく、等間隔に分割するのが一般的なやり方です。その場合、各矩形の幅は Δx = ( b - a ) / n と表すことができます。

関数 f(x) が単調に増加または減少しているなら、矩形の端で高さを決めるよりも中央で高さを決めた方が精度がよくなります。分割位置の間のちょうど中央での値 f( [ (x_k) + f(x_k-1) ] / 2 ) を各矩形の高さとして面積の総和を求める方法を「中点則 ( または中点公式 ; Midpoint Rule )」といいます。

また、両端の値の平均値 [ f(x_k) + f(x_k-1) ] / 2 ( k = 1, 2, ... n ) を矩形の高さとした場合、これに Δx を掛けたものは、上底と下底が f(x_k), f(x_k-1)、高さが Δx の台形の面積を表すことになります。この台形の面積の総和から近似値を求める方法は「台形則 ( または台形公式 ; Trapezoidal Rule )」といいます。台形則を展開すると、

のように、非常にシンプルな式になります。

中点則と台形則を使って定積分の値を求めるサンプル・プログラムを以下に示します。

/**
   Integ_Midpoint : 中点則 ( Midpoint Rule ) による定積分計算用クラス

   関数 Integral へ引数として渡して利用する。

   例) Integral( Integ_Midpoint(), f, 0.0, 1.0, 10000 );
**/
struct Integ_Midpoint
{
  // 積分値の初期値設定
  //
  // f : 積分する関数
  // a, b : 積分する範囲
  // h : 矩形の幅
  template< class T, class F >
  T initValue( F f, T a, T b, T h )
  { return( f( a + h / 2 ) ); }

  // x の初期値設定
  //
  // a : 積分範囲の開始位置
  // h : 矩形の幅
  template< class T >
  T initX( T a, T h )
  { return( a + 3 * h / 2 ); }

  // 加算していく積分値の計算
  //
  // f : 積分する関数
  // x : f(x) を計算する x の値
  // h : 矩形の幅
  template< class T, class F >
  T operator()( F f, T x, T h )
  { return( f( x ) ); }

  // 和を計算した後の補正
  //
  // val : 補正前の積分値
  // h : 矩形の幅
  template< class T >
  T cleanUp( T val, T h )
  { return( val * h ); }
};

/**
   Integ_Trapezoidal : 台形則 ( Trapezoidal Rule ) による定積分計算用クラス

   関数 Integral へ引数として渡して利用する。

   例) Integral( Integ_Trapezoidal(), f, 0.0, 1.0, 10000 );
**/
struct Integ_Trapezoidal
{
  // 積分値の初期値設定
  //
  // f : 積分する関数
  // a, b : 積分する範囲
  // h : 矩形の幅
  template< class T, class F >
  T initValue( F f, T a, T b, T h )
  { return( ( f( a ) + f( b ) ) / 2 ); }

  // x の初期値設定
  //
  // a : 積分範囲の開始位置
  // h : 矩形の幅
  template< class T >
  T initX( T a, T h )
  { return( a + h ); }

  // 加算していく積分値の計算
  //
  // f : 積分する関数
  // x : f(x) を計算する x の値
  // h : 矩形の幅
  template< class T, class F >
  T operator()( F f, T x, T h )
  { return( f( x ) ); }

  // 和を計算した後の補正
  //
  // val : 補正前の積分値
  // h : 矩形の幅
  template< class T >
  T cleanUp( T val, T h )
  { return( val * h ); }
};

/*
  Integral : 区分求積法による積分値の計算

  計算手法 ( 中点則・台形則など ) は引数 op にて指定する。
  積分範囲の大小関係が逆転している ( a > b の ) 場合、積分値の符号が逆転することに注意。

  op : 積分値を計算する手法
  f : 積分値を計算する対象の関数
  a,b : 積分範囲
  n : 分割数 ( 大きいほど精度が上がるが、処理時間は長くなる )
*/
template< class T, class F, class Op >
  T Integral( Op op, F f, T a, T b, unsigned int n )
{
  assert( n > 0 );
  bool invSign; // 符号反転するか

  // a > b なら swap して後で符号を反転する
  if ( a > b ) {
    std::swap( a, b );
    invSign = true;
  } else {
    invSign = false;
  }

  T h( ( b - a ) / n ); // 矩形の幅

  // 計算開始

  T ans = op.initValue( f, a, b, h ); // 積分値の初期状態
  T x = op.initX( a, h );       // x の初期値
  for ( unsigned int i = 1 ; i < n ; ++i ) {
    ans += op( f, x, h ); // f(x) を求めて加算する
    x += h;
  }

  // 符号反転
  if ( invSign ) ans *= -1;

  // 正しい積分値にする
  return( op.cleanUp( ans, h ) );
}

計算処理は非常に単純で、関数 f(x) の値を求めながら変数 ans に総和を求めているだけです。その処理を行うのが関数 Integral です。前処理として、分割数がゼロだった場合は処理を中断するようにし ( assert を使っているためプログラム自体が停止します。その他に例外を投げたり、NaN などの無効値を返すやり方もあります )、また a > b である場合は変数を入れ替えると同時に、後で総和の符号を反転する処理を行えるように変数 invSign のフラグを立てています。これは、

という定義によるものです。ところで、関数 f(x) の値は負値にもなり得るので、矩形の面積も負になる可能性があります。これも定積分の定義により、x 軸よりも下側の面積は負値として扱うため、符号は変えずにそのまま総和を計算します。

処理の方式は引数 op で切り替えられるようにしています。その引数として渡すのが中点則の Integ_Midpoint クラスと台形則の Integ_Trapezoidal クラスです。これらのクラスは 4 つのメンバ関数 initValue, initX, operator(), cleanUp を共通で持っています。initValue は積分値の初期値を計算するためのもので、最初の矩形の面積はこの関数で計算します。次に x 座標を初期化しますが、これは二番目の矩形に対する x 座標を計算することになるのに注意してください。ループ処理では、二番目の矩形から面積を算出します。このとき使うのが operator() です。最後の cleanUp は求めた和に定数等を作用させるために利用します。和の計算時に共通となる部分は、最後の cleanUp 時に乗除させた方が精度が向上します。

いくつかの関数を使い、実際に計算を行った結果を以下に示します。計算した定積分の式は次の通りです。

上側の式の値は 1 になり、下側の ( 1 - x² )^1/2 は半径 1 の円の第一象限における方程式を表しているので、定積分の値は π になります。

結果を見ると、少なくとも 1000 等分、できれば 10000 等分に分割しなければ充分な精度は得られないことが分かります。関数が複雑になればそれだけ処理時間も長くなるため、総和を求める時の分割数を少なくしても充分な精度が得られるような方法が必要となります。

2) シンプソン則 ( Simpson's Rule )

台形則では、区間内の定積分値を台形の面積として近似計算したわけですが、これは、区間内の関数 f(x) を直線で近似したことと同じ意味になります。

とすれば、各区間での関数 y = f(x) を直線 y = a_k x + b_k に置き換えたことになり、以前紹介した「画像のサンプル補間」における「線形補間法」と同様な手法を利用しています。サンプル補間では直線近似をしてピクセル間の値を決めていたのに対して、ここでは面積を求めるために関数を直線で置き換えていることになります。

図 2-1. 台形則が面積として計算する部分

「画像のサンプル補間」の中で紹介した「多項式補間」では、ピクセル間の値を多項式を使って補間しています。定積分においても任意の関数を多項式で近似すれば、多項式の定積分は簡単にできるので、より正確な近似計算が少ない分割数で実現できそうです。

図 2-2. 多項式による関数の近似

まずは、二次式による近似を検討してみます。二次式 P(x) = ax² + bx + c の各未知数 a, b, c を求めるためには三点を決めればよいので、両端の座標 ( x_k-1, f(x_k-1) ), ( x_k, f(x_k) ) とその中点の座標 ( x_m, f(x_m) ) ( 但し、x_m = ( x_k-1 + x_k ) / 2 ) を使って連立方程式を解きます。

として、解きやすくするために P(x) = a( x - x_M )² + b( x - x_M ) + c と表すと、連立方程式は

となります。ところが、

なので、x_R - x_L = h とおくと、

と、非常に簡単になります。この連立方程式を解くと、

P(x) の積分値を求めると、

よって、

という近似式ができます。これを「シンプソン則 ( Simpson's Rule )」といいます。積分する区間を偶数個に分割し、( x_2k, f(x_2k) ) から ( x_2k+2, f(x_2k+2) ) までの区間にシンプソン則を適用すると、各区間の幅を h としたとき

なので、x₀ から x_2n までの和を取れば

となり、分割数を多くすればそれだけ精度を上げることもできます。これを「合成シンプソン則 ( Composite Simpson's Rule )」といいます。

図 2-3. シンプソン則による区間分割

シンプソン則を使って定積分の値を求めるサンプル・プログラムを以下に示します。

/**
   Integ_Simpson : シンプソン則 ( Composite Simpson's Rule ) による定積分計算用クラス

   関数 Integral へ引数として渡して利用する。

   例) Integral( Integ_Simpson(), f, 0.0, 1.0, 10000 );
**/
struct Integ_Simpson
{
  // 積分値の初期値設定
  //
  // f : 積分する関数
  // a, b : 積分する範囲
  // h : 矩形の幅
  template< class T, class F >
  T initValue( F f, T a, T b, T h )
  { return( f( a ) + 4 * f( b - h / 2 ) + f( b ) ); }

  // x の初期値設定
  //
  // a : 積分範囲の開始位置
  // h : 矩形の幅
  template< class T >
  T initX( T a, T h )
  { return( a ); }

  // 加算していく積分値の計算
  //
  // f : 積分する関数
  // x : f(x) を計算する x の値
  // h : 矩形の幅
  template< class T, class F >
  T operator()( F f, T x, T h )
  { return( 4 * f( x + h / 2 ) + 2 * f( x + h ) ); }

  // 和を計算した後の補正
  //
  // val : 補正前の積分値
  // h : 矩形の幅
  template< class T >
  T cleanUp( T val, T h )
  { return( val * h / 6 ); }
};

中点則、台形則を含めて、以下の定積分値を求めた結果は次のようになります。

上側は中点則・台形則のテストで利用したものと同じ式です。下側の定積分値は 1 になります。

cos θ の積分は 100 等分程度でかなりの精度まで一致していることが分かります。π の値を求める上側の式においても、中点則や台形則に比べると精度が向上しているように見えますが、cos θ の場合と比較すると劇的な変化とはいえないようです。この件については次の章でもう少し詳しく調べてみたいと思います。

3) ニュートン・コーツの公式 ( Newton-Cotes Formulas )

台形則は一次式、シンプソン則は二次式で関数を近似することによって定積分値を求めていました。よって、より高次の多項式で近似をすれば、精度を向上させることが期待できそうです。しかし、高次になるほど未知数となる係数の数は多くなるため、連立方程式の計算に手間がかかるようになります。そこで、「多項式補間」で利用した「ラグランジュの多項式 ( Lagrange Polynomial )」を使います。

関数 f(x) 上の任意の点を ( x_i, f( x_i ) ) としたとき、

但し、

とすれば、

となるので f( x_i ) = L_N( x_i ) を満たします。そこで、f(x) の定積分値を L_N(x) で近似して、

∫{a→b} l_i( x ) dx = w_i、f( x_i ) = f_i として、

これを「ニュートン・コーツの公式 ( Newton-Cotes Formulas )」といいます。N = 1 であれば台形則、N = 2 であればシンプソン則を表し、ニュートン・コーツの公式はそれらを一般化したものと考えることができます。N = 3 の場合、三等分したときの各区間を h とすれば

ここで l₀( x ) に対して x - x₂ = X とおくと、

なので

になります。また、l₁( x ) に対して x - x₂ = X とおくと、

となって、

と求めることができます。w₂、w₃ も同様な方法によってそれぞれ 9h / 8、3h / 8 と求められるので、

となります。主要な公式について、その分割数と式をまとめたものを以下に示します。

シンプソンの 3/8 則とブール則によって定積分の値を求めるサンプル・プログラムを以下に示します。

/**
   Integ_Simpson38 : シンプソンの 3/8 則 ( Simpson's 3/8 Rule ) による定積分計算用クラス

   関数 Integral へ引数として渡して利用する。

   例) Integral( Integ_Simpson38(), f, 0.0, 1.0, 10000 );
**/
struct Integ_Simpson38
{
  // 積分値の初期値設定
  //
  // f : 積分する関数
  // a, b : 積分する範囲
  // h : 矩形の幅
  template< class T, class F >
  T initValue( F f, T a, T b, T h )
  { return( f( a ) + 3 * ( f( b - h / 3 ) + f( b - 2 * h / 3 ) ) + f( b ) ); }

  // x の初期値設定
  //
  // a : 積分範囲の開始位置
  // h : 矩形の幅
  template< class T >
  T initX( T a, T h )
  { return( a ); }

  // 加算していく積分値の計算
  //
  // f : 積分する関数
  // x : f(x) を計算する x の値
  // h : 矩形の幅
  template< class T, class F >
  T operator()( F f, T x, T h )
  { return( 3 * ( f( x + h / 3 ) + f( x + 2 * h / 3 ) ) + 2 * f( x + h ) ); }

  // 和を計算した後の補正
  //
  // val : 補正前の積分値
  // h : 矩形の幅
  template< class T >
  T cleanUp( T val, T h )
  { return( val * h / 8 ); }
};

/**
   Integ_Boole : ブール則 ( Boole's Rule ) による定積分計算用クラス

   関数 Integral へ引数として渡して利用する。

   例) Integral( Integ_Boole(), f, 0.0, 1.0, 10000 );
**/
struct Integ_Boole
{
  // 積分値の初期値設定
  //
  // f : 積分する関数
  // a, b : 積分する範囲
  // h : 矩形の幅
  template< class T, class F >
  T initValue( F f, T a, T b, T h )
  { return( 7 * f( a ) + 32 * ( f( b - 3 * h / 4 ) + f( b - h / 4 ) ) + 12 * f( b - h / 2 ) + 7 * f( b ) ); }

  // x の初期値設定
  //
  // a : 積分範囲の開始位置
  // h : 矩形の幅
  template< class T >
  T initX( T a, T h )
  { return( a ); }

  // 加算していく積分値の計算
  //
  // f : 積分する関数
  // x : f(x) を計算する x の値
  // h : 矩形の幅
  template< class T, class F >
  T operator()( F f, T x, T h )
  { return( 32 * ( f( x + h / 4 ) + f( x + 3 * h / 4 ) ) + 12 * f( x + h / 2 ) + 14 * f( x + h ) ); }

  // 和を計算した後の補正
  //
  // val : 補正前の積分値
  // h : 矩形の幅
  template< class T >
  T cleanUp( T val, T h )
  { return( val * h / 90 ); }
};

シンプソン則、シンプソンの 3/8 則、ブール則を使った以下の定積分値の計算結果は次のようになります。

上式はどちらもシンプソン則のテストで利用したものと同じ式です。

∫{0→1} 4( 1 - x² )^1/2 dx と ∫{0→π/2} cos θ dθ を台形則から Boole 則まで、分割数を 1, 10, 100 として計算したときの誤差は次のようになります。

cos θ の定積分は、高次の多項式を利用するほど劇的に精度が向上する反面、π を求める上側の定積分は台形則を除くと大きな精度の向上は得られていないことが分かります。注意点として、高次の多項式を利用している方法では区間内をさらに分割して f(x) の値を求めているので、f(x) を取得する回数は倍々に増えています。例えば分割数が 10 の場合、台形則は f(x) の演算を初期化時に 2 回、ループ内で 9 回の計 11 回行っているのに対し、Simpson 則は初期化で 3 回、ループ内で 18回の計 21 回、3/8 則では初期化で 4 回、ループ内で 27 回の計 31 回、最後の Boole 則では初期化で 5 回、ループ内で 36 回の計 41 回 f(x) を計算しています。それを考慮すると、上側の式に対しては次数を増やすことによる効果はほとんど見られないことになります。

( 1 - x² )^1/2 の定積分において精度が向上しない理由は、x = 1 の付近で関数の値が急激に減少し、特に x = 1 では導関数の値が発散するためなので、x = 1 付近での定積分が行われないような工夫をすることで誤差を小さくすることができます。例えば、直線 y = x によって円弧は等分されるので、円弧と直線によって囲まれた部分の面積を定積分法で求めてそれを 8 倍することで、π の近似値を求めることができます。円弧と直線との交点は ( √2 / 2, √2 / 2 ) であり、面積は ( 1 - x² )^1/2 - x の定積分によって求められるので、計算式は

になります。

図 3-1. 円弧と直線で囲まれた部分の定積分

 

表 3-6. ∫{0→√2/2} 8[ ( 1 - x² )^1/2 - x ] dx の計算結果

	台形則		Simpson則		3/8則		Boole則
n	計算結果	誤差	計算結果	誤差	計算結果	誤差	計算結果	誤差
1	2.82842712475	-0.31316552884	3.13714412300	-0.00444853059	3.13948828738	-0.00210436621	3.14145643965	-0.00013621394
10	3.13826261523	-0.00333003836	3.14159183222	-0.00000082137	3.14159228795	-0.00000036564	3.14159265300	-0.00000000059
100	3.14155932059	-0.00003333300	3.14159265351	-0.00000000008	3.14159265355	-0.00000000004	3.14159265359	0.00000000000
1000	3.14159232026	-0.00000033333	3.14159265359	-0.00000000000	3.14159265359	-0.00000000000	3.14159265359	-0.00000000000
10000	3.14159265026	-0.00000000333	3.14159265359	0.00000000000	3.14159265359	0.00000000000	3.14159265359	-0.00000000000
100000	3.14159265355	-0.00000000004	3.14159265359	-0.00000000000	3.14159265359	-0.00000000000	3.14159265359	-0.00000000000

 

表 3-7. ∫{0→√2/2} 8[ ( 1 - x² )^1/2 - x ] dx の各公式における誤差

定積分法	分割数 「()内は台形則を1.00としたときの比率」
	1	10	100
台形則	-0.31316552884 (1.00)	-0.00333003836 (1.00)	-0.00003333300 (1.00)
Simpson則	-0.00444853059 (1.42E-2)	-0.00000082137 (2.47E-4)	-0.00000000008 (2E-6)
3/8則	-0.00210436621 (6.72E-3)	-0.00000036564 (1.10E-4)	-0.00000000004 (1E-6)
Boole則	-0.00013621394 (4.35E-4)	-0.00000000059 (1.8E-7)	0.00000000000 (-)

今度は高い精度で結果が得られ、高次の多項式を利用することによる効果も見られます。下に示したのは、√2/2 から 1 までの範囲の定積分値を求め、それに三角形部分の面積を加算した場合の算出結果です。こちらは精度が悪くなることが予想され、また結果も予想通りのものとなっています。

この結果から、数値積分法では関数の性質によって専用の処理方法を用意することがある程度必要で、汎用に使える手法というのは一般的にはないということになります。

これまでは、関数を高次多項式で近似して積分値を求める手法について紹介してきましたが、次数を大きくするほど公式は複雑になるので、できるだけ高次にすればよいというわけでもなく、また Simpson 則でも工夫次第で充分な精度が得られるので、次は精度を上げるための別の手法を紹介したいと思います。その手法を説明する前に、その中で利用されている「ベルヌーイ数」を先に紹介しておきます。

4) ベルヌーイ数 ( Bernoulli Number )

1 から n までの和は「三角数 ( Triangular Number )」と呼ばれ、点を三角形の形に並べたときの点の個数になっています。この数に関して、数学者のガウスが小学生時代、先生から出された課題をあっという間に解いてしまったという逸話があります。その時の問題は「1 から 100 までの数の和を求める」という、三角数を求める問題と同じものでした。n = 4 のときの三角数は、下図の中の黒点の数になります。これを上下反転したものを右側に追加すると(白点の部分)、横に五個、縦に四個の点が並んだ平行四辺形の形になるので、三角数は 5 x 4 / 2 = 10 と計算できます。

図 4-1. 三角数 ( n = 4 の場合 )

この考えかたを応用すると、1 から n までの和を計算するための公式が得られます。まず、1 から n までの数を横に並べ、その下に n から 1 まで数を逆に並べます。

1   2   3   4 ... n-1   n
n n-1 n-2 n-3 ...   2   1

上下の数どうしをそれぞれ足し合わせていきます。

    1   2   3   4 ... n-1   n
    n n-1 n-2 n-3 ...   2   1
+)______________________________
  n+1 n+1 n+1 n+1 ... n+1 n+1

それぞれの和は n + 1 になり、それが n 個できるので、その和は n( n + 1 ) になります。しかし、この中で 1 から n までの数が二回ずつ加算されているので、2 で割ることで求める三角数 n( n + 1 ) / 2 が得られます。

次に、n 個の数の和ではなく n 個の平方数の和を求めることを検討してみます。この和を S₂(n) と表すと、

になります。ここでは三角数を求める方法を利用することはできないので、別の手法を考える必要があります。そこで、次のような等式を利用します。

中辺の各項は互いに打ち消しあって左辺の形になります。このような和の形は「望遠鏡和 ( Telescoping Sum )」と呼ばれ、典型的な例としては

というものがあります。先ほどの変形は ( n + 1 )³ に対して無理矢理、望遠鏡和の形に表したと考えることができます。

( k + 1 )³ = k³ + 3k² + 3k + 1 なので、先ほどの式を変形すると

ここで、Σ_k{1→n}( k² ) は求める和 S₂(n) になります。また、Σ_k{1→n}( k ) は三角数の和を表し、Σ_k{1→n}( 1 ) = n なので

これを S₂(n) について解くと、

これで、S₂(n) を解くことができました。同様の手法を使って、一般解として

を求めてみたいと思います。まず、( n + 1 )^m+1 を望遠鏡和の形に表して

「二項定理」より、二項係数を C( n, r ) で表すと

を代入すると、

Σ_k{1→n}( k^r ) = S_r(n) なので、

総和の中から r = m の部分だけを抽出してその部分について解くと

C( m + 1, m ) = C( m + 1, 1 ) = m + 1 なので

これで、S_m(n) を計算するための漸化式が得られたことになります。三角数の和は S₁(n) = n( n + 1 ) / 2 でした。これを上記の漸化式で求めてみると

S₀(n) はゼロ乗の総和を意味するので n とすれば、

になります。S₂(n) は

で、先ほど計算した結果と一致します。さらに、

より、三乗の総和を計算する公式も得ることができます。

S_m(k) に対しては次の公式が成り立ちます。

k は自然数なので S_m(k) は離散値を取りますが、明らかに多項式の形をとるため実数にまで範囲を広げて微分しても問題はありません。そこで、上式を微分してみます。

この式について、0 ≤ k ≤ n - 1 の間で総和を計算すると、

よって、

ここで、先ほど得られた結果から S_m(n) は次数の最大が m + 1 になることがわかるので

とすると、S_m'( 0 ) は次数 1 の係数を表すので b_m1 (定数) になります。従って、

これをさらに微分していくと、

S_m(n) の漸化式を見ると、定数項は必ずゼロになることが分かります。従って、S_m(0) = 0 となります。よって、

以上から、S_m(n) の「マクローリン展開」は

ここで、k = m - j + 1 とすると j = m + 1 - k で、j が 1 から m + 1 まで変化するとき k は m から 0 まで変化するので

よって、S_m( n ) の m + 1 - k 次の係数が { C( m + 1, k ) / ( m + 1 ) } b_k1 という値で表されることになります。ここで、b_k1 は S_k(n) の一次項の係数になります。この数は「ベルヌーイ数 ( Bernoulli Number )」と呼ばれ、その導き方から分かるように全て有理数となります。また、べき乗の総和をベルヌーイ数を使った多項式で表したものは「ファウルハーバーの公式 ( Faulhaber's Formula )」または「ベルヌーイの公式 ( Bernoulli's Formula )」と呼ばれています。

b_k1 を新たに b_k と表すことにして、得られた漸化式に n = 1 を代入すると

これを使って b_k を求めてみます。

ここで、C( 2, 1 ) b₁ = 1, C( 3, 1 ) b₁ = 3 / 2, ... C( m + 1, 1 ) b₁ = ( m + 1 ) / 2 なので、右辺の m + 1 を減算することで必ず - ( m + 1 ) / 2 となります。これを利用して、

と変形しても、b₁ = -1 / 2 となるだけでそれ以降は上記と同じ値が得られます。ベルヌーイ数を計算する場合はこちらの漸化式が利用され、b₁ = -1 / 2 で表します。

数列 { a_k } = a₀, a₁, ... を係数とするべき級数

は、数列 { a_k } の「母関数 ( Generating Function ; または生成関数 )」といいます。数列は離散値なのに対し、母関数で表すことで連続的な値をとる関数として扱うことができるようになるので、関数の持つ性質 ( 例えば微積分や極限など ) を利用して数列の性質を解明する目的で母関数は用いられます。上記形式の母関数は「通常型母関数 ( Ordinary Generating Function )」と呼ばれ、他にも次のような種類のものがあります。

通常型母関数において { a_k } = { 1, 1, 1, ... } とすれば、G(x) = 1 + x + x² + ... となるので、これは初項 1、公比 x の等比級数を表します。第 n 項までの和は

これに x を掛けると

両式の差分を計算すると

よって、G(x,n) = ( 1 + xⁿ⁺¹ ) / ( 1 - x ) となり、|x| < 1 ならば右辺は 1 / ( 1 - x ) ( n → ∞ ) となるので、

という公式が得られます。、G(x) の両辺を微分すると、

両辺に x を掛ければ

これは、自然数を数列とした場合の母関数に対する公式になります。同様に微分を続けていけば、平方数・立方数を数列とした母関数の公式も求めることができます。ところで、二項定理から

なので、これは二項係数を数列とした場合の母関数に対する公式でもあります。

指数型母関数を見ると、{ a_k } = { 1, 1, 1, ... } とすれば EG(x) = 1 + x + x² / 2! + ... + x^k / k! + ... = Σ_k{0→∞}( x^k / k! ) となり、この右辺は指数関数 e^x のマクローリン展開です。「指数型」と名付けられたのはこの性質のためです。

ベルヌーイ数 b_k を求める漸化式

の左辺を 0 から m - 1 までの和に置き換えた式からなる数列 { a_m } = { Σ_k{0→m-1}( C( m, k ) b_k ) } の指数型母関数として、

を定義します。総和は 1 から始まるので、母関数は定数項を持ちません。このとき、

なので、

となります。一方、

母関数の積の公式 ( 補足 2 ) から上式は

こうして得られた二式から

で表すことができます。この式から、ベルヌーイ数 b_k は関数 x / ( e^x - 1 ) のマクローリン展開における、次数 k の項の係数を表すことになります。

ベルヌーイ数を求めるためのサンプル・プログラムを以下に示します。

/*
  AccurateCombination : 組み合わせ ( n 個の要素から r 個を選ぶ場合の数 ) の計算

  nCr = n! / r!( n - r )! を正確に計算する。
  計算可能な範囲は戻り値の型 Res に依存し、Res の取りうる値の範囲を超えた場合は正しい結果は得られない。

  n : 要素数
  r : 要素から選択する数
  戻り値 : 組み合わせの計算結果
*/
template< class Res, class T >
  Res AccurateCombination( T n, T r )
{
  if ( n < r ) return( Res() );   // n < r なら nCr = 0
  if ( r < 0 ) return( Res() );   // nC-r = 0

  if ( r > n - r ) r = n - r; // nCn-r = nCr

  if ( r == 0 ) return( Res( 1 ) );  // nC0 = 1
  if ( r == 1 ) return( Res( n ) );  // nC1 = n

  // まずは順列 n( n - 1 )...( n - r + 1 ) を求める
  Res ans = n;
  for ( T t = n - r + 1 ; t < n ; ++t )
    ans *= t;

  // 階乗 r! を求める
  Res fact = r;
  for ( --r ; r > 1 ; --r )
    fact *= r;

  ans /= fact;

  return( ans );
}

/*
  AccurateBernoulliNum : ベルヌーイ数 Bn を求める

  ベルヌーイ数の計算に利用する組み合わせ(二項係数)は AccurateCombination を使って正確に求めている。

  n : 何番目のベルヌーイ数か？
  戻り値 : ベルヌーイ数 Bn
*/
template< class Res, class I >
  Res AccurateBernoulliNum( I n )
{
  if ( n < I() ) return( Res() );
  if ( n == I() ) return( Res( 1 ) );

  I np( n + 1 );
  //Res b = Res( n + 1 );
  Res b = Res();
  for ( I i = I() ; i < n ; ++i )
    b -= Res( AccurateCombination< Res, I >( np, i ) ) * AccurateBernoulliNum< Res, I >( i );

  b /= Res( np );

  return( b );
}

ベルヌーイ数の計算は、関数 AccurateBernoulliNum の中で行います。また、組み合わせの計算が必要となるため、AccurateCombination 関数を実装しています。どちらも「正確に」値を計算しますが、利用する型はテンプレート引数で指定するようになっているので、例えば浮動小数点数を利用した場合は近似値として結果が得られることになります。AccurateBernoulliNum においてテンプレート引数 Res は戻り値 ( ベルヌーイ数 ) の型を表し、有理数や浮動小数点数の利用を想定しています。また、I は引数 ( 何番目のベルヌーイ数か ) を表す型で、整数型の利用を想定しています。AccurateCombination については戻り値 ( 組み合わせ ) と引数 ( 要素数と選択数 ) それぞれで型を分けていますが、どちらも整数型であることを想定しています。

ベルヌーイ数は有理数なので、「数値演算法 (5) 有理数と無理数の演算」の中の「2) 有理数の表現方法」で作成した BigNum::Rational クラスを利用することができます。このクラスを使って計算した結果は次のようになります。

n = 1 を除き、奇数番目のベルヌーイ数は全てゼロになっています。また、偶数番目のベルヌーイ数は正負が交互に入れ替わっています。

ところで、サンプル・プログラムで、結果を格納するための変数 b の初期値をゼロにしている記述の上に n + 1 で初期化する記述があって、こちらはコメントアウトされています。ここで、n + 1 で初期化する方を有効にした場合、b₁ の値の正負が逆転するだけで、他の値は全て同じになります。これは先述した通り、漸化式

において、C( m + 1, 1 ) b₁ = ( m + 1 ) / 2 となることから、両辺から m + 1 を減算した上で b₁ の符号を反転しても漸化式が成り立つためです。ここで、

とした場合の { a_m } を係数とした指数型母関数を求めると、

となるので、

但し、この時は b₁ = 1 / 2 になります。

ベルヌーイ数の名前は、ヨーロッパにおいて少なくとも八人の優れた数学者を輩出したとされる数学者一族ベルヌーイ家が由来となっています。ベルヌーイ家の一人「ヤコブ・ベルヌーイ ( Jakob Bernoulli )」が、べき乗の総和を求める問題の一般解の中でベルヌーイ数を使い、その内容は、彼の死後に出版された「推測法 ( Ars Conjectandi )」によって広く知られるようになりました。しかし、ほぼ同時期に日本の数学者「関考和」も独自に同じ結果を得て、その内容は「括要算法」の中で発表されています。

5) オイラー・マクローリンの公式 ( Euler-Maclaurin Formula )

ベルヌーイ数を係数に含んだ次のような多項式を「ベルヌーイ多項式 ( Bernoulli Polynomials )」といいます。

B_m(x) の定数項は b_m になるので、B_m(0) = b_m が成り立ちます。m の値をいくつか決めて多項式の実際の式を求めると次のようになります。

ここで、ベルヌーイ多項式に関するいくつかの性質を証明しておきます。まず、B_m(x) を微分すると、

という関係式が得られます ( べき乗の総和について得られた関係式 S_m'(n) = mS_m-1(n) + b_m とよく似た形となっています )。また、ベルヌーイ多項式に対する指数型母関数は、

ベルヌーイ数 b_k を求める漸化式 Σ_k{0→m}( C( m + 1, k ) b_k ) = 0 の指数型母関数から Σ_m{0→∞}( b_m t^m / m! ) = t / ( e^t - 1 ) であることはすでに前節で示しました。また、e^x のマクローリン展開から Σ_m{0→∞}( (tx)^m / m! ) = e^tx となるので、

と求めることができます。これを利用すると、

なので、係数を比較することで、m = 1 を除いて B_m(0) = B_m(1) が成り立つことが分かります ( m = 1 の場合は B₁(1) - B₁(0) = 1, B₁(0) = b₁ = -1/2 より B₁(0) = - B₁(1) = - 1/2 になります )。B_m(0) = b_m なので、これらはベルヌーイ数と等しいということも意味しています。また、

なので、m = 1 以外の奇数次の項は B_m(0) = B_m(1) であると同時に B_m(0) = - B_m(1) でもあり、従って B_m(0) = B_m(1) = 0 が成り立つことになります。つまり、奇数番目のベルヌーイ数は m = 1 を除いて全てゼロになります。

以上、ベルヌーイ多項式の性質をまとめると、

ある関数 f(x) の 0 から 1 までの定積分 ∫{0→1} f(x) dx において、B₀(x) = 1 なので ∫{0→1} B₀(x) f(x) dx と強引に代入します。ここで、微分に関する関係式から B'₁(x) = B₀(x) が成り立つので、部分積分法 ∫f'(x)g(x)dx = f(x)g(x) - ∫f(x)g'(x)dx を利用して

上式において、f^(m)(x) は m 階導関数を表すものとします。B₁(0) = - B₁(1) = -1/2 なので、

この操作は何度も繰り返すことができて、その度に二つの式に分解されていきます。∫{0→1} B_m(x) f^(m)(x) dx に対しては B'_m+1(x) = ( m + 1 )B_m(x) が成り立つので、

上記操作の繰り返しの中で、奇数番号のベルヌーイ数が 1 を除いて全てゼロであることから偶数番号のベルヌーイ数だけが残ることになり、以下のような式が得られることになります。

定積分を行う区間を任意の範囲 [ a, a + h ] としたとき、y = ( x - a ) / h とすれば dx = h・dy で、y は [ 0, 1 ] の区間に変換されます。また、x から y への変数変換によって、f(x) の y による微分を行うと h が係数として現れる ( df/dy = (df/dx)・(dx/dy) = h・(df/dx) ) ので、

この操作を繰り返すと

になります。区間を、[ a, a + h ], [ a + h, a + 2h ], ... , [ a + ( n - 1 )h, a + nh ] と細かく分割してた上で、上記操作を各区間ごとに行ってその総和を取ったとき、

ここで、

これは台形則による定積分値の求め方と同じなので、T_n とおきます。また、

になります。残る剰余項については、

と表します。但し [ ( x - a ) / h ] は ( x - a ) / h を超えない最大の整数を表します。これは、ベルヌーイ多項式の積分範囲が常に 0 〜 1 の範囲になることを表しています。

以上の内容から、

項を置き換えると

右辺は、台形則で求めた近似値と実際の定積分値との差を表すことになり、これによって誤差を評価することができることを意味しています。この公式を「オイラー・マクローリンの公式 ( Euler-Maclaurin Formula )」といいます。この公式は、「レオンハルト・オイラー ( Leonhard Euler )」と「コリン・マクローリン ( Colin Maclaurin )」がそれぞれ独立に発見したことからその名が付けられています。マクローリンは、例の「マクローリン展開」の名でも有名な数学者です。

ベルヌーイ数を含む総和部分を T_n から減算すれば精度が上がることになり、m を大きくすればそれだけ誤差も小さく抑えることができます。m = 1, 2 の場合の式は、次のようになります。

残念ながら、上式を利用するためには関数 f(x) の導関数を求めておく必要があるので、任意の関数に対して適用することは簡単にはできません。あらかじめ関数が決まっていて、その高階導関数に対する値が簡単に求められれば、少ない分割数でかなり精度の高い近似式が得られるようになります。例として、指数関数 e^x に対して適用した結果を以下に示します。e^x は何回微分しても関数の形が変化しないため、非常に簡単に結果を調べることができます。

「補正1」は一階導関数、「補正2」は三階導関数の値まで利用した補正を行った結果です。「補正2」の結果を見ると、10 個程度の分割でかなり精度のよい近似値が得られているので、台形則のみで計算した場合と比較すると圧倒的に早く収束させることができることが分かります。

6) スターリングの公式 (Stirling's Formula)

階乗は、順列や組み合わせの計算などでよく利用されます。階乗の計算は 1 から順番に自然数を掛け合わせていくことで計算することができますが、すぐに巨大な数値になるため正確な値を計算するためには「多倍長整数」による演算が必要になります。例えば、100 個の要素を並べる場合の数は 100! で、

100! =  93,326,215,443,944,152,681,699,238,856,266,700,490,715,968,264,381,621,468,592,
       963,895,217,599,993,229,915,608,941,463,976,156,518,286,253,697,920,827,223,758,
       251,185,210,916,864,000,000,000,000,000,000,000,000

という 158 桁の巨大な数になります。コンピュータなどなかった時代は階乗を手計算で求める必要があったので、このような乗算を毎回繰り返すようなことはしたくないわけで、必然的に近似値を求めるための手法が編み出されました。このとき、数値積分法が利用されています。

まず、階乗の対数をとります。

一方、

上記は、f(x) = x、g(x) = ln(x) として、部分積分法 ∫ f'(x)g(x) dx = f(x)g(x) - ∫ f(x)g'(x) dx を利用しています。

前者は区分求積法によって ∫{1→n} ln(x) dx の値を近似計算したものと考えれば、両者の値は近いことが推測できます。よって、

より

∫{1→n} ln(x) dx に対して台形則を利用した場合は、

となるので、

大雑把な近似式として、無視できる部分を全て除外した

という式が用いられることもあります。n が増大するにしたがって分割数も増大するので、近似値は実際の値に近づいていきます。つまり、n の増大に伴い、ln(n!) と n・ln(n) - n は値が近づいていき、n → ∞ においては両者の差が一定の値に収束することが予想され、実際にそうなります。以下で、その証明方法について説明します。

オイラー・マクローリンの公式を ∫{1→n} ln(x) dx に適用してみます。

よって、

上式の左辺および右辺の式を R(n) とします。このとき、

となるので、両式を辺々引いて

両辺を整理すると

になるので、この式に対して n → ∞ のときの極限を考えます。まず、R(n) について、R(n) - R(n + 1) を計算すると、

ネイピア数 e の定義から lim{n→∞}( 1 + 1/n )ⁿ = e なので、上記の値は n → ∞ でゼロに収束します。つまり、R(n) はある値に収束することになるので、その値を C とします。また、

となるため、

になります。左辺は、以下に示す「ウォリス積 ( 補足 3 )」を使って正確な値を求めることができます ( 先程 ln( 2・4・...・2n ) と ln( ( 2n + 1 )! ) が唐突に現れたのはこれを利用するためです )。

ウォリス積を利用すると、

となるので、

(左辺) = (右辺) なので、C の値を求めることができて、

よって、

となって、

つまり、n が充分に大きければ、以下の近似式が成り立つことになります。

この式は、「スターリングの公式 ( Stirling's Formula )」として知られる有名な近似式です。しかし、公式を最初に発見したのは「ド・モアブルの定理」で有名な「アブラーム・ド・モアブル ( Abraham de Moivre )」で、最初は

であることまでが示されていました。スターリングの名は、その後「ジェイムズ・スターリング ( James Stirling )」が C = ( 2π )^1/2 であることを示したことから付けられたようです。

ところで、

でもあり、これを使って R(n) - R(M) (但し n < M) を計算すると、

M → ∞ とすれば、R(M) → ( 2π )^1/2 であり、右辺の M が入った項は無視することができるので、

と表すことができます。右辺の総和部分のいくつかを展開してみると

になります。よって、スターリングの公式は次のように表すこともできます。

この式は、n がそれほど大きくない場合でもある程度高い精度が得られる近似式になります。

実際に計算を行った結果を以下に示します。対象の値は 10!, 100!, 200!, 500!, 1000! の 5 種類です。

各表の一番上が正確な値 ( 但し有効桁数は 15 桁です ) を表しており、近似計算結果の右側には正確な値との比率が示されています。数が大きくなるほど、どの計算式についても誤差が小さくなっていることが分かります。上側にある近似式 ( n / e )ⁿ と e ( n / e )ⁿ は大雑把すぎてあまり正確とはいえませんが、他の式に関してはかなり近い値となっています。特に、下側の二つは 10! ですでにほぼ一致しているため、比較的小さい数に対しても精度が高いことがわかります。

ところで、上記検証のために利用した変数は long double 型で、この型によって表現可能な有効桁数は 4900 桁強なので、計算可能な階乗の数は 1700! 程度になります。10000! にもなると、桁数が 35000 を超えるほど巨大な数になります。そのままでは扱いづらいので、通常は対数で表現することになります。

今回は、区分求積法による定積分値の計算法からはじまって、スターリングの公式が区分求積法の考え方を逆利用していることを最後に紹介しました。実は、精度の高い公式は他にもあり、またほとんどの関数に対して収束を早めることのできる万能型の変数変換法なども存在します。さらに広義積分や重積分に対する手法など、なかなか奥深いテーマでもあるので、次の機会にまた取り上げてみたいと思います。

補足 1) 定積分の定義

f(x) は実数値関数で、a ≤ x ≤ b ( つまり閉区間 I = [ a, b ] ) の範囲で連続であるとします。また、この区間では f(x) が有限な範囲の値 M 以下になる( |f(x)| ≤ M )とします。

区間 I 上に有限個の点 a = x₀ < x₁ < ... < x_n-1 < x_n = b をとり、その点によって区間 I を分割します。この点を「分点」といいます。分割した各区間の中から代表値を決めて、区間 [ x_k-1, x_k ] 上の点を ξ_k ( k = 1, 2, ... n ) としたとき、

は、区間 I の範囲の f(x) の面積に対する近似値となります。これを「リーマン和 ( Riemann Sum )」といいます。分割の仕方も代表値の決め方も任意であり、等間隔に分割する必要も決められた方法で代表値を選ぶ必要もないので、その方法によってリーマン和の値は変化します。

分割の仕方を Δ で表し、各区間での上限を M_k、下限を m_k として、これらを代表値としたときのリーマン和を

とすると、-M ≤ m_k ≤ M_k ≤ M なので

が成り立ちます。

分割の仕方が Δ₁ と Δ₂ の二つあって、それぞれに対する上で定義したリーマン和を R_M(Δ₁), R_m(Δ₁) と R_M(Δ₂), R_m(Δ₂) としたとき、二つの分割に対する分点を全て含む分割 Δ₃ のリーマン和 R_M(Δ₃), R_m(Δ₃) との大小関係はどのようになるかを考えます。Δ₁ を基準に考えると、Δ₁ の分点の間に Δ₂ の分点が存在することになります。Δ₁ のある区間 [ x_1,m, x_1,m+1 ] の間に Δ₂ の分点 x_2,n, x_2,n+1, ... x_2,n+k があって、区間 [ x_1,m, x_1,m+1 ] の上限と下限はそれぞれ M₁, m₁、区間 [ x_2,n+j-1, x_2,n+j ] のそれらは M_2j, m_2j ( j = 1, 2, ... k ) とすると、

が成り立ちます。また、逆に Δ₂ を基準とすれば、同様に考えて

となります。全区間の総和をとると、中辺は上側と下側のどちらの式についても R_m(Δ₃) ≤ R_M(Δ₃) であることを表しているので、

つまり、任意の分割 Δ, Δ' に対して R_m(Δ) ≤ R_M(Δ') となります。簡単に言えば、どのように分割の仕方を変えてみても、総和の上限と下限について大小関係が逆になるような分割の組み合わせは存在しないということです。

図 1-1. リーマン和の大小関係

f(x) は区間内で有限な値をとるとしているので、分割をいろいろと変化させる中で必ず R_M(Δ) の中の下限と R_m(Δ) の中の上限が存在します。それぞれ inf( R_M ), sup( R_m ) ( inf は 下限 ( infimum )、sup は 上限 ( supremum ) の数学用語です ) とすると、上記内容から

になります。ここで等号が成り立つ時、f(x) は「積分可能」または「可積分」であるといって、その共通の値が f(x) の閉区間 [ a, b ] 上の「定積分」になり、

 b
∫ f(x) dx ( このドキュメントの中では ∫{a→b} f(x) dx と表現しています )
 a

で表されます。例えば f(x) = K (定数) ならば、閉区間 [ a, b ] 上の R_Δ はどのような分割の仕方に対しても常に K( b - a ) になるので等号が成り立って可積分であることになります。

分割 Δ の中で、最大となる区間の幅を「Δ の幅」と呼んで、||Δ|| で表します。f(x) が可積分であれば、分割の仕方や代表値の選び方に依存せず、R_Δ は ||Δ|| → 0 において一定の値に収束し、逆に lim{||Δ||→0}( R_Δ ) = α ならば、f(x) は可積分で α = ∫{a→b} f(x) dx が成り立ちます。||Δ|| を小さくすると、inf( R_M ), sup( R_m ) と ∫{a→b} f(x) dx の差がいくらでも小さくなることを厳密に証明した定理は「ダルブーの定理 ( Darboux's Theorem )」と呼ばれ、R_M, R_m はそれぞれ「上ダルブー和 ( Upper Darboux Sum )」、「下ダルブー和 ( Lower Darboux Sum )」、さらに inf( R_M ), sup( R_m ) が「上ダルブー積分 ( Upper Darboux Integral )」、「下ダルブー積分 ( Lower Darboux Integral )」と名付けられています。リーマン和 R_Δ は sup( R_m ) と inf( R_M ) の間の値をとるので、lim{||Δ||→0}( R_Δ ) = ∫{a→b} f(x) dx が成り立つことになるわけです。

積分可能な関数の条件や、上記説明の厳密な証明法については、解析学の入門書に書かれており、また他のホームページなどでも公開されているので、気になる方は探してみて下さい。参考にした書籍やホームページは、このページの最後に記載しています。

補足 2) 母関数の積

二つの母関数 f(x) = Σ_n{0→∞}( a_nxⁿ ) と g(x) = Σ_n{0→∞}( b_nxⁿ ) の積は

よって、

が成り立ちます。係数の連番と次数は一致するので、積の中の項の係数は全て、a_i と b_j の連番 i, j の和がその次数と等しくなります。よって、上記のように展開することができるわけです。

補足 3) ウォリス積 ( Wallis Product )

ウォリス積は「数値演算法 (5) 有理数と無理数の演算」の中でも紹介しており、sin x の乗積展開を使った証明を示してあります ( 証明法 )。ここでは、別の証明法を示します。

とおくと、部分積分法を利用して

となるので、漸化式

が得られます。よって、

より

S_2n+1 / S_2n の、n → ∞ での極限を考えると、0 < x < π/2 では 0 < sinx < 1 なので 0 < S_2n+1 < S_2n < S_2n-1 が成り立って、

よって、S_2n+1 / S_2n → 1 ( n → ∞ ) になります。最後に S₀ = ∫{0→π/2} 1 dx = π/2, S₁ = ∫{0→π/2} sinx dx = 1 なので、

これで、ウォリス積が証明されました。ちなみに、この一連の証明は、その一部が頻繁に大学入試などで出題されているそうです。

「ウォリス積」の名はイギリスの数学者「ジョン・ウォリス ( John Wallis )」から付けられたものです。ウォリスは、円の面積を求める定積分値 4∫{0→1} ( 1 - x² )^1/2 dx が、二つの定積分値 ∫{0→1} ( 1 - x² )⁰ dx と ∫{0→1} ( 1 - x² )¹ dx の間を補間する形で得られるのではないかと類推しました。前者の値は 1、後者の値は 2/3 となり、その幾何平均は [ 1 x ( 2 / 3 ) ]^1/2 = ( 2 / 3 )^1/2 ≅ 0.8165 で、π / 4 ≅ 0.7854 なので、それぞれ近い値となります。この考え方を発展させて、ウォリス積が誕生したそうです。

1. 「数値計算の常識」 伊理正夫・藤野和建 著 (共立出版)
2. 「基礎課程 解析入門」 野本久夫・岸正倫 著 (サイエンス社)
3. 「はじめての数論」 ジョセフ・H・シルヴァーマン 著 (ピアソン・エデュケーション)
4. 「アトムの物理ノート」(http://letsphysics.blog17.fc2.com/) : オイラー・マクローリンの総和公式 (証明)
5. 「数学のページ」(http://www.igaris.com/math/) : べき乗和の公式(pdf)
6. MIT COMPUTER SCIENCE AND ARTIFICIAL INTELLIGENCE LABORATORY - Kuat Yessenov
   Kuat's course projects - [18.704] A short expository paper on Euler-Maclaurin formula(pdf)
7. Wikipedia

前に戻る

タイトルに戻る