11. scipyの基本と応用

scipyは科学技術計算を行うライブラリで、numpyで作成したデータに基づいて 様々な数値解析を高速かつ容易に行うことができます。なお、numpy.linalgで 利用できる関数は、同じものがscipy.linalgでも利用できます。scipy.linalg で呼び出す関数の方が、numpy.linalgより高速になっている場合があります。 scipyをimportした場合の線形代数計算は、scipy.linalgを使う方が良いです。

numpy同様、scipyは巨大なライブラリです。信号解析からデータ圧縮まで非常 に多くの機能を有しますが、化学と関係のありそうなものだけピックアップし て紹介します。

https://docs.scipy.org/doc/scipy/reference/index.html

今回の学習目標は以下のとおりです。

  • scipyを使ったパラメータの最適化法を理解する。

  • scipyで数値積分を行う。

  • scipyで電子軌道と関係する特殊関数を呼び出してみる。

最適化とパラメータフィッテング

https://docs.scipy.org/doc/scipy/reference/optimize.html

fsolve

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

x軸との交点を求める。

\[f(x) = 2x^2 - 20x - 15\]
example-01
 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy.optimize import fsolve
 4
 5x = np.linspace(-10, 10, 1000)
 6def func(x):
 7    f = 3 * x**2 - 4 * x - 80
 8    return f
 9
10v1 = fsolve(func, 0)
11v2 = fsolve(func, 2.5)
12fig, ax = plt.subplots()
13ax.plot(x, func(x))
14ax.plot(x, np.zeros(x.shape))
15ax.plot(v1, func(v1), "o")
16ax.plot(v2, func(v2), "o")
17fig.savefig("fsolve.png")
18plt.show()
../_images/fsolve.png

2つの関数の交点を求める。

\[ \begin{align}\begin{aligned}&x^2 + y^2 = 2\\&y = x^2 - 4x + 1\end{aligned}\end{align} \]
example-02
 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy.optimize import fsolve
 4
 5x = np.linspace(-10, 10, 1000)
 6def func(p):
 7    x, y = p[0], p[1]
 8    f = np.zeros(2)
 9    f[0] = x**2 + y**2 - 2
10    f[1] = x**2 - 4*x + 1 - y
11    return f
12
13o1 = fsolve(func, [1, 1])
14o2 = fsolve(func, [-1, -1])
15# x**2 + y**2 = 2
16theta = np.linspace(0, 2*np.pi, 1000)
17r = np.sqrt(2)
18x = np.sqrt(2) * np.cos(theta)
19y = np.sqrt(2) * np.sin(theta)
20# x**2 - 4*x2 + 1
21x2 = np.linspace(-2, 2, 1000)
22y2 = x2**2 - 4*x2 + 1
23fig, ax = plt.subplots()
24ax.plot(x, y)
25ax.plot(x2, y2)
26ax.plot(o1[0], o1[1], "o")
27ax.plot(o2[0], o2[1], "o")
28ax.set_xlim(-2, 2)
29ax.set_ylim(-2, 2)
30fig.savefig("fsolve2.png")
31plt.show()
../_images/fsolve2.png

curve_fit (leastsq)

最小自乗法でパラメータ \(a, b, c\) の最適値を決める。

\[f(x) = a\exp(bx) + cx^2\]

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

example-03
 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy.optimize import curve_fit
 4
 5x = np.linspace(0, 5, 20)
 6y = 3*np.exp(-2*x) + 0.11*x**2
 7y = y + (np.random.random(20) - 1)*0.5
 8fig, ax = plt.subplots()
 9
10func = lambda x, a, b, c: a*np.exp(b*x) + c*x**2
11
12p0 = [0, 0, 0] # 初期値
13p, cov = curve_fit(func, x, y, p0)
14print(p)   # [ 2.81148394 -2.62913599  0.09389399]
15print(cov) # [[  3.15666320e-02  -2.55900954e-02   7.83101093e-06] 
16           #  [ -2.55900954e-02   1.03468885e-01  -4.74614449e-05] 
17           #  [  7.83101093e-06  -4.74614449e-05   1.25025906e-05]]
18
19ax.plot(x, y, "ro", label="exp")
20ax.plot(x, func(x, p[0], p[1], p[2]), "r-", label="fit")
21fig.savefig("curve_fit.png")
22plt.show()
../_images/curve_fit.png

