2018-04-26

sympyの虚数単位を通常に戻してnumpyが使えるようにする

sympy moduleを一部使ってゴニョゴニョ計算させて、それを元にnumpy moduleで処理（例えば逆行列計算）させようとした時に虚数単位が引っかかって怒られることがあった。
sympyの代数表記を数値表記に戻す方法として"N()"や"~.evalf()"があるが、これではsympyの虚数単位"I"は通常python上で扱われる"j"に戻らない。
戻すためには、"complex()"を使う必要がある。
ちなみに、Clebsh-Gordan係数のような関数丸ごとがobjectとしてsympy上で扱われている場合、"~.doit()"で一度崩してあげる必要がある。
つまり、"complex( ~.doit() )"のような感じ。

Sympyのウェブマニュアルには一箇所だけ"complex()"が出てきて、後は"N()"とか別のを使っているので、流し読みしていると見つからなかった。
Numerical evaluation — SymPy 1.2 documentation
検索機能って本当大事ですね。

追記
例えば"N()"や"evalf()"で数値化したものが実数だったとして、それに"+1j"で虚数を足すと、sympyの虚数単位"I"に変換されてしまう。
そのため、実数だと分かっている場合には、"float()"で通常の実数にし、そのあとで"+1j"を足すことが必要。

2018-03-25

平方根の行列表現

数学

平方根を評価するのに、級数展開以外のもので何か無いか考えた時に、虚数の行列表現を思い出した。

とりあえず $\sqrt{2}$ に対応する行列を求めたい。
つまり、二乗したら単位行列に2をかけたものを返す行列を考える。

$\displaystyle \begin{align} \left( \begin{array}{cc} 2 & 0 \\ 0 & 2 \end{array} \right) = \left( \begin{array}{cc} a & b \\ c & d \end{array} \right)^2 = \left( \begin{array}{cc} a^2 + bc & b ( a + d ) \\ c ( a + d ) & d^2 + bc \end{array} \right) \end{align}$

ここで非対角成分を消すために $b = c = 0$ としてしまうと、 $a^2=d^2=\sqrt{2}$ となって全く面白くない。
そのため $a + d=0$ の時の行列を探す。
この時、 $bc = 2 - a^2$ という条件しか出て来ないため、一意には決まらないことがわかる。
とりあえず、 $a = 0, c = 1$ を課せば、

$\displaystyle \begin{align} \left( \begin{array}{cc} 0 & 2 \\ 1 & 0 \end{array} \right) \end{align}$

が $\sqrt{2}$ に対応した行列になる。
$b$ と $c$ を入れ替えた転置の表式でも二乗して同じ値になる。
もし仮に $a = 1$ とした時には、

$\displaystyle \begin{align} \left( \begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array} \right) = \left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array} \right) + \left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array} \right) = \sigma_3 + \sigma_1 \end{align}$

となり、Pauli行列の和で思わず書き直したくなる形になる。

$a$ が二乗の形で $bc$ が担う寄与を引いてくれているので、例えば $\sqrt{ x^2+x+1 }$ は

$\displaystyle \begin{align} \left( \begin{array}{cc} x & x+1 \\ 1 & -x \end{array} \right) \end{align}$

と、全ての項を一次の形で表すことが出来る。

級数展開と行列表現が行き来出来るようになると楽しそうに思う。

2018-03-20

ローレンツブーストとローレンツ収縮

力学

名前的に非常にややこしい。ブーストするのか収縮するのかどっちかにして欲しいところである。

ざっくり言えば、

ローレンツブースト：観測者（静止系）に対して、運動している系の座標軸は運動方向に伸びる（時間軸も同時に影響を受ける）。
ローレンツ収縮：（ローレンツブーストの結果）静止系から見て、運動している物体は縮んで見える。

であろう。

ポイントは、

与えられた（固定する）パラメータが、観測者（静止系）と運動系のどちらに属するものか

と思われる。

一次元系に話を絞る。
静止系 $( x(t), t )$ と速度 $v$ で運動している運動系 $( x'( t' ), t' )$ を考える。
ローレンツ因子を $\gamma = 1 / \sqrt{ 1 - v^2/c^2 }$ と定義しておく（したがって $\gamma > 1$ ）。

観測者を静止系に取る（つまり、あんまり何も考えなくて良い場合）。
ロケット（運動系）の $x'$ 座標は、観測者（静止系）のパラメータ $x, t$ を使って（基準にして）、

