3次方程式の解を求めたとき、生じないはずの現象が起こるので解決したいです。

Question

前提・実現したいこと

Pythonを始めたばかりで、初歩的な質問かもしれません。

(x,a,b,c)を変数として3次方程式x^{3}-ax^{2}-bx-c=0をxについて解いたとき、虚数単位Iが含まれないか解βが1つ現れます。このβに対し、後で(a,b,c)に具体的な整数値を代入したとき、特定の(a,b,c)に対しては何故か虚数単位Iが出現します(例えば(a,b,c)=(4,-1-1)など)。

もしこの現象がsympyではしばしば起こってしまうものでしたら、あきらめようかと思いますが、そうでなければ、こうした現象が起こらないようにしたいです。ソースコードに問題があるのなら、ご指摘いただきたいと思います。

追記1：『生じないはずの現象』ですが、(a,b,c)=(4,-1,-1)のような実数解しか持たない場合でも常に虚数単位Iが出現します。この点についての質問でした。いくつか回答を頂きましたが、私の書き方が悪いために意図した内容とは違う回答でしたので、この場で謝罪いたします。申し訳ありません。

追記2：『生じないはずの現象』について追記です。実は(a,b,c)=(4,-1,-1)は1より大きい解と絶対値が1未満の解を2つもちます(Pisot数と呼ばれる代数的整数の性質です)。そこで(a,b,c)=(4,-1,-1)として

sympy.solve(sympy.Eq(f,0),x)

で解を出力すると、すべて虚数単位を有する解を出力します。

発生している問題・エラーメッセージ

エラーは発生していません。

該当のソースコード

import sympy
import math
import numpy

print("3次の代数的整数βの最小多項式をx^3-ax^2-bx-cとする.")

sympy.var("x,a,b,c")
f=x**3-a*x**2-b*x-c

sol=sympy.solve(sympy.Eq(f,0),x)

print("このとき, x^3-ax^2-bx-c=0の解は")
print("x1=",sympy.simplify(sol[0]))
print("x2=",sol[1])
print("x3=",sol[2])

print("今, a,b,cを次のように定める.")
a1=int(input("aを入力せよ."))
b1=int(input("bを入力せよ."))
c1=int(input("cを入力せよ."))

s0=sol[0].subs([(a,a1),(b,b1),(c,c1)])
s1=sol[1].subs([(a,a1),(b,b1),(c,c1)])
s2=sol[2].subs([(a,a1),(b,b1),(c,c1)])

beta=s0

print("このとき, β=",beta)

Comments

[can110]

「**」がマークアップ文字と誤認されているためか
「NameError: name 'x3' is not defined」が発生します。
コードはコードブロック(```)で囲んでください。

[TTATK]

ご指摘ありがとうございます。ソースコードの入力には、コード挿入用のツールがあるのですね。x3はもともとx**3と記載していたのですが、何故か'x3'で認識されたので、このサイトの仕様かと勘違いしていました。

Answer 1

wolframalpha
計算すると実数値なんでしょう。

近似値が解の一つとちゃんと一致してますし。

Comments

[TTATK]

一瞬？と思いましたが、なるほど、きちんと計算すると虚数単位が消えるということですね。そこは盲点でした。回答ありがとうございます。他の方法も模索しつつ、虚数単位が消える方向にも考えてみようかと思います。

Answer 2

解

((2*a - 2**(2/3)*(-2*a**3 - 9*a*b - 27*c + sqrt(-4*(a**2 + 3*b)**3 + (2*a**3 + 9*a*b + 27*c)**2))**(1/3))*(-2*a**3 - 9*a*b - 27*c + sqrt(-4*(a**2 + 3*b)**3 + (2*a**3 + 9*a*b + 27*c)**2))**(1/3) - 2*2**(1/3)*(a**2 + 3*b))/(6*(-2*a**3 - 9*a*b - 27*c + sqrt(-4*(a**2 + 3*b)**3 + (2*a**3 + 9*a*b + 27*c)**2))**(1/3))

には虚数部は含まれませんがa=4, b=-1, c=-1において

sqrt(-4*(a**2 + 3*b)**3 + (2*a**3 + 9*a*b + 27*c)**2)

の中の

-4*(a**2 + 3*b)**3 + (2*a**3 + 9*a*b + 27*c)**2

は-4563になるので、結果として虚数部を含みます。

Comments

[TTATK]

わざわざ計算までしていただき、回答して下さりありがとうございます。追記させていただきましたが、実は生じないはずの現象とは『実数解しか持たないはずなのに虚数単位が出現する』点でした。お手数をおかけしまして、大変申し訳ありません。もしこの現象について心当たりがありましたら、再度ご回答いただけたら幸いです。

Answer 3

sol[0] の解に平方根が含まれているので、a, b, c の値によっては平方根の中身が負になって、その結果、虚数単位が表れます。

x ** 3 - a * x ** 2 - b * x - c = 0, a=4,b=-1, c=-1 - Wolfram|Alpha

追記

evalf() で数値解として評価すると、虚数解が微小ですがでてきますね。

sols = sympy.solve(sympy.Eq(f, 0), x)

print("このとき, x^3-ax^2-bx-c=0の解は")
print("x1=", sols[0])
print("x2=", sols[1])
print("x3=", sols[2])

print("今, a,b,cを次のように定める.")
a1 = 4
b1 = -1
c1 = -1

for sol in sols:
    print(sol.subs([(a, a1), (b, b1), (c, c1)]).evalf())

# 0.726109445035782 + 0.e-20*I
# -0.377202853972958 + 0.e-22*I
# 3.65109340893718 - 0.e-22*I

Wolframe Alpha で a=4,b=c=-1 のときの厳密解は以下のようになっています。

Comments

[TTATK]

回答していただきありがとうございます。おっしゃる通り解の中に平方根が含まれているので、虚数単位が現れることはあるかもとは思ったのですが、(a,b,c)=(4,-1,-1)の場合は実数解しか持たず虚数単位が現れることは本来あり得ない現象なのです。お手数をおかけしまして大変申し訳ありませんが、こういった現象にもし心当たりがありましたら、ご回答いただけたら幸いです。

[tiitoi]

前提となる数式ですが
x^3-ax^-bx-c=0
この -bx は x^{-bx} ですか？
累乗がどこにかかっているかがわかりづらいので、tex 記法で記載していただけますか？

[tiitoi]

^2 が抜けてただけっぽいですね
失礼しました

[TTATK]

ご指摘ありがとうございます。質問文にいろいろと抜けがあり大変申し訳ありません。今、訂正を致しました。その他、不明な点等ありましたら遠慮なく申し付けてください。

[tiitoi]

微小な値ですが、解析解に虚部が含まれているようですね。グラフを描画すると x 軸と3点とも交わっているようにも見えるので、なんで虚数解が出てくるのか謎ですね
ただ、wolframe alpha など他の数式処理ソフトでも同様の結果になったので、なにか理由があるのだと思います。

[tiitoi]

最終的にほしいのが数値解なのであれば、scipy などのモジュールを使えば数値解として得られると思います。

[TTATK]

細かいところまで計算していただき、誠にありがとうございます。

他の数式処理でも同様でしたか。ということは、もしかしたら数式処理全般に起こり得る仕様かもしれませんね。私の研究では、この段階の数値解の誤差はちょっと許容できないのでどうにかならないかと足搔いていたのですが……

また別の方法を考えてみようかと思います。付き合って下さり誠にありがとうございました。