台形公式

この章の目標

台形則と複合台形則の導出、誤差 $O(h^2)$ の評価、リチャードソン外挿による精度向上を理解する。

前提知識

目次

1. 基本台形公式

関数 $f(x)$ の区間 $[a,b]$ 上の定積分を、両端点を結ぶ直線(台形)の面積で近似する。

基本台形公式

$$\displaystyle\int_a^b f(x)\,dx \approx \dfrac{b-a}{2}\bigl[f(a) + f(b)\bigr]$$

これは $f$ を区間 $[a,b]$ 上で1次多項式(直線)で補間し、その下の面積を求めることに対応する。1次の精度を持ち、$f$ が1次関数ならば正確な値を与える。

下のアニメーションでは、区間 $[0,1]$ を $n$ 等分し、各小区間を台形で近似してそれらを合計した値を近似値 $T_n$ と書く(区間をいくつの台形に分けたかが $n$)。$n$ を大きく(分割を細かく)するほど $T_n$ は真の積分値に近づく。この $T_n$ の厳密な式は第2節(複合台形公式)で与える。

\(\displaystyle I=\int_0^1 e^x\,dx=e-1=1.718282\)
0 / 0
図1. 台形公式の幾何学的意味(アニメーション)。曲線 $y = f(x) = e^x$ の下の面積(真の積分 $\displaystyle\int_0^1 e^x\,dx = e-1 \approx 1.71828$)を台形で近似する。「▶ 自動再生」または「次へ/戻る」で1段ずつ進み、分割数 $n$ が $1 \to 2 \to 4 \to 8 \to 16 \to 32$ と倍になるたびに台形が細かくなって推定値 $T_n$ が真値へ収束する。薄赤の領域が真の面積、青い台形が近似であり、ログには誤差が各段でおよそ $1/4$ に減る($O(h^2)$)様子が表示される。

2. 複合台形公式

区間 $[a,b]$ を $n$ 等分し、各小区間で台形公式を適用して合計する。刻み幅 $h = (b-a)/n$、分点 $x_k = a + kh$ とすると、

複合台形公式

$$\displaystyle\int_a^b f(x)\,dx \approx T_n = h\!\left[\dfrac{f(a)}{2} + \displaystyle\sum_{k=1}^{n-1} f(x_k) + \dfrac{f(b)}{2}\right]$$

両端点は重み $1/2$、内部の分点は重み $1$ で合計する。

3. 誤差評価

定理(複合台形公式の誤差)

$f \in C^2[a,b]$ のとき、

$$\displaystyle\int_a^b f(x)\,dx - T_n = -\dfrac{(b-a)h^2}{12}\,f''(\xi), \qquad \xi \in (a,b)$$

すなわち誤差は $O(h^2)$ である。

刻み幅 $h$ を半分にすると誤差はおよそ $1/4$ になる。この性質はリチャードソン外挿ロンバーグ積分の基礎となる(第4節も参照)。

log h log |誤差| 台形 O(h²) Simpson O(h⁴) 傾き 2 傾き 4
図2. 両対数スケールでの刻み幅 $h$ と誤差の関係。誤差 $\sim C h^p$ は直線になり、その傾きが収束次数 $p$ を表す。台形公式は傾き 2($O(h^2)$)、シンプソン公式は傾き 4($O(h^4)$)。$h$ を小さくする(左へ進む)ほど両者とも誤差が減り、傾きの急なシンプソン公式の方が速く減少する。

4. 端点補正

オイラー=マクローリン公式により、複合台形公式の誤差を $h$ のべき級数で展開できる:

$$T_n - I = \displaystyle\sum_{k=1}^{p} \dfrac{B_{2k}}{(2k)!}\,h^{2k}\bigl[f^{(2k-1)}(b) - f^{(2k-1)}(a)\bigr] + O(h^{2p+2})$$

ここで $B_{2k}$ はベルヌーイ数である($B_2 = \tfrac{1}{6},\ B_4 = -\tfrac{1}{30},\ \ldots$)。端点での導関数値 $f^{(2k-1)}(a),\,f^{(2k-1)}(b)$ が既知であれば、この展開の補正項を $T_n$ から差し引いて精度を上げられる。これを端点補正(勾配補正)と呼ぶ。

一次の端点補正

