読者です 読者をやめる 読者になる 読者になる

計算物理屋の研究備忘録

Linux, Cygwin, Mac, zsh, emacs等の使い方、設定などのメモ

Scipyとmatplotlibで球ベッセル関数、球ノイマン関数をプロット

Python

目次

必要なもの

  • python
  • scipy
  • numpy
  • matplotlib

scipyのインストール

Cygwinにscipyをpipでインストールするのは難しいので、素直にWindowsネイティブのpythonを使う。Anacondaでインストールすればよい。Anacondaならデフォルトでscipy, numpy, matplotlibが入っているはず。minicondaの場合はcondaコマンドでインストール。

python スクリプト

ipython notebookを使った。1行目はnotebookでinline表示するためのコマンド。

%matplotlib inline
import numpy as np
from scipy.special import sph_jn, sph_yn
import matplotlib.pyplot as plt

xmin = 0
xmax = 20
lmax = 2
x = np.linspace(xmin, xmax)
jn = []
yn = []
for i in x:
    if jn == []:
        jn = sph_jn(lmax, i)[0]
    else:
        jn = np.vstack([jn, sph_jn(lmax, i)[0]])
# j0 --> jn[:, 0]
# j1 --> jn[:, 1]
# j2 --> jn[:, 2]
    if yn == []:
        yn = sph_yn(lmax, i)[0]
    else:
        yn = np.vstack([yn, sph_yn(lmax, i)[0]])
# y0 --> yn[:, 0]
# y1 --> yn[:, 1]
# y2 --> yn[:, 2]

#----- plot
plt.figure(figsize=(10, 4))

plt.subplot(121)
plt.hlines(0, xmin, xmax, linestyle='dashed')
plt.ylim(-1.2, 1.2)
plt.plot(x, jn[:, 0], label='$j_0$')
plt.plot(x, jn[:, 1], label='$j_1$')
plt.plot(x, jn[:, 2], label='$j_2$')
plt.xlabel('$x$', fontsize=18)
plt.legend(loc='best', fontsize=14)

plt.subplot(122)
plt.hlines(0, xmin, xmax, linestyle='dashed')
plt.ylim(-1.2, 1.2)
plt.plot(x, yn[:, 0], label='$n_0$')
plt.plot(x, yn[:, 1], label='$n_1$')
plt.plot(x, yn[:, 2], label='$n_2$')
plt.xlabel('$x$', fontsize=18)
plt.legend(loc='best', fontsize=14)

plt.savefig("spherical_Bessel.png", bbox_inches="tight")

プログラムメモ

from scipy.special import sph_jn, sph_yn

球ベッセル関数(sph_jn)と球ノイマン関数(sph_yn)をインポート
マニュアルにはこう書いてある

sph_jn(n, z) Compute the spherical Bessel function jn(z) and its derivative for all orders up to and including n.
sph_yn(n, z) Compute the spherical Bessel function yn(z) and its derivative for all orders up to and including n.

例えば、sph_jn(2, 0)はこんな結果になる

In [1]: sph_jn(2, 0)
Out[1]: (array([ 1.,  0.,  0.]), array([ 0.        ,  0.33333333,  0.        ]))

sph_jn(2, 0)[0]が[j0(0), j1(0), j2(0)]で、sph_jn(2, 0)[1]がその微分を返すというちょっと変な関数になっている。ここでは微分は必要ない。

x = np.linspace(xmin, xmax)

x軸に使用するarrayを生成。

Signature: np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)

np.linspace(0, 20)なら0から20の区間で(20も含む)50点サンプルをとる。np.linspace(0, 20, 100)と第3引数を指定すれば100点になる。

jn = np.vstack([jn, sph_jn(lmax, i)[0]])

上に書いたように、sph_jn(2, 0)[0]が[j0(0), j1(0), j2(0)]というリストになるので、for文でxを変えながら縦にarrayをスタックしている。
- j0 --> jn[:, 0]
- j1 --> jn[:, 1]
- j2 --> jn[:, 2]
- y0 --> yn[:, 0]
- y1 --> yn[:, 1]
- y2 --> yn[:, 2]

plt.savefig("spherical_Bessel.png", bbox_inches="tight")

bbox_inches="tight"と書くことで図の余白を小さくできる。

グラフ

f:id:keisanbutsuriya:20151028144723p:plain