$\int_0^1 e^{-x^2}\,dx$ の値はいくつでしょうか。この積分は初等関数で表せないため、不定積分を求めて代入するという教科書的な方法は使えません。しかし、数値的には任意の精度で計算できます。このような「手で解けない積分」は、理論研究でもエンジニアリングの実務でも頻繁に現れます。たとえば正規分布の累積分布関数 $\Phi(x) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-t^2/2}\,dt$ も同種の積分であり、統計学のあらゆる場面で数値的に計算されています。
数値積分（numerical integration, quadrature）は、定積分の値を近似的に計算する手法です。”quadrature” という言葉は「正方形化」を意味し、古代ギリシャで曲線で囲まれた面積を正方形に等しい面積に変換する問題に由来します。現代では、解析的に求まらない積分、高次元の積分、被積分関数が離散データとしてのみ与えられる場合など、実用上不可欠な技術です。
数値積分の考え方は非常にシンプルです。関数の値をいくつかの点でサンプリングし、それらの値を「うまく」重み付けして足し合わせることで積分値を近似します。日常的なイメージでいえば、不規則な形の土地の面積を求めるために、測量点を何か所か設定して長方形や台形で近似するようなものです。ここで「うまく」の部分に数学的な工夫が詰まっており、それが台形則やシンプソン則やガウス求積法といった手法の違いとなって現れます。
数値積分を理解すると、以下のような場面で活用できます。
- 物理シミュレーション: 軌道計算、電磁場計算での面積分・体積積分。たとえばロケットの推力プロファイルから得られる速度変化量 $\Delta v = \int_0^T \frac{F(t)}{m(t)}\,dt$ は一般に解析的に解けず、数値積分が不可欠です
- 統計学: 確率分布の期待値計算、ベイズ推論での周辺化。事後分布の正規化定数 $\int p(D|\theta)p(\theta)\,d\theta$ の計算は多くの場合、数値積分に頼ります
- 有限要素法: 要素剛性行列 $\bm{K}^e = \int_{\Omega^e} \bm{B}^T \bm{D}\bm{B}\,d\Omega$ の計算にガウス求積法が標準的に使われています
- 信号処理: スペクトル解析での数値フーリエ変換。連続信号のフーリエ変換 $\hat{f}(\omega) = \int_{-\infty}^{\infty} f(t)e^{-i\omega t}\,dt$ も離散近似です
本記事の内容
- 数値積分の基本的な考え方
- 台形則の導出と誤差解析
- シンプソン則の導出と誤差解析
- ガウス求積法の理論
- 適応型求積法の考え方
- Pythonでの実装と精度比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
数値積分の基本的な考え方
求積公式の一般形
数値積分の根底にある発想は、「関数の値をいくつかの点で調べて、それらを足し合わせれば積分値がわかるはずだ」というものです。積分とは結局のところ面積ですから、関数の高さ $f(x_i)$ に何らかの「幅」を表す重み $w_i$ を掛けて合計すれば、面積の近似値が得られます。
この考え方を定式化すると、ほぼすべての数値積分公式は次の形で表せます。
$$ \int_a^b f(x)\,dx \approx \sum_{i=0}^n w_i f(x_i) $$
$x_i$ は求積点（nodes）、$w_i$ は重み（weights）です。手法の違いは、$x_i$ と $w_i$ をどう選ぶかの違いにほかなりません。求積点をどこに配置するか、各点にどれだけの重みを与えるかで、近似の精度が劇的に変わります。
この一般的な枠組みから、なぜ「重みの選び方」が重要なのかを直感的に理解しておきましょう。たとえば、すべての点に等しい重みを与える方法（矩形則）は最も素朴ですが、関数の端の情報を正しく反映できません。端点を半分の重みにするだけで精度が大きく向上するのが台形則であり、このような「たった少しの工夫」が数値積分の研究を奥深いものにしています。
ニュートン=コーツ型とガウス型
数値積分の手法は、大きく分けて2つのアプローチに分類できます。
- ニュートン=コーツ型: 等間隔の求積点を使う（台形則、シンプソン則）。実装が容易で、離散データ（実験データなど）にもそのまま適用できるのが利点です
- ガウス型: 求積点と重みの両方を最適化する（ガウス=ルジャンドル求積）。関数を自由に評価できる場合、同じ点数でニュートン=コーツ型より格段に高い精度を達成できます
ニュートン=コーツ型は「測定点が等間隔に決まっている」場合に自然な選択であり、ガウス型は「測定点を自由に選べる」場合に最適な選択です。この使い分けを理解しておくことが、実際の問題に適切な手法を選ぶための鍵となります。
まずは最も基本的で直感的な台形則から始めましょう。
台形則
直感的な理解
台形則の考え方は、グラフ用紙に描かれた曲線の下の面積を、定規だけで測ろうとする場面をイメージするとわかりやすいでしょう。隣り合う2点を直線でつなぎ、その直線と $x$ 軸で囲まれた台形の面積を合計すれば、曲線の下の面積の近似値が得られます。
曲線が直線に近いほど近似は正確になりますから、細かく分割するほど各小区間での曲線は直線に近づき、精度が上がることが直感的にわかります。問題は「どのくらい細かくすれば十分か」であり、それを教えてくれるのが誤差解析です。
導出
区間 $[a, b]$ を $n$ 等分し、$h = (b-a)/n$, $x_k = a + kh$ とします。各小区間 $[x_k, x_{k+1}]$ で $f(x)$ を1次関数（直線）で近似すると
$$ \int_{x_k}^{x_{k+1}} f(x)\,dx \approx \frac{h}{2}[f(x_k) + f(x_{k+1})] $$
これは台形の面積です。全区間で足し合わせると複合台形則が得られます。
$$ \begin{equation} \int_a^b f(x)\,dx \approx T_n = \frac{h}{2}\left[f(x_0) + 2\sum_{k=1}^{n-1} f(x_k) + f(x_n)\right] \end{equation} $$
端点 $f(x_0)$ と $f(x_n)$ の重みが $1/2$、内点の重みが1であることに注目してください。この重み構造を直感的に理解するには、各内点が左右2つの台形に共有されていることを考えるとよいでしょう。左の台形からの寄与 $1/2$ と右の台形からの寄与 $1/2$ を合わせると、内点の重みは1になります。一方、端点は1つの台形にしか属さないため、重みは $1/2$ のままです。
この重みの構造は、台形則が「端点を特別扱い」するという特徴を示しています。実は、この端点の扱い方の工夫が、単純な矩形則（すべての重みが等しい）よりも高い精度をもたらしているのです。
誤差解析
台形則がどの程度正確なのかを定量的に評価してみましょう。各小区間 $[x_k, x_{k+1}]$ での台形則の誤差は、テイラー展開から
$$ \int_{x_k}^{x_{k+1}} f(x)\,dx – \frac{h}{2}[f(x_k) + f(x_{k+1})] = -\frac{h^3}{12}f”(\xi_k) $$
ここで $\xi_k \in [x_k, x_{k+1}]$ です。この誤差が $f”$（2階微分）に比例していることには意味があります。1次関数近似では曲がり具合（曲率）が無視されるため、関数がどれだけ曲がっているかを表す $f”$ が誤差に直結するのです。直線に近い関数（$f”$ が小さい関数）ほど台形則の近似は正確であるということを、この式は定量的に示しています。
$n$ 個の小区間の誤差を合計すると
$$ E_T = \int_a^b f(x)\,dx – T_n = -\frac{(b-a)h^2}{12}f”(\xi) $$
ここで $\xi \in [a, b]$ は中間値の定理によって存在が保証される点です。合計の際に $nh = b – a$ の関係を使っていることに注意してください。
重要な結論: 台形則の誤差は $O(h^2)$ です。刻み幅を半分にすると誤差は約1/4になります。たとえば $h = 0.1$ で $10^{-4}$ 程度の誤差がある場合、$h = 0.01$ にすれば $10^{-6}$ 程度の誤差になることが期待できます。この予測可能な振る舞いは、必要な精度を達成するために何点の評価が必要かを事前に見積もるうえで非常に重要です。
オイラー=マクローリン公式
台形則の誤差のより詳しい構造は、オイラー=マクローリン公式で明らかになります。
$$ T_n = \int_a^b f(x)\,dx + \sum_{k=1}^{p} \frac{B_{2k}}{(2k)!}h^{2k}[f^{(2k-1)}(b) – f^{(2k-1)}(a)] + O(h^{2p+2}) $$
ここで $B_{2k}$ はベルヌーイ数（$B_2 = 1/6$, $B_4 = -1/30$, …）です。誤差が $h$ の偶数乗のみを含むことが重要で、これはリチャードソン外挿の理論的基盤になります。
リチャードソン外挿の基本的なアイデアは次の通りです。刻み幅 $h$ での台形則の値を $T(h)$、$h/2$ での値を $T(h/2)$ とすると、オイラー=マクローリン公式から
$$ T(h) = I + c_1 h^2 + c_2 h^4 + \cdots $$
$$ T(h/2) = I + c_1 \frac{h^2}{4} + c_2 \frac{h^4}{16} + \cdots $$
が成り立ちます。ここで $I$ は厳密な積分値です。$h^2$ の項を消去するように線形結合をとると
$$ \frac{4T(h/2) – T(h)}{3} = I + O(h^4) $$
となり、精度が $O(h^2)$ から $O(h^4)$ に向上します。実はこの結果はシンプソン則に一致しており、シンプソン則は台形則のリチャードソン外挿という解釈もできるのです。
台形則は単純ですが、誤差の構造を理解することで精度を体系的に改善できることがわかりました。次に、最初から2次関数近似を用いることでより高い精度を達成するシンプソン則を見てみましょう。
シンプソン則
直感的な理解
台形則は関数を直線で近似しましたが、曲がった関数には曲がった近似を使うほうが自然です。シンプソン則は、隣り合う3点を通る放物線（2次関数）で関数を近似し、その放物線の下の面積を計算します。
たとえば、山の断面を紙に描いた場面をイメージしてください。定規で直線的に近似するよりも、柔軟な曲線定規を当てたほうが山の形状をよく捉えられます。シンプソン則はまさにこの「曲線定規」に相当する方法です。
導出
台形則が1次関数近似だったのに対し、シンプソン則は2次関数（放物線）近似です。3点 $x_0, x_1 = (x_0 + x_2)/2, x_2$ を通る放物線を積分すると
$$ \int_{x_0}^{x_2} f(x)\,dx \approx \frac{h}{3}[f(x_0) + 4f(x_1) + f(x_2)] $$
ここで $h = (x_2 – x_0)/2$ です。重みが $1:4:1$ の比であることに注意してください。中点 $x_1$ に最も大きな重み4がかかっているのは、放物線を通す際に中点の情報が曲率を決定するのに最も重要だからです。
この公式の名前は、18世紀のイギリスの数学者トーマス・シンプソンに由来しますが、同様の公式はニュートンやケプラーによっても独立に発見されていました。実際、ドイツ語圏では「ケプラーの樽の公式」（Keplersche Fassregel）とも呼ばれます。ケプラーはワイン樽の体積を計算するためにこの公式を使ったと伝えられています。
区間全体を $n$（偶数）等分した複合シンプソン則は
$$ \begin{equation} S_n = \frac{h}{3}\left[f(x_0) + 4\sum_{k=0}^{n/2-1} f(x_{2k+1}) + 2\sum_{k=1}^{n/2-1} f(x_{2k}) + f(x_n)\right] \end{equation} $$
奇数番目の点に重み4、偶数番目の点に重み2がかかっています。
誤差解析
2次関数近似なので誤差は $O(h^3)$ と思うかもしれませんが、実は1次高い $O(h^4)$ になります。これは偶然の一致ではなく、シンプソン則が対称な公式であることに起因しています。
この「ボーナス精度」を理解するために、3次の単項式 $f(x) = x^3$ で確認してみましょう。$[-h, h]$ でシンプソン則を適用すると
$$ \frac{h}{3}[f(-h) + 4f(0) + f(h)] = \frac{h}{3}[-h^3 + 0 + h^3] = 0 $$
一方、厳密な積分も $\int_{-h}^{h} x^3\,dx = 0$ です。奇関数を対称区間で積分するとゼロになりますが、シンプソン則の対称な重み構造がこの性質を自然に保存しているのです。このため、2次関数近似でありながら3次の多項式も厳密に積分でき、結果として誤差は4階微分 $f^{(4)}$ に依存します。
$$ E_S = -\frac{(b-a)h^4}{180}f^{(4)}(\xi) $$
シンプソン則の誤差は $O(h^4)$ です。台形則に比べて2次高い精度が得られます。具体的な数値で感覚をつかむと、$h = 0.1$ のとき台形則の誤差が $10^{-4}$ のオーダーなら、シンプソン則の誤差は $10^{-8}$ のオーダーまで小さくなります。同じ精度を達成するために必要な分割数が大幅に減るため、計算コストの面でもシンプソン則は台形則より優れています。
シンプソンの3/8則
4点を使うシンプソンの3/8則もあります。
$$ \int_{x_0}^{x_3} f(x)\,dx \approx \frac{3h}{8}[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)] $$
これは3次ラグランジュ補間を積分したもので、精度は通常のシンプソン則と同じ $O(h^4)$ です。3/8則の名前は、先頭の係数 $3h/8$ に由来しています。通常のシンプソン則（1/3則）が $n$ が偶数のときにしか適用できないのに対し、3/8則と組み合わせることで奇数分割にも対応できる点が実用上の利点です。
ここまでで、ニュートン=コーツ型の2大手法を見てきました。台形則（1次近似、$O(h^2)$）からシンプソン則（2次近似、$O(h^4)$）へと精度が向上することがわかりました。しかし、さらに高次のニュートン=コーツ公式（4次、5次の多項式近似）を使えば際限なく精度が上がるのでしょうか？
実は、ニュートン=コーツ型には限界があります。7点以上の等間隔ニュートン=コーツ公式では負の重みが現れ、数値的に不安定になることが知られています（ルンゲ現象の積分版です）。この問題を根本的に解決するのが、求積点の位置自体を最適化するガウス求積法です。
ガウス求積法
基本的な考え方
ニュートン=コーツ型では求積点の位置が等間隔に固定されていたため、自由に調整できるのは重み $w_i$ だけでした。しかし、「どこで関数を評価するか」も自由に選べるなら、その追加の自由度を活かしてより高い精度を達成できるはずです。これがガウス求積法の出発点です。
$n+1$ 個の点で関数を評価するとき、ニュートン=コーツ型は $2n-1$ 次（台形）から $2n$ 次（シンプソン）の多項式に対して厳密です。しかし、求積点の位置も自由に選べるなら、$2n+1$ 次の多項式に対して厳密にできます。
$$ \int_{-1}^{1} f(x)\,dx \approx \sum_{i=0}^{n} w_i f(x_i) $$
求積点 $x_i$ と重み $w_i$ の合計 $2(n+1)$ 個のパラメータを、$f(x) = 1, x, x^2, \ldots, x^{2n+1}$ の $2(n+1)$ 個の条件から決定します。
具体例で見てみましょう。$n = 0$（1点求積）の場合、パラメータは $x_0$ と $w_0$ の2つです。$f(x) = 1$ と $f(x) = x$ の2条件から
$$ w_0 \cdot 1 = \int_{-1}^{1} 1\,dx = 2, \quad w_0 \cdot x_0 = \int_{-1}^{1} x\,dx = 0 $$
と求まり、$w_0 = 2$, $x_0 = 0$ が得られます。つまり1点ガウス求積は「区間の中点で関数を評価し、区間の長さを掛ける」という中点則に一致します。
$n = 1$（2点求積）の場合、$x_0, x_1, w_0, w_1$ の4つのパラメータを $f(x) = 1, x, x^2, x^3$ の4条件で決定します。計算すると $x_{0,1} = \pm 1/\sqrt{3}$, $w_{0,1} = 1$ が得られます。たった2点で3次多項式まで厳密に積分できるのは驚くべきことです。
ガウス=ルジャンドル求積の理論
具体的にパラメータを求めるには連立方程式を解く必要がありますが、驚くべきことに、最適な求積点はルジャンドル多項式 $P_{n+1}(x)$ の零点に一致するという美しい結果が知られています。
$n+1$ 点のガウス=ルジャンドル求積は $2n+1$ 次までの多項式を厳密に積分します。
| 点数 $n+1$ | 精度（次数） | 求積点 | 重み |
|---|---|---|---|
| 1 | 1次 | $0$ | $2$ |
| 2 | 3次 | $\pm 1/\sqrt{3}$ | $1, 1$ |
| 3 | 5次 | $0, \pm\sqrt{3/5}$ | $8/9, 5/9, 5/9$ |
なぜルジャンドル多項式なのか
求積公式が $2n+1$ 次の多項式に対して厳密であるための条件は
$$ \int_{-1}^{1} q(x)\omega(x)\,dx = 0, \quad \forall q \in \mathcal{P}_n $$
ここで $\omega(x) = \prod_{i=0}^n (x – x_i)$ は節多項式です。
これは $\omega$ が次数 $n$ 以下のすべての多項式と直交することを意味します。$[-1, 1]$ 上の重み1に対する直交多項式はルジャンドル多項式 $P_{n+1}$ なので、$\omega(x) \propto P_{n+1}(x)$ となり、求積点は $P_{n+1}$ の零点です。
重みの計算
求積点 $x_i$ が決まれば、重みはラグランジュ補間から
$$ w_i = \int_{-1}^{1} \prod_{j \neq i} \frac{x – x_j}{x_i – x_j}\,dx $$
で計算できます。すべての重みは正であることが示せ、これは数値的安定性にとって重要な性質です。重みがすべて正であるということは、関数値の誤差が打ち消し合うのではなく加算されるため、丸め誤差に対してロバストであることを意味します。
一方、高次のニュートン=コーツ公式では負の重みが現れ得るため、大きな値と小さな値の差を計算する桁落ちが起こりやすくなります。ガウス求積法ではこの問題が原理的に発生しない点が、実用上の大きな利点です。
一般区間への変換
ガウス=ルジャンドル求積は $[-1, 1]$ 上で定義されますが、変数変換
$$ x = \frac{b-a}{2}t + \frac{b+a}{2}, \quad dx = \frac{b-a}{2}dt $$
により一般の区間 $[a, b]$ に適用できます。
$$ \int_a^b f(x)\,dx = \frac{b-a}{2}\int_{-1}^{1} f\left(\frac{b-a}{2}t + \frac{b+a}{2}\right)\,dt $$
ガウス求積の拡張
重み関数付きの積分にも対応できます。
| 名称 | 積分 | 重み関数 | 直交多項式 |
|---|---|---|---|
| ガウス=ルジャンドル | $\int_{-1}^{1} f(x)\,dx$ | $1$ | ルジャンドル |
| ガウス=チェビシェフ | $\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}}\,dx$ | $(1-x^2)^{-1/2}$ | チェビシェフ |
| ガウス=ラゲール | $\int_0^{\infty} e^{-x}f(x)\,dx$ | $e^{-x}$ | ラゲール |
| ガウス=エルミート | $\int_{-\infty}^{\infty} e^{-x^2}f(x)\,dx$ | $e^{-x^2}$ | エルミート |
半無限・無限区間の積分にも自然に対応できるのがガウス求積の大きな利点です。たとえば、量子力学での水素原子の波動関数の規格化条件 $\int_0^{\infty} |R_{nl}(r)|^2 r^2\,dr = 1$ では半無限区間の積分が必要ですが、ガウス=ラゲール求積を使えば効率的に計算できます。また、統計力学での状態密度の計算 $\int_{-\infty}^{\infty} g(\epsilon) f(\epsilon)\,d\epsilon$ にはガウス=エルミート求積が適しています。
ガウス求積法の理論的な美しさは、直交多項式・近似理論・線形代数が見事に融合している点にあります。ここまでの理論を踏まえたうえで、各手法の理論が揃ったところで、Pythonで精度を比較してみましょう。
適応型求積法
動機
ここまでの手法では、区間全体で一様な刻み幅を使っていました。しかし、実際の被積分関数は区間内で一様に振る舞うとは限りません。ある領域では非常に滑らかで粗い分割で十分なのに、別の領域では急激に変化して細かい分割が必要になることがあります。
たとえば $\int_0^{10} \frac{\sin(100x)}{x + 0.01}\,dx$ のような積分では、$x$ が小さい領域で被積分関数が急激に振動しますが、$x$ が大きい領域では振動が緩やかです。一様な刻み幅で全区間を細かく分割するのは計算の無駄であり、振動が激しい領域だけ細かくすべきです。適応型求積法は、この「必要な場所だけ細かくする」戦略を自動的に実現します。
アルゴリズム
適応型求積法の基本的な考え方は、「2つの精度の異なる近似を比較して、誤差が十分小さいかを判定する」というものです。具体的な手順は以下の通りです。
- 区間 $[a, b]$ 全体でシンプソン則 $S_1$ を計算する
- 区間を2分割し $[a, (a+b)/2]$ と $[(a+b)/2, b]$ の各部分でシンプソン則を計算して合計 $S_2$ を得る
- $|S_1 – S_2|$ が十分小さければ（$|S_1 – S_2| < 15\varepsilon$）、$S_2$ を採用する。ここで係数 $15$ はリチャードソン外挿から導かれるもので、$S_2 + (S_2 - S_1)/15$ がより高精度な推定値になることに基づいています
- そうでなければ、各部分区間に対して再帰的に同じ手順を繰り返す
この再帰的な分割は二分木として理解できます。被積分関数が急変する領域では木が深くなり（細かい分割）、滑らかな領域では浅いまま（粗い分割）になります。結果として、計算資源が「本当に必要な場所」に集中的に配分されます。
実用上の注意点
適応型求積法を使う際には、いくつかの注意点があります。
- 特異点の処理: 被積分関数が端点や内部に特異点を持つ場合、適応分割が過度に深くなる可能性があります。この場合、特異点を除去する変数変換を事前に行うのが効果的です
- 振動関数: 高周波で振動する関数では、振動の周期より細かい分割が必要になるため、計算コストが大きくなります。このような場合にはフィロン求積法（Filon quadrature）など、振動を考慮した専用の手法が適しています
- 再帰深度の制限: 無限ループを避けるため、最大再帰深度を設定するのが一般的です
SciPyの scipy.integrate.quad はこの考え方に基づく高精度な適応型求積法を実装しており、内部では QUADPACK ライブラリの Fortran コードが使われています。多くの実用的な問題では、手法の選択に迷ったらまず quad を試すのがよいでしょう。
ここまでの理論的な議論を踏まえて、いよいよ各手法をPythonで実装し、精度を実際に比較してみましょう。
Pythonによる実装
以下のコードでは、台形則・シンプソン則・ガウス=ルジャンドル求積の3手法を実装し、異なる性質を持つテスト関数に対して精度を比較します。特に、関数の滑らかさが収束率にどう影響するかに注目してください。
import numpy as np
import matplotlib.pyplot as plt
# 数値積分の実装
def trapezoidal(f, a, b, n):
"""複合台形則"""
x = np.linspace(a, b, n + 1)
h = (b - a) / n
return h * (0.5 * f(x[0]) + np.sum(f(x[1:-1])) + 0.5 * f(x[-1]))
def simpson(f, a, b, n):
"""複合シンプソン則（nは偶数）"""
if n % 2 != 0:
n += 1
x = np.linspace(a, b, n + 1)
h = (b - a) / n
return h / 3 * (f(x[0]) + 4 * np.sum(f(x[1::2])) +
2 * np.sum(f(x[2:-1:2])) + f(x[-1]))
def gauss_legendre(f, a, b, n_points):
"""ガウス=ルジャンドル求積"""
# numpy.polynomial.legendreから求積点と重みを取得
nodes, weights = np.polynomial.legendre.leggauss(n_points)
# [-1, 1] → [a, b] に変換
x = 0.5 * (b - a) * nodes + 0.5 * (b + a)
w = 0.5 * (b - a) * weights
return np.sum(w * f(x))
# テスト関数
def f1(x):
return np.exp(-x**2) # ガウス関数
def f2(x):
return 1 / (1 + 25 * x**2) # ルンゲ関数
def f3(x):
return np.sqrt(np.abs(x - 0.5)) # 微分不連続点あり
# 厳密値（高精度数値積分）
from scipy.integrate import quad
exact1 = quad(f1, 0, 1)[0]
exact2 = quad(f2, -1, 1)[0]
exact3 = quad(f3, 0, 1)[0]
fig, axes = plt.subplots(2, 2, figsize=(14, 11))
# --- (a) 収束の比較（滑らかな関数） ---
ax = axes[0, 0]
ns = np.array([2, 4, 8, 16, 32, 64, 128, 256, 512])
errors_trap = [abs(trapezoidal(f1, 0, 1, n) - exact1) for n in ns]
errors_simp = [abs(simpson(f1, 0, 1, n) - exact1) for n in ns if n >= 2]
ns_gauss = np.array([1, 2, 3, 4, 5, 6, 8, 10, 15, 20])
errors_gauss = [abs(gauss_legendre(f1, 0, 1, n) - exact1) for n in ns_gauss]
ax.loglog(ns, errors_trap, "bo-", linewidth=2, markersize=6, label="Trapezoidal")
ax.loglog(ns, errors_simp, "rs-", linewidth=2, markersize=6, label="Simpson")
ax.loglog(ns_gauss, errors_gauss, "g^-", linewidth=2, markersize=6, label="Gauss-Legendre")
# 理論的収束率
ax.loglog(ns, 0.3 * ns**(-2.0), "b:", linewidth=1, alpha=0.5, label="$O(n^{-2})$")
ax.loglog(ns, 0.1 * ns**(-4.0), "r:", linewidth=1, alpha=0.5, label="$O(n^{-4})$")
ax.set_xlabel("Number of points $n$", fontsize=12)
ax.set_ylabel("Absolute error", fontsize=12)
ax.set_title("(a) Convergence: $\\int_0^1 e^{-x^2}\\,dx$", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")
ax.set_ylim(1e-16, 1)
# --- (b) ルンゲ関数での比較 ---
ax = axes[0, 1]
errors_trap2 = [abs(trapezoidal(f2, -1, 1, n) - exact2) for n in ns]
errors_simp2 = [abs(simpson(f2, -1, 1, n) - exact2) for n in ns]
errors_gauss2 = [abs(gauss_legendre(f2, -1, 1, n) - exact2) for n in ns_gauss]
ax.loglog(ns, errors_trap2, "bo-", linewidth=2, markersize=6, label="Trapezoidal")
ax.loglog(ns, errors_simp2, "rs-", linewidth=2, markersize=6, label="Simpson")
ax.loglog(ns_gauss, errors_gauss2, "g^-", linewidth=2, markersize=6, label="Gauss-Legendre")
ax.set_xlabel("Number of points $n$", fontsize=12)
ax.set_ylabel("Absolute error", fontsize=12)
ax.set_title("(b) Convergence: $\\int_{-1}^1 \\frac{1}{1+25x^2}\\,dx$", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")
ax.set_ylim(1e-16, 1)
# --- (c) 被積分関数の可視化 ---
ax = axes[1, 0]
x = np.linspace(0, 1, 500)
x2 = np.linspace(-1, 1, 500)
ax.plot(x, f1(x), "b-", linewidth=2, label="$e^{-x^2}$ (smooth)")
ax.plot(x2, f2(x2), "r-", linewidth=2, label="$1/(1+25x^2)$ (peaked)")
ax.plot(x, f3(x), "g-", linewidth=2, label="$\\sqrt{|x-0.5|}$ (non-smooth)")
ax.set_xlabel("$x$", fontsize=12)
ax.set_ylabel("$f(x)$", fontsize=12)
ax.set_title("(c) Test Functions", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# --- (d) ガウス求積の求積点 ---
ax = axes[1, 1]
for n_pts in [2, 3, 5, 7, 10]:
nodes, weights = np.polynomial.legendre.leggauss(n_pts)
ax.scatter(nodes, [n_pts] * n_pts, s=weights * 200, c="blue", alpha=0.6)
ax.scatter(nodes, [n_pts] * n_pts, s=5, c="black")
# 等間隔点との比較
for n_pts in [2, 3, 5, 7, 10]:
eq_nodes = np.linspace(-1, 1, n_pts)
ax.scatter(eq_nodes, [n_pts + 0.3] * n_pts, s=30, c="red", alpha=0.6,
marker="x")
ax.set_xlabel("Node position $x_i$", fontsize=12)
ax.set_ylabel("Number of points", fontsize=12)
ax.set_title("(d) Gauss-Legendre Nodes (blue) vs\nEqui-spaced (red x)", fontsize=12)
ax.set_xlim(-1.2, 1.2)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("numerical_integration.png", dpi=150, bbox_inches="tight")
plt.show()
この4つの図から、数値積分手法の特性が明確に読み取れます。
-
左上（滑らかな関数での収束）: 台形則は $O(n^{-2})$、シンプソン則は $O(n^{-4})$ の理論的収束率に沿って誤差が減少しています。ガウス=ルジャンドル求積はわずか5-6点で機械精度に到達しており、滑らかな関数に対する圧倒的な効率がわかります
-
右上（ルンゲ関数での収束）: 尖った関数に対してもガウス求積は優れた収束を示しています。ただし台形則・シンプソン則の収束率はやや悪化しており、関数の滑らかさが精度に直接影響することが確認できます
-
左下（テスト関数）: 3つの関数の形状の違いが見えます。ガウス関数は非常に滑らか、ルンゲ関数は中央に鋭いピーク、$\sqrt{|x-0.5|}$ は $x = 0.5$ で微分不連続です。このような性質の違いが数値積分の精度に大きく影響します
-
右下（求積点の配置）: ガウス=ルジャンドル求積点（青丸、サイズが重みに比例）は区間の端に集中する傾向があり、等間隔点（赤×）とは明らかに異なります。端点付近に求積点を密に配置することで、ルンゲ現象を回避しつつ多項式近似の精度を最大化しています。この最適な配置がガウス求積の高い精度の源であり、チェビシェフ節点の分布と類似した構造を持っていることも興味深い点です
以上の結果は、手法選択の実践的な指針を与えてくれます。滑らかな関数で高精度が必要な場合はガウス求積を、等間隔データや中程度の精度で十分な場合はシンプソン則を、実装の簡潔さを優先する場合は台形則を選ぶのが合理的です。
手法選択のガイドライン
ここまでの理論と実験を踏まえて、実際の問題に対して数値積分手法をどう選ぶかの指針をまとめておきます。
| 状況 | 推奨手法 | 理由 |
|---|---|---|
| 等間隔の離散データ | 台形則 or シンプソン則 | 求積点を自由に選べないため |
| 滑らかな関数、高精度が必要 | ガウス=ルジャンドル求積 | 少ない点数で機械精度に到達 |
| 特異点や急変のある関数 | 適応型求積法（scipy.integrate.quad） |
自動的に問題領域に分割を集中 |
| 重み関数付きの積分 | 対応するガウス求積（ラゲール、エルミート等） | 重み関数を直交多項式で吸収 |
| 高次元の積分 | モンテカルロ法 | 次元の呪いを回避（本記事の範囲外） |
特に注意すべきは、関数の滑らかさが手法選択において最も重要な要因だということです。どんなに高度な手法を使っても、被積分関数に微分不連続点や特異点がある場合は、理論的な収束率が達成されません。そのような場合は、まず変数変換や区間分割で問題の性質を改善してから数値積分を適用するアプローチが有効です。
まとめ
本記事では、主要な数値積分手法の理論と精度を体系的に解説しました。
- 台形則は1次関数近似に基づき精度 $O(h^2)$、シンプルだがオイラー=マクローリン公式による誤差の構造が知られており、リチャードソン外挿の基盤となる
- シンプソン則は2次関数近似に基づくが対称性により精度 $O(h^4)$ を達成する。台形則のリチャードソン外挿としても導出できる
- ガウス=ルジャンドル求積は求積点の配置も最適化し、$n+1$ 点で $2n+1$ 次の多項式を厳密に積分できる。滑らかな関数に対して少ない点数で機械精度に到達する圧倒的な効率を持つ
- ガウス求積の求積点はルジャンドル多項式の零点であり、直交多項式の理論に基づいている。すべての重みが正であることが数値的安定性を保証する
- 適応型求積法は局所的な誤差推定に基づいて分割を自動調整し、被積分関数が一様でない場合に計算資源を効率的に配分する
数値積分は「積分を数値的に計算する」という単純な動機から出発しますが、その背後には近似理論、直交多項式、誤差解析といった豊かな数学的構造が存在します。特にガウス求積法は、直交多項式の零点が最適な求積点になるという美しい定理に支えられており、数値解析と純粋数学の接点を示す好例といえます。
次のステップとして、以下の記事も参考にしてください。
- 補間と近似の理論 — 数値積分の基礎
- 常微分方程式の数値解法 — 数値積分の応用
- 級数解法とべき級数展開 — 特殊関数の求積