展開の $k=1$ の項だけを使う。$\dfrac{B_2}{2!} = \dfrac{1}{12}$ なので、誤差の主要項は $\dfrac{h^2}{12}\bigl[f'(b) - f'(a)\bigr]$ である。これを $T_n$ から引くと、

端点補正した台形公式

$$T_n^{\text{corr}} = T_n - \dfrac{h^2}{12}\bigl[f'(b) - f'(a)\bigr]$$

$h^2$ の項が打ち消され、$T_n^{\text{corr}}$ の誤差は $O(h^4)$ に改善される。端点の1階導関数 $f'(a),\,f'(b)$ さえ分かれば、関数評価を増やさずに精度が2次向上する点が利点である。

二次の端点補正

さらに $k=2$ の項($f^{(3)}$ を用いる)まで補正すれば $O(h^6)$ になる。$\dfrac{B_4}{4!} = \dfrac{-1/30}{24} = -\dfrac{1}{720}$ なので、補正項の符号はになり、

二次まで端点補正した台形公式

$$T_n^{\text{corr2}} = T_n - \dfrac{h^2}{12}\bigl[f'(b) - f'(a)\bigr] + \dfrac{h^4}{720}\bigl[f^{(3)}(b) - f^{(3)}(a)\bigr]$$

と書ける。$h^2$ と $h^4$ の項が両方打ち消され、$T_n^{\text{corr2}}$ の誤差は $O(h^6)$ となる。一般に $k$ を増やすほど高次の端点導関数 $f^{(2k-1)}(a),\,f^{(2k-1)}(b)$ を要するが、各補正で精度が2次ずつ上がっていく。

計算例:補正の効果

$\displaystyle\int_0^1 e^x\,dx$ を例にとる。$f(x)=e^x$ では $f'(x)=f^{(3)}(x)=e^x$ なので、$f'(1)-f'(0) = f^{(3)}(1)-f^{(3)}(0) = e - 1$ となる。図1の $n=4$ の値 $T_4 = 1.727222$(誤差 $8.94\times10^{-3}$)を、$h=\tfrac14$ として補正していく。

一次補正($O(h^4)$):

$$T_4^{\text{corr}} = T_4 - \dfrac{h^2}{12}\,(e-1) = 1.727222 - 0.008949 = 1.718273$$

真値 $e-1 = 1.718282$ との誤差は $9.3\times10^{-6}$ で、補正前の約 $1/1000$ に縮む。

二次補正($O(h^6)$):$f^{(3)}$ の項をさらに加えると、$\dfrac{h^4}{720}(e-1) = 0.0000093$ なので

$$T_4^{\text{corr2}} = T_4 - \dfrac{h^2}{12}(e-1) + \dfrac{h^4}{720}(e-1) = 1.727222 - 0.008949 + 0.0000093 = 1.718282$$

誤差は $1.4\times10^{-8}$ まで下がり、一次補正よりさらに約 $1/700$ に縮む(補正前のおよそ60万分の1)。

近似真値との誤差
$T_4$(補正なし)1.727221905$8.94\times10^{-3}$
$T_4^{\text{corr}}$(一次・$O(h^4)$)1.718272520$9.31\times10^{-6}$
$T_4^{\text{corr2}}$(二次・$O(h^6)$)1.718281842$1.4\times10^{-8}$
真値 $e-1$1.718281828

端点の導関数を1つずつ足すだけで、$n$ を増やさずに精度が桁違いに向上していくことが分かる。

導関数を使わない別ルート

端点補正は $f'(a),\,f'(b)$ を必要とするが、同じオイラー=マクローリン展開($h^2,\,h^4,\,\ldots$ の偶数べきのみ)を使い、導関数を一切使わずに精度を上げる方法もある。刻み幅の異なる複数の $T_n$ を組み合わせて $h^{2k}$ の項を順に消すリチャードソン外挿と、それを表として繰り返し適用するロンバーグ積分である。導関数が手に入らない一般の場合は、こちらが実用的な選択肢となる。

5. 周期関数での超収束

$f$ が周期 $T = b - a$ の滑らかな周期関数である場合、$f^{(k)}(a) = f^{(k)}(b)$ が全ての $k$ で成り立つため、オイラー=マクローリン展開の補正項が全て消える。結果として、台形公式は指数関数的に収束する。

これは数値的に極めて重要であり、周期関数の積分には台形公式が最適な手法の一つである。

例1:三角多項式は厳密に積分できる

