速度分布関数についてsympyで計算し可視化するまで

はじめに

大学院の講義で出た課題で計算した結果をグラフにまとめる必要がありました.
せっかくなので手計算ではなくpythonのライブラリであるsympyで計算させ, 続けてmatplotlibで可視化することにしました.

やること

陽子の速度分布関数をプロットし, 特徴的な物理量をマークする.

sympy(微分・積分の計算)

• 3次元速度分布関数を変数空間全体で積分し1に規格化する.
• 規格化した速度分布関数に速度を掛けて積分することで平均速度$\bar{v}$をだす.
• 同様のやり方で根二乗平均速度$\sqrt{\bar{v^2}}$をだす.
• 微分して速度分布関数がピークをとるときの速度をだす.

matplotlib(可視化)

• 規格化した速度分布関数をプロットする.
• 平均速度, 根二乗平均速度, ピーク時の速度の部分に印をつける.

背景

3次元速度分布関数$f(v)$は以下のように表されます。
$$
f(v) = 4\pi Cv^2\exp(-mv^2/2k_\mathrm{B}T)
$$
ここで$v$は速さ, $m$は質量, $k_\mathrm{B}$はボルツマン定数, $T$は温度, $C$は規格化定数です.
$C$は注目したい物理量によって規格化し, 速度分布関数から密度の分布を得たい場合は,
$$
\int^\infty_0 f(v)\mathrm{d}v=n
$$
となり, 確率を得たい場合は1に規格化します. (確率密度関数になる)
$$
\int^\infty_0f(v)\mathrm{d}v=1
$$
ここでは1に規格化することとすれば平均速度$\bar{v}$, 根二乗平均速度$\sqrt{\bar{v^2}}$, ピークをとる速度$v_\mathrm{p}$は以下の式を解くことで得ることができます.
$$
\bar{v}=\int^\infty_0vf(v)\mathrm{d}v
$$
$$
\sqrt{\bar{v^2}}=\sqrt{\int^\infty_0v^2f(v)\mathrm{d}v}
$$
$$
\left.\frac{\mathrm{d}f}{\mathrm{d}v}\right|_{v=v_p}=0
$$
そこで, 以上の計算を手計算ではなくsympyを使って計算していきます.

やってみる

インポート

matplotlibとnumpyとsympyをインポートします

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp
import sympy as sm

plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.major.width'] = 1.2
plt.rcParams['ytick.major.width'] = 1.2
plt.rcParams['axes.linewidth'] = 1.2
plt.rcParams['axes.grid']=True
plt.rcParams['grid.linestyle']='--'
plt.rcParams['grid.linewidth'] = 0.3
plt.rcParams["font.size"] = 14
plt.rc('text', usetex=True)

定義

速度分布関数をsympyで扱うため, シンボルを定義します.

k = 1.38064852E-23 # ボルツマン定数
m = 1.67262190E-27 # 陽子の質量(kg)
T = 1000000        # 温度(K)
a = m/(2*k*T)      # 簡単化のための係数
    
# 変数の定義
v = sm.Symbol('v')
f = sm.Function('f')
F = sm.Function('F')
C = sm.Symbol('C')
    
# 速度分布関数(三次元)
f = 4*C*sm.pi*v**2*sm.exp(-v**2*a)

この時点で$f$を表示させると以下のようになります.
$$
f(v) = \displaystyle 4 \pi C v^{2} e^{- 6.05737765901491 \cdot 10^{-11} v^{2}}
$$

規格化

次にこの速度分布関数を規格化します.

sm.integrate(f, (v, 0, sm.S.Infinity))

この結果はこのようになります.
$$
\int^\infty_0f(v)\mathrm{d}v=\displaystyle 5.82306221790302 \cdot 10^{-35} \pi^{\frac{3}{2}} C
$$
1に規格化するため, (右辺)=1すなわち, (右辺)-1=0を$v$について解きます.

sm.solve(sm.integrate(f, (v, 0, sm.S.Infinity))-1)[0])

結果はタプルで返ってきます.
[8.46646598960099e-17]
この値を$C$に代入し, $f$を規格化します.

F = f.subs([(C,sm.solve(sm.integrate(f, (v, 0, sm.S.Infinity))-1)[0])])

規格化定数$C$に具体的な数字を代入した$f$は以下のようになります.

$$
F(v)=\displaystyle 3.38658639584039 \cdot 10^{-16} \pi v^{2} e^{- 6.05737765901491 \cdot 10^{-11} v^{2}}
$$
この$F$について色々みていこうと思います.

諸物理量

平均速度

$vF$を積分すると平均速度$\bar{v}$を出すことができます.

vmean = sm.integrate(v*F,(v, 0, sm.S.Infinity))

$$
\bar{v}=\displaystyle 46149.0601591185 \pi=144981.5483659602
$$
100万度下において、陽子の平均速度はおよそ144 km/sだということがわかります.

根二乗平均

$v^2F$を積分し平方根を出すと根二乗平均速度$\sqrt{\bar{v^2}}$を出すことができます.

vsquare = sm.sqrt(sm.integrate(v**2*F,(v, 0, sm.S.Infinity)))

$$
\sqrt{\bar{v^2}}=\displaystyle 66686.9568088195 \pi^{\frac{3}{4}}=157363.24542811327
$$

ピーク時の速度

