LoginSignup
23
21

More than 5 years have passed since last update.

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

Last updated at Posted at 2017-07-20

scipy.integrateのcumtrapzメソッド(台形則)およびsimpsメソッド(シンプソン則)を利用して離散データの数値積分を行う。
例として,$\int_0^1 \frac{4}{1+x^2} dx = \pi $を考える。

内容

(1.A)scipyを利用したコード。急いでいるときはコレ。
(2) 台形則とシンプソン則の計算精度に関する補遺。
(3) 台形則とシンプソン則のpythonコードによる簡易な実装。計算手続きを知りたい場合に有用。


(1) scipyを利用したコード

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

t.png

図の直線の傾きが収束次数の$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()
23
21
0

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
23
21