$\sin^2 x = \dfrac{1 - \cos 2x}{2}$ は周波数2までの三角多項式である。1周期 $[0,2\pi]$ での真値は

$$\int_0^{2\pi}\sin^2 x\,dx = \pi \approx 3.1415927$$
0 / 0
図3. 被積分関数 $f(x)=\sin^2 x$ を区間 $[0,2\pi]$ で台形公式により積分するアニメーション。薄い赤が真の面積 $\pi$、青い台形が近似 $T_n$。「▶ 自動再生」または「次へ/戻る」で分割数 $n$ を $1,2,3,4,8$ と変えられる。$\sin^2 x$ は周波数2の三角多項式なので $n\ge 3$ で $T_n$ は厳密に $\pi$ に一致する($n=1,2$ は標本がすべて零点に当たり $T_n=0$)。

複合台形公式($h = 2\pi/n$)の値は次のようになる。

$n$$T_n$真値 $\pi$ との誤差
10$3.1\times10^{0}$
20$3.1\times10^{0}$
33.14159270(厳密)
43.14159270(厳密)
83.14159270(厳密)

$n \ge 3$ で厳密に真値 $\pi$ に一致する。1周期上の台形公式は、分割数 $n$ が被積分関数の最高周波数を超えれば(ここでは周波数2なので $n \ge 3$)三角多項式を誤差ゼロで積分するからである。$n = 1, 2$ では周波数2の成分が平均値へ折り返して(エイリアシング)誤差が出る。なお、最も単純な $\displaystyle\int_0^{2\pi}\sin x\,dx = 0$ も同様に、すべての $n$ で厳密に $0$ が得られる。

ただしこの厳密性は区間が完全な1周期のときの特権である。半周期 $\displaystyle\int_0^{\pi}\sin^2 x\,dx = \dfrac{\pi}{2}$ を台形公式で求めると周期性が効かず、通常の $O(h^2)$ 誤差に戻る。

例2:三角多項式でない場合の指数関数的収束

被積分関数が三角多項式でなければ厳密にはならないが、滑らかで周期的なら真値へ指数関数的に収束する。$e^{\cos x}$ はフーリエ級数が無限に続くため有限の $n$ では厳密にならないが、真値は第1種変形ベッセル関数 $I_0$ を用いて次のように書ける。

$$\int_0^{2\pi}e^{\cos x}\,dx = 2\pi I_0(1) \approx 7.9549265$$
0 / 0
図4. 被積分関数 $f(x)=e^{\cos x}$ を区間 $[0,2\pi]$ で台形公式により積分するアニメーション。薄い赤が真の面積 $2\pi I_0(1)\approx 7.955$、青い台形が近似 $T_n$。「▶ 自動再生」または「次へ/戻る」で分割数 $n$ を $1,2,4,8,16$ と倍にできる。三角多項式ではないため厳密にはならないが、滑らかな周期関数ゆえ誤差が指数関数的に減り、$T_n$ が急速に真値へ収束する。

複合台形公式の値 $T_n$ と誤差は次のようになる。

$n$$T_n$真値との誤差
29.6953563$1.7\times10^{0}$
47.9893277$3.4\times10^{-2}$
67.9552091$2.8\times10^{-4}$
87.9549278$1.3\times10^{-6}$
107.9549265$3.5\times10^{-9}$

分割数を $4 \to 8$ にしただけで誤差は $3.4\times10^{-2} \to 1.3\times10^{-6}$ と約2万分の1になる。$O(h^2)$($n$ を倍にして誤差 $1/4$)をはるかに上回り、誤差が $n$ に対して幾何級数的に減る超収束である。被積分関数が滑らかで周期的なほど、オイラー=マクローリン展開の補正項が次々に消え、この高速収束が得られる。

なぜ両端の補正だけで効くのか

「両端の台形だけが特別」なわけではない。各小区間の局所誤差はどれも曲率で決まり、$-\dfrac{h^3}{12}f''(\xi_k)$ で同じオーダーである。特別なのは、足し合わせたときに内部の寄与が打ち消し合うことである。

オイラー=マクローリンの形では、1つの区間が出す誤差は、その区間の両端での導関数の差として書ける。これを全区間で足すと、各内部ノード $x_k$ には左隣の区間からの $+f^{(2k-1)}(x_k)$ と右隣の区間からの $-f^{(2k-1)}(x_k)$ が現れ、符号が逆なので完全に相殺する(テレスコープ)。残るのは一番外側の $-f^{(2k-1)}(a)$ と $+f^{(2k-1)}(b)$ だけで、これが補正項に端点しか現れない理由である。