パラメータの境界制限を指定してフィッテングを行う。

example-04
 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy.optimize import curve_fit
 4
 5x = np.linspace(0, 50, 100)
 6y = 10*np.exp(-0.2*x) - 1*np.exp(-0.05*x) 
 7y = y + (np.random.random(x.size) - 0.5)*0.1
 8
 9f = lambda x, a, b: a*np.exp(b*x)
10func = lambda x, a, b, c, d: f(x, a, b) + f(x, c, d)
11
12# 境界制限なし
13p0 = [0, 0, 0, 0]
14p, cov = curve_fit(func, x, y, p0)
15print(p)
16fig, ax = plt.subplots(ncols=1, nrows=2)
17ax[0].plot(x, y, "ro", label="exp")
18ax[0].plot(x, func(x, *p), "r-", label="fit") # *p -> p[0], p[1], p[2], p[3]
19ax[0].plot(x, f(x, p[0], p[1]), "-", label="exp1")
20ax[0].plot(x, f(x, p[2], p[3]), "-", label="exp2")
21ax[0].legend()
22ax[0].set_title("non constraint")
23
24# 境界制限付き
25b = [(-np.inf, -np.inf, -np.inf, -np.inf), (np.inf, np.inf, 100, np.inf)]
26p, cov = curve_fit(func, x, y, p0, bounds=b)
27print(p)
28ax[1].set_title("constraint")
29ax[1].plot(x, y, "ro", label="exp")
30ax[1].plot(x, func(x, *p), "r-", label="fit") # *p -> p[0], p[1], p[2], p[3]
31ax[1].plot(x, f(x, p[0], p[1]), "-", label="exp1")
32ax[1].plot(x, f(x, p[2], p[3]), "-", label="exp2")
33ax[1].legend()
34fig.savefig("curve_fit_const.png")
35plt.show()
../_images/curve_fit_const.png

minimize

パラメータ同士の関係またはパラメータの範囲に拘束条件を課して最適値を求める。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

\[\begin{split}&-3x^3 - x^2y + \frac{1}{y} = 3z - 20\\ &x + y \gt z\\ &x^2 + y^2 = -20\\ &-10 \le x \le 30\\ &-10 \le y \le 30\\ &-10 \le z \le 30\\\end{split}\]
example-05
 1import numpy as np
 2from scipy.optimize import minimize
 3
 4
 5def func(p):
 6    x, y, z = p
 7    res = -3 * x**3 - x**2 * y  - 1 / (x*y) - 3 * z + 20
 8    return res
 9
10
11def cond_ineq(p):
12    x, y, z = p
13    return x + y - z
14
15
16def cond_eq(p):
17    x, y, z = p
18    return x**2 + y**2 + 20
19
20bounds = ((-10, 30), (-10, 30), (-10, 30))
21conds = ({'type': 'ineq', 'fun': cond_ineq},
22         {'type': 'eq', 'fun': cond_eq})
23
24p0 = [5, 5, 5]  # 初期値
25o = minimize(func, p0, constraints=conds, bounds=bounds)
26print("func", func(o.x))
27print("cond_ineq", cond_ineq(o.x))
28print("cond_eq", cond_eq(o.x))

数値積分

以下の関数について、数値解をscipy.integrate.quad()を使って定積分を行います。

https://docs.scipy.org/doc/scipy/reference/integrate.html

\[f(x) = \int_0^2 (2x^2 + 3x - 1)dx\]
example-06
1from scipy.integrate import quad
2
3def func(x):
4    return 2*x**2 + 3*x - 1
5
6v, err = quad(func, 0, 1)
7print(v)

積分範囲に \(\infty\) を含む定積分も可能です。関数い追加の引数を与 えたい場合は、quad()にargsで渡します。gauss積分をscipyで計算して解析解 と比較してみます。