$\displaystyle x' = \gamma( x - v t )$

と与えられる。
座標が伸びていることをわかり易くする為、静止系から見て $v_0$ で運動しているガンダムを別に考える（ $v_0$ は静止系に属するパラメータ）。
観測者、ロケット、そしてガンダムは $t = t' = 0$ の時に原点に居たとする。
つまり、t秒後には観測者（静止系）から見てガンダムは $( v_0 t, t )$ に居ることになる。
この時に、ロケット（運動系）からガンダムはどう見えるかと、

$\displaystyle x' = \gamma( v_0 t - v t ) = \gamma ( v_0 - v ) t$

であり、 $\gamma > 1$ であるから、非相対論（ガリレイ変換）に比べて差が大きくなっていることがわかる。
これはローレンツブーストによって座標が伸びた帰結である。

（ちなみに、「観測者（ロケット）から見たら」というのは $x (x')$ の値がどうなるか、ということを言っている。）

で、これだけ見ると何でもかんでも（ロケットから見ると）長くなるように思ってしまうが、ロケット上でのパラメータが基準となるものは逆になる。
「観測者のパラメータを基準に取った時、ロケットから見ると長さが長くなる」ということは、例えば、「ロケット上から見たロケットの長さ $L'$ は、観測者から見ると短い長さ $L( < L' )$ に見える」ということと同義である。
（ロケットが地上から宇宙に「加速して」飛び立つということを考え始めると一般相対論の範疇なので、ロケットが速度を持つ過程は省略するが、）元々のロケットの長さは作成時に $L'$ と決まっていて、ロケットが飛んでいる時でも「ロケット上から見たら」長さは $L'$ である。
一方で、止まっている観測者から見ると、（ユニットを分離したり、先端が凹んだりしない限り）変わらないはずのロケットの長さが $L$ に縮んで見えてしまう。
これが「ローレンツ収縮」と言っていることである。
この場合には、そもそも観測者からは縮んで見えているので、そこからロケットでは伸びて見えると言っても、結局は元の長さで見えるというだけの話である。

静止系と運動系、どちらのパラメータを固定するべきかで見え方が変わるので、非常にややこしい。。。

2018-03-19

Sympyを使って平方完成をする。

プログラミング数学

Sympyはpythonモジュールの一つで、Mathematicaっぽく数式記号そのままで式を扱えるものである。
しかし、どうも式変形をゴリゴリ進めるノート的な使い方をするよりかは、ややこしい項が出て来た時に展開が合ってるかとか手計算が合っているかサッと電卓的に確かめるような立ち位置な気がしてならない。

ここでは、式変形の練習として、平方完成をしてみる。
https://minireference.com/static/tutorials/sympy_tutorial.pdf

from sympy import *

x, h, k = symbols( 'x h k' )

f = x**2 - 4*x + 7
g = ( x - h )**2 + k

print( solve( f - g, [ h, k ] ) )

出力は

[(2, 3)]

と返ってくるはず。
リストかつタプルなことに注意すれば、solve(...)[0][0 or 1]から h, k の値を直接持ってくることが出来る。

個人的には、専門書の式の間を埋めるのに上手く使えないかなぁと思っているが、Sympyで式変形を自在に扱うには工夫が要るように思える。

2018-03-17

ローレンツ変換行列要素における反変・共変の関係

力学量子力学数学

反変ベクトルに対する変換行列 $a^{\mu}_{ \,\, \nu }$ と共変ベクトルに対する変換行列 $a_{\mu}^{ \,\, \nu }$ は、計量テンソルで以下のように結び付けられる。

$\displaystyle a_{ \mu }^{ \,\, \nu} = g_{\mu \sigma }a^{\sigma}_{ \,\, \rho} g^{\rho \nu}$

しかし、やっぱりわかったようなわからないような、という感じに苛まれる。
ここで $k,k' = 1,2,3$ と定義すると、以下のような関係であることがわかる。

$\displaystyle \begin{align} a_0^{ \,\, 0} &= a^0_{ \,\, 0} \\ a_k^{ \,\, 0} &= - a^k_{ \,\, 0} \\ a_0^{ \,\, k} &= - a^0_{ \,\, k} \\ a_k^{ \,\, k'} &= a^k_{ \,\, k'} \\ \end{align}$

つまり、ローレンツ変換では時間と座標が混ざるが、共変ベクトルの場合には、反変ベクトルに比べてその「時間と座標の混ざり」の部分が逆符号（逆方向とまで言って良いかどうかは自信がないです）になる。
時間から時間、座標から座標の変換は反変と同じ。
これで多少見通しが良くなった気がする。

