scipy.integrateのcumtrapzメソッド(台形則)およびsimpsメソッド(シンプソン則)を利用して離散データの数値積分を行う。
例として，$\int_0^1 \frac{4}{1+x^2} dx = \pi $を考える。
内容
(1.A)scipyを利用したコード。急いでいるときはコレ。
(2) 台形則とシンプソン則の計算精度に関する補遺。
(3) 台形則とシンプソン則のpythonコードによる簡易な実装。計算手続きを知りたい場合に有用。
(1) scipyを利用したコード
（１）scipy利用コード。簡潔。
from scipy import integrate
import numpy as np
x=np.linspace(0,1,5) # [0,1]を5等分したグリッドをxに格納
y=4/(1+x**2) # 数値積分の被積分関数
y_integrate_trape = integrate.cumtrapz(y, x) #台形則による数値積分計算
y_integrate_simps = integrate.simps(y, x) #シンプソン則による数値積分計算
print(y_integrate_trape[-1]) #結果の表示
print(y_integrate_simps) # 結果の表示
結果(1):scipy利用
3.13117647059 #台形則
3.14156862745 #シンプソン則
3.14159265358979... 厳密値(π)
(2)[補遺] 台形則とシンプソン則の精度・収束速度の比較
よく知られているように， 十分小さな刻み幅hに対して，
台形則による積分誤差は$O(h^2)$，シンプソン則のそれは$O(h^3)$となる。
グリッド数をNとすると，これらはそれぞれ$O(N^{-2})$，$O(N^{-3})$となる。
これを本例題を通じて直接検証することは教育的である。
以下にそれを確かめるためのコードを記す。
数値積分による相対誤差をy軸に， グリッド数Nを横軸にとり，両対数プロットしている。
上述の関係が成り立つ領域で$log y \propto n logN$となり， $n=-2$が台形則， $n=-3$がシンプソン則に相当する。
(2)計算精度の比較
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
err_trape=[]
err_simps=[]
nn=[4,8,16,32,64,128,256,512,1024,2048] # 4から2048までのグリッド数をリストnnに格納
for j in nn:
x=np.linspace(0,1, j)
y=4/(1+x**2)
y_integrate_trape = integrate.cumtrapz(y, x) #台形則による数値積分計算
y_integrate_simps = integrate.simps(y, x) #シンプソン則による数値積分計算
err_trape.append(abs(np.pi-y_integrate_trape[-1])/np.pi) # 台形則による数値積分の相対誤差評価
err_simps.append(abs(np.pi-y_integrate_simps)/np.pi) # シンプソン則による数値積分の相対誤差評価
# for plot
ax = plt.gca()
ax.set_yscale('log') # y軸をlogスケールで描く
ax.set_xscale('log') # x軸をlogスケールで描く
plt.plot(nn,err_trape,"-", color='blue', label='Trapezoid rule')
plt.plot(nn,err_simps,"-", color='red', label='Simpson rule')
plt.xlabel('Number of grids',fontsize=18)
plt.ylabel('Error (%)',fontsize=18)
plt.grid(which="both") # グリッド表示。"both"はxy軸両方にグリッドを描く。
plt.legend(loc='upper right')
plt.show()
図の直線の傾きが収束次数の$n$に相当する。予想されたとおり，台形則(青線)は$n\sim-2$, シンプソン則(赤線)は$n\sim-3$となっている。
(3) 台形則とシンプソン則のpythonコードによる簡易な実装。
利用法: XYデータファイル('xy.dat')を用意して，
python3 fuga.py xy.dat
とすると台形則・シンプソン則による数値計算結果を表示する。
fuga.py
import os, sys, math
def integrate_trape(f_inp): #台形則による数値積分の関数
x_lis=[]; y_lis=[]
# f_inp.readline()
for line in f_inp:
x_lis.append(float(line.split()[0]))
y_lis.append(float(line.split()[1]))
sum_part_ylis=sum(y_lis[1:-2])
del_x=(x_lis[1]-x_lis[0])
integrated_value = 0.5*del_x*(y_lis[0]+y_lis[-1]+2*sum_part_ylis)
print("solution_use_trapezoid_rule: ", integrated_value)
def integrate_simpson(f_inp): #シンプソン則による数値積分の関数
x_lis=[]; y_lis=[]
n_counter = 0
for line in f_inp:
x_lis.append(float(line.split()[0]))
if n_counter % 2 == 0:
y_lis.append(2*float(line.split()[1]))
else:
y_lis.append(4*float(line.split()[1]))
n_counter += 1
sum_part_ylis=sum(y_lis[1:-2]) # sum from y_1 to y_n-1
del_x=(x_lis[1]-x_lis[0])
integrated_value = del_x*(y_lis[0]+y_lis[-1]+sum_part_ylis)/3
print("solution_use_Simpson_formula: ", integrated_value)
##
def main():
fname=sys.argv[1]
f_inp=open(fname,"r")
integrate_trape(f_inp)
f_inp.close()
f_inp=open(fname,"r")
integrate_simpson(f_inp)
f_inp.close()
if __name__ == '__main__':
main()