nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

pythonでヒストグラムを作成し、その半値全幅を算出する。

参考サイト
NumPyでヒストグラムを作るnp.histogram関数の使い方 - DeepAge
scipy.signal.peak_widths — SciPy v1.12.0 Manual

import numpy as np
from scipy.signal import find_peak, peak_widths
from matplotlib import pyplot as plt

# data file: data.dat
# data.dat has 2 raw data, and second raw data is to be plotted as a histogram

# read data from data file
data = np.loadtxt('data.dat').T

# check histogram by plotting
plt.hist(data[1])
plt.show()

# get half maximum of full width of the histogram
histogram = np.histgram(data[1])
# histogram[0]: intensity (y-axis)
# histogram[1]: values (x-axis)
peaks = find_peaks(histogram[0])
results_half = peak_widths(histogram[0], peaks, rel_height=0.5)
half_widths = results_half[0] * (histgram[1][1] - histogram[1][0])

# check peak and widths by plotting
plt.plot(histogram[0])
plt.plot(peaks, histogram[0][peaks], "x")
plt.hlines(*results_half[1:])
plt.show()

sedで、不要な行を消去しながら、範囲を指定して置換し、結果を上書きする

参考サイト
sed で指定した範囲の行を置換する - まくまくsed/awkノート
sed | テキストの置換処理を得意とするスクリプト言語
sedでこういう時はどう書く? #Linux - Qiita

  • text.txt内の1~10行目を無視したい。
sed -e '1,10d' text.txt
  • 無視した後の残された行を、頭から1行目と数え直して、100行目から1000行目の間にある全てのhogeHOGEに置換したい。
sed -e '1,10d' -e '100,1000 s/hoge/HOGE/g' text.txt
  • 上記の内容を、同じtext.txtに反映させたい。
sed -i -e '1,10d' -e '100,1000 s/hoge/HOGE/g' text.txt

※giにすると、大文字と小文字を区別しなくなる。

  • #を含む行を消したい。
sed -e '/#/ d' text.txt
  • 行頭に#がある行を消したい。
sed -e '/^#/ d' text.txt
  • 行頭の#を消したい。
sed -e 's/^#//' text.txt
  • 行頭に#がある行のhogeだけ全てHOGEに置換したい。
sed -e '/^#/ s/hoge/HOGE/g' text.txt

- 行頭に#がない行のhogeだけ全てHOGEに置換したい。
>|bash|
sed -e '/^#/! s/hoge/HOGE/g' text.txt

オイラー・ラグランジュ方程式の導出におけるラグランジュアンの偏微分

ラグランジュアンの微分について、佐久間さんの呟きを見つけた。

昔の自分の記事で考察したことがあったのを思い出した。
koideforest.hatenadiary.com
改めて見直すと、汎関数微分を全微分で考えようとしていたり、不正確なところもあり、反省した。
その上で、数学と物理の考え方の違いが自分の中で見えた気がする。

ラグランジュアン L(q(t),\dot q(t),t)の考え方に、二つあると見た。

  • 物理:最終的には時間の関数と思う。
  • 数学:独立な2変数(もしくは3変数)を引数に持つ関数と思う。

 q(t) \dot q(t)は時間の関数だから、 L(q(t),\dot q(t),t) q(t) \dot q(t)の具体的な形を入れると、結局は時間変数 tだけで表されることになる。
荒っぽく言えば、時間だけの関数を都合の良いように (q(t),\dot q(t),(t))で分けて、それぞれを独立だと思って偏微分すると、それらしいものが得られる。

一方で、本来、 Lを関数として考えるならば、定義域をしっかり規定する必要がある。 L(q,\dot q,(t))の形にする時点で、独立な変数の組 (w,v,\tau)の全領域(もしくは物理として必要な領域)で L(w,v,(\tau))が定義されていて、その微分導関数)も(物理として使う上で)定義されている必要がある。
この導関数は、独立な変数の組 (w,v,\tau)に対する微分操作によって得られる。
物理として L(q(t),\dot q(t),(t))やその導関数を使う場合には、独立な変数の組 (w,v,\tau)値として (q(t),\dot q(t),(t))を代入したと考えるべきなのだろう。

