.. toctree:: ===================== 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軸との交点を求める。 .. math:: f(x) = 2x^2 - 20x - 15 .. literalinclude:: fsolve.py :caption: example-01 :linenos: .. image:: fsolve.png :width: 400px :align: center 2つの関数の交点を求める。 .. math:: &x^2 + y^2 = 2 &y = x^2 - 4x + 1 .. literalinclude:: fsolve2.py :caption: example-02 :linenos: .. image:: fsolve2.png :width: 400px :align: center curve_fit (leastsq) ------------------- 最小自乗法でパラメータ :math:`a, b, c` の最適値を決める。 .. math:: f(x) = a\exp(bx) + cx^2 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html .. literalinclude:: curve_fit.py :caption: example-03 .. image:: curve_fit.png :width: 400px :align: center パラメータの境界制限を指定してフィッテングを行う。 .. literalinclude:: curve_fit_const.py :caption: example-04 .. image:: curve_fit_const.png :width: 400px :align: center minimize -------- パラメータ同士の関係またはパラメータの範囲に拘束条件を課して最適値を求める。 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html .. math:: &-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\\ .. literalinclude:: minimize.py :caption: example-05 数値積分 ======== 以下の関数について、数値解をscipy.integrate.quad()を使って定積分を行います。 https://docs.scipy.org/doc/scipy/reference/integrate.html .. math:: f(x) = \int_0^2 (2x^2 + 3x - 1)dx .. literalinclude:: quad1.py :caption: example-06 積分範囲に :math:`\infty` を含む定積分も可能です。関数い追加の引数を与 えたい場合は、quad()にargsで渡します。gauss積分をscipyで計算して解析解 と比較してみます。 .. math:: \int_{-\infty}^\infty \exp(-a(x-b)^2)dx .. literalinclude:: quad2.py :caption: example-07 .. note:: 記号を定義して代数計算(因数分解や微分積分、微分方程式を解くことまで!) 行う方法は、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 .. math:: &\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. ルジャンドル(Legendre)陪多項式 ------------------------------ 球面調和関数 :math:`Y^l_m (m=0)` の結果に等しいルジャンドルの多項式で す。 https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.legendre.html 直行多項式であるため量子力学でたまに出てきます。以下の微分方程式を満たす :math:`P_n(x)` です。 .. math:: 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 :math:`x=\cos(\theta)` で :math:`P_n(x)` を計算してみます。 .. literalinclude:: legendre.py :caption: example-08 :linenos: 水素の波動関数の関数系である :math:`x = \cos(\theta)` として、lとmのパターンをいくつか計算してプロットしてみます。 .. image:: legendre.png :width: 400px :align: center ラゲール(Laguerre)多項式 ------------------------ 水素原子の波動関数の動径方向の減衰を決める関数です。 https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.laguerre.html .. image:: laguerre.png :width: 400px :align: center 球面調和関数(Spherical harmonics) --------------------------------- 球の表面形状を記述するための直行多項式です。軌道の形を決めために使われま す。degree(方位量指数lに相当)とorder(磁気量子数mに相当)、極座標の :math:`\theta` と :math:`\phi` を与えると、複素数として :math:`Y^m_n` が得られます。 :math:`P^m_n(\cos(\theta)` はルジャンドルの多項式です。 .. math:: 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。 :math:`Y^m_n` の :math:`m` * degree: 磁気量子数l。:math:`Y^m_n` の :math:`n` * theta: 極座標の偏角。一般に記号として :math:`\phi [0, 2\pi]` を使うので注意 * phi: 極座標の偏角。一般に記号として :math:`\theta [0, \pi]` を使うので注意 m, n, theta, phiを与えると球の表面までの距離が得られます。matplotlibで3次元プロットしてみます。実空間では実数のみ扱うので、realパートだけプロットします。軌道の形(:math:`\theta` と :math:`\phi` 依存性)は、球面調和関数で表現されます。 直行座標で可視化するためには、 :math:`|Y^m_n|` の先端がなぞる面を表示します。 :math:`Y^m_n` の正負は、色を変えて表現することとします。 .. literalinclude:: sph_harm.py :caption: example-09 :linenos: プログラムを実行するとp, d, f軌道に相当する球面調和関数の図を作成します。画像をならべた結果は以下のとおりです。 .. image:: sph_harm.png :width: 800px :align: center 水素の波動関数 -------------- :math:`n, l, m` を与えて :math:`|\psi^2|` を計算します。 matplotlibは、:math:`|\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は位相が変わるだけで軌道の形は同じです。 .. literalinclude:: hydrogen.py :caption: example-10 :linenos: プログラムを実行して生成されるcubefileです。vestaで開いて確認してみてください。 * :download:`nml650.cube (n,l,m) = (6,5,0) ` * :download:`nml651.cube (n,l,m) = (6,5,1) ` * :download:`nml652.cube (n,l,m) = (6,5,2) ` * :download:`nml653.cube (n,l,m) = (6,5,3) ` * :download:`nml654.cube (n,l,m) = (6,5,4) ` * :download:`nml655.cube (n,l,m) = (6,5,5) ` vestaで可視化して作成した図をならべた結果です。 .. image:: hydrogen.png :width: 800px :align: center cubeファイルは少なくとも原子座標を1つ含まなければvestaで可視化できない ようです。そのためH原子を1つだけ中心に置いています。上の図ではShow modelを offにしてH原子を表示していません。 vesta上でboundaryを指定すれば、任意の面でスライスすることもできます。(n,l,m)=(6,5,5)のデータをzmin=0.5からzmax=0.5でスライスした場合です。 .. image:: nml655-zslice.png :width: 300px :align: center クイズ ====== Q1 -- 実験データ(:download:`xy_fit.dat `)を2つのガウス関数でピーク分 離することを考える。 .. math:: y = p_0\exp\left(-\frac{x-p_1}{p_2}\right) + p_3\exp\left(-\frac{x-p_4}{p_5}\right) ただし2つのガウス関数の面積は等しいとして以下の拘束条件を満たすように :math:`p` を求めます。ガウス関数の面積はquad()で求めます。 .. math:: \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で拘束条件付きでフィッテングします。 結果は以下のとおりです。 .. image:: xy_fit.png :width: 400px :align: center .. code-block:: none 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 .. note:: quadで積分するのはコストがかかります。面積で規格化された下式 のガウス関数を使えば、わざわざ数値積分する必要がないので、効率的にフィッ テングできるでしょう。 .. math:: &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 Q2 -- 量子数が :math:`(n, l, m) = (3, 2, 0)` の場合、原子核から3.0 Ang.以内の距離に電子が存在する確率を計算する。微小体積の確率を 求めて足し合わせで求めてください。3.0 Ang.以内の領域を分割する分割数 :math:`N` を10, 50, 100, 200と増やして2桁の精度が確認できれば十分です。 .. code-block:: none N=10 7.77% N=50 8.85% N=100 8.88% N=200 8.89% 答え 8.9%