\[\int_{-\infty}^\infty \exp(-a(x-b)^2)dx\]
example-07
 1from scipy.integrate import quad
 2import numpy as np
 3
 4
 5def func(x, a, b):
 6    return np.exp(-a * (x-b)**2)
 7
 8
 9a, b = 3, 2
10v, err = quad(func, -np.infty, np.infty, args=(a, b))
11print("scipy integral: {:.6f}".format(v))
12print("gauss integral: {:.6f}".format(2 * np.sqrt(np.pi / a) / 2))

注釈

記号を定義して代数計算(因数分解や微分積分、微分方程式を解くことまで!) 行う方法は、Sympy(scipy.Symbolとは別パッケージ)でまとめて扱います。 Sympyが使えれば大学1年くらいまでの数学演習は全部pythonで解けるかも?

特殊関数

面倒な特殊関数もscipyを使えば簡単に計算してくれます。

https://docs.scipy.org/doc/scipy/reference/special.html

特殊関数で解析解が得られる水素原子の波動関数を計算するために、ルジャン ドル(Legendre)陪多項式、ラゲール(Laguerre)多項式、球面調和関数 (Spherical harmonics)を計算してみます。

https://ja.wikipedia.org/wiki/%E6%B0%B4%E7%B4%A0%E5%8E%9F%E5%AD%90

