7.1. 数値積分の方法

7.1.1. 概観

この章では1変数の積分

(1)I=\int_{a}^{b} f(x) dx

を計算する方法を解説します。

以下の2つの場合で最適なアルゴリズムは変わります。

  1. f(x) が任意の点 x で計算できる → 求積法

  2. f(x) が離散的な点 x_i 上でのみ与えられている → 離散積分公式

最初のケースは f(x) が解析的に与えられている場合です。この場合は原理的には積分値を任意の精度まで求められます。 第二のケースは f(x_i) が何らかの数値計算を通して得られる場合、もしくは実験データの場合です。 積分の精度には限界があります。また、f(x_i) が誤差を含む場合は、高次の（複雑な）公式を使うほど精度が上がるわけではありません。

さて、pythonでは数値積分は scipy.integrate モジュールを使って実行できます。 数値積分のアルゴリズムの説明に移る前に、scipy.integrate に含まれる関数を見ておきましょう。

Methods for Integrating Functions given function object.

  quad          -- General purpose integration.
  dblquad       -- General purpose double integration.
  tplquad       -- General purpose triple integration.
  fixed_quad    -- Integrate func(x) using Gaussian quadrature of order n.
  quadrature    -- Integrate with given tolerance using Gaussian quadrature.
  romberg       -- Integrate func using Romberg integration.

Methods for Integrating Functions given fixed samples.

  trapz         -- Use trapezoidal rule to compute integral from samples.
  cumtrapz      -- Use trapezoidal rule to cumulatively compute integral.
  simps         -- Use Simpson's rule to compute integral from samples.
  romb          -- Use Romberg Integration to compute integral from
                   (2**k + 1) evenly-spaced samples.

通常はquad関数（quadrature, 求積法）を使えば問題ありません。 離散的なデータしかない持っていない場合は、trapz（trapezoidal rule, 台形公式）やsimps（simpson's rule, シンプソンの公式）を使うことになります。

7.1.2. 離散積分公式

離散的な点 x_i でデータが与えられている場合には、多項式で内挿して積分します。 内挿式を1次式、2次式、3次式・・・のように高次のものにすることで、誤差は小さくなっていきます。 1次式を使った場合、つまり単に f(x_i) を直線で結んだ場合の近似式は台形公式と呼ばれます。 x_i が等間隔に区切られている場合、台形公式は以下で与えられます：

I \simeq \frac{h}{2} \left[ f(x_0) + 2f(x_1) + 2f(x_2) + \cdots + 2f(x_{N-1}) + f(x_N) \right] + \mathcal{O}(h^2)

ここで、h = x_{i+1} - x_{i} です。 次に、f(x_{i-1}), f(x_{i}), f(x_{i+1}) （i は奇数）の3点を使って2次式で内挿して積分するとシンプソンの公式が得られます：

I \simeq \frac{h}{3} \left[ f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + 2f(x_{N-2}) + 4f(x_{N-1}) + f(x_N) \right] + \mathcal{O}(h^4)

このように、公式ごとに f(x_i) を足し合わせる重みが変わります。

同様にして、より高次の公式も導けますが、3次式以上の内挿公式を使うのであれば、スプライン等で内挿して積分を求めた方が簡単です。 スプラインは各区間を3次式で近似し、隣同士が滑らかに接続するように係数を決めたものです。 各区間は3次式なので、積分値は解析的に求まります。 内挿については scipy.interpolate モジュールの公式ドキュメントを参照してください。

7.1.3. ガウス求積法

ガウス求積法（Gaussian quadrature） は、関数 f(x) が任意の点で計算できる場合に利用できる、非常に効率の良い数値積分法です。 (1) 式の積分を

I \approx \sum_{i=1}^n w_i f(x_i)

の形で近似することを考えます。上述の離散積分公式の場合、座標点 x_i は等間隔です。台形公式とシンプソン公式では重み w_i が異なり、それによって精度が変わります。ここでは、さらに座標点 x_i も変えられるとします。 最適な座標点 x_i と重み w_i の組はどのようなものでしょうか。

関数 f(x) が多項式の場合に、少ない座標点で厳密な結果を与えるのがガウス・ルジャンドルの公式です。 より正確には、f(x) が 2n-1 次の多項式の場合に、x_i として、n 次のルジャンドル多項式 P_n(x) のゼロ点（n 個ある）を取り、w_i を適切に選ぶと積分を厳密に評価できます。

どのくらい効率が良いかを理解するために、例えば、被積分関数 f(x) が3次の多項式の場合を考えます。 関数 f(x) の各項の係数を決定するには4つの点が必要です。したがって、4点あれば（座標はどこでもよい）、問題なく積分を厳密に計算できます。 ガウス・ルジャンドルの公式は、座標点を適切に選べば、たった2点だけで厳密な積分値が得られると言っています。 2つの点からは直線しか決まらないので、3次多項式を直線で近似しているように思いますが、この情報だけで積分を厳密に評価できてしまうというのは驚きです。

これに関する直感的な説明は、英語版の Wikipadia - Gaussian quadrature の図が分かりやすいです。

7.1.4. scipy.integrateの使い方

積分はSciPyの scipy.integrate モジュールを使って計算できます。 1変数積分の求積法は scipy.integrate.quad 関数を使います。 公式ドキュメントによるとquad関数の引数は次のようになってます：

scipy.integrate.quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08,
                     epsrel=1.49e-08, limit=50, points=None, weight=None,
                     wvar=None, wopts=None, maxp1=50, limlst=50)

引数がたくさんありますが、とりあえず最初の4つだけ見ておけば問題ありません。 funcは被積分関数、aとbはそれぞれ積分の上限と下限です。 funcに与える関数は以下のように、x を受け取って f(x) を返す関数です：

def integrand(x):
    処理
    return f

f(x) がパラメーターを取る場合（param1, param2とする）には

def integrand(x, param1, parama2):
    処理
    return f

と定義します。 パラメーターの値は args=(0.1, 2,5) というようにしてquad関数に渡します。

特異性がある場合

被積分関数に積分可能な特異性がある場合（\log (x-x_0) や 1/(x-x_0)）には注意が必要です。 事前に、変数変換で特異性を除いておくのが確実です。

scipy.integrate.quad では weight 引数を使って、特異性のある関数形を指定することで積分できる場合があります。例えば、f(x)=g(x)/(x-a) の場合、weight='cauchy' として、関数 g(x) をquad関数に与えます。 また、points 引数で特異性や不連続のある座標を与えることで、特異性を適切に扱ってくれるようです。 詳細は公式ドキュメント scipy.integrate.quad を参照してください。