「数式処理」と「アインシュタイン縮約記法」を使用してローレンツ変換してみた（SymPy、einsum）

# はじめに
本記事は2024/2/3（土）秋葉原ロボット部 「理論グーループ勉強会」 14:00〜17:00 にて、話した内容をまとめたものです。

理論グループ 勉強会 https://akbrobot.connpass.com/event/307848/
「座標系をローレンツ変換してみた」  使用したつぎのPythonライブラリは良さそう。
・numpy.einsum: アインシュタインの縮約記法 
・SymPy: 数式処理システム 
ローレンツ変換の参考本「一般相対性理論を一歩一歩数式で理解する」 P320、P324

# 本記事のポイント
- デバッグ出力が数式表現 「$\beta = V/c$」で見やすい。（sympy.display, sympy.Eq, sympy.subs）
- はじめは数式世界で、途中から数値世界へ移行する。（sympy.evalf または sympy.N）
- Python プログラム行数が短い、グラフ生成を含めて４０行。
   アインシュタイン縮約記法（numpy.einsum）により多重ループが１行で書ける。走る添字数分のループを書かなくてよい。

 <!-- 
# 発見
- 基底ベクトルと成分の２つの変換
- どの座標系で何のグラフを書きたいの？
--> 

# デバッグ出力 （sympy.display、ローレンツ変換行列）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/5a01fce5-cade-85fd-5995-c34474a2589f.png)

# グラフ出力　（1/2光速、基底ベクトルのローレンツ変換）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/78212/ef46f95d-4f6f-063b-0e65-7d5124491c38.png)

# Python プログラム （SymPy、numpy.einsum）

```py
import numpy as np 
import sympy
import matplotlib.pyplot as plt
from mpl_toolkits.axisartist.axislines import AxesZero
X = Y = np.arange(0, 1, 0.02)
sympy.var('V,c,beta,gamma,Lambda,e')
l = sympy.Matrix([[gamma, -gamma*beta],[-gamma*beta, gamma]])
b = sympy.Eq(beta,V/c)
g= sympy.Eq(gamma,1/sympy.sqrt(1-(V**2/c**2)))
display("Lambda=",l,g,b)
l = l.subs([(b.lhs, b.rhs), (g.lhs, g.rhs)]); display("Lambda=",l)
l0 = np.array(l.subs(V, 0).evalf(3), dtype='float64'); display("V=0",l0)
l1 = np.array(l.subs(V, c/2).evalf(3), dtype='float64'); display("V=c/2",l1)
l1inv= np.linalg.inv(l1); display("inv",l1inv) # basis vector transformation
fig = plt.figure()
ax = fig.add_subplot(1,1,1,axes_class=AxesZero)
for direction in ["xzero", "yzero"]:ax.axis[direction].set_axisline_style("-|>");ax.axis[direction].set_visible(True)
for direction in ["left", "right", "bottom", "top"]: ax.axis[direction].set_visible(False)
ax.text(-0.06, -0.02, "O")
ax.axes.xaxis.set_ticks([]);ax.axes.yaxis.set_ticks([])
# ax.autoscale(enable=False)
for color in  ['black', 'red']:
  if color == 'black':
    A = l0
    mark="x","ct"
  else:
    A = l1inv
    mark="x'","ct'"
  x = y = np.array([0] *X.shape[0])
  for i in range(Axis:=3):
    Xy = np.array([X,y]) # x axis
    XyP = np.einsum("ij,jk->ik",A,Xy)
    plt.plot(XyP[0], XyP[1], ".",color=color)
    if i==0 :plt.text(0.05+(XyP[0])[X.shape[0]-1], -0.01+ (XyP[1])[X.shape[0]-1], mark[0])
    xY = np.array([x,Y]) # ct axis
    xYP = np.einsum("ij,jk->ik",A,xY)
    plt.plot(xYP[0], xYP[1], ".",color=color)
    if i==0 :plt.text(0.01+(xYP[0])[X.shape[0]-1], 0.05+ (xYP[1])[X.shape[0]-1], mark[1])
    x = y = x+(1/(Axis-1))
```

# おわりに
つぎの２つのPythonライブラリは、使えそう。
・numpy.einsum: アインシュタインの縮約記法
・SymPy: 数式処理システム

数式レベルでのデバッグ出力は楽しい。行列を含んだ等式（sympy.Eq）は使用できない？
einsum は強力でプログラムを短くできる。

相対論や量子論の分野のPythonライブラリが多数ある。機会を見つけて使用したい。

# 参考資料
一般相対性理論を一歩一歩数式で理解する 2017 石井俊全

アインシュタインの縮約記法ってナンだ？
https://qiita.com/tshimizu8/items/a56795bdc34aaeaa2594
　
https://docs.sympy.org/dev/tutorials/intro-tutorial/intro.html#what-is-symbolic-computation

# 参考 Wikipedia
　
https://ja.wikipedia.org/wiki/SymPy
https://ja.wikipedia.org/wiki/数式処理システム
https://ja.wikipedia.org/wiki/ローレンツ変換

# 付録
## 使用環境
- Google Colaboratory 

```py
!python -V
! pip list | grep -E "numpy" 
! pip list | grep -E "sympy" 
! pip list | grep -E "matplot" 
!cat /etc/os-release | head -n 2 | tr "\n" " "; echo
!uname -a
# !nvidia-smi
# !cat /proc/cpuinfo
# !cat /proc/meminfo
```

```py
Python 3.10.12
numpy                            1.23.5
sympy                            1.12
matplotlib                       3.7.1
matplotlib-inline                0.1.6
matplotlib-venn                  0.11.10
PRETTY_NAME="Ubuntu 22.04.3 LTS" NAME="Ubuntu" 
Linux e96af331b621 6.1.58+ #1 SMP PREEMPT_DYNAMIC Sat Nov 18 15:31:17 UTC 2023 x86_64 x86_64 x86_64 GNU/Linux
```
## 相対論と量子論のPythonライブラリ情報（ChatGPT）
### 相対論



- EinsteinPy: 一般相対性理論と宇宙物理学のためのPythonライブラリ、軌道計算や時空の可視化が可能。
- NRPy+ (Numerical Relativity in Python): 数値相対論のためのPythonパッケージ、一般相対性理論の方程式から数値コードを自動生成。
- PyCBC: 重力波信号の解析と検出に特化したPythonライブラリ、LIGOやVirgoのデータ処理に使用。

### 量子論
- QuTiP (Quantum Toolbox in Python): 量子力学と量子情報学のためのツールボックス、量子システムのダイナミクスシミュレーションをサポート。
- QuantumPy: 量子力学のシミュレーションに特化したPythonライブラリ、基本的な量子計算を行うための関数を提供。
- Qutip-qip: QuTiPの一部であり、量子情報処理の実験をシミュレートするためのモジュール。

### 量子化学

- Psi4: 高度な量子化学計算を行うためのオープンソースソフトウェアパッケージです。エネルギー計算、勾配計算、その他多くの量子化学的特性の計算をサポートしています。
- PySCF: 軽量で柔軟性が高く、拡張しやすい量子化学計算のためのPythonライブラリです。密度汎関数理論（DFT）、結合クラスター理論など多岐にわたる方法をサポートしています。
- RDKit: 化学情報学および機械学習アプリケーション向けのオープンソースのツールキットです。分子の生成、操作、分析、可視化などが可能で、量子化学計算と組み合わせて使用されることがあります。