\[\begin{split}&\Psi(r, \theta, \phi) = \sqrt{ \left(\frac{2}{na_0}\right)^3 \frac{(n-l-1)!}{2n(n+l)!} } \exp\left(-\frac{\rho}{2}\right) \rho^l L^{2l+1}_{n-l-1}(\rho) Y_{lm}(\theta, \phi)\\ &\rho = \frac{2r}{na_0}\\ &Y_{lm}(\theta, \phi) = \left\{ \begin{array}{ll} \frac{i}{\sqrt{2}}\left(Y_l^m - (-1)^mY_l^{-m} \right) & ({\rm if}\ m<0) \\ Y_l^m & ({\rm if}\ m=0)\\ \frac{1}{\sqrt{2}}\left(Y_l^{-m} - (-1)^mY_l^{m} \right) & ({\rm if}\ m>0) \\ \end{array} \right.\end{split}\]

ルジャンドル(Legendre)陪多項式

球面調和関数 \(Y^l_m (m=0)\) の結果に等しいルジャンドルの多項式で す。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.legendre.html

直行多項式であるため量子力学でたまに出てきます。以下の微分方程式を満たす \(P_n(x)\) です。

\[ \begin{align}\begin{aligned}P_v^m = (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}P_v(x)\\P_v = \sum_{k=0}^{\infty}\frac{(-v)_k(v+1)_k}{(k!)^2}\left(\frac{1-x}{2}\right)^k\end{aligned}\end{align} \]

\(x=\cos(\theta)\)\(P_n(x)\) を計算してみます。

example-08
 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy.special import lpmv
 4
 5# scipy.special.lpmv(m, v, x)
 6# order m(磁気量子数; m), degree v(方位量子数; l)
 7theta = np.linspace(0, 2*np.pi, 1000)
 8x = np.cos(theta)
 9fig, ax = plt.subplots()
10ax.plot(theta * 180/np.pi, lpmv(0, 1, x), label="$l=1, m=0$")
11ax.plot(theta * 180/np.pi, lpmv(1, 1, x), label="$l=1, m=1$")
12ax.plot(theta * 180/np.pi, lpmv(0, 2, x), label="$l=2, m=0$")
13ax.plot(theta * 180/np.pi, lpmv(1, 2, x), label="$l=2, m=1$")
14ax.plot(theta * 180/np.pi, lpmv(2, 2, x), label="$l=2, m=2$")
15ax.legend()
16ax.set_xlabel("angle/deg.")
17ax.set_ylabel("$P^m_l$")
18fig.savefig("legendre.png")
19plt.show()
20

水素の波動関数の関数系である \(x = \cos(\theta)\) として、lとmのパターンをいくつか計算してプロットしてみます。

../_images/legendre.png

ラゲール(Laguerre)多項式

水素原子の波動関数の動径方向の減衰を決める関数です。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.laguerre.html

../_images/laguerre.png

球面調和関数(Spherical harmonics)

球の表面形状を記述するための直行多項式です。軌道の形を決めために使われま す。degree(方位量指数lに相当)とorder(磁気量子数mに相当)、極座標の \(\theta\)\(\phi\) を与えると、複素数として \(Y^m_n\) が得られます。 \(P^m_n(\cos(\theta)\) はルジャンドルの多項式です。

\[Y^m_n(\theta, \phi) = \sqrt{ \frac{2n+1}{4\pi} \frac{(n-m)!}{(n+m)!} } \exp(i\theta) P_n^m(\cos(\phi))\]

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html

scipy.special.sph_harm(degree, order, theta, phi)を呼び出せば計算してくれます。引き数は以下のとおり。

  • order: 方位量子数m。 \(Y^m_n\)\(m\)

  • degree: 磁気量子数l。\(Y^m_n\)\(n\)

  • theta: 極座標の偏角。一般に記号として \(\phi [0, 2\pi]\) を使うので注意

  • phi: 極座標の偏角。一般に記号として \(\theta [0, \pi]\) を使うので注意

    m, n, theta, phiを与えると球の表面までの距離が得られます。matplotlibで3次元プロットしてみます。実空間では実数のみ扱うので、realパートだけプロットします。軌道の形(\(\theta\)\(\phi\) 依存性)は、球面調和関数で表現されます。

直行座標で可視化するためには、 \(|Y^m_n|\) の先端がなぞる面を表示します。 \(Y^m_n\) の正負は、色を変えて表現することとします。

example-09
 1import numpy as np
 2import matplotlib.pyplot as plt
 3from scipy.special import sph_harm
 4from mpl_toolkits.mplot3d import Axes3D
 5
 6def PlotData(Ylm, theta, phi, lm, outfile=None):
 7    """
 8    データのプロット
 9    """
10    l, m = lm[0], lm[1]
11    # |Ylm|がなぞるベクトルの先端を計算
12    r = np.abs(Ylm.real)
13    x = r * np.sin(theta) * np.cos(phi)
14    y = r * np.sin(theta) * np.sin(phi)
15    z = r * np.cos(theta)
16    color = np.ones(r.shape)  # 正負の色つけのため   
17    color[Ylm < 0] = -1       # 正負の色つけのため
18    fig = plt.figure()
19    ax = fig.gca(projection='3d')
20    ax.plot_surface(x, y, z, facecolors=plt.cm.jet(color),
21                    rstride=1, cstride=1)
22    plt.title("$Y^{{m={}}}_{{l={}}}$".format(m, l), size=20)
23    length = 0.5
24    ax.set_xlim(-length, length)
25    ax.set_ylim(-length, length)
26    ax.set_zlim(-length, length)
27    ax.set_xlabel('x')
28    ax.set_ylabel('y')
29    ax.set_zlabel('z')
30    plt.setp(ax.get_xticklabels(), visible=False)
31    plt.setp(ax.get_yticklabels(), visible=False)
32    plt.setp(ax.get_zticklabels(), visible=False)
33    if outfile is not None:
34        fig.savefig(outfile)
35        print(outfile, "was created.")
36
37def CalcYlm(lm):
38    l, m = lm[0], lm[1]
39    # order=m degree=lで球面調和関数を計算
40    t = np.linspace(0, np.pi, 50)
41    p = np.linspace(0, 2 * np.pi, 100)
42    theta, phi = np.meshgrid(t, p)
43    if m == 0:
44        Y = sph_harm(m, l, phi, theta)
45    if m > 0:
46        Ya = sph_harm(-m, l, phi, theta)
47        Yb = sph_harm(m, l, phi, theta)
48        Y = 1/np.sqrt(2) * (Ya + (-1)**m*Yb)
49    if m < 0:
50        Ya = sph_harm(m, l, phi, theta)
51        Yb = sph_harm(-m, l, phi, theta)
52        Y = 1j/np.sqrt(2) * (Ya - (-1)**m*Yb)
53    return Y, theta, phi
54
55# (l, m)の組み合わせを好きなだけ
56lmset = [(1, -1), (1, 0), (1, 1),
57         (2, -2), (2, -1), (2, 0), (2, 1), (2, 2),
58         (3, -3), (3, -2), (3, -1), (3, 0), (3, 1), (3, 2), (3, 3)]
59
60for lm in lmset:
61    Y, theta, phi = CalcYlm(lm)
62    # lとmから出力ファイル名を決める
63    outfile = "l{}_m{:+d}.png".format(lm[0], lm[1])
64    PlotData(Y, theta, phi, lm, outfile)
65# plt.show()

プログラムを実行するとp, d, f軌道に相当する球面調和関数の図を作成します。画像をならべた結果は以下のとおりです。

../_images/sph_harm.png

水素の波動関数

\(n, l, m\) を与えて \(|\psi^2|\) を計算します。 matplotlibは、\(|\psi^2(x, y, z)|\) のような3次元の密度データをプロットできないので、可視化はgaussian cubeファイルを作成してvestaで行うこととします。cubeファイルは、小さな体積に空間を分割(voxel)して各voxelに密度データを割り当てるグリッドデータです。空間(voxel)を定義するためのベクトルと分割数を必要とします。

あまり教科書にでてこないg軌道(n=6, l=5, m=0,1,2,3,4,5)のcubeファイルを作成 してみます。m=-1,-2,-3,-4,-5は位相が変わるだけで軌道の形は同じです。

example-10
 1import numpy as np
 2from scipy.special import sph_harm
 3from scipy.special import factorial
 4from scipy.special import assoc_laguerre
 5
 6a0 = 0.529177  # global, Distance: Bohr unit a0 = 0.529177 Ang.
 7
 8def OutputCube(nlm, L, psi2, outfile):
 9    """
10    gaussian cube fileで書き出し
11    a, b, c vectorはBohr単位で指定する
12    """
13    d = (L[1] - L[0]) / L[2] / a0 # Bohr単位
14    p = (L[1] - L[0])/a0 * 0.5
15    o = open(outfile, "w")
16    o.write("hydrogen orbital cube file\n")
17    o.write("(n, l, m) = {} {} {}".format(nlm[0], nlm[2], nlm[2]))
18    o.write("\n")
19    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(1, 0, 0, 0)) # 原子数、origin
20    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(L[2], d, 0, 0)) # a vector bohr
21    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(L[2], 0, d, 0)) # b vector bohr
22    o.write("{:6d} {:8f} {:8f} {:8f}\n".format(L[2], 0, 0, d)) # c vector bohr
23    o.write("{:6d} {:8f} {:8f} {:8f} {:8f}\n".format(1, 0, p, p, p)) # Hの座標
24    psi2 = psi2.reshape(-1) * 100
25    for i, rho in enumerate(psi2):
26        o.write("{:12g} ".format(rho))
27        if i % 6 == 5:
28            o.write("\n")
29    print(outfile, "was created.")
30
31def CalcPsi(nlm, L):
32    n, l, m = nlm[0], nlm[1], nlm[2]
33    # voxcelデータの生成
34    xi = np.linspace(L[0], L[1], L[2])
35    yi = np.linspace(L[0], L[1], L[2])
36    zi = np.linspace(L[0], L[1], L[2])
37    x, y, z = np.meshgrid(xi, yi, zi)
38    # 直行座標 -> 極座標
39    r = np.sqrt(x**2 + y**2 + z**2)
40    theta = np.arccos(z / r)
41    phi = np.arctan2(y, x)
42    # 係数
43    rho = 2*r/a0/n
44    coef1 = np.sqrt((2 / (n * a0))**3 * factorial(n - l - 1) /
45                 (2 * n * factorial(n + l)))
46    coef2 = np.exp(-rho/2) * rho**l
47    # ラゲール(Laguerre)
48    L = assoc_laguerre(rho, n - l - 1, 2 * l + 1)
49    # 球面調和関数(spherical harmonic function)
50    if m == 0:
51        Y = sph_harm(m, l, phi, theta)
52    if m > 0:
53        Ya = sph_harm(-m, l, phi, theta)
54        Yb = sph_harm(m, l, phi, theta)
55        Y = 1/np.sqrt(2) * (Ya + (-1)**m*Yb)
56    if m < 0:
57        Ya = sph_harm(m, l, phi, theta)
58        Yb = sph_harm(-m, l, phi, theta)
59        Y = 1j/np.sqrt(2) * (Ya - (-1)**m*Yb)
60    psi =coef1 * coef2 * L * Y
61    psi2 = psi.real * np.conj(psi.real)
62    # 色つけのための符号
63    s = np.ones(psi.shape)
64    s[Y < 0] = -1
65    psi2 = psi2 * s
66    return psi2
67
68# 主量子数, 方位量子数, 磁気量子数
69nlmset = [(6, 5, 0), (6, 5, 1), (6, 5, 2), (6, 5, 3),  (6, 5, 4), (6, 5, 5)]
70nlmset = [(6, 5, -5), (6, 5, -4), (6, 5, -3), (6, 5, -2),  (6, 5, -1)]
71# xyz軸の min(Ang.), max(Ang.), voxcelの分割数
72L = (-30, 30, 100)
73for nlm in nlmset:
74    outfile = "nml{}{}{}.cube".format(nlm[0], nlm[1], nlm[2])
75    psi2 = CalcPsi(nlm, L)
76    OutputCube(nlm, L, psi2, outfile)