$F$を$v$について微分し0となる時の$v$を求めます.

F.diff().doit() # 微分を実行

$$
\displaystyle - 4.10276655489749 \cdot 10^{-26} \pi v^{3} e^{- 6.05737765901491 \cdot 10^{-11} v^{2}} + 6.77317279168079 \cdot 10^{-16} \pi v e^{- 6.05737765901491 \cdot 10^{-11} v^{2}}
$$

以下のようにすると, $F$が極値をとる時の$v$が得られます.

sm.solve(F.diff().doit())
vpeak = sm.solve(F.diff().doit())[2]

[-128486.551855745, 0.0, 128486.551855745]
ただし, $v>0$のところしか意味を持たないので, 3番目をピークをとる速度として採用します.

可視化

sympyにもプロットする機能(matplotlibをラッパーしているらしい)がありますが扱いにくいので, matplotlibに直接プロットしたいデータを投げたいと思います.
このまま$F$を投げてもプロットできないのでnumpyで扱える形(=ラムダ式)にしてやります.

Flam = sm.lambdify(v,F,'numpy')

このFlamにはnumpyの配列を入れることができます.
次のようにプロットします

xlim = 400000
ylim = 1E-5
Flam = sm.lambdify(v,F,'numpy')
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111)

x = np.linspace(0,xlim,10000)
ax.plot(x, Flam(x), label='$T$=%d K' % T)

ax.legend()
ax.set_xlim(0,xlim)
ax.set_ylim(0,ylim)
ax.set_xlabel('velocity (m/s)')
ax.set_ylabel('rate ((m/s)${}^{-1}$)')
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_major_formatter(mp.ticker.ScalarFormatter(useMathText=False))
ax.yaxis.set_major_formatter(mp.ticker.ScalarFormatter(useMathText=False))
ax.ticklabel_format(style="sci", axis="x", scilimits=(0,0))
ax.ticklabel_format(style="sci", axis="y", scilimits=(0,0))

ax.vlines(vpeak,0, Flam(float(vpeak)), linestyles='dashed')
ax.vlines(vmean,0, Flam(float(vmean)), linestyles='dashed')
ax.vlines(float(vsquare),0, Flam(float(vsquare)), linestyles='dashed')
ax.text(vpeak,1.2*Flam(float(vsquare)),"$v_\mathrm{p}$")
ax.text(vmean,1.1*Flam(float(vsquare)),"$\overline{v}$")
ax.text(float(vsquare),1*Flam(float(vsquare)),"$\sqrt{\overline{v^2}}$")

fig.savefig('hoge.png')

可視化するとこのようになります.

全コード

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp
import sympy as sm

# プロットの全体設定
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.major.width'] = 1.2
plt.rcParams['ytick.major.width'] = 1.2
plt.rcParams['axes.linewidth'] = 1.2
plt.rcParams['axes.grid']=True
plt.rcParams['grid.linestyle']='--'
plt.rcParams['grid.linewidth'] = 0.3
plt.rcParams["font.size"] = 14
plt.rc('text', usetex=True)

k = 1.38064852E-23 # ボルツマン定数
m = 1.67262190E-27 # 陽子の質量(kg)
T = 1000000        # 温度(K)
a = m/(2*k*T)      # 簡単化のための係数
    
# 変数の定義
v = sm.Symbol('v')
f = sm.Function('f')
F = sm.Function('F')
C = sm.Symbol('C')
    
# 速度分布関数(三次元)
f = 4*C*sm.pi*v**2*sm.exp(-v**2*a)


F = f.subs([(C,sm.solve(sm.integrate(f, (v, 0, sm.S.Infinity))-1)[0])])

vmean = float(sm.integrate(v*F,(v, 0, sm.S.Infinity)))
vsquare = float(sm.sqrt(sm.integrate(v**2*F,(v, 0, sm.S.Infinity))))
vpeak = sm.solve(F.diff().doit())[2]

xlim = 400000
ylim = 1E-5
Flam = sm.lambdify(v,F,'numpy')
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111)

x = np.linspace(0,xlim,10000)
ax.plot(x, Flam(x), label='$T$=%d K' % T)

ax.legend()
ax.set_xlim(0,xlim)
ax.set_ylim(0,ylim)
ax.set_xlabel('velocity (m/s)')
ax.set_ylabel('rate ((m/s)${}^{-1}$)')
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_major_formatter(mp.ticker.ScalarFormatter(useMathText=False))
ax.yaxis.set_major_formatter(mp.ticker.ScalarFormatter(useMathText=False))
ax.ticklabel_format(style="sci", axis="x", scilimits=(0,0))
ax.ticklabel_format(style="sci", axis="y", scilimits=(0,0))

ax.vlines(vpeak,0, Flam(float(vpeak)), linestyles='dashed')
ax.vlines(vmean,0, Flam(float(vmean)), linestyles='dashed')
ax.vlines(float(vsquare),0, Flam(float(vsquare)), linestyles='dashed')
ax.text(vpeak,1.2*Flam(float(vsquare)),"$v_\mathrm{p}$")
ax.text(vmean,1.1*Flam(float(vsquare)),"$\overline{v}$")
ax.text(float(vsquare),1*Flam(float(vsquare)),"$\sqrt{\overline{v^2}}$")

fig.savefig('hoge.png')

最後に

sympyは文字式を計算できるのでとても便利です!