数値積分

まずは簡単な区間が一定の場合で考えてみよう。 各区間の長さを $h$ とする。

すると、

• 最初の台形の面積 $S_1$ は、 $$ S_1 = \frac{\left(f(x_0) + f(x_1)\right) h}{2} $$

• 次の台形の面積 $S_2$ は、 $$ S_2 = \frac{\left(f(x_1) + f(x_2)\right) h}{2} $$

• …

• 最後の台形の面積 $S_n$ は、 $$ S_n = \frac{\left(f(x_{n-1}) + f(x_n)\right) h}{2} $$

となる。

あとはこれらを足せば良いだけである。 すなわち、求める積分の近似値 $S'$ は

$$\begin{align*} S' & = \frac{\left(f(x_0) + f(x_1)\right) h}{2} + \frac{\left(f(x_1) + f(x_2)\right) h}{2} + \cdots + \frac{\left(f(x_{n-1}) + f(x_n)\right) h}{2} \end{align*}$$

となる。 これをプログラムとして書けばよい。

と言いたいところだが、ちょっとだけ立ち止まってよく考えてみよう。

まず一般論として、

このことを踏まえて、上の式を見直してみよう。

1. $S_1$ と $S_2$ を良く見て欲しい。 どちらにも、$\frac{f(x_1) h}{2}$ が含まれている。 半分にしているものを二つ足しているのだから、そもそも 2 で割る必要がない。

2. 同様に、$f(x_1)$ を二回計算しているが、これも一回で済ませられる。

3. また、全ての $S_i \; (i = 1, 2, ..., n)$ には $h$ が掛けられているが、毎回掛け算せずとも最後に一回だけ掛ければ良い。

一般的には $f(x)$ の計算には非常に時間がかかる (時間がかからないような場合は、そもそも数値積分するまでもなく計算できることがほとんど)。 例えば $f(x) = \mathrm{e}^{-x^2}$ だったとすると、$x^2$ の計算に時間がかかり、$\exp$ の計算でさらに時間がかかる。 従って $f(x)$ の計算回数もできるだけ減らすべきである。

結局のところ、上の式は次のように書きかえることができる。 \[ S' = \left( \frac{f(x_0)}{2} + f(x_1) + f(x_2) + ... + f(x_{n-1}) + \frac{f(x_n)}{2} \right) h \] あるいは、さらにこだわれば、 \[ S' = \left( \frac{f(x_0) + f(x_n)}{2} + f(x_1) + f(x_2) + ... + f(x_{n-1}) \right) h \] となる。

掛け算の回数が $n$ 回から 1 回 (!) になっており、割り算の回数も $n$ 回から 2 ないし 1 回になっている。 また、最初の式では $f(x)$ は $2n$ 回出てきているが、後の式では $n + 1$ 回でほぼ半分になっている。 これらは数学的には全く同一だが、コンピューターでの計算と考えると異なる。

(積分に限らず) ある計算をさせるプログラムを書く手順は、概ね、

• 計算方法の原理を考える

• その為に必要な計算式を導出する

• その式を検討し、簡略化等の余地がないか考える

• プログラムを書く

となる。

このことから分かるのは、プログラムを書くよりも、 その前段階 (準備) の方が多く、重要である ということである。 プログラムを書く前の準備を適当にやってしまうと、とてつもなく遅いプログラムになったり、誤差が大きくなったりするおそれがある。

もちろん特に最初のうちはとりあえず手を動かしてみるというのも悪くない。 ただしその場合も「あとで検討しなおすことが必要で、その際プログラムが大きく変更になる可能性が高い」ことを認識しておくことが重要である。

最も素朴な方法は、最初に $h$ を求めておいて、それを $x_0$ ($= x_\mathrm{l}$) に足していく方法である。

x = xl
h = (xu - xl) / n

for i in range(...):
    ...
    x += h
    ...

しかし $h$ の値は (コンピューター内部では) 誤差を含みうるので、この方法では誤差が蓄積していく可能性がある (数値計算で注意すべきこと: 誤差入門 も参照せよ)。

もう一つの方法は、

$$x_i = x_\mathrm{l} + \frac{\left(x_\mathrm{u} - x_\mathrm{l}\right) \times i}{n}$$ と書けることから、これを毎回計算する方法である。 こちらは $h$ の誤差に起因する誤差を抑えることができるが、掛け算・割り算の回数が増えてしまう。

これらの方法を誤差の観点から比較してみよう。 いずれの方法でも数学的には同一の値が得られるはずなのに、なぜ違いが出る (場合がある) のか疑問に思うかもしれない。 一言で言うと、「コンピューターは無限桁を扱えない」ことに起因する。

分かりやすくするために十進数で考えてみよう。 区間 [0, 1] を 3分割し、その間の点を求めたいとする。 ここでは有効数字 3 桁としよう。

第一の方法では、$h = (1 - 0) / 3 = 1 / 3 = 0.333$ となる。 従って $x_i$ は次のようになる: $$\begin{align*} x_0 &= 0 \\ x_1 &= x_0 + h = 0.333 \\ x_2 &= x_1 + h = 0.666 \\ x_3 &= x_2 + h = 0.999 \end{align*}$$ $x_3 = 1$ であるはずなのだが、$h$ の誤差に起因した誤差が生じている。

第二の方法では、 $$\begin{align*} x_0 &= 0 \\ x_1 &= x_0 + \frac{(1 - 0) \times 1}{3} = 0.333 \\ x_2 &= x_0 + \frac{(1 - 0) \times 2}{3} = 0.667 \\ x_3 &= x_0 + \frac{(1 - 0) \times 3}{3} = 1 \end{align*}$$ となり、第一の方法より誤差が抑えられていることが分かる。 (第二の方法でも、$x_1$, $x_2$ は厳密な値ではないことに注意しよう。実際にはこれらは無限小数である。)

ただし、第一の方法で必ず誤差が発生する訳ではない。 $h$ が誤差を含まなければ (例えば十進数の場合であれば $h = 0.2$ などの場合)、それによる誤差は発生しない。