SageMath/SageManifoldsの紹介¶

by YAMAMOTO, Noboru

SageMathとは？¶

SageMathは数式を記号のまま取り扱う数式処理システムの一つです。

  • 数式処理システムとしてはMathematicaが著名ですが、SageMathはOpenSourceの複数の数式処理エンジンを統合して利用できるシステムです。
  • NumPy, SciPy, matplotlib, Sympy, Maxima（老舗数式処理システム）, GAP(a System for Computational Discrete Algebra), FLINT( number theory), R(統計処理) などのプログラムを統合した環境を提供します。

SageManifoldsとは？¶

  • 微分幾何学とテンソルを取り扱うための様々な機能をSageMathのパッケージとして実現しています。
  • 一般相対性理論のさまざまな数式を取り扱うために開発されました。

NoteBook¶

  • SageMathの入力ファイルはIPythonで使われるノートブック形式ファイル(.ipynb)です。
  • Jupyterが現在の既定のインターフェイスです。
  • JupyterLabでもSageMath/SageManifoldsを実行できます。(要設定作業)

Jupyter Notebook　と Jupyterlab¶

Jupyter Notebook¶

image.png

Jupyter Notebook　と Jupyterlab¶

Jupyterlab¶

image.png

SageMath: 数式処理とは？¶

  • Mathematica/SageMath/sympyなどのプログラムは数式処理プログラムと呼ばれます。
  • その名が示すように、数式を数値に置き換えて処理するのではなく、記号のまま 処理します。
  • Jupyter/Jupyter Labでは出力結果をLaTexの数式表現を利用して、 綺麗に出力することもできます。
In [1]:
%display plain 
var('a b c x y t')
solve(a*x**2 + b*x+c == 0,x)
Out[1]:
[x == -1/2*(b + sqrt(b^2 - 4*a*c))/a, x == -1/2*(b - sqrt(b^2 - 4*a*c))/a]

出力をLaTex形式を使って整形します。HTML/MathJaxを使っています。

In [2]:
%display latex 
solve(a*x**2 + b*x+c == 0,x)
Out[2]:
\(\displaystyle \left[x = -\frac{b + \sqrt{b^{2} - 4 \, a c}}{2 \, a}, x = -\frac{b - \sqrt{b^{2} - 4 \, a c}}{2 \, a}\right]\)

因数分解¶

数式の因数分解も、整数の素因数分解も簡単です。

In [3]:
factor(x**2 + (a + b)*x + a*b)
Out[3]:
\(\displaystyle {\left(a + x\right)} {\left(b + x\right)}\)
In [4]:
factor(x^3-1,algorithm="magma")
Out[4]:
\(\displaystyle {\left(x^{2} + x + 1\right)} {\left(x - 1\right)}\)
In [5]:
factor(2**31-1)
Out[5]:
\(\displaystyle 2147483647\)

式の展開と整理¶

In [58]:
var('x y')
ex=(x+y)**4
show(ex)
ex=ex.expand()
show(ex)
\(\displaystyle {\left(x + y\right)}^{4}\)
\(\displaystyle x^{4} + 4 \, x^{3} y + 6 \, x^{2} y^{2} + 4 \, x y^{3} + y^{4}\)
In [59]:
show(ex.factor())
\(\displaystyle {\left(x + y\right)}^{4}\)

simplify_***　メソッド : 式の単純化¶

式を等価な単縦化された表現に変換します。

  • simplify, simplify_full
  • simplify_factorial, simplify_log
  • simplify_trig, simplify_hypergeometric
  • simplify_rational
  • simplify_real
In [60]:
ex=(cos(x)+1)**4 + (cos(x)-1)**4
show(ex)
show(ex.expand())
\(\displaystyle {\left(\cos\left(x\right) + 1\right)}^{4} + {\left(\cos\left(x\right) - 1\right)}^{4}\)
\(\displaystyle 2 \, \cos\left(x\right)^{4} + 12 \, \cos\left(x\right)^{2} + 2\)
In [61]:
show(ex.expand().trig_reduce())
\(\displaystyle \frac{1}{4} \, \cos\left(4 \, x\right) + 7 \, \cos\left(2 \, x\right) + \frac{35}{4}\)
In [62]:
show((x**6-1).factor())
(x**6-1).factor_list()
\(\displaystyle {\left(x^{2} + x + 1\right)} {\left(x^{2} - x + 1\right)} {\left(x + 1\right)} {\left(x - 1\right)}\)
Out[62]:
\(\displaystyle \left[\left(x^{2} + x + 1, 1\right), \left(x^{2} - x + 1, 1\right), \left(x + 1, 1\right), \left(x - 1, 1\right)\right]\)
In [63]:
var('x y')
show((log(x)+log(y)).simplify_full())
show((log(x)+log(y)).simplify_log())
\(\displaystyle \log\left(x\right) + \log\left(y\right)\)
\(\displaystyle \log\left(x y\right)\)