熱力学の偏微分も、同じようなことが起きていると思われる。

単純な(熱力学)エントロピーのモデル

単純なエントロピーのモデルを弄ってみる。
参考文献:

  • 清水明、『熱力学の基礎(第2版)』東京大学出版会
  • H.B. Callen, Thermodynamics and an introduction to thermostatistics, 2nd edition (Wiley, 1985)

エントロピー表示

\displaystyle
S(U,V,N) = (UVN)^{1/3} \, (U \ge 0)

  • 同次性


\displaystyle
S(\lambda U,\lambda V,\lambda N) = (\lambda^3 UVN)^{1/3} = \lambda (UVN)^{1/3} = \lambda S(U,V, N)

  • 凸性(上に凸)


\displaystyle
\frac{ \partial^2 S}{\partial U^2} = -\frac{2}{9} \left(\frac{VN}{U^5}\right)^{1/3} \le 0

エネルギー表示

\displaystyle
U(S,V,N) = \frac{S^3}{VN}

  • 同次性


\displaystyle
U(\lambda S,\lambda V,\lambda N) = \frac{\lambda^3}{\lambda^2} \frac{S^3}{VN}= \lambda U(S,V, N)

  • 凸性(下に凸)


\displaystyle
\frac{ \partial^2 U}{\partial S^2} = 6 \frac{S}{VN} \ge 0
\\
\displaystyle
\frac{ \partial^2 U}{\partial V^2} = 2 \frac{S}{V^3N} \ge 0

温度

\displaystyle
\frac{ \partial U}{\partial S} = 3 \frac{S^2}{VN} = \frac{3}{S} U = T
\\
\displaystyle
U(S, V, N) =\frac{S}{3} T(S, V, N)
\\
\displaystyle
T(S, V, N) =3\frac{S^2}{VN} \,\, \therefore S = \left(\frac{T}{3} V N\right)^{1/2}

ヘルムホルツの自由エネルギー

\displaystyle
F(T,V,N) = U(T,V,N) - T S(T,V,N) = -\frac{2}{3} T S(T, V, N)
\\
\displaystyle
\quad = -\frac{2}{3} \left(\frac{T^3}{3} V N\right)^{1/2}

  • 凸性(温度に対して上に凸)


\displaystyle
\frac{ \partial^2 F}{\partial T^2} = -\frac{1}{2} \left(\frac{1}{3T} V N\right)^{1/2} \le 0
温度一定の条件を課せば、凸性は内部エネルギーと一致する。
このことは、凸性がルジャンドル変換で逆になっていても、熱浴下による温度一定条件であれば、平衡状態は  U と同様に  F を最小することで得られることに対応している。

(静的)構造因子の正値性について

光・X線中性子回折/散乱において、構造因子が重要である。
X線も光であるが、一般的な光散乱とX線回折で使い方が異なるのでここでは分けた)

時間依存性の無い場合には静的構造因子と呼ばれる。
(静的)構造因子は、密度関数のフーリエ変換を指す場合と、密度自己相関関数のフーリエ変換を指す場合がある。
前者はX線回折で、後者は各種散乱実験で使われる。

ここでは、後者の密度自己相関関数のフーリエ変換の正値性について考察する。

\displaystyle
\rho(\vec r) \equiv \sum_{i=1}^N \delta( \vec r - \vec r_i )
  = \int d\vec k\, e^{i \vec k \cdot \vec r} \rho(\vec k)
\\
\displaystyle
\rho(\vec k) = \int \frac{d \vec r}{2\pi}\, e^{-i \vec k \cdot \vec r} \rho(\vec r)
  = \frac{1}{2\pi} \sum_{i=1}^N e^{-i \vec k \cdot \vec r_i}
  = \rho^*(-\vec k)
