Open1ヶ月前にコメント追加9

統計物理学について

統計物理学

普遍的な物理学に用いられている統計的な考え方を習得して、自然現象のモデル化の一助にしたい。
まずは、統計物理学の文献がある程度読めるように大学レベルの統計力学の基礎知識を習得したい。
そこで「統計力学キャンパス・ゼミ 前編 Kindle版　馬場敬之 (著) 」

https://amzn.asia/d/0WcWWtB

をまずは読んで、習得した知識を整理したい

覚えるべきこと：

・ニュートンの運動方程式とラジグランジュの運動方程式またはハミルトンの正準方程式は等価

・統計力学で重要な役割を担うのはハミルトンの正準方程式

・ハミルトンの正準方程式：
 \frac{dq_j}{dt} = \frac{\partial H}{\partial p_j} , \frac{dp_j}{dt} = - \frac{\partial H}{\partial q_j} 覚えるために考えたこと

きちんと理解してイメージで覚えるために、落ち着いてハミルトンの正準方程式についてよく考えてみる。

以下は（現時点で）式を眺めながら考えたことを記録したもの。

・ニュートンの運動方程式と等価ということは、粒子に対する運動方程式である。

・左の式はある粒子の速度は系全体の量であるハミルトニアンをその粒子の運動量で偏微分したもの

・右の式はある粒子の運動量の時間変化はハミルトニアンをその粒子の位置で偏微分したもの

・ハミルトニアンさえ計算できれば、粒子の速度、運動量の時間変化を計算できる

（そのうち、ハミルトニアンの詳細がわかれば実際に具体例で計算してみるべきだろう）
頭の中で色々イメージをしてみたのを何とか表現しようとした図

 １次元の調和振動子の運動をハミルトンの正準方程式で計算してみる
 正準方程式の計算ハミルトニアンはH = \frac{p^2}{2m}+\frac{1}{2}kq^2で与えられる。

この正準方程式を計算すると

\dot{q}=\frac{\partial H}{\partial p} = p/m = v

\dot{p}=-\frac{\partial H}{\partial q} = -kq \equiv -m\omega^2q

と確かにニュートンの運動方程式となっている。

 位相空間とトラジェクトリを求める(q, p)を座標軸とする空間を位相空間という

時刻tを進展させたときの粒子の位相空間上の軌跡をトラジェクトリという
t=0でq=A, p=0とすると

H = K + \Phi = \frac{p^2}{2m} + \frac{1}{2}m\omega^2q^2 = \frac{1}{2}m\omega^2A^2

これを右辺で割ると楕円の方程式

\frac{q^2}{A^2} + \frac{p^2}{m^2\omega^2 A^2}= 1

が得られる

位相空間のトラジェクトリを可視化したい

とりあえずmatplotlibでアニメーションを作る方法を調べて作成することを目標として調べていくことにする

 matplotlibでトラジェクトリを描けるようになるまでのメモhttps://qiita.com/osanshouo/items/3c66781f41884694838b

でアニメーション作成手順を学習中

 覚えるべきポイント（コード）anim = FuncAnimation(fig, update, frames=range(8), interval=1000)
anim.save("c01.gif", writer="pillow") 
それぞれの引数の型とそれを使って何をやってそうなのかを説明できるようになること

推測だが

１行目はアニメーション作成用オブジェクトを生成している

２行目はそれをつかってアニメーション画像ファイルを生成している
引数の型について

fig : <class 'matplotlib.figure.Figure'>

update : <class 'function'>

range(8) : <class 'range'> （多分イテラブルなオブジェクトなら何でも良い）
FuncAnimationインスタンス化の前に決めること
アニメーションのfps→interval
フレーム数→framesのイテラブルオブジェクトの要素数
フレームごとの描画に必要なデータ→framesのイテラブルオブジェクトの要素に1フレーム描画のためのデータを含める
イテラブルオブジェクトの要素を一つ取り出して描画する関数→update（引数にはイテラブルオブジェクトの要素を想定すること）
matplot.figure.Figureのインスタンス

 FuncAnimationを使うために決めなければならないこと１フレーム１フレームをどう描画していくか
１フレームを描画するために必要なデータはなにか
１フレームを描画するために必要なデータをどうイテラブルオブジェクトとしてまとめるか

 トラジェクトリをプロットするどうもハミルトニアンからプロットするのはあんまり筋がよろしくない気がするので、

正準方程式をといてq, pの時間発展をプロットするのが良い気がした。
調和振動子の正準方程式の解は

q(0) = q_0, p(0) = 0 とすると

q(t) = q_0 \, \mathrm{cos}(\omega t)

p(t) = -q_0 m\omega \, \mathrm{sin}(\omega t)

 トラジェクトリプロットで考えたステップ１．パラメータである。質量と振幅と振動数を決める

q0 = 1, m = 1, ω=2π

２．１周期で十分なのでt = 0 ~ 1とする

３．tを何分割するか決める＝サンプル数を決める＝フレーム数を決める

frames = 100

４．tのリストを作成する

５．mapping用の関数を作成する

６．ｔのリストでマップする（q,p）のタプルのリストを作るようにする

７．タプルのリストでプロットする関数を作成する（２次元プロット）

８．FuncAnimationはframesにtのリスト、上記(q,p)を出力する関数とその点をプロットする関数を作る

 特に何も考えずに作成したコードimport matplotlib.pyplot as plt
