pythonの強みの一つは世界中の開発者によって開発された豊富なライブラリ(モジュール)です。
数値処理には、基本的な機能を提供し、その他のモジュールのベースとして広く使われているnumpy
、グラフ生成のmatplotlib
,
それらの基盤の上に構築された高機能なpandas
などがよく使われます。また、機械学習の分野ではTensorFlow
,PyTorch
, Keras
などが知られています。
Pythonには数値計算だけでなく, Mathematica
やmaxima
のように数式を記号のまま取り扱うモジュール, sympy
も存在しています。
つまり、sympy
では、次の例のように(数値ではなく)数式そのものを取り扱うことができます。
sympyでの処理の例をみてみましょう。
from sympy import *
var('x')
expr=(x+1)*(x-1)
expr
数式を展開することもできます。
expr.expand()
逆に因数分解することも、
_.factor()
三角関数を含む式の変換だって、
expr=sin(x)*cos(x)
expr.simplify()
sympy
でできること¶sympy
には次のような機能があります(sympy チュートリアル)。
* 数式処理:数式の展開(expand), 整理(simplify)、因数分解(factor)
* 方程式の解を求める。 (solve)
* 式のままの微分、積分 (diff, integrate)
* 行列計算 (Matrix)
* プロット (plot, plot_implicit, plot_parametric,...)
* その他 (physics, control, ....)
というように数式そのものを計算対象として取り扱えます。
pythonをベースにした数式処理システムとしてSageMath
もよく知られています。
SageMath
は:
Jupyter
を使ったMathematica
ライクなNotebook
インタフェースを備えたシステムです。sage
コマンドでコマンドラインからの利用もできます。SageMath
はsympy
やmaxima
,R
など複数のライブラリを統合して、一つのシステムとしてまとめ上げています。 python
を基本としていますが、数式処理に合わせた拡張(変更)が行われた独自の言語となっています。今回は sympy
を使ったpythonでの数式処理だけをご紹介します。
sympy
を利用するpythonプログラムでは、sympy
をインポートしておきます。
import sympy
from sympy import var, solve, dsolve, plot, Matrix # from sympy import * でもよい
from sympy import plot, plot_implicit, plot_parametric
モジュールのバージョンは、そのモジュールの "__version__" で確認します。 (Python処理系のバージョンは、sys.version_infoなどで確認します。)
インターネット上のドキュメントを利用する際には、そのドキュメントが対象とするバージョンと利用中のモジュールのバージョン が一致していることを、確認しましょう。
print("sympy version:", sympy.__version__)
import sys
sympy.__version__, sys.version_info, sys.version
sympy version: 1.9
('1.9', sys.version_info(major=3, minor=9, micro=10, releaselevel='final', serial=0), '3.9.10 (v3.9.10:f2f3f53782, Jan 13 2022, 16:55:46) \n[Clang 13.0.0 (clang-1300.0.29.30)]')
sympy
を使ってみよう。¶sympy
のシンボル(数式中の変数を表す記号)¶sympy
で数式を操作する際には、数式の中に現れる記号(シンボル)の名前 (例えば、$x^2 - y^2$ の $x$ や $y$ ) を宣言しておく必要があります。 これには通常var
関数が使われます。var
関数の引数は数式で使われる変数名(Symbol
)を並べた文字列です。var()
を実行すると、文字列中の変数名のシンボルが作成され、同名のpythonの変数に割り当てられます。#数式で使用するシンボルを定義します。
from sympy import var,symbols
var('x y z t')
(x, y, z, t)
type(x),x.name
(sympy.core.symbol.Symbol, 'x')
var('a b c d')
print( a*b + c*d)
a*b + c*d # jupyter+sympy では出力をLaTex形式を使って表示します。
a*b + c*d
var
関数呼び出しは、<name> = symbols("<name>")
と等価です。
var('a')
asym=symbols('a')
aSym=Symbol('a')
a+asym+aSym
pythonの変数 a
には名前が 'a'であるSymbolオブジェクトが割り当てられていることを確認してみましょう。
type(a), a.name
(sympy.core.symbol.Symbol, 'a')
symbols()
はvar()
と同じく複数の記号名を登録できます。
Symbol()
は一つのシンボルだけを登録します。
a,b,c,d=symbols('a b c d')
abcd=Symbol('a b c d')
((a+a*b*c*d)**2+a*abcd).factor()
a
とSymbol
'a'の関係は固定されたものではないことに注意しましょう。a
はいつでも別のオブジェクトに再割り当てすることができます。var('a b c')
u=a; v=b
print(type(u),u.name)
u*(a+v)
<class 'sympy.core.symbol.Symbol'> a
sympy
のlatex
関数を使って、数式のLaTex表現の文字列を入手できます。from sympy import latex
var('a b c d')
print((sqrt(a*b + c*d)))
print(latex(sqrt(a*b + c*d)))
sqrt(a*b + c*d)
sqrt(a*b + c*d) \sqrt{a b + c d}
名前がギリシャ文字名のシンボルは、jupyter Notebookではギリシャ文字で表示されます。
$\lambda$(\lambda
)を使う際にはpythonの予約語lambda
との衝突を避けるため、lambda_
とするなどの注意が必要です。
var('x y z t u v varphi alpha phi ell rho lambda_ Alpha Beta Gamma_mu ') # lambda_は予約語lambdaとの衝突をさけるため
sin(varphi)**2 * alpha * cos(phi) + ell * tan(rho) * sqrt(lambda_) + Alpha*Beta*Gamma_mu
Symbol
関数を使うことで、変数名とシンボルとして印刷されるときの名前を変えておくこともできます。
lambda_=Symbol("\lambda")
lambda_**2
LaTeXあるいはMathJaxで結果を利用したい場合には latex()
関数でLaTeX表記の文字列を入手しましょう。
print(latex(sin(varphi)**2 * alpha * cos(phi) + ell * tan(rho) * sqrt(lambda_) ))
\sqrt{\lambda} \ell \tan{\left(\rho \right)} + \alpha \sin^{2}{\left(\varphi \right)} \cos{\left(\phi \right)}
余談::
ギリシャ文字あるいはアルファベットの一文字を名前とするシンボルはsympy.abc
モジュールに定義済みです。 sympy.abc._clash
に含まれる名前は,sympy
で定義されている変数とsympy.abc
で定義されている名前で競合が発生するので、注意が必要です。
import sympy.abc
sympy.abc._clash.keys()
dict_keys(['O', 'Q', 'N', 'I', 'E', 'S', 'beta', 'zeta', 'gamma', 'pi'])
expand
)、整理(simplify
)、因数分解(factor
)¶記号として $x$, $y$ を使った数式を考えてみます。
var('x y')
expr=(x**2 + x*y + y**2)*(x-y)
expr
このように、sympyでは数式自体が計算結果となります。
type()
関数で調べてみると、expr
の中身は数値ではなく、sympy
のオブジェクトであることがわかります(ここでは、sympy.core.mul.Mul
オブジェクトなっています)。
type(expr)
sympy.core.mul.Mul
この数式を展開(.expand()
)してみましょう。
expr=expr.expand()
expr
print(type(expr))
<class 'sympy.core.add.Add'>
この式はsympy.core.add.Add
オブジェクトになっています。多項式は、積(sympy.core.mul.Mul
)で作られた単項を足し合わせ(sympy.core.add.Add
)て作られていることがわかります。
此の式を因数分解(.factor()
)してみると、
expr=expr.factor()
print(type(expr))
expr
<class 'sympy.core.mul.Mul'>
となります。なお、factor
, expand
は関数としても利用できます。
expand(expr)
factor(expand(expr))
式の簡素化に.simplify()
メソッドが用意されています。
expr.simplify()
式の変数を別の値で置き換えた式をつくるには、.subs
メソッドを使います。
式中の変数$x$を数値$123$に置き換えてみましょう。
expr= x**2 + 3*x + 2
expr.subs(x, 123)
式中の変数$x$を別の数式で置き換えることも可能です。式に現れる一部の項だけを置き換えることも可能です。
expr.subs({x:2+y}) #置換ルールを辞書型データで指定しています。
expr.subs({x:2+y, y:4})
expr.subs({x**2:y,y:(x+1)**3})
expr.subs(3*x+2, y )
三角関数を含む数式の整理には、trigsimp()
やexpand_trig()
関数を使うことができます。
expr=cos(2*x) + sin(2*x)
expr
三角関数専用の簡素化を行う関数trigsimp()
や展開関数expand_trig
が用意されています。
trigsimp(expr)
expand_trig(expr)
expr=cos(2*x) + 2*sin(x)*cos(x)
expr.simplify()
expand_trig(expr)
expand_trig(expr).simplify()
expr.rewrite(cos)
expr.as_terms()
([(2*sin(x)*cos(x), ((2.0, 0.0), (1, 1, 0), ())), (cos(2*x), ((1.0, 0.0), (0, 0, 1), ()))], [sin(x), cos(x), cos(2*x)])
2次方程式 $a x^2 + b x + c = 0$ の解をsympy
を使って求めてみましょう。
まず、$x$についての2次の多項式を作ります。
var('a b c x')
expr=a*x**2 + b* x + c
expr
solve()
関数は与えらてた数式を0
とする解があれば、それらの解のリストを返します。
sols=solve(expr, x)
sols
[(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]
Eq()
クラスを使って、方程式 $expr = 0$ を明示的に表現することもできます。
Eq(expr,0)
solve(Eq(expr,0),x)
[(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]
解を辞書型や集合型のデータとして求めることも可能です。
sols_dict=solve(expr, x, dict=True)
sols_dict
[{x: (-b + sqrt(-4*a*c + b**2))/(2*a)}, {x: -(b + sqrt(-4*a*c + b**2))/(2*a)}]
sols_set=solve(expr, x, set=True)
sols_set
([x], {((-b + sqrt(-4*a*c + b**2))/(2*a),), (-(b + sqrt(-4*a*c + b**2))/(2*a),)})
解(sols
)の一つを、方程式に代入(substitute)してみます。
expr.subs(x,sols[0])
式を整理してみると、方程式を満足していることがわかります。
expr.subs(x,sols[1]).simplify()
辞書型で求めた解の一つを代入する時には、次のように書きます。
expr.subs(sols_dict[0])
expr.subs(sols_dict[0]).simplify()
二つの解がどちらも方程式を満たしていることを,リスト内包表記 [<expr> for i in <list>]
をつかって 確認してみます。
[expr.subs(sol).simplify() == 0 for sol in sols_dict]
[True, True]
等式(Eq
)を使った場合はこちら。
[Eq(expr, 0).subs(sol).simplify() for sol in sols_dict]
[True, True]
set型の戻り値を使った場合は、このようになるでしょう。
[expr.subs(sols_set[0][0],sol[0]).simplify() == 0 for sol in sols_set[1]]
[True, True]
sympy
では var()
やsymbols()
に与える文字列の中にスライス記法をつかって一連の変数名をまとめて宣言できます。
Symbolとしての名前は、LaTeXのルールに従うので、"_"を含む名前は下付き文字を使って表現されます。
var("a_:3 b_:2(:2)") # sympyのsymbols/varでは複数のシンボルをスライス記法を用いて指示できます。
(a_0, a_1, a_2, b_00, b_01, b_10, b_11)
print(latex(a_0+a_1*x+ a_2*x**2))
a_0+a_1*x+ a_2*x**2 +b_11 # jupyter環境では`_`は下付き指標に整形されます。
a_{0} + a_{1} x + a_{2} x^{2}
LaTeXで上付き文字を指定する ^
はpythonの変数名の中で使うことができません。この場合には、sympy()
のsymbols()
を使うことになります。 (c^0
などを直接 pythonの変数名とすることはできません)
c0,c1,c2,d00,d01,d10,d11=symbols("c^:3 d^:2(:2)") #
c0.name
'c^0'
a_0*c0+a_1*c1
Matirx
)¶Matrix
オブジェクトを作ることで、行列計算も行えます。
行列オブジェクトの和、積を求めてみましょう。
var("a_:2(:2) b_:2(:2)")
m=Matrix(
((a_00,a_01),
(a_10,a_11)
)
)
n=Matrix(((b_00,b_01),(b_10,b_11)))
m+n
m*n
m.det() #行列式
m.transpose() # 転置行列
(m**(-1)).subs(det(m), 1) # 行列式が1の行列の逆行列
simplify(m*m**(-1))
m.det()
m.trace()
各種の共役行列も求めることができます。
m.conjugate() #複素共役
m.adjoint() #随伴行列/エルミート共役
simplify(m* m.adjoint())
sympyの数式が表す数値の近似値を.evalf()
メソッドを用いて求めることができます。引数で、有効桁数を指定します。
.evalf()
メソッドの短縮形として、.n()
メソッドを使うこともできます。N(expr, n)
関数はexpr.evalf(n)
と等価です。
sin(1).evalf(10), sin(1).evalf(20), sin(1).evalf(30)
(0.8414709848, 0.84147098480789650665, 0.841470984807896506652502321630)
cos(1).n(20),tan(1).n(40)
(0.54030230586813971740, 1.557407724654902230506974807458360173087)
pi.n(10), E.evalf(20) # piは円周率, E=exp(1)=E1
(3.141592654, 2.7182818284590452354)
Rational
クラス)¶sympy
はPythonの文法に従いますので、1/3
などの式はfloat
として評価されてしまいます。
有理数として1/3
を取り扱うためのRational
クラスが用意されています。また、Singleton リポジトリ関数S()
を使って、
有理数を作りだすこともできます。
1/3, Rational(1,3), 1/3 == Rational(1,3)
(0.3333333333333333, 1/3, False)
Rational(1,3) + Rational(1/2) + Rational(1,3)/Rational(1/2)
Integer(1)/Integer(3) == Rational(1,3)
True
$\lim_{x\rightarrow \infty} f(x)$のような極限値を求めることも可能です。
expr=(1+1/x)**x
expr.limit(x,oo) #sympy.core.numbers.Exp1
sin(x).nseries(x,0,n=10)
無限大を表す記号として、oo
を使います。oo
はsympy.core.numbers.Infinity
に割り当てられています。
数式(多項式)の因数分解だけでなく、整数の素因数分解(factorint
)の機能も使えます。
factorint(5120,verbose=True,multiple=False)
Factoring 5120 Trial division with ints [2 ... 32768] and fail_max=600 2 ** 10 Factorization is complete.
{2: 10, 5: 1}
factorint(5120, visual=True)
import random
ri=random.randint(1,100000)
print(ri)
factorint(ri)
15187
{15187: 1}
gcd(447550424579,304874938368000)
lcm(447_550_424_579, 304874_938_368_000) == 447550424579*304874938368000
True
isprime(447550424579), isprime(2850639647)
(False, True)
factorint(447550424579), factorint(304874938368000)
({157: 1, 2850639647: 1}, {2: 16, 3: 7, 5: 3, 7: 1, 11: 1, 13: 1, 17: 1})
鶴亀(蟻)算を例に、連立一次方程式を解いてみます。 問題は:
という時の鶴、亀、蟻の数を求めよ
というものです。eqs
の各要素が0になる解を求めます。
var('鶴 亀 蟻')
eqs=[ 2*鶴+4*亀+6*蟻 - 34, 鶴+亀+蟻 - 10, 蟻 + 1 - 亀 ]
sol=solve( eqs, (鶴, 亀, 蟻))
sol
{鶴: 5, 亀: 3, 蟻: 2}
Eq()
を使う¶等式を表現するEq()
関数を使って、解くべき方程式を表現することも可能です。
Eq(2*鶴+4*亀+6*蟻 , 34)
eqs=[Eq( 2*鶴 + 4*亀 + 6*蟻 , 34), #足が全部で34本
Eq( 鶴 + 亀 + 蟻 , 10), #全部で10体
Eq( 蟻 + 1 , 亀) ] # 亀は蟻より1多い
sol=solve( eqs, (鶴, 亀, 蟻))
sol
{鶴: 5, 亀: 3, 蟻: 2}
解sol
を方程式に代入して、確かに解になっていることを確認します。
Eq(2*鶴 + 4*亀 + 6*蟻 , 34).subs(sol)
[eq.subs(sol) for eq in eqs]
[True, True, True]
expr=x**2+ sin(x)**3
expr.diff(x)
diff()
関数、Derivative
クラスも使えます
diff(expr,x)
Derivative
クラスは式に対する微分の操作を表現します。
Derivative(expr,x)
.doit
メソッドで微分操作を実行します。
Derivative(expr,x).doit()
.integrate()
メソッドで積分を実行します。
expr=(x**2-2*x-2)/(x**3-1)
expr.integrate(x)
integrate()
関数をつかうこともできます。
integrate(expr,x)
被積分関数expr
を部分分数展開すると、この積分が正しいことは一目瞭然です。
expr.apart()
定積分には、変数範囲をタプルで指定します( (<var>,<lower range>,<upper range>)
。sympy
ではoo
は無限大を表します。
integrate(exp(-x**2), (x, -oo, oo) )
(exp(-x**2)).integrate((x,-oo,oo))
Integral
クラスは積分の数式を表現します。(Eq
が等式を表現するクラスであったのと同じです)
expr=Integral(exp(-x**2), (x, -oo, oo))
expr
積分クラスの積分を実行した結果を求めるには.doit()
メソッドを使います。
expr=Integral(exp(-x**2), (x, -oo, oo))
expr.doit()
Function
クラスは関数を表すオブジェクトを作成します。
f=Function('f') # `f`は関数を表すオブジェクト
f(x)
f(x).diff(x)
f(x,y).diff(x,y,x)
f(x,y).diff(x,2,y) #これも可能
# n階の微分をあ羅和すには?
var('n m',integer=True)
f(x,y).diff(x,n,y,m)
f=Function('f') # `f`は関数を表すオブジェクト
g=Function('g')
f(g(t)).diff(t)
( f(x)**2 ).diff(f(x)) # `f(x)**2.diff(f(x))` は ❌
( f(x)**n ).diff(x).simplify() # `f(x)**2.diff(x)` は ❌
dsolve
関数で微分方程式を解くこともできます。
f=Function('f')
g=Symbol('g')
m=Symbol('m')
dsolve(Eq(m*f(t).diff(t,2), -g) , f(t))
$C_1, C_2$は初期条件で決まる定数を表現しています。初期条件$\left( f(0) = x_0, \frac{d f}{d t}(0) = v_0\right)$ を与えてみましょう。
var('x_0 v_0')
dsolve(Eq(m*f(t).diff(t,2), -g) , f(t), ics={f(0):x_0,f(t).diff(t).subs(t,0):v_0})
var('omega', positive=True) #パラメータ`omega`は正の数であることを仮定します。
eq=Eq(Derivative(f(t), t, t) , -(omega)**2*f(t) )
dsolve(eq, f(t))
var('omega',complex=True) # `omegaは複素数と仮定します
eq=Eq(Derivative(f(t), t, 2) , -(omega)**2*f(t) )
dsolve(eq, f(t))
var('omega',complex=True) # `omegaは複素数と仮定します
eq=Eq(f(t).diff(t,2) , -(omega)**2*f(t) )
dsolve(eq, f(t))
2階の微分方程式を、1階の連立微分方程式に書き直して解いてみましょう。
var('omega')
var('f g', cls=Function)
dsolve([Eq(diff(g(t),t), - omega**2* f(t)),
Eq(diff(f(t),t), g(t))], [f(t),g(t)])
[Eq(f(t), I*C1*exp(-I*omega*t)/omega - I*C2*exp(I*omega*t)/omega), Eq(g(t), C1*exp(-I*omega*t) + C2*exp(I*omega*t))]
var('omega',real=True) # aは実数であるとの条件をつける。
g=Function('g')
dsolve([Eq(diff(g(t),t), - omega**2*f(t)),
Eq(diff(f(t),t), g(t))], [f(t),g(t)])
[Eq(f(t), C1*sin(omega*t)/omega + C2*cos(omega*t)/omega), Eq(g(t), C1*cos(omega*t) - C2*sin(omega*t))]
matplotlib
がインストールされた環境では、sympy
はmatplotlib
を利用してグラフを描画します。
matplotlib
が無い場合でも、文字を使ってグラフを表示します。
var('x y')
P=plot(sin(x),(x,-3*pi,3*pi),line_color="red")
先ほどのグラフに次のグラフを合成したグラフを表示してみましょう。
Q=plot(cos(x),(x,-3*pi,3*pi),line_color="green")
P.extend(Q)
P.show()
陰関数で定義される曲線を印刷することもできます。
Q=plot_implicit(x**2+y**2-1, (x,-1,1),(y,-1,1),line_color="purple")
P=plot_parametric((cos(x),sin(x)), (x,-pi,pi),line_color="red")
作成したグラフを保存します。
P.save('foo.png')
関数のTaylor級数展開を求めることも可能です。
sin(x).series(x,0,8) # x=0の周りで8次までの級数展開を求める。
異なる点($x=\frac{\pi}{2}$)の回りで展開してみましょう。
sin(x).series(x,x0=pi/2,n=8)
級数どうしの積を求めてみます。
(sin(x).series(x,0,8)*cos(x).series(x,0,10)).expand()
展開前の関数の積を展開しても同じ級数となります。
(sin(x)*cos(x)).series(x,0,8)
(sin(2*x)/2).series(x,0,8)
無限級数の和を計算してみましょう。oo
は無限大を表しています(sympy.core.numbers.Infinity
)。
var('n x')
f=lambda n:1/n**n
S=Sum(f(n), (n, 1, oo))
S.evalf()
有理関数を部分分数展開してみます。
var('z')
expr=1/(z**2 + 3*z +2)
apart(expr)
.apart()
メソッドも使えます。
expr.apart()
この様な展開はLaplace逆変換を求める際に便利です。
sympy.integrals.transforms.inverse_laplace_transform(expr, z, t)
var('omega omega_0')
sympy.integrals.transforms.inverse_fourier_transform(1/(omega-omega_0), omega,t)
式の中の特定の変数について整理するときは.collect()
を使います。その変数について特定の次数の係数を.coeff()
メソッド
で取り出します。
expr=x**2-2*x+a*x*y - (y **3 - 3*x*y**2)
expr.simplify()
expr.collect(x)
collect(expr,y)
collect(expr,(y,x)).simplify().collect(x)
expr.coeff(x,1)
import sympy.physics as Physics
if sympy.__version__ > "1.9":
import sympy.pysics.control as control
これでおしまい。
以下はメモ書き
(classes_sympy.png)
次の連立方程式。 $ \left\lbrace \begin{eqnarray*} y^{2} = x^{3} - 3 x^{2} + 2 x, \\ y^{3} - 3 y^{2} = x^{2} - 2 x \\ \end{eqnarray*} \right. $を解いてみる。
まずは、solve
を使ってみる。
と解はでるが、どうやって解いたのかは解らないので一歩づつ進めてみる。 まず、第一式を第2式に代入し、yについて解いてみる。 方程式を因数分解してみやすい形にしておく。
eqs=factor(eqs)
eqs[0]
eqs[1]
第一式を第二式に代入する。
factor(eqs[1].subs(eqs[0].lhs,eqs[0].rhs))
これを解くと次のようになる。
sols_y=solve(factor(eqs[1].subs(eqs[0].lhs,eqs[0].rhs)),(x,y),dict=True)
sols_y
これらの解を第一式に代入すると、
[eqs[0].subs(sol) for sol in sols_y]
[eqs[1].subs(sol) for sol in sols_y]
それぞれの場合に、これらの方程式の解を求める。
[solve(eqs[0].subs(sol),(x,y),dict=True) for sol in sols_y]
$x=2$ および $x=2$ の場合はいずれも$y=0$が解となることがわかる。 $y=\frac{3 x -2}{x-1}$の場合には、$x$として五個の解がある。これをもう少し調べてみる。 $x\ne0$かつ$x\ne2$の時の$x,y$の関係を求める。
これを方程式の第1式に代入して、整理してみる。
eq0=eqs[0].subs(sols_y[-1]).simplify()
expr0=eq0.lhs-eq0.rhs
num,den=expr0.as_numer_denom()
expr0 *=den
expr0=expr0.factor()
print(latex(expr0))
expr0
第2式に代入しても同じ式にたどり着く。
eq1=eqs[1].subs(sols_y[-1]).simplify()
expr1=eq1.lhs-eq1.rhs
num,den=expr1.as_numer_denom()
expr1 *=den
expr1=expr1.factor()
expr1
$\left(x^{2} - 4 x + 2\right) = 0$ および $\left(x^{3} - x^{2} + 3 x - 2\right)=0$ をそれぞれ解いて、 $y$を$y=\frac{3 x -2}{x-1}$で求めれば、全ての解が確定する。
.args
をつかって、因数の各項を取り出せる。これらが0となる解をそれぞれ求める。
expr1.args
var('x')
sol1=solve(Eq(expr1.args[0],0),x)
sol1
これらの解を、$y=(3 x -2)/(x-1)$に代入して、元の方程式の解 $ x, y$を求める。
sol1_y=[{x:x1, y:((3*x1-2)/(x1-1)).simplify()} for x1 in sol1]
sol1_y
expr1
の因数の後半から、
var('x')
sol2=solve(Eq(expr1.args[1],0),x)
sol2
$x,\,y$のペアを作成
sol2_y=[{x:x1, y:((3*x1-2)/(x1-1)).simplify()} for x1 in sol2]
sol2_y
solutions=sol1_y+sol2_y
solutions
解が方程式を満たしていることを確認する。
check0=eqs[0].subs(sols_y[2]).simplify()
check0
[check0.subs(r, simultaneous=True).simplify() for r in solutions]
check1=eqs[1].subs(sols_y[2]).simplify()
check1
[check1.subs(r, simultaneous=True).simplify() for r in solutions]
式$x^3 - x^2 + 3 x - 2 = 0$は$ x\rightarrow x+\frac{1}{3}$を行うことで、 カルダノの公式の標準形に変換できる。
(x**3 - x**2 + 3*x - 2).subs(x,x+Rational(1,3)).expand().collect(x)
solutions[-1][x].expand().simplify()
p=Rational(8,3);q=-Rational(29,27)
u3=-q/2+sqrt((q/2)**2+(p/3)**3)
v3=-q/2-sqrt((q/2)**2+(p/3)**3)
(-q/2+sqrt((q/2)**2+(p/3)**3)).simplify()
(u3+v3).expand().simplify(),factor((u3*v3-(-p/3)**3).expand().simplify())
(u3**(Rational(1,3)).simplify() + ((-p/3)*u3**Rational(-1,3)).simplify())
factorint(512),factorint(729)
solutions[-1][x].as_real_imag()
var('omega')
sols[2][x].as_real_imag()
これらの解は数値的には0と見做せる小さな違いを除いて一致している。
N(sols[2][x]-solutions[-1][x])
N(sols[3][x]-solutions[-2][x])
N(sols[4][x]-solutions[2][x])
factorint(14331349955611230)
factorint(799890941846430)
((29+3*sqrt(321))*(3*sqrt(321)-29)).expand().simplify()
160646206240*3/8953569504
omega=exp(2*I*pi/3)
omega**2
(1+omega+omega**2).simplify()
expr=(x+y+z)*(x+omega*y+omega**2*z)*(x+omega**2*y+omega*z)
print(expr.expand().simplify())
print(expr.expand().simplify().collect(x))
(x**3 + y**3 + z**3 - 3*x*y*z).factor()
var('omega x y z p q')
expr=(x+omega*y+z*omega**2)*(x+y*omega**2+omega*z);
expr.expand().subs({omega**4:omega,omega**3:1,}).subs({omega**2:-1-omega}).simplify()
solve([Eq((-p/3)**3,y**3*z**3),Eq(q,y**3+z**3)],(y**3,z**3),dict=True)
solve(Eq(y**3+(-p/3)**3/y**3,q),y**3,dict=True)
$y^3= \left(q/2+\sqrt{12 p^3+81 q^2}/18\right)$は1の複素3乗根を$\omega$として、 $y = y_0, \omega y_0, \omega^2 y_0$ ただし、$y_0=(q/2-\sqrt(12 p^3+81 q^2)/18)^{1/3}$. 此の時、$z$は$z=(-p/3)/y_0, (-p/3)/y_0 \omega^2, (-p/3)/y_0 \omega$ということになる。
$z_0= (-p/3)/y_0=
y_0=(q/2-sqrt(12*p**3+81*q**2)/18)**(Integer(1)/Integer(3))
z_0=(q/2+sqrt(12*p**3+81*q**2)/18)**(Integer(1)/Integer(3))
((y_0**3*z_0**3).expand())**(Integer(1)/Integer(3))
help(z_0.expand)
(z_0/(q/2+sqrt(12*p**3+81*q**2)/18)**(Integer(1)/Integer(3))).expand(deep=True)
expr.expand().simplify()
from sympy import *
var('x k y')
expr=(k+1)*x**2+k*x-2*k
sols=solve(expr,x,dict=True)
solve(Eq( x.subs(sols[0]), x.subs(sols[1])),k,dict=True)
from sympy import *
var('x k y')
expr=(k+1)*x**2+k*x-2*k
sols=solve(expr, x)
solve(Eq(sols[0], sols[1]), k, dict=True)
((29 +3*sqrt(321))**5).expand().simplify()
factorint(23631496),factorint(108537648)
factorint(3730)
((29 +3*sqrt(321))**3).expand().simplify()
factorint(275732), factorint(16236)