実行制御¶

  • SageMathはpythonをベースに作成されています。
  • Pythonの文法を使って、高度なプログラムの作成が可能です。
In [12]:
for i in range(3,32,8):
  print(2**i-1, "=", factor(2**i-1))

7 = 7
2047 = 23 * 89
524287 = 524287
134217727 = 7 * 73 * 262657

show()関数¶

show 関数は引数として与えれたオブジェクトが"rich output"を持っている場合、 それを利用して結果を表示します。図形オブジェクトではその図形が表示されます。

In [13]:
for i in range(3,32,8):
  show( factor(2**i-1))
\(\displaystyle 7\)
\(\displaystyle 23 \cdot 89\)
\(\displaystyle 524287\)
\(\displaystyle 7 \cdot 73 \cdot 262657\)

マジックコマンド¶

  • %display latexなどのような%で始まる行は(ipythonの）マジックコマンドです。
  • マジックコマンドには　セルマジック(%%)とラインマジック(%)があります。
In [14]:
%%bash
date
echo $EPICS_ROOT

2023年 6月28日 水曜日 16時36分34秒 JST
In [15]:
%display plain
%lsmagic
Out[15]:
<IPython.core.magics.basic.MagicsDisplay object at 0x11a6b5310>
In [16]:
%lsmagic?

Docstring: List currently available magic functions.
File:      /private/var/tmp/sage-10.0-current/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/IPython/core/magics/basic.py
In [17]:
var('x y')
solve([x**2+y**2 == 25, y-x==1], [x,y])
Out[17]:
[[x == -4, y == -3], [x == 3, y == 4]]
In [18]:
sols=solve([x**2+y**2 == 25, y-x==1], [x,y])
for sol in sols:
  show(sol," ", (x**2+y**2 == 25).substitute(sol), ":",(y-x ==1).subs(sol))

[x == -4, y == -3] ' ' 25 == 25 ':' 1 == 1
[x == 3, y == 4] ' ' 25 == 25 ':' 1 == 1

SageManifolds¶

  • SageManifoldsは[SageMath]の上に実装された微分可能多様体を取り扱うためのモジュールです。
  • SageManifoldsはSageMathとは独立したプロジェクトですが、現在はSageMathの正式配布版に含まれており、SageMath をインストールすれば直ちに利用可能となっています。
  • プロジェクトを開始したÉric Gourgoulhon氏は、LUTE - パリ天文台、宇宙理論研究室のDirecteur de rechercheです。一般相対性理論、ブラックホールなどを研究されています。
  • チュートリアルの日本語版がwebで閲覧可能です。 (山本が翻訳）

多様体　$M$ を定義する。¶

In [64]:
%display latex
M = Manifold(4, 'M', r'\mathcal{M_0}',  structure='Lorentzian')
print(M)

4-dimensional Lorentzian manifold M

座標系 $LC$を導入します。 $LC$の座標軸のラベルはt,x,y,z とします。

In [20]:
LC.<t,x,y,z> = M.chart()
show(LC)
\(\displaystyle \left(\mathcal{M_0},(t, x, y, z)\right)\)

多様体の計量$g$ を 導入します。

In [21]:
g = M.lorentzian_metric('g',signature="negative")
g[0,0],g[1,1],g[2,2],g[3,3]=1,-1,-1,-1
g.display(LC)
Out[21]:
\(\displaystyle g = \mathrm{d} t\otimes \mathrm{d} t-\mathrm{d} x\otimes \mathrm{d} x-\mathrm{d} y\otimes \mathrm{d} y-\mathrm{d} z\otimes \mathrm{d} z\)

別の座標系(極座標）を導入します。極座標ではそれぞれの座標の取りうる範囲が異なります。 (+ooは+Infinityの別の表記方法です)

In [22]:
XS.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
XS
Out[22]:
\(\displaystyle \left(\mathcal{M_0},(t, r, {\theta}, {\phi})\right)\)

座標系の軸は配列のように取り出せます。

In [23]:
XS[:]
Out[23]:
\(\displaystyle \left(t, r, {\theta}, {\phi}\right)\)

二つの座標系の関係を変換関数を使って定義します。

In [55]:
XS_to_LC= XS.transition_map(
  LC, [t, r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th) ]
)
XS_to_LC.display()
Out[55]:
\(\displaystyle \left\{\begin{array}{lcl} t & = & t \\ x & = & r \cos\left({\phi}\right) \sin\left({\theta}\right) \\ y & = & r \sin\left({\phi}\right) \sin\left({\theta}\right) \\ z & = & r \cos\left({\theta}\right) \end{array}\right.\)
In [25]:
#XS_to_LC= XS.transition_map(LC, 
#                              [XS[0],
#                               XS[1]*sin(XS[2])*cos(XS[3]),
#                               XS[1]*sin(XS[2])*sin(XS[3]),
#                               XS[1]*cos(XS[2])
#                             ]
#                  )
#XS_to_LC.display()

この座標系での計量の表示がどうなるかみてみましょう。

In [26]:
g.display(XS)
Out[26]:
\(\displaystyle g = \mathrm{d} t\otimes \mathrm{d} t-\mathrm{d} r\otimes \mathrm{d} r -r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} -r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}\)

逆変換を求めてみましょう。

In [27]:
try:
  XS_to_LC.inverse().dislay()
except Exception as m:
  print(m)

no solution found; use set_inverse() to set the inverse manually

残念ながらこの場合、SageManifoldsは逆変換を求めることができませんでした.

別途計算した逆変換を設定します。

In [28]:
XS_to_LC.set_inverse(t,
                     sqrt(x**2+y**2+z**2),
                     arctan2(sqrt(x**2+y**2),z),
                     arctan2(y,x))
XS_to_LC.inverse().display()

Check of the inverse coordinate transformation:
  t == t  *passed*
  r == r  *passed*
  th == arctan2(r*sin(th), r*cos(th))  **failed**
  ph == arctan2(r*sin(ph)*sin(th), r*cos(ph)*sin(th))  **failed**
  t == t  *passed*
  x == x  *passed*
  y == y  *passed*
  z == z  *passed*
NB: a failed report can reflect a mere lack of simplification.
Out[28]:
\(\displaystyle \left\{\begin{array}{lcl} t & = & t \\ r & = & \sqrt{x^{2} + y^{2} + z^{2}} \\ {\theta} & = & \arctan\left(\sqrt{x^{2} + y^{2}}, z\right) \\ {\phi} & = & \arctan\left(y, x\right) \end{array}\right.\)
In [29]:
(M.coord_change(LC,XS)*M.coord_change(XS,LC)).display()
Out[29]:
\(\displaystyle \left\{\begin{array}{lcl} t & = & t \\ r & = & r \\ {\theta} & = & \arctan\left(r \sin\left({\theta}\right), r \cos\left({\theta}\right)\right) \\ {\phi} & = & \arctan\left(r \sin\left({\phi}\right) \sin\left({\theta}\right), r \cos\left({\phi}\right) \sin\left({\theta}\right)\right) \end{array}\right.\)
In [30]:
(M.coord_change(XS,LC)*M.coord_change(LC,XS)).display()
Out[30]:
\(\displaystyle \left\{\begin{array}{lcl} t & = & t \\ x & = & x \\ y & = & y \\ z & = & z \end{array}\right.\)
In [57]:
LC_to_XS= LC.transition_map(
  XS, [t, sqrt(x**2+y**2+z**2), arctan2(sqrt(x**2+y**2),z), arctan2(y,x) ]
)
show(LC_to_XS.display())
show(LC_to_XS.inverse().display())
\(\displaystyle \left\{\begin{array}{lcl} t & = & t \\ r & = & \sqrt{x^{2} + y^{2} + z^{2}} \\ {\theta} & = & \arctan\left(\sqrt{x^{2} + y^{2}}, z\right) \\ {\phi} & = & \arctan\left(y, x\right) \end{array}\right.\) 
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [57], line 5
      1 LC_to_XS= LC.transition_map(
      2   XS, [t, sqrt(x**Integer(2)+y**Integer(2)+z**Integer(2)), arctan2(sqrt(x**Integer(2)+y**Integer(2)),z), arctan2(y,x) ]
      3 )
      4 show(LC_to_XS.display())
----> 5 show(LC_to_XS.inverse().display())

File /private/var/tmp/sage-10.0-current/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/sage/manifolds/chart.py:3608, in CoordChange.inverse(self)
   3606         list_x2_to_x1.append(x2_to_x1)
   3607 if len(list_x2_to_x1) == 0:
-> 3608     raise ValueError("no solution found; use set_inverse() to " +
   3609                      "set the inverse manually")
   3610 if len(list_x2_to_x1) > 1:
   3611     print("Multiple solutions found: ")

ValueError: no solution found; use set_inverse() to set the inverse manually

さらに一定速度で運動する系$MC$を考えてみましょう。 静止系 $LC$との変換はローレンツ変換に成ります。 この場合、逆変換はSageManifoldsが計算してくれます。

In [31]:
MC.<t_o,x_o,y_o,z_o> = M.chart()
show(MC)
\(\displaystyle \left(\mathcal{M_0},(t_{o}, x_{o}, y_{o}, z_{o})\right)\)
In [46]:
var('v')
assume(abs(v) < 1)
var('gamma')
#gamma=1/sqrt(1-v**2)
MC_to_LC= MC.transition_map(LC, 
                              [gamma*(t_o- v*x_o),
                               gamma*(x_o- v*t_o),
                               y_o,
                               z_o
                             ]
                  )
show(MC_to_LC.display())
show(MC_to_LC.inverse().disp())
\(\displaystyle \left\{\begin{array}{lcl} t & = & -{\left(v x_{o} - t_{o}\right)} \gamma \\ x & = & -{\left(t_{o} v - x_{o}\right)} \gamma \\ y & = & y_{o} \\ z & = & z_{o} \end{array}\right.\)
\(\displaystyle \left\{\begin{array}{lcl} t_{o} & = & -\frac{v x + t}{\gamma v^{2} - \gamma} \\ x_{o} & = & -\frac{t v + x}{\gamma v^{2} - \gamma} \\ y_{o} & = & y \\ z_{o} & = & z \end{array}\right.\)

運動系での計量を確認してみましょう。

In [47]:
g.display(MC)
Out[47]:
\(\displaystyle g = \mathrm{d} t_{o}\otimes \mathrm{d} t_{o}-\mathrm{d} x_{o}\otimes \mathrm{d} x_{o}-\mathrm{d} y_{o}\otimes \mathrm{d} y_{o}-\mathrm{d} z_{o}\otimes \mathrm{d} z_{o}\)
In [49]:
(M.coord_change(LC,MC)*M.coord_change(MC,LC)).display()
Out[49]:
\(\displaystyle \left\{\begin{array}{lcl} t_{o} & = & t_{o} \\ x_{o} & = & x_{o} \\ y_{o} & = & y_{o} \\ z_{o} & = & z_{o} \end{array}\right.\)
In [50]:
M.atlas()
Out[50]:
\(\displaystyle \left[\left(\mathcal{M_0},(t, x, y, z)\right), \left(\mathcal{M_0},(t, r, {\theta}, {\phi})\right), \left(\mathcal{M_0},(t_{o}, x_{o}, y_{o}, z_{o})\right), \left(\mathcal{M_0},(t_{1}, x_{1}, y_{1}, z_{1})\right)\right]\)
In [51]:
M.coord_change(XS,LC)
Out[51]:
\(\displaystyle \left(\mathcal{M_0},(t, r, {\theta}, {\phi})\right) \rightarrow \left(\mathcal{M_0},(t, x, y, z)\right)\)
In [52]:
M.coord_change(MC,LC)
Out[52]:
\(\displaystyle \left(\mathcal{M_0},(t_{o}, x_{o}, y_{o}, z_{o})\right) \rightarrow \left(\mathcal{M_0},(t, x, y, z)\right)\)
In [53]:
graph = XS.plot(LC, ambient_coords=(x,y,z), fixed_coords={t:0}, 
                number_values={r:2,th:5,ph:9}, plot_points=100, color={r:'green',th:'red',ph:'blue'}, 
                 thickness=1.5)
show(graph)

座標変換の合成¶

座標変換の合成を行ってみます。 運動系 $MC$ から　極座標系 $XS$への変換関数を求めてみます。 多様体 $M$のメソッド.atlas()は　＄M$に定義されている座標系の一覧を与えます。

In [39]:
M.atlas()
Out[39]:
\(\displaystyle \left[\left(\mathcal{M_0},(t, x, y, z)\right), \left(\mathcal{M_0},(t, r, {\theta}, {\phi})\right), \left(\mathcal{M_0},(t_{o}, x_{o}, y_{o}, z_{o})\right)\right]\)
In [40]:
MC_to_XS = M.coord_change(LC,XS) * M.coord_change(MC, LC)
show(MC_to_XS)
MC_to_XS.display()
\(\displaystyle \left(\mathcal{M_0},(t_{o}, x_{o}, y_{o}, z_{o})\right) \rightarrow \left(\mathcal{M_0},(t, r, {\theta}, {\phi})\right)\)
Out[40]:
\(\displaystyle \left\{\begin{array}{lcl} t & = & -\frac{v x_{o} - t_{o}}{\sqrt{v + 1} \sqrt{-v + 1}} \\ r & = & \frac{\sqrt{t_{o}^{2} v^{2} - 2 \, t_{o} v x_{o} - {\left(v^{2} - 1\right)} y_{o}^{2} - {\left(v^{2} - 1\right)} z_{o}^{2} + x_{o}^{2}}}{\sqrt{v + 1} \sqrt{-v + 1}} \\ {\theta} & = & \arctan\left(\frac{\sqrt{t_{o}^{2} v^{2} - 2 \, t_{o} v x_{o} - {\left(v^{2} - 1\right)} y_{o}^{2} + x_{o}^{2}}}{\sqrt{v + 1} \sqrt{-v + 1}}, z_{o}\right) \\ {\phi} & = & \arctan\left(y_{o}, -\frac{t_{o} v - x_{o}}{\sqrt{v + 1} \sqrt{-v + 1}}\right) \end{array}\right.\)
In [41]:
var('v_1')
assume(abs(v_1) < 1)
MC1.<t_1,x_1,y_1,z_1> = M.chart()
MC1_to_LC= MC1.transition_map(
  LC,
  [ (t_1- v_1*x_1)/sqrt(1-v_1**2),
    (x_1- v_1*t_1)/sqrt(1-v_1**2),
   y_1,z_1
]
                  )
show(MC1_to_LC.display())
show(MC1_to_LC.inverse().disp())
\(\displaystyle \left\{\begin{array}{lcl} t & = & -\frac{v_{1} x_{1} - t_{1}}{\sqrt{-v_{1}^{2} + 1}} \\ x & = & -\frac{t_{1} v_{1} - x_{1}}{\sqrt{-v_{1}^{2} + 1}} \\ y & = & y_{1} \\ z & = & z_{1} \end{array}\right.\)
\(\displaystyle \left\{\begin{array}{lcl} t_{1} & = & -\frac{{\left(v_{1} x + t\right)} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{v_{1}^{2} - 1} \\ x_{1} & = & -\frac{{\left(t v_{1} + x\right)} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{v_{1}^{2} - 1} \\ y_{1} & = & y \\ z_{1} & = & z \end{array}\right.\)
In [42]:
(M.coord_change(LC,MC1) * M.coord_change(MC, LC)).inverse().display()
Out[42]:
\(\displaystyle \left\{\begin{array}{lcl} t_{o} & = & -\frac{{\left(t_{1} v v_{1} - {\left(v - v_{1}\right)} x_{1} - t_{1}\right)} \sqrt{v + 1} \sqrt{-v + 1} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{{\left(v^{2} - 1\right)} v_{1}^{2} - v^{2} + 1} \\ x_{o} & = & \frac{{\left(t_{1} v - t_{1} v_{1} - {\left(v v_{1} - 1\right)} x_{1}\right)} \sqrt{v + 1} \sqrt{-v + 1} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{{\left(v^{2} - 1\right)} v_{1}^{2} - v^{2} + 1} \\ y_{o} & = & y_{1} \\ z_{o} & = & z_{1} \end{array}\right.\)
In [43]:
M.coord_change(LC,MC1).display()
Out[43]:
\(\displaystyle \left\{\begin{array}{lcl} t_{1} & = & -\frac{{\left(v_{1} x + t\right)} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{v_{1}^{2} - 1} \\ x_{1} & = & -\frac{{\left(t v_{1} + x\right)} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{v_{1}^{2} - 1} \\ y_{1} & = & y \\ z_{1} & = & z \end{array}\right.\)
In [44]:
M.coord_change(MC,MC1).display()
Out[44]:
\(\displaystyle \left\{\begin{array}{lcl} t_{1} & = & -\frac{{\left(t_{o} v v_{1} + {\left(v - v_{1}\right)} x_{o} - t_{o}\right)} \sqrt{v + 1} \sqrt{-v + 1} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{{\left(v^{2} - 1\right)} v_{1}^{2} - v^{2} + 1} \\ x_{1} & = & -\frac{{\left(t_{o} v - t_{o} v_{1} + {\left(v v_{1} - 1\right)} x_{o}\right)} \sqrt{v + 1} \sqrt{-v + 1} \sqrt{v_{1} + 1} \sqrt{-v_{1} + 1}}{{\left(v^{2} - 1\right)} v_{1}^{2} - v^{2} + 1} \\ y_{1} & = & y_{o} \\ z_{1} & = & z_{o} \end{array}\right.\)