[Pythonによる科学・技術計算] 数値積分, 台形則・シンプソン則，数値計算, scipy

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()