2018-05-01

行列要素を視覚的にチェックする方法

行列要素をなんとか可視化したくて、三次元プロット散布図で表現してみた。
Pythonでデータを可視化するmatplotlib初級（3D散布図編） - しぷぜん

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

N = 20
row = range( N )
column = range( N )

X, Y = np.meshgrid( row, column )
X = X.reshape( N * N )
Y = Y.reshape( N * N )

fig = plt.figure()
ax = fig.add_subplot( 111, projection='3d' )
for i in range(3):
    matrix_reshape_1d = np.random.rand( N * N )
    # default marker size = 20
    ax.scatter( X, Y, [i] * ( N * N ), s = matrix_reshape_1d * 20. )
plt.show()

f:id:koideforest:20180501190738p:plain

例えば行列全体が時間に依存していて、それの時間発展を見たい時には良さそう。
複素数の場合には、色で位相を表現することになると思うが、色は配列で渡せないので、一点一点をいちいちプロットしないと行けない気がする。行列サイズが大きいと、おそらく途中で限界が来るだろう。

2018-04-29

Legendre多項式の直交性を python (scipy) で確かめる

プログラミング数学

いつもnumpyの多項式の使い方を忘れるので、Legendre多項式の直交関係式で練習する。
ルジャンドル多項式 - Wikipedia

$\displaystyle \int^{1}_{-1} P_{ l' }( x ) \, P_{ l }( x ) \, dx = \frac{ 2 }{ 2 l + 1 } \delta_{ l' l }$

scipy.special.legendreでLegendre多項式を作ると、poly1dというタイプの関数が帰ってくるので、そのままxを代入できる（下記リンク or サンプルスクリプト参照）。
scipy.special.legendre — SciPy v1.0.0 Reference Guide
numpy.poly1d — NumPy v1.14 Manual

積分は一次元積分のquadを使用。
Python SciPy : SciPy の積分関数の基本的使い方 | org-技術

from scipy.special import legendre
from scipy.integrate import quad

def Plx( l, x ):
    # legendre( l ) >>> poly1d([ ... ])
    # poly1d([ ... ])( x ) >>> float
    return legendre( l )( x )

def orthogonal_relation( lp, l ):
    # quad >>> ( result, error )
    return quad( lambda x: Plx( lp, x ) * Plx( l, x ), -1., 1. ) [0]

if __name__ == '__main__'
    for lp in range( 5 ):
        for l in range( 5 ):
            print( 'lp, l = ', lp, l )
            print ( 'numerical', orthogonal_relation( lp, l ) )
            if lp == l:
                print( 'analytical', 2. / ( 2 * l + 1. ) )
            else:
                print( 'analytical', 0. )

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}$

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

nano_exit

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