各区間は両端に ∓f′ を出す ── 内部で相殺、両端だけ残る 区間1 区間2 区間3 区間4 + + + + = 0 = 0 = 0 −f′(a) +f′(b) x₀ = a x₁ x₂ x₃ x₄ = b
図5. 端点補正がなぜ両端だけで効くかの模式図。各小区間(区間 $i$)は、誤差として左端に $-f'(x_{i-1})$、右端に $+f'(x_i)$ を出す。内部のノードでは、左隣の区間が出す $+f'(x_k)$ と右隣の区間が出す $-f'(x_k)$ が打ち消し合って $0$ になる。相殺されずに残るのは最も左端の $-f'(a)$ と右端の $+f'(b)$ だけであり、これが誤差の主要項 $(h^2/12)[f'(b)-f'(a)]$ の正体である。

言葉だけでは半信半疑かもしれないので、和が実際に崩れる様子を書き下してみる。区間 $i$($x_{i-1}$ から $x_i$)の局所誤差 $I_i - T_i = -\dfrac{h^3}{12}f''(\xi_i)$ で、曲率 $f''$ を傾きの変化に読み替えるのが鍵である。$\displaystyle\int_{x_{i-1}}^{x_i} f''\,dx = f'(x_i) - f'(x_{i-1})$、すなわち $h\,f''(\xi_i) \approx f'(x_i) - f'(x_{i-1})$ なので、各区間の誤差はその両端の傾きの差で書ける($I_i$ は区間 $i$ の真の積分、$T_i$ はその台形値)。

$$I_i - T_i \;\approx\; -\frac{h^2}{12}\,\bigl[\,f'(x_i) - f'(x_{i-1})\,\bigr]$$

これを全区間で足すと、隣り合う項が次々に消えるテレスコープ和(望遠鏡和)になる。

$$\begin{aligned} I - T_n \;&=\; -\frac{h^2}{12}\sum_{i=1}^{n}\bigl[\,f'(x_i) - f'(x_{i-1})\,\bigr] \\[2pt] &=\; -\frac{h^2}{12}\Bigl[\,\bigl(f'(x_1)-f'(x_0)\bigr) + \bigl(f'(x_2)-f'(x_1)\bigr) + \cdots + \bigl(f'(x_n)-f'(x_{n-1})\bigr)\,\Bigr] \\[2pt] &=\; -\frac{h^2}{12}\,\bigl[\,f'(x_n) - f'(x_0)\,\bigr] \;=\; -\frac{h^2}{12}\,\bigl[\,f'(b) - f'(a)\,\bigr] \end{aligned}$$

中央のどのノード $x_k$($1 \le k \le n-1$)の傾き $f'(x_k)$ も、隣り合う2つの括弧に $+f'(x_k)$ と $-f'(x_k)$ として1回ずつ現れる。左隣の区間が右端で出す $+f'(x_k)$ と、右隣の区間が左端で出す $-f'(x_k)$ が、符号が逆なのでちょうど打ち消し合うのである。生き残るのは最も外側の $-f'(x_0)=-f'(a)$ と $+f'(x_n)=+f'(b)$ の2つだけ。台形値の誤差で書けば $T_n - I = \dfrac{h^2}{12}\bigl[f'(b)-f'(a)\bigr]$ となり、これが端点補正で差し引く項そのものである ── 「内部は相殺し、両端だけが残る」ことの数式的な正体である。

最も本質的には、台形公式は滑らかな周期関数を1周期で積分すると誤差ゼロであり、有限区間 $[a,b]$ の誤差は、$f$ を周期化したときの両端の継ぎ目の食い違い $f^{(2k-1)}(b) - f^{(2k-1)}(a)$ だけから生じる。内部のノードには継ぎ目がないので誤差を生まない。$f$ が周期関数ならこの食い違いが全次数で消え、第5節の超収束になる。端点補正と超収束は、「誤差は境界量である」という同じ事実の表裏なのである。

6. 端点で特異な場合

複合台形公式 $T_n = h\left[\dfrac{f(a)}{2} + \displaystyle\sum_{k=1}^{n-1} f(x_k) + \dfrac{f(b)}{2}\right]$ は端点の値 $f(a),\,f(b)$ を直接含む。したがって端点で $f$ が定義されない(値が存在しない、または発散する)場合、そのままでは計算できない。これは広義積分(improper integral)の扱いに関わる問題で、状況を3つに分けて整理する。

