nano_exit

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

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

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

時間演算子については、様々な議論がなされている。
以下の文献では、ハミルトニアンと非可換な時間演算子について数学的に精密な調査をしている。
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)

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

(古典)解析力学:ラグランジュアンとハミルトニアンとリウビリアン

ラグランジュアンから出発して、ハミルトニアン、そしてリウビリアンまで概説する。

ラグランジュアン L(q,\dot q, t)が満たすラグランジュ方程式

\displaystyle
\frac{\partial L}{\partial q} - \frac{d}{d t} \frac{\partial L}{\partial \dot q} = 0

 L(q,\dot q, t)の全微分

\displaystyle
dL = \frac{\partial L}{\partial q} dq + \frac{\partial L}{\partial \dot q} d\dot q + \frac{\partial L}{\partial t} dt

ルジャンドル変換により、ハミルトニアン H(q,p,t)を定義。

\displaystyle
p(\dot q) \equiv \frac{\partial L}{\partial \dot q}
\\
\displaystyle
dL = \frac{\partial L}{\partial q} dq + p (\dot q) d\dot q + \frac{\partial L}{\partial t} dt
\\
\displaystyle
H(q,p,t) = \dot q(p) p - L(q, \dot q(p), t)
 p(\dot q) p \dot qの関数である。
 \dot q(p) \dot q pの関数として表し直したものである。

ルジャンドル変換の性質より

\displaystyle
\frac{\partial H}{\partial p} = \dot q

また、 Hの定義とラグランジュ方程式から、

\displaystyle
\frac{\partial H}{\partial q} = - \frac{\partial L}{\partial q} = - \frac{d}{d t}  \frac{\partial L}{\partial \dot q} = -  \frac{d}{d t} p = - \dot p

これで正準方程式が得られた。

位相空間の座標ベクトルを \Gamma = (q,p)とする。
位相空間の密度分布 \rho(\Gamma, t)の時間発展を流体として扱う。圧縮・非圧縮はまだ仮定していないが、後で正準方程式を満たす運動系は位相空間で非圧縮流体的に振る舞うことが分かる(リウビルの定理)。
保存流は (\rho, \rho \dot \Gamma)であるため、連続の方程式は、

\displaystyle
\frac{\partial \rho}{\partial t} + \nabla_\Gamma \cdot (\rho \dot \Gamma) = 0
\\
\displaystyle
\frac{\partial \rho}{\partial t} + (\nabla_\Gamma \rho) \cdot \dot \Gamma + \rho_\Gamma (\nabla \cdot \dot \Gamma) = 0

非圧縮流体の場合には、 \nabla_\Gamma \cdot \dot \Gamma=0となるが、これが正準方程式により満たされることを示す。

\displaystyle
\nabla_\Gamma \cdot \dot \Gamma = \frac{\partial}{\partial q} \dot q + \frac{\partial}{\partial p} \dot p = \frac{\partial}{\partial q} \frac{\partial H}{\partial p} - \frac{\partial}{\partial p} \frac{\partial H}{\partial q} = 0

これにより、 \rhoが満たすべき方程式が得られる。

\displaystyle
\frac{\partial \rho}{\partial t} + (\nabla_\Gamma \rho) \cdot \dot \Gamma = 0
\\
\displaystyle
\frac{\partial \rho}{\partial t} = - \left(\frac{\partial \rho}{\partial q}, \frac{\partial \rho}{\partial p}\right) \cdot (\frac{\partial H}{\partial p}, -\frac{\partial H}{\partial q}) = \frac{\partial H}{\partial q} \frac{\partial \rho}{\partial p} - \frac{\partial H}{\partial p} \frac{\partial \rho}{\partial q} \equiv \left\{H, \rho \right\}_{PB}

したがって、古典リウビル方程式は、

\displaystyle
\frac{\partial \rho}{\partial t} = \left\{H, \rho \right\}_{PB} \equiv \mathcal{L} \rho

2次元ヒストグラムにおける最適輸送距離

異なる2次元ヒストグラム間での最適輸送距離を計算してみた。

import numpy as np
import itertools
import ot
from matplotlib import pyplot as plt

n = 100
bins_  = 10
range_ = ([0,1],[0,1])
d_TF   = True

# input 2-dim. histogram-1
r = np.random.rand(n,2).T
f_,x_,y_ = np.histogram2d(r[0],r[1],bins=bins_,range=range_,density=d_TF)
plt.hist2d(r[0],r[1],bins=bins_,range=range_,density=d_TF)
plt.show()
a = f_.reshape(bins_*bins_)  # input 2-dim. histogram converted into 1-dim. array.

# input 2-dim. histogram-2
r = np.random.rand(n,2).T
f_,x_,y_ = np.histogram2d(r[0],r[1],bins=bins_,range=range_,density=d_TF)
plt.hist2d(r[0],r[1],bins=bins_,range=range_,density=d_TF)
plt.show()
b = f_.reshape(bins_*bins_)  # input 2-dim. histogram converted into 1-dim. array.

# input cost matrix (distance between each coordinate)
v_ = itertools.product(x_[:-1],y_[:-1])
v  = np.array(list(v_))
c_ = v[np.newaxis]-v[:,np.newaxis]
c  = np.linalg.norm(c_,axis=2)

P  = ot.emd(a,b,c)  # optimal transfer matrix
OT = ot.emd2(a,b,c)  # optimal transfer distance
print(OT)

二次元ヒストグラムを一次元配列に変換(行列の要素を一行に並べ直した)し、コスト行列を各座標間の距離で定義。
座標行列はitertools.productでx,yの直積で表現。

三次元以上の空間にも容易に拡張できると思われる。

