by YAMAMOTO, Noboru
SageMathは数式を記号のまま取り扱う数式処理システムの一つです。
IPython
で使われるノートブック形式ファイル(.ipynb
)です。%display plain
var('a b c x y t')
solve(a*x**2 + b*x+c == 0,x)
[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を使っています。
%display latex
solve(a*x**2 + b*x+c == 0,x)
数式の因数分解も、整数の素因数分解も簡単です。
factor(x**2 + (a + b)*x + a*b)
factor(x^3-1,algorithm="magma")
factor(2**31-1)
var('x y')
ex=(x+y)**4
show(ex)
ex=ex.expand()
show(ex)
show(ex.factor())
式を等価な単縦化された表現に変換します。
simplify
, simplify_full
simplify_factorial
, simplify_log
simplify_trig
, simplify_hypergeometric
simplify_rational
simplify_real
ex=(cos(x)+1)**4 + (cos(x)-1)**4
show(ex)
show(ex.expand())
show(ex.expand().trig_reduce())
show((x**6-1).factor())
(x**6-1).factor_list()
var('x y')
show((log(x)+log(y)).simplify_full())
show((log(x)+log(y)).simplify_log())
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"を持っている場合、
それを利用して結果を表示します。図形オブジェクトではその図形が表示されます。
for i in range(3,32,8):
show( factor(2**i-1))
%display latex
などのような%
で始まる行は(ipythonの)マジックコマンドです。%%bash
date
echo $EPICS_ROOT
2023年 6月28日 水曜日 16時36分34秒 JST
%display plain
%lsmagic
<IPython.core.magics.basic.MagicsDisplay object at 0x11a6b5310>
%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
var('x y')
solve([x**2+y**2 == 25, y-x==1], [x,y])
[[x == -4, y == -3], [x == 3, y == 4]]
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
は[SageMath
]の上に実装された微分可能多様体を取り扱うためのモジュールです。SageManifolds
はSageMath
とは独立したプロジェクトですが、現在はSageMath
の正式配布版に含まれており、SageMath
をインストールすれば直ちに利用可能となっています。%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
とします。
LC.<t,x,y,z> = M.chart()
show(LC)
多様体の計量$g$ を 導入します。
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)
別の座標系(極座標)を導入します。極座標ではそれぞれの座標の取りうる範囲が異なります。
(+oo
は+Infinity
の別の表記方法です)
XS.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
XS
座標系の軸は配列のように取り出せます。
XS[:]
二つの座標系の関係を変換関数を使って定義します。
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()
#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()
この座標系での計量の表示がどうなるかみてみましょう。
g.display(XS)
逆変換を求めてみましょう。
try:
XS_to_LC.inverse().dislay()
except Exception as m:
print(m)
no solution found; use set_inverse() to set the inverse manually
残念ながらこの場合、SageManifoldsは逆変換を求めることができませんでした.
別途計算した逆変換を設定します。
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.
(M.coord_change(LC,XS)*M.coord_change(XS,LC)).display()
(M.coord_change(XS,LC)*M.coord_change(LC,XS)).display()
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())
--------------------------------------------------------------------------- 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が計算してくれます。
MC.<t_o,x_o,y_o,z_o> = M.chart()
show(MC)
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())
運動系での計量を確認してみましょう。
g.display(MC)
(M.coord_change(LC,MC)*M.coord_change(MC,LC)).display()
M.atlas()
M.coord_change(XS,LC)
M.coord_change(MC,LC)
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$に定義されている座標系の一覧を与えます。
M.atlas()
MC_to_XS = M.coord_change(LC,XS) * M.coord_change(MC, LC)
show(MC_to_XS)
MC_to_XS.display()
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())
(M.coord_change(LC,MC1) * M.coord_change(MC, LC)).inverse().display()
M.coord_change(LC,MC1).display()
M.coord_change(MC,MC1).display()