問題4.19の解答(台形公式、シンプソンの公式)

概説微分積分の問題4.19の解答。

台形公式は定積分を台形で近似するもの。台形の公式を使って面積を出して整理すればすぐ出てくる。

def trapezoidal(f, range_, n):
    a, b = range_
    result = f(a) + f(b)
    for i in xrange(1, n):
        result += f(a + (b - a) * i / n) * 2
    return result * (b - a) / 2 / n

シンプソンの公式は定積分を放物線で近似するもの。これも愚直に計算して整理すれば出てくる。

def simpsons(f, range_, n):
    a, b = range_
    result = f(a) + f(b)
    for i in xrange(1, n):
        result += f(a + (b - a) * (2 * i) / 2 / n) * 2
    for i in xrange(0, n):
        result += f(a + (b - a) * (2 * i + 1) / 2 / n) * 4
    return result * (b - a) / 6 / n

実際に利用して誤差を見てみる。sympy を使って実際の値を求め、近似値と比較する。

import math
import sympy as sym
from sympy.utilities import lambdify
from sympy.abc import x

# compute the real answer
f = 1 / (x**2 + 1)
F = sym.integrate(f, x)
answer = F.subs(x, 1) - F.subs(x, 0)
print "f: %s, F: %s, answer=%s" % (f, F, answer)

# compute approximations
integral_range = (0., 1.)
n = 2
f_lambda = lambdify(x, f)
evaled_answer = answer.evalf()
trapezoidal_result = trapezoidal(f_lambda, integral_range, n * 2)
trapezoidal_error = abs(trapezoidal_result - evaled_answer)
print "trapezoidal: %f, error=%f" % (trapezoidal_result, trapezoidal_error)

simpsons_result = simpsons(f_lambda, integral_range, n)
simpsons_error = abs(simpsons_result - evaled_answer)
print "simpsons: %f, error=%f" % (simpsons_result, simpsons_error)

f: 1/(x**2 + 1), F: atan(x), answer=pi/4
trapezoidal: 0.782794, error=0.002604
simpsons: 0.785392, error=0.000006

この例のように滑らかな曲線だと、シンプソンの公式の方が有利に働くことが多いが、必ずしもこうなるとは限らない。特に滑らかではない曲線の場合は、台形公式の方が有利らしい。

これらの近似の誤差を評価するには、台形公式では二階微分の絶対値の上限、シンプソンの公式では四階微分の絶対値の上限が必要となる。