p1_ex03

数値積分

モジュールのインポート

import numpy as np
from scipy import integrate # 数値積分
import matplotlib.pyplot as plt

import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])

1重積分

1変数の被積分関数は，

def func3(x): # 1変数の被積分関数
    y = np.sqrt(4.0-x**2)
    return y

関数 integrate

準備として，まず変数および被積分関数のサンプル点での値を配列に入れる．

a = 0.0
b = 2.0
nx = 151
xx = np.linspace(a,b,nx)
ff = func3(xx)

数値積分は，
https://docs.scipy.org/doc/scipy/reference/integrate.html
を用いて，台形公式 (trapezoidal rule) の場合，

y4 = integrate.trapz(ff, xx) # 台形公式
print("integrate.trapz: ",y4,y4-np.pi) # 積分値と誤差

integrate.trapz:  3.1409526607810467 -0.0006399928087463813

また，シンプソンの公式 (Simpson's rule) の場合，

y5 = integrate.simps(ff, xx) # シンプソンの公式
print("integrate.simps: ",y5,y5-np.pi) # 積分値と誤差

integrate.simps:  3.1413425976445577 -0.00025005594523541674

Numpyの関数
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad
では，Fortran のライブラリ QUADPACK を使ったガウスの数値積分 (Gaussian quadrature) があり，

y3 = integrate.quad(func3, a, b) # 1次元積分(適応積分)
print("integrate.quad :",y3) # 積分値と推定誤差

integrate.quad : (3.1415926535897922, 3.533564552071766e-10)

プロット

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,ff)
plt.axis('square')
fig.savefig('p1_ex03_1.pdf')

2重積分

2変数の被積分関数は，

def func33(y,x):# 2変数の被積分関数
    f = np.sqrt(x**2+y**2)
    return f

2重積分を求める関数は，
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html

x_min, x_max = 0.0, 1.0
y33 = integrate.dblquad(func33, x_min, x_max, 0.0, lambda x: np.sqrt(1-x**2)) # 2重積分
print("integrate.dblquad:",y33) # 積分値と推定誤差
print("sol:",1.0/6.0*np.pi)

integrate.dblquad: (0.5235987755981053, 1.4792290955080084e-08)
sol: 0.5235987755982988