プログラムを実行して生成されるcubefileです。vestaで開いて確認してみてください。

vestaで可視化して作成した図をならべた結果です。

../_images/hydrogen.png

cubeファイルは少なくとも原子座標を1つ含まなければvestaで可視化できない ようです。そのためH原子を1つだけ中心に置いています。上の図ではShow modelを offにしてH原子を表示していません。

vesta上でboundaryを指定すれば、任意の面でスライスすることもできます。(n,l,m)=(6,5,5)のデータをzmin=0.5からzmax=0.5でスライスした場合です。

../_images/nml655-zslice.png

クイズ

Q1

実験データ(xy_fit.dat)を2つのガウス関数でピーク分 離することを考える。

\[y = p_0\exp\left(-\frac{x-p_1}{p_2}\right) + p_3\exp\left(-\frac{x-p_4}{p_5}\right)\]

ただし2つのガウス関数の面積は等しいとして以下の拘束条件を満たすように \(p\) を求めます。ガウス関数の面積はquad()で求めます。

\[\int_{-\infty}^{\infty} p_0\exp\left(-\frac{x-p_1}{p_2}\right) = \int_{-\infty}^{\infty} p_3\exp\left(-\frac{x-p_4}{p_5}\right)\]

minimizeは初期値が最適値から遠いとうまく最適値を求めてくれません。そこ で、最初にleast_squareで拘束条件なしでパラメータを求め、求めたパラメー タを初期値としてminimizeで拘束条件付きでフィッテングします。

