2023-01-21

（古典）解析力学：ラグランジュアンとハミルトニアンとリウビリアン

力学統計力学

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

ラグランジュアンからハミルトニアンまで

ラグランジュアン $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$

2023-01-18

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の直積で表現。

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

2023-01-18

二次元ヒストグラムの計算およびプロット（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で呼び出す関数名が異なることに注意。

2023-01-18

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を取って、残りの次元は残す」という風に理解した。
２次元の行列を例に取ると、

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

また、ndarrayの使用で、要素数が１の次元はあたかもスカラーのように加減乗除できる。
したがって、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次元配列で差を取ったものが並んでいる」とでも言えば良いだろうか。

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

2022-10-02

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）。

2021-08-12

多面体の内側に点があるかどうかを判定する。

プログラミング

多面体の頂点の座標が分かっているとする。
多面体の内側の点を原点に取り、対象とする点が多面体の内側に入っているかどうかを判断したい。

もし対象点が多面体の外側にあるとき、「多面体の頂点ベクトルを（原点から見た）対象点の方向に射影させると、射影成分は常に（原点からの）対象点の距離より小さい」が成り立つ。
逆に言えば、各頂点について調べていき、一つでも射影成分の方が大きければ、多面体の内側にあることが分かる。
訂正（2021/0812）
多面体の頂点との内積では上手く行かないことが分かった。
原点から各面に垂線を下し、面と垂線の交点と垂点と定義する。
原点から対象点へ向かうベクトルを各垂点方向へ射影したときに、射影成分が垂線の長さよりも短ければ、多面体の内側にいると判断できる。

以下、テストコード。

import numpy as np
import itertools as it

def make_vertical_points( vertices ):
    vps = []
    for v1, v2, v3 in it.combinations( vertices, 3 ):
        vv = np.cross( v2 - v1, v3 - v1 )
        ev = vv / np.linalg.norm( vv )
        pv = np.dot( v1, ev )
        nv = abs( pv )
        if np.all( nv >= np.array( [ np.dot( vertex, ev ) for vertex in vertices ] ) ) and \
           not np.any( [ np.all( pv * ev == vp ) for vp in vps ] ):
            vps.append( pv * ev )
    return vps

def check_inside_polyhedron( point, vps ):
    return np.all( [ np.dot( point, vp/np.linalg.norm(vp) ) <= np.linalg.norm( vp ) for vp in vps ] )

# example
vertices = np.array( list( it.product( [-1,1], repeat=3 ) ) )
vps = make_vertical_points( vertices )
check_inside_polyhedron( np.array( [0.9, 0, 0] ), vps )

2021-08-12

pythonでLebedev quadratureを使う（更新）

プログラミング数学

以前にquadpyについて紹介した。
koideforest.hatenadiary.com
しかし、記事が古くなっていて、仕様がいろいろ変更されたようなので、更新する。
詳しい内容は公式のGitHubページを参照。
GitHub - nschloe/quadpy: Numerical integration (quadrature, cubature) in Python

テストコード。

import numpy as np
import quadpy

def f_cart( xyz ):
    # np.shape( xyz ) -> ( 3, number of points )
    n_points = np.shape( xyz )[1]
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

def f_sphere( theta_phi ):
    # np.shape( theta_phi ) -> ( 2, number of points )
    # theta: polar : theta_phi[0]
    # phi: azimuth : theta_phi[1]
    n_points = np.size( theta_phi[0] )
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

#scheme = quadpy.u3._lebedev.lebedev_019()
scheme = quadpy.u3.get_good_scheme(19)

# \int^\pi_0 d\theta \sin\theta \int^{2\pi}_0 d\phi = 4 \pi

xyz_center = ( 0., 0., 0. )
radius = 1.
val_cart = scheme.integrate( f_cart , xyz_center, radius )

val_sphere  = scheme.integrate_spherical( f_sphere )

print( val_cart, val_sphere, 4*np.pi )

nano_exit

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

（古典）解析力学：ラグランジュアンとハミルトニアンとリウビリアン

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

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

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

FFTでフーリエ変換

多面体の内側に点があるかどうかを判定する。

pythonでLebedev quadratureを使う（更新）