台形法とシンプソン法を用いて以下の数値積分を行う。

\int_0^2 {x^5}dx

厳密解

\begin{eqnarray}
\int_0^2 {x^5}dx &=& [\frac{1}{6} x^6]^2_0 \\
&=& 10.6666666666... 
\end{eqnarray}

台形法

被積関数f(x)を積分区間でN個に分割し、i+1番目の区間[x_i, x_i+1]を台形で近似し足し合わせていく方法。

各台形の面積は(上底+下底)×高さ÷2より

S_i=\frac{(f(x_{i+1})+f(x_{i}))×(x_{i+1}-x_{i})}{2}

区間x=[a,b]での和は

\begin{eqnarray}
S &=& \int_a^b {f(x)}dx \\
&=& \sum_{i=0}^{N-1} \int_{x_{i}}^{x_{i+1}} {f(x)}dx \\
&≃& \sum_{i=0}^{N-1}S_i \\
&=& \frac{1}{2}\sum_{i=0}^{N-1}(f(x_{i+1})+f(x_{i}))(x_{i+1}-x_{i}) \\
&=& \frac{1}{2}\Biggl(f(x_0)(x_1-x_0)+\sum_{i=1}^{N-1}f(x_i)(x_{i+1}-x_{i-1}+f(x_{N})(x_{N}-x_{N-1})\Biggr)
\end{eqnarray}

となる。
N等分の場合は、h = (x_i+1 - x_i) = (b - a)/Nとすると

S = \frac{h}{2}(f(x_0)+2\sum_{i=1}^{N-1}f(x_i)+f(x_N))...①

誤差はO(h^2)

・計算コード

N=10000

def f(x):
    return x**5.0 # x^5

N = 10000

a = 0
b = 2
h = (b-a)/N

S = (h/2) * (f(a) + 2*sum(f(h*i) for i in range(1,N-1)) + f(b)) # ①の計算

print(S)

・結果

>>> print(S)
10.660270132693393

シンプソン法

3点間を2次のラグランジュ補間を行い、各3点間に対して積分を行い和を計算する方法。

3点間x_i、x_i+1、x_i+2の2次のラグランジュ補間関数は以下

\begin{eqnarray}
f^{(i)}(x)&=&\frac{(x-x_{i+1})(x-x_{i+2})}{(x_i-x_{i+1})(x_i-x_{i+2})}f(x_i) + \frac{(x-x_{i})(x-x_{i+2})}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})}f(x_{i+1}) + \frac{(x-x_{i})(x-x_{i+1})}{(x_{i+2}-x_{i})(x_{i+2}-x_{i+1})}f(x_{i+2})
\end{eqnarray}

補間関数を区間[x_i,x_i+2]で積分すると

\begin{eqnarray}
S_i &=& \int_{x_i}^{x_{i+2}}f^{(i)}(x)dx\\
&=&\frac{1}{6}[\frac{x_{i+2}-x_i}{x_i-x_{i+1}}(2x_i+x_{i+2}-3x_{i+1}f(x_i) + \frac{(x_i-x_{i+2}^3)}{(x_{i+1}-x_i)(x_{i+2}-x_{i+1})}f(x_{i+1}) + \frac{x_{i+2}-x_i}{x_{i+2}-x_{i+2}}(2x_{i+2}+x_i-3x_{i+1})f(x_{i+2})]
\end{eqnarray}

N等分の場合は、h = (x_i+1 - x_i) = (b - a)/Nとすると

\begin{eqnarray}
S_i &=& \frac{h}{3}(f(x_i)+4f(x_{i+1})+f(x_{i+2}))
\end{eqnarray}

区間x=[a,b]での和は

\begin{eqnarray}
S &≃& \sum_{i=0}^{\frac{N}{2}-1}S_{2i}
\end{eqnarray}

誤差はO(h^4)

・計算コード

N=10000

def f(x):
    return x**5.0

N = 10000

a = 0
b = 2
h = (b-a)/N

S = (h/3) * sum((f(h*i) + 4*f(h*(i+1)) + f(h*(i+2))) for i in range(0,N-1, 2))

print(S)

・結果

>>> print(S)
10.666666666666648