[Python에 의한 과학·기술 계산] 수치 적분, 사다리꼴 법칙·심슨 법칙, 수치 계산, scipy

scipy.integrate의 cumtrapz 메소드(사다리꼴 법칙)와 simps 메소드(심슨 법칙)를 이용하여 이산 데이터의 수치 적분을 실시한다.
예를 들어, $\int_0^1\frac{4}{1+x^2} dx =\pi $를 고려한다.

내용

(1.A) scipy를 이용한 코드. 서둘러 있을 때는 이거.
(2) 사다리꼴 법칙과 심슨 법칙의 계산 정밀도에 관한 부록.
(3) 사다리꼴 법칙과 심슨 법칙의 파이썬 코드에 의한 간단한 구현. 계산 절차를 알고 싶을 때 유용합니다.

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

그림의 직선 기울기는 수렴 차수의 $ n $에 해당합니다. 예상대로 사다리꼴 법칙(청색 선)은 $n\sim-2$, 심슨 법칙(적선)은 $n\sim-3$이다.

(3) 사다리꼴 법칙과 심슨 법칙의 파이썬 코드에 의한 간단한 구현.

사용법 : 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()