#はじめに
Pythonの使い方に慣れるため，色々な記事を参考にしながら，Pythonを使った簡単な数値積分の方法をまとめてみました．
今回は 例として，$\int_{a}^{b}x^2dx\quad a=0，b=3\hspace{5pt}$の積分を行います．
#目次
1. SciPy，Sympyを使用する方法
・定積分
・不定積分
2. ライブラリを使用せずに求める方法
・区分求積法
・台形公式
・シンプソン法
#SciPy，Sympyでの数値積分
Scipyでの積分方法の例
| 種類 | |
|---|---|
| quad | 求積法(1変数) |
| dblquad | 2重積分 |
| tplquad | 3重積分 |
| trampz | 台形公式 |
| simps | シンプソンの公式 |
###定積分
Sympyでもできますが，今回はScipyを使ってみます．
SciPyを用いるとintegrate.quadによって積分近似値と推定誤差の2つの値が得られます．
# scipyのintegrateをインポート
from scipy import integrate
# 積分する関数の定義(今回はx^2)
def f(x)
return x**2
# 上限値と下限値の入力
a = int(input("下限"))
b = int(inpput("上限"))
ans,err = integrate.quad(f,a,b)
# 積分近似値の表示
print(ans)
# 推定誤差の表示
print(err)
▼入力
a=0 b=3
▼出力
積分近似値： 9.000000000000002
推定誤差： 9.992007221626411e-14
###不定積分
SciPyでは (多分) できないのでsympyを使用します．
詳しくは"sympyによるシンボリック変数での微分積分"を参考に．
# Sympyのインポート
from sympy import *
# 変数を定義
x = Symbol('x')
# 関数を定義
f = x**2
# 積分
a = integrate(f,x)
# 結果の表示
print(a)
▼出力
x**3/3
#ライブラリを使用せずに求める方法
以下では積分方法の公式をScipy，Sympyなどのライブラリを用いずにpythonで行う方法を記述していますが，Numpyを使用した方が少し綺麗に書けます．(参考：数値積分-台形公式とシンプソンの公式の導出)
###区分求積法
短冊状に区切って長方形の面積を足していって面積を求める方法です．
分割数nが大きい程精度はあがります．
a = int(input("下限"))
b = int(input("上限"))
n = int(input("分割数"))
# 短冊の幅Δx
dx = (b-a)/n
# 積分する関数の定義
def f(x):
return x**2
# 面積の総和
s = 0
for i in range(1,n):
# x,yの値
x1 = a + dx*i
f1 = f(x1)
# 面積
s += dx*f1
# 結果の表示(小数点以下10桁)
print("{:.10f}".format(s))
▼入力
a=0 b=3 n=100
▼出力
8.8654500000
###台形公式
台形公式はある関数をn個に分割し，x=iとi+1の区間を近似し，足し合わせていく方法です．
a = int(input("下限"))
b = int(input("上限"))
n = int(input("分割数"))
# 短冊の幅Δx
dx = (b-a)/n
# 積分する関数の定義
def f(x):
return x**2
# 面積の総和
s = 0
for i in range(1,n-1):
x1 = a + dx*i
x2 = a + dx*(i+1)
f1 = f(x1) # 上底
f2 = f(x2) # 下底
# 面積
s += dx * (f1+f2)/2
# 結果の表示(小数点以下10桁)
print("{:.10f}".format(s))
▼入力
a=0 b=3 n=100
▼出力
8.7331230000
###シンプソン法
シンプソン法は台形公式とは異なり，二次関数に近似し，変化率の大きい関数でも精度よく近似できる方法です．
\int_{a}^{b}f(x)dx≅\frac{b-a}{6}\Biggl(f(a)+4f(\frac{a+b}{2})+f(b)\Biggl)
a = int(input("下限"))
b = int(input("上限"))
m = int(input("分割数"))
n = int(m/2)
# 積分する関数の定義
def f(x):
return x**2
def Simpson(f,a,b,m):
df = (b-a)/(2*n)
h = (b-a)/m
sum = 0
for i in range(1,n+1):
x1 = a + df*(2*i -2)
x2 = a + df*(2*i -1)
x3 = a + df*(2*i)
f1 = f(x1)
f2 = f(x2)
f3 = f(x3)
sum += h*(f1 + 4*f2 + f3)/3
return sum
# 結果の表示
print(Simpson(f,a,b,m))
▼入力
a=0 b=3 m=100
▼出力
8.999999999999998
##参考
Pythonでざっくりと数値解析の基礎をしてみる(3_数値積分編) -Hatena Blog
Pythonで学ぶ計算物理 7.1.数値積分の方法
コンピューター処理 11.scipyの基本と応用