結果は以下のとおりです。

../_images/xy_fit.png
least_squares fitting:
p[0] = 2.01
p[1] = -2.00
p[2] = 1.00
p[3] = 1.00
p[4] = 0.99
p[5] = 3.00
area ratio gauss1/gauss2: 1.16
minimize fitting:
p[0] = 1.97
p[1] = -2.03
p[2] = 0.92
p[3] = 0.99
p[4] = 0.90
p[5] = 3.63
area ratio gauss1/gauss2: 1.00

注釈

quadで積分するのはコストがかかります。面積で規格化された下式 のガウス関数を使えば、わざわざ数値積分する必要がないので、効率的にフィッ テングできるでしょう。

\[\begin{split}&f(x) = p_0\sqrt{\frac{\log(2)}{\pi}}\frac{1}{p_2}\exp\left[ -\log(2)\left\{\frac{x - p_1}{p_2} \right\}^2 \right]\\ &\int_{-\infty}^{\infty} f(x) dx = p_0\end{split}\]

Q2

量子数が \((n, l, m) = (3, 2, 0)\) の場合、原子核から3.0 Ang.以内の距離に電子が存在する確率を計算する。微小体積の確率を 求めて足し合わせで求めてください。3.0 Ang.以内の領域を分割する分割数 \(N\) を10, 50, 100, 200と増やして2桁の精度が確認できれば十分です。

N=10 7.77%
N=50 8.85%
N=100 8.88%
N=200 8.89%

答え 8.9%