数値積分

0.はじめに

学校の課題です。

1.概説

2.説明

数値積分という言葉自体は微分方程式の解法や原関数を求める手法とされることが多いが、本レポートでは数値積分は「定積分の値を数値的に求める求積法」のことである。数値積分はまず、乱数の使用の有無で分類されるが、本手法は乱数を使用しない手法である。また、今回の手法は「台形公式」と「シンプソンの公式」の二つだが、これらは数値積分において「『ニュートン・コーツの公式』を用いた複数の手法の内のひとつである」という位置づけである。

2.1 ニュートン・コーツの公式

ニュートン・コーツの公式とは被積分関数に対して等間隔の点の値を用いて、定積分の値を求める手法の総称である。名前は互いに協力し研究し、本公式などの積分法等などに業績を遺したアイザック・ニュートンとロジャー・コーツに由来する。なお、この方式以外の手法での求積法はガウス求積やクレンショー・カーティス法などが存在する。

2.2台形公式

台形公式は1次関数を用いて、定積分を近似的に求める手法である。一定の間隔に積分範囲を分割し、その各範囲を台形として計算し近似する。この場合の計算方法は以下のようになる。

$ \int_{a}^{b}f(x)dx = \int_{a_0}^{a_1}f(x)dx + \int_{a_1}^{a_2}f(x)dx ・・・ \int_{a_n-1}^{a_n}f(x)dx ≒ \sum_{k=1}^{n}(a_k - a_{k-1})\frac{f(a_{k-1})+f(a_k)}{2}

https://upload.wikimedia.org/wikipedia/commons/thumb/4/42/Composite_trapezoidal_rule_illustration.png/330px-Composite_trapezoidal_rule_illustration.png

図1.台形公式のイメージ

2.3シンプソンの公式

シンプソンの公式は2次関数を用いて、定積分を近似的に求める手法である。台形公式よりも、より近い形に近似するので精度が良い。名称はイギリスの数学者のトーマス・シンプソンからだが、本手法自体はシンプソン以前より知られていた。この場合の計算方法は以下のようになる。

3.精度

台形公式、シンプソンの公式、どちらも分割量を大きくすれば一定以上の精度は出るが、その分の計算量が増えてしまう。下記に示すサンプルプログラムでは100の分割で実用には問題ないと思われる程度の精度は出ている。以下にそれぞれの分割数と近似した値の結果の相関グラフを示す。(横軸を分割数、縦軸を出力値とする)

https://gyazo.com/b52b2e7d56cfbb0e814b9001a14510b5

図.2 台形公式(赤)とシンプソンの公式(青)の近似への収束図(全体)

https://gyazo.com/af7a5ef559108ed8738a862fa07d0548

図.3 台形公式(赤)とシンプソンの公式(青)の近似への収束図(収束地点)

4.参考文献

森北出版 数値計算法 第2版 新装版

5.サンプルプログラム

以下に台形公式、シンプソンの公式の順でサンプルプログラムを添付する。

code:daikei.py

import math

def function(x):

y = 4*(x**3)-10*(x**2)+4*x+5

return y

x_1 = int(input("Please input Range min.\n"))

x_2 = int(input("Please input Range max.\n"))

split = int(input("Please input Split num.\n"))

split = (x_2 - x_1)/split

result = function(x_1)/2 + function(x_2)/2

while x_1 <= x_2:

x_1 += split

result += function(x_1)

result *= split

print(result)

code:simpson.py

import math

def function(x):

y = 4*(x**3)-10*(x**2)+4*x+5

return y

x_1 = int(input("Please input Range min.\n"))

x_2 = int(input("Please input Range max.\n"))

split = int(input("Please input Split num.\n"))

split = (x_2 - x_1)/split

result = function(x_1) + function(x_2)

flag = -1

multi=4

while x_1 <= x_2:

x_1 += split

result += function(x_1)*multi

multi += 2*flag

flag *= -1

result *= (split/3)

print(result)