12
14

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

Python - 数値積分 台形法 & シンプソン法

Last updated at Posted at 2017-06-10

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

\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
12
14
4

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
12
14

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?