場合1:端点で極限が有限(除去可能な特異点)

例として $f(x)=\dfrac{\sin x}{x}$ は $x=0$ で未定義だが $\displaystyle\lim_{x\to 0}\frac{\sin x}{x}=1$ である。このような除去可能な特異点では、$f(a)$ を極限値で定義し直せば($f(0):=1$)関数は連続になり、台形公式はそのまま適用できる。$f$ が端点まで2回連続微分可能な状態が回復すれば、誤差評価 $O(h^2)$ も成り立つ。

場合2:端点で発散するが積分は収束する

例:$\displaystyle\int_0^1 \frac{1}{\sqrt{x}}\,dx = 2$ は収束するが $f(0)=+\infty$ であり、$f(a)$ を書けない。これは台形公式が最も苦手とする状況で、次の対処がある。

  • 開いた求積公式を使う:中点則のように端点を評価しない公式に切り替え、特異点を踏まないようにする。
  • 変数変換で特異性を消す:適切な置換で被積分関数を滑らかにする。端点特異性には二重指数関数型公式(DE公式・tanh-sinh 変換、高橋・森)が特に強力で、実務で広く使われる。
  • 特異部分を解析的に分離する:$f(x)=\dfrac{g(x)}{\sqrt{x}}$ のように特異な因子を取り出し、その積分は手計算し、残りの滑らかな部分だけを数値積分する。
  • 重み付き求積を使う:端点の特異性を重み関数として取り込むガウス・ヤコビ求積などを用いる。

例:変数変換で特異性を除く

$\displaystyle\int_0^1 \frac{1}{\sqrt{x}}\,dx$ で $x=t^2$($dx = 2t\,dt$)とおくと、

$$\int_0^1 \frac{1}{\sqrt{x}}\,dx = \int_0^1 \frac{1}{t}\cdot 2t\,dt = \int_0^1 2\,dt = 2$$

被積分関数が定数 $2$(無限回微分可能)に変わり、特異性が消える。変換後であれば、台形公式は刻み幅によらず厳密値 $2$ を返す。

場合3:積分そのものが発散する

例:$\displaystyle\int_0^1 \frac{1}{x}\,dx = \infty$。積分値が存在しないため、どの数値積分法でも有限の正しい値は得られない(刻みを細かくするほど値が増え続ける)。この場合は、数値的に発散を検知する、という使い方になる。

3つの場合のまとめ

ここまでの場合1〜3を一覧にすると次のようになる。

場合端点の状況台形公式対処
1極限が有限(除去可能)値を補えば使える$f(a)$ を極限値で定義
2発散するが積分は収束そのままでは不可中点則/変数変換(DE公式)/特異部分の分離
3積分が発散不可解析的に発散と判断

要点は、台形公式の誤差 $O(h^2)$ が$f \in C^2[a,b]$(端点まで2回連続微分可能)を前提にしていることである。端点で $f$ やその導関数が発散すると、たとえ値を補って計算できても収束は極端に遅くなり、精度の保証は失われる。したがって端点特異性は、変数変換や特異部分の分離で滑らかにしてから台形公式を適用するのが基本方針となる。

7. よくある質問

Q1. 台形公式とは何か

定積分を被積分関数の両端の値を結ぶ台形の面積で近似する数値積分法である。複合台形公式では区間を $n$ 等分し、各小区間で台形近似を行い合計する。

Q2. 台形公式の誤差はどのくらいか

複合台形公式の誤差は $O(h^2)$ である。刻み幅 $h$ を半分にすると誤差はおよそ $1/4$ になる。端点補正を用いればさらに高精度化できる。

Q3. 台形公式が特に高精度になる場合はあるか

周期関数を1周期にわたって積分する場合、台形公式は超収束(指数関数的収束)を示す。これはオイラー=マクローリン展開の補正項が全て消えるためである。

8. 参考資料

  • Wikipedia「台形公式」(日本語版)
  • Wikipedia「Trapezoidal rule」(英語版)
  • L. N. Trefethen & J. A. C. Weideman, "The Exponentially Convergent Trapezoidal Rule," SIAM Review, 56(3), 2014.
  • R. L. Burden & J. D. Faires, Numerical Analysis, 10th ed., Cengage, 2016.