2018-03-17

特殊相対論の計量テンソルと4元ベクトルの内積

力学量子力学

計量テンソル $g_{ \mu \nu }$ を以下のように定義。

$\displaystyle \begin{align} g_{ \mu \nu } = \left\{ \begin{array}{rl} 1, & \mu = \nu = 0 \\ -1, & \mu = \nu = 1,2,3 \\ 0, & \mu \neq \nu \end{array} \right. \end{align}$

それで以下の4元ベクトルの内積が、ローレンツ変換で不変であることが特殊相対論の特徴であった。

$\displaystyle \begin{align} g_{ \mu \nu } x'{}^\mu x'{}^\nu = g_{ \mu \nu } x^\mu x^\nu \end{align}$

で、これはアインシュタインの規約を用いているから実際には和である。
しかし、頭ではわかっているのだが、心がどうも胡散臭く感じているので、全部展開してみた、というのがこの記事の目的。

$\begin{align} g_{ \mu \nu } x^\mu x^\nu = g_{ 0 0 } x^0 x^0 + g_{ 1 0 } x^1 x^0 + g_{ 2 0 } x^2 x^0 + g_{ 3 0 } x^3 x^0 \\ + g_{ 0 1 } x^0 x^1 + g_{ 1 1 } x^1 x^1 + g_{ 2 1 } x^2 x^1 + g_{ 3 1 } x^3 x^1 \\ + g_{ 0 2 } x^0 x^2 + g_{ 1 2 } x^1 x^2 + g_{ 2 2 } x^2 x^2 + g_{ 3 2 } x^3 x^2 \\ + g_{ 0 3 } x^0 x^3 + g_{ 1 3 } x^1 x^3 + g_{ 2 3 } x^2 x^3 + g_{ 3 3 } x^3 x^3 \\ = g_{ 0 0 } x^0 x^0 + g_{ 1 1 } x^1 x^1 + g_{ 2 2 } x^2 x^2 + g_{ 3 3 } x^3 x^3 \\ = x^0 x^0 - x^1 x^1 - x^2 x^2 - x^3 x^3 \\ = (ct)^2 - (x)^2 - (y)^2 -(z)^2 \end{align}$

まぁ、やる前から当たり前だが、やった後だともっと当たり前に感じることが出来るようになった。
要は、

$\displaystyle \begin{align} g_{ \mu \nu } x^\mu x^\nu = g_{ \mu \mu } x^\mu x^\mu \end{align}$

であることが既に明らかであるから、愚直に全部展開することまでする必要はなかった。

一般性を保つならば、記号を分けておいた方が良いことは良いが、一度は（必要とあれば何度でも）アホなことをするのは大事だと思った次第である。

2018-02-03

ASE: 構造クラスターを得るスクリプト

プログラミング固体物理量子力学

半径Rの構造クラスターをASEで作る。
例として単体アルミニウムを使う。

import math
from ase.spacegroup import crystal
from ase.visualize import view

# parameters
LC = 4.05
R = 7.
NSL = math.ceil( 2. * R / LC )

# make bulk structure
aluminium = crystal( 'Al', [ ( 0, 0, 0 ) ], spacegroup=225,
                     cellpar=[ LC, LC, LC, 90, 90, 90 ] )
al_sc = aluminium.repeat( ( NSL, NSL, NSL ) ) # super cell

# cut to make sphere
v_com = al_sc.get_center_of_mass()
al_sc.translate( v_com )

al_sc.wrap() # take outside atoms into unit cell by periodic boundary condition

distances = al_sc.get_all_distances()[0] # distances from the atom at the origin
del al_sc[ [ i for i, dis in enumerate( distances ) if dis > R ] ]

al_sc.translate( - v_com )

view(al_sc)

ここでは原点の原子を中心として構造クラスターを得るようにしている。
unit cell の重心は atoms.get_center_of_mass() ですぐ求まるし、距離もatoms.get_all_distances()で丸ごと求まる。
atoms.wrap()は unit cell の外側にある原子を周期的境界条件で内側に入れ込むというもの。
atoms.translate( xyz ) だけでは、ただ並進移動するだけで、unit cell の外側に原子がはみ出たままになる。
この状態でwrap をしないと、del で不要な原子を除いた時に形がおかしくなる。
最後は一応、元の原点に戻るようにしている。

nano_exit

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