二次元ヒストグラムの計算およびプロット(python)

いつも忘れるので個人用のメモ

参考ページ
matplotlib - hist2d で2次元ヒストグラムを作成する方法 - pystyle
[Python]Matplotlibで2次元ヒストグラムを作成する方法 - Qiita
NumPyでヒストグラムを作るnp.histogram関数の使い方 - DeepAge

ヒストグラムの棒の数がbinsで、bin"s"なのはbinの複数系なのではなく、binsでビンと読むらしい。
ヒストグラム - Wikipedia
コードを書くときにいつもsを付け忘れてエラーになり、その度に引数をググっていた。

適当に二次元配列の乱数から二次元ヒストグラムを作る。

import numpy as np
from matplotlib import pyplot as plt

# setting for distribution
n = 100
r = np.random.rand(n,2).T

# setting for histogram
bins_ = 5  # if you want to specify different bins between x and y, you can use tuple or list of (bins_x, bins_y).  
range_ = ([0,1],[0,1])
d_TF = True  # True: histogram will be normalized

# more fine tuning for bins
# min_ = 0
# max_ = 1
# interval_ = 0.2
# bins_ = int((max_-min_)/interval_)

# plot histogram
plt.hist2d(r[0],r[1],bins=bins_,range=range_,density=d_TF)
plt.show()

# get histogram data
frequency, x_axis, y_axis = np.histogram2d(r[0],r[1],bins=bins_,range=range_,density=d_TF)

numpyとmatplotlib.pyplotで呼び出す関数名が異なることに注意。

numpyで異なる離散分布における各点間のノルムを総当たりに得る方法

最適輸送距離について調べていた時に、何をやっているのかパッとわからないノルムの取り方をしているコードがあった。
最適輸送入門

n = 100 # 点群サイズ
mu = np.random.randn(n, 2) # 入力分布 1
nu = np.random.randn(n, 2) + 1 # 入力分布 2
C = np.linalg.norm(nu[np.newaxis] - mu[:,np.newaxis], axis=2) # コスト行列

mu, nuは二次元分布を表しており、Cがmuの各点とnuの各点の間の距離を表す行列になっている。

Cの計算(定義)の部分は、以下の愚直な内容と同じになっている。

C = np.zeros([n,n])
for i1, mu_ in enumerate(mu):
    for i2, nu_ in enumerate(nu):
        C[i1][i2] = np.linalg.norm(mu_-nu_)

np.linalg.normのaxisは以下で説明されている。
normのordとaxisの役割がわからない

個人的には、normのaxisの意味は「指定したaxisの自由度を潰す様にnormを取って、残りの次元は残す」という風に理解した。
2次元の行列を例に取ると、

  • axis=0: 行を潰す様にノルムを取るため、残るのは列方向の横ベクトル
  • axis=1: 列を潰す様にノルムを取るため、残るのは行方向の縦ベクトル(実際の出力では普通に1次元配列(横ベクトル))

また、ndarrayの使用で、要素数が1の次元はあたかもスカラーのように加減乗除できる。
したがって、np.newaxisで増やした次元部分は要素数が1しかないため、v1, v2が要素数nの1次元配列だとすると、

  • v1[np.newaxis]: (1,n)の2次元配列(行列)
  • v2[:, np.newaxis]: (n,1)の2次元配列(行列)

 \therefore v1[np.newaxis] - v2[:, np.newaxis]  \Rightarrow (n,n)の2次元配列となる。
v1[np.newaxis] - v2[:, np.newaxis]の各列は、v1のある1要素に対して、v2の各要素で差を取ったものが列の要素として並ぶ。

C = np.linalg.norm(nu[np.newaxis] - mu[:,np.newaxis], axis=2)に対し、

  • nu[np.newaxis]: (1,n,2)の3次元配列
  • mu[:, np.newaxis]: (n,1,2)の3次元配列
  • nu[np.newaxis] - mu[:,np.newaxis]: (n,n,2)の3次元配列

であるため、np.linalg.norm(nu[np.newaxis] - mu[:,np.newaxis], axis=2)で要素数2の次元(x,y座標に対応するベクトル部分)でノルムを取り、その次元を潰した(n,n)の2次元配列になる。
説明が難しいが、「nu[np.newaxis]のaxis=1の要素(要素数2の1次元配列)を固定したときに、mu[:, np.newaxis]のaxis=0における(要素数2の)各1次元配列で差を取ったものが並んでいる」とでも言えば良いだろうか。

そういう訳で、一見何をしているかわからなかったが、結果は非常に単純なものが得られていることがわかった。

FFTでフーリエ変換

昔の記事で、FFTフーリエ変換する内容を紹介した。
koideforest.hatenadiary.com
koideforest.hatenadiary.com

そこから、フーリエ変換をパッと計算出来るヘルパー関数をまとめていなかったので、自分用にまとめる。

from scipy.fftpack import fft, fftfreq, fftshift, ifft

def calc_FT(x,y):
    N=len(x)
    dx=x[1]-x[0]
    tilde_y=fftshift(fft(y)*dx)
    k      =fftshift(fftfreq(N,dx)*2*np.pi)
    return k, tilde_y

def calc_IFT(k,tilde_y):
    N=len(k)
    dk=k[1]-k[0]
    # timing of fftshift is different between x and y
    y=ifft(fftshift(tilde_y))
    x=fftshift(fftfreq(N,dk/2/np.pi))
    y/=x[1]-x[0]
    return x, y

def calc_ILFT(k,tilde_y):
    x,y = calc_IFT(k,tilde_y)
    return x-x[0], y

calc_ILFTは横軸の原点が0始まりの場合の逆フーリエ変換(LはLaplaceのL)。