import numpy as np

# 時間とそのフレーム数から時間のリストを作成
t_start = 0
t_end = 1
frame_num = 100
t_list = np.linspace(t_start, t_end, frame_num)

# 解を具体的に計算するためのパラメータを決める ただし、t=0でp=0とする
m = 1
q0 = 1
omega = 2*np.pi

# tからその時刻の解のq, p のタプルを出力する関数作成
def get_qp(t):
    q = q0 * np.cos(omega*t)
    p = q0 * m * omega * np.sin(omega*t)
    return q, p

# qのリストとpのリストをスタックした２次元配列作成
qp_list = np.apply_along_axis(get_qp, 0, t_list)

print(qp_list.shape)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(qp_list[0],qp_list[1])


from matplotlib.animation import FuncAnimation


def update(frame):
    q, p = get_qp(frame)
    ax.plot(q, p, "o")
    
anim = FuncAnimation(fig, update, frames=t_list, interval=10)
anim.save("c01.gif", writer="pillow")
粒子が動いているようには見えないので要改善

すこし改善したトラジェクトリとその粒子起動

import matplotlib.pyplot as plt
import numpy as np

# 時間とそのフレーム数から時間のリストを作成
t_start = 0
t_end = 1
frame_num = 100
t_list = np.linspace(t_start, t_end, frame_num)

# 解を具体的に計算するためのパラメータを決める ただし、t=0でp=0とする
m = 1
q0 = 1
omega = 2*np.pi

# tからその時刻の解のq, p のタプルを出力する関数作成
def get_qp(t):
    q = q0 * np.cos(omega*t)
    p = q0 * m * omega * np.sin(omega*t)
    return q, p

# qのリストとpのリストをスタックした２次元配列作成
qp_list = np.apply_along_axis(get_qp, 0, t_list)

print(qp_list.shape)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(qp_list[0],qp_list[1])


from matplotlib.animation import FuncAnimation


def update(frame):
    ax.cla()
    ax.plot(qp_list[0],qp_list[1])
    q, p = get_qp(frame)
    ax.plot(q, p, "o")
    
anim = FuncAnimation(fig, update, frames=t_list, interval=25)
anim.save("c01.gif", writer="pillow")

from IPython.display import HTML
HTML(anim.to_jshtml())

横軸と縦軸にラベルがついていないのが問題

 放物運動の正準方程式とトラジェクトリハミルトニアンは

H = K + \Phi = \frac{1}{2m}(p_1^2+p_2^2) + mgq_1
ここから 正準方程式は

\dot{q_1} = \frac{\partial H}{\partial p_1} = p_1/m

\dot{p_1} = -\frac{\partial H}{\partial q_1} =-mg
\dot{q_2} = \frac{\partial H}{\partial p_2} = p_2/m

\dot{p_2} = -\frac{\partial H}{\partial q_2} =0
正準方程式をのp1, p2解くと

p_1 = -mgt,  \quad p_2 = c(定数)
これをもとにq1, q2を解くと

q_1 = -mgt^2 + d(定数) , \quad q_2 = ct + e(定数)
初期条件q_1=q_2=p_1=p_2=0とすると

着目すべきは(q1, p1)の位相空間のみ
この方程式を元にトラジェクトリをプロットするコードを以下に記載する
import matplotlib.pyplot as plt
import numpy as np

# 時間とそのフレーム数から時間のリストを作成
t_start = 0
t_end = 10
frame_num = 100
t_list = np.linspace(t_start, t_end, frame_num)

# 解を具体的に計算するためのパラメータを決める ただし、t=0でp=0とする
m = 1
g = 9.8

# tからその時刻の解のq, p のタプルを出力する関数作成
def get_qp(t):
    q = -m*g*(t**2)
    p = -m*g*t
    return q, p

# qのリストとpのリストをスタックした２次元配列作成
qp_list = np.apply_along_axis(get_qp, 0, t_list)

print(qp_list.shape)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(qp_list[0],qp_list[1])
ax.set_xlabel("q")
ax.set_ylabel("p")


from matplotlib.animation import FuncAnimation


def update(frame):
    ax.cla()
    ax.plot(qp_list[0],qp_list[1])
    q, p = get_qp(frame)
    ax.plot(q, p, "o")
    ax.set_xlabel("q")
    ax.set_ylabel("p")
    
anim = FuncAnimation(fig, update, frames=t_list, interval=25)
anim.save("c01.gif", writer="pillow")

from IPython.display import HTML
HTML(anim.to_jshtml())


放物運動のトラジェクトリ

 今後の展開についていくための準備
 以下のガウス積分を覚えておく\int_{-\infty}^{\infty} e^{-ax^2}dx = \sqrt{\frac{\pi}{a}}

\int_{-\infty}^{\infty}x^2 e^{-ax^2}dx = \sqrt{\frac{\pi}{2}}a^{-3/2}

そらでかけるようになる.かつ.どういう関数形かふんわりとしたイメージが言語化できるまで、一時的にウェイト

 以下のガンマ関数の定義を覚えておく\Gamma (\alpha) = \int_0^{\infty} t^{\alpha - 1}e^{-t}dt \quad (\alpha > 0)

そらでかけるようになる.かつ.どういう関数形かふんわりとしたイメージが言語化できるまで、一時的にウェイト