\\
\displaystyle
S(\vec r)
  \equiv \frac{1}{N} \int d\vec r'\rho(\vec r') \rho(\vec r' + \vec r)
\\
\displaystyle
\quad
  = \frac{1}{N}\int d\vec r' \int d\vec k_1 \int d\vec k_2\, e^{i \vec k_1 \cdot \vec r'} e^{i \vec k_2 \cdot (\vec r' + \vec r)} \rho(\vec k_1) \rho(\vec k_2)
\\
\displaystyle
\quad
  = \frac{1}{N} \int d\vec k_1 \int d\vec k_2\, (2\pi \delta(\vec k_1 + \vec k_2)) e^{i \vec k_2 \cdot \vec r} \rho(\vec k_1) \rho(\vec k_2)
\\
\displaystyle
\quad
  = \frac{2\pi}{N} \int d\vec k\, e^{i \vec k \cdot \vec r} \rho(\vec k) \rho(-\vec k)
\\
\displaystyle
S(\vec k) = \int \frac{d\vec r}{2\pi}\, e^{-i \vec k \cdot \vec r } S(\vec r)
  = \frac{1}{N} \rho(\vec k) \rho(-\vec k)
  = \frac{1}{N} \left| \rho(\vec k) \right|^2
 S(\vec k)の正値性は、 \left| \rho(\vec k) \right|^2より明らかではあるが、具体的に中身を展開したときに一見すると正値性が保たれているか判別付かない。

\displaystyle
N S(\vec k)
  = \rho(\vec k) \rho(-\vec k)
  = \sum_{i=1}^N \sum_{j=1}^N e^{-i \vec k \cdot (\vec r_i - \vec r_j ) }
  = N + \sum_{i \neq j} e^{-i \vec k \cdot (\vec r_i - \vec r_j ) }
\\
\displaystyle
\sum_{i \neq j} e^{-i \vec k \cdot (\vec r_i - \vec r_j ) }
  = \sum_{i=2}^N \sum_{j=1}^{i-1} \left( e^{-i \vec k \cdot (\vec r_i - \vec r_j ) } + e^{i \vec k \cdot (\vec r_i - \vec r_j ) } \right)
\\
\displaystyle
\quad
  = 2 \sum_{i=2}^N \sum_{j=1}^{i-1} \cos \left( \vec k \cdot (\vec r_i - \vec r_j ) \right)
この cos関数の和が -N/2より小さくなると、 S(\vec k)は負になる。
ただし、 S(\vec k) \propto |\rho(\vec k)|^2であるため、それはあり得ないということになるが、本当にそうなのかパッとわからない。
そこで、簡単のため一次元の場合にこれを考察してみる。

簡略化として k = 1とし、 cos関数がなるべく負の値を持つことを考え、 r_{i+1} - r_{i} = \piの場合を考える。
これは粒子が波の逆位相の位置に整列していることに対応する。
そして、よく考えると、一個飛ばした先の粒子とは同位相になる。すなわち r_{i+2} - r_i = 2\pi
そのため、cos関数の値はプラスとマイナスが混合することになる。
和を計算して評価すると

\displaystyle
\sum_{i=2}^N \sum_{j=1}^{i-1} \cos \left( \vec k \cdot (\vec r_i - \vec r_j ) \right)
  = \sum_{i=2}^N \sum_{j=1}^{i-1} (-1)^{i-j}
  = \sum_{i=2}^N \sum_{j=1}^{i-1} (-1)^j
  = - \sum_{i=2}^N \chi_{i}^e
ここで、 \chi_i^e iが偶数の時に1、奇数のときにゼロとなる特性関数である。
同様に、 \chi_i^o iが奇数の時に1、偶数のときにゼロとなる特性関数として定義しておく。
よって、

\displaystyle
  \sum_{i=2}^N \chi_{i}^e
  = \sum_{i=1}^{N-1} \chi_{i-1}^e
  = \sum_{i=1}^{N-1} \chi_{i}^o
  \le N/2
したがって、今の設定において

\displaystyle
\sum_{i=2}^N \sum_{j=1}^{i-1} \cos \left( \vec k \cdot (\vec r_i - \vec r_j ) \right) \ge - N/2
が示せた。不等式で書いているが、今の場合にはほとんどイコールである。つまり、 S(\vec k)全体としてほぼゼロになる。
逆に言えば、 S(\vec k)全体としてほぼゼロになるのは、逆位相に粒子が配列しているときということになる。

逆位相ではなく、同位相に粒子が配列するときには、cos関数が常に+1となるため、 S(\vec k)が大きな値を持つ。
一方、粒子が波数kの持つ周期に対して等間隔、もしくは乱雑に位置しているとき、cos関数の和の部分は打ち消しあってゼロになり、 S(\vec k)\approx 1となる。

このように、 S(\vec k)が正値を持つことを、具体的な場合を考えて考察することができた。

リウビリアンから時間演算子を構築する妄想

プリゴジンの『存在から発展へ』では、パイこね変換の回数を「時間」として見做し、パイこね変換を基にして時間演算子を定義した後にリウビリアンとの非可換性を論じている。
存在から発展へ【新装版】 | みすず書房

時間演算子については、様々な議論がなされている。
以下の文献では、ハミルトニアンと非可換な時間演算子について数学的に精密な調査をしている。
https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1266-4.pdf

ここでは、位置と運動量の非可換性のアナロジーから、リウビリアンから時間演算子を構築できないか妄想する。

  • (空間)並進操作

位置演算子が先にあるとして、運動量演算子を定めていく過程をおさらいする。(参照:JJサクライ)
 Tが空間並進操作の演算子として、以下の性質を満たすと「定義」する。

\displaystyle
T(dx) \left| x \right\rangle
  = \left| x + dx\right\rangle
  = \left( 1 + dx \frac{\partial}{\partial x} \right) \left| x \right\rangle 
\\
\displaystyle
T(dx) \left| \psi \right\rangle
  = T(dx) \int dx' \left| x' \right\rangle \left\langle x' | \psi \right\rangle
  = \int dx' \left| x' + dx \right\rangle \left\langle x' | \psi \right\rangle
\\
\displaystyle
\quad
  = \int dx'' \left| x'' \right\rangle \left\langle x''- dx | \psi \right\rangle
  = \int dx'' \left| x'' \right\rangle \left( 1 - dx \frac{\partial}{\partial x} \right) \left\langle x'' | \psi \right\rangle

並進操作として、更に以下の性質を要請する。

  1. ユニタリー性(並進操作の前後でノルムが保存): T(dx)^\dagger T(dx) = 1
  2. 加法性: T(dx_1) T(dx_2) = T(dx_1 + dx_2)
  3. 恒等操作: T(0)=1
  4. 逆元:  T(-dx) = T^{-1}(dx)

2~4により、並進操作はアーベル群を為す。
これらを満たす T(dx)の表現として、エルミート演算子 O(=O^\dagger)を用いて以下のように表せる。

\displaystyle
T(dx) = 1 - i O dx

 Oを用いた表現と、元々の T(dx)の性質を比べると、 O微分演算を含むことがわかる。

\displaystyle
T(dx) \left| \psi \right\rangle
  = \int dx'' \left| x'' \right\rangle \left( 1 - i dx O \right) \left\langle x'' | \psi \right\rangle
\\
\displaystyle
\therefore
O = \frac{1}{i} \frac{\partial}{\partial x}

 Oの意味を考えるために、 x, T(dx)の交換関係に着目すると、

\displaystyle
[x, T(dx)] \left| x \right\rangle
  = (x+dx) \left| x + dx \right\rangle - x \left| x + dx \right\rangle
  = dx \left| x + dx \right\rangle
\\
\displaystyle
\quad
  = dx \left| x \right\rangle
\\
\displaystyle
\quad
  = - i dx [x, O] \left| x \right\rangle
\\
\displaystyle
\therefore
[x, O] = i

正準交換関係に合わせるために、次元の調整でプランク定数(を 2\piで割ったディラック定数 \hbar)を用いることで、運動量演算子 pが導入される。

\displaystyle
O = \frac{p}{\hbar}
\\
\displaystyle
T(dx) = 1 - i \frac{p}{\hbar} dx
\\
\displaystyle
\therefore
[x, p] = i \hbar

また、無限小ではなく有限の並進操作では、

\displaystyle
T(\Delta x+dx) = T(dx) T(\Delta x) = \left(1 - i \frac{p}{\hbar} dx\right) T(\Delta x)
\\
\displaystyle
\frac{T(\Delta x+dx) - T(\Delta x)}{dx} = - i \frac{p}{\hbar} T(\Delta x)
\\
\displaystyle
\frac{\partial}{\partial x} T(\Delta x) = - i \frac{p}{\hbar} T(\Delta x)
\\
\displaystyle
\therefore
T(\Delta x) = e^{- i \frac{p}{\hbar} \Delta x}

  • リウビリアンからの時間演算子(妄想)

リウビリアン \mathcal{L}は密度行列 \rhoに対して、(量子)リウビル方程式を満たす。

\displaystyle
i \frac{\partial}{\partial t} \rho = [ H, \rho ] = \mathcal{L} \rho

時間演算子 \mathcal{T}に満たして欲しい性質として、 \mathcal{L}との非可換性を課す。

\displaystyle
[\mathcal{L}, \mathcal{T}] \neq 0
なぜなら、位置と運動量の交換関係から非可換性には微分演算が重要であり、 \mathcal{L}には既に(絶対時間・外部時間としての)時間微分の意味が含まれている。
そのため、(絶対時間・外部時間と対応するような内部時間を測る) \mathcal{T} \mathcal{L}と非可換であって欲しい。
また、非可換な \mathcal{T}が作れれば、 \mathcal{T}の減少関数を介してリアプーノフ関数が作れるらしい(by プリゴジン)。

 [\mathcal{L}, \mathcal{T}]の非可換性から話を進めるためには、 \mathcal{L}固有値・固有状態があることが要請される。

\displaystyle
\mathcal{L} \rho(z) = z \rho(z)
 \mathcal{L}を対角化する基底を選択した場合、逆に \mathcal{T} \sim \frac{\partial}{\partial z}となることが予想される。これは、運動量基底を選択した際に、位置演算子が運動量の微分演算になることに対応する。
 \mathcal{T} z微分演算になるということは、 zを並進移動させる操作の母関数になるということである。
そのため、 zを並進移動させる意味がわかれば、一般的な物理系に対して \mathcal{T}を構築できる可能性があるように思える。

、、、というのが妄想話である。

そもそも、 \mathcal{L}固有値・固有状態がよくわからない。
開放系量子系において、 \mathcal{L}の固有状態は共鳴・反共鳴状態に対応するらしいが、詳しくは存じ上げない。
https://yamadazaidan.jp/wordpress/wp-content/uploads/2021/09/2012s_Hatano.pdf

また、 zの無限小変換が定義されている必要があるため、必然的に zは連続的でなければならない。
これについては、離散系では [\mathcal{L}, \mathcal{T}]が可換になるという話や、参考文献で基本的に散乱のような連続状態において議論されていることと対応しているように思える。

まずは(量子)リウビル方程式ともっと仲良くなる必要がありそうだ。

摩擦のある直線運動のラグランジュアン、ハミルトニアン、リウビリアン

前回、ラグランジュアンからハミルトニアン、リウビリアンまでの流れをまとめた。
koideforest.hatenadiary.com

今回は、具体的な問題として、摩擦のある直線運動を扱う。


\displaystyle
m \ddot q = - m \gamma \dot q

普通に解くと、

\displaystyle
\log \dot q(t) = - \gamma t + C
\\
\displaystyle
\dot q(t) = v_0 e^{- \gamma t} \quad (\dot q(0) = e^C \equiv v_0)
\\
\displaystyle
q(t) = q_0 + \frac{v_0}{\gamma} \left( 1 - e^{- \gamma t} \right)

天下り的だが、ニュートン方程式の解を満足するようにラグランジュアンを構築する。

\displaystyle
T = \frac{m}{2} \dot q^2 e^{\gamma t}, \, V = 0
\\
\displaystyle
L = T - V = \frac{m}{2} \dot q^2 e^{\gamma t}
\\
\displaystyle
\frac{\partial L}{\partial q} - \frac{d}{dt} \frac{\partial L}{\partial \dot q} = 0 = \frac{d}{dt} \frac{\partial T}{\partial \dot q} = \frac{d}{dt} \left( m \dot q e^{\gamma t} \right) = m\ddot q e^{\gamma t} + m \gamma \dot q e^{\gamma t} = 0
\\
\displaystyle
\therefore
m \ddot q = - m \gamma \dot q
構築したラグランジュアンが、元のニュートン方程式を再現することが確認できた。


\displaystyle
p(\dot q) = \frac{\partial L}{\partial \dot q} =  \frac{\partial T}{\partial \dot q} = m \dot q e^{\gamma t}
\\
\displaystyle
\dot q (p) = \frac{p}{m} e^{-\gamma t}
\\
\displaystyle
L(q, \dot q(p), t) = \frac{m}{2} \dot q^2(p) e^{\gamma t} = \frac{m}{2} e^{\gamma t} \left( \frac{p}{m} e^{-\gamma t} \right)^2 = \frac{p^2}{2m} e^{-\gamma t}
\\
\displaystyle
H = \dot q (p) p - L(q,\dot q(p), t) =  \frac{p^2}{m} e^{-\gamma t} - \frac{p^2}{2m} e^{-\gamma t} = \frac{p^2}{2m} e^{-\gamma t}

  • リウビリアン


\displaystyle
\frac{\partial \rho}{\partial t}
  = \frac{\partial H}{\partial q} \frac{\partial \rho}{\partial p} - \frac{\partial H}{\partial p} \frac{\partial \rho}{\partial q}
  = - \frac{p}{m} e^{-\gamma t} \frac{\partial}{\partial q} \rho
  \equiv \mathcal{L} \rho
\\
\displaystyle
\mathcal{L} = - \frac{p}{m} e^{-\gamma t} \frac{\partial}{\partial q}

リウビル方程式から軌跡を復元する。
以下の偏微分方程式を考える。

\displaystyle
\left(\frac{\partial}{\partial t} + \frac{p}{m} e^{-\gamma t} \frac{\partial}{\partial q}\right) \rho = 0

一階の偏微分方程式は以下の方法で解ける。
自然科学のための数学2014年度第29講
以下の関数 f(q,p,t)は、対象の偏微分方程式を満足する。

\displaystyle
f(q,p,t) = \frac{p}{m} \frac{1}{\gamma} e^{-\gamma t} + q
\\
\displaystyle
\left(\frac{\partial}{\partial t} + \frac{p}{m} e^{\gamma t} \frac{\partial}{\partial q}\right) f(q,p,t) = 0
この fを引数とする任意の関数 F(f)偏微分方程式を満足する。
よって、 f(q,p,t)=constとなる q,p位相空間上での運動の軌跡に対応する。

\displaystyle
\frac{p}{m} \frac{1}{\gamma} e^{-\gamma t} + q = C
\\
\displaystyle
q = C -\frac{p}{m} \frac{1}{\gamma} e^{-\gamma t}
\\
\displaystyle
q(t=0) = q_0 = C - \frac{p}{m} \frac{1}{\gamma}, \quad C = q_0 + \frac{p}{m} \frac{1}{\gamma}
\\
\displaystyle
q(t) = q_0 + \frac{p}{m} \frac{1}{\gamma} \left( 1 - e^{-\gamma t} \right)
  \equiv q_0 + \frac{v_0}{\gamma} \left( 1 - e^{-\gamma t} \right), \quad (p/m = v_0)

よって、摩擦項に由来する時間依存性が入ったハミルトニアンの場合でも、リウビル方程式は元々のニュートン方程式がもたらす運動の軌跡を再現することがわかった。