python sympyを使って方程式を解く

Question

前提・実現したいこと

pythonを使って方程式を解きたいです。
sympyモジュールを使って次の方程式を解きます:

a * x**3 - b*x - c = 0

係数a, b, cは動的に与えられるものとしたいのですが、どう書くのがルールなのでしょうか。

ちなみに、係数があらかじめ決まっている場合は次のように書くことはわかっています。

x = Symbol('x', positive = True)
f = 4*x**3 -2 *x - 1
x0 = solve(f, x)

よろしくおねがいします。

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

エラーメッセージ

該当のソースコード

ソースコード

試したこと

ここに問題に対して試したことを記載してください。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Answer 1

次の二つのコードの動作は同じです。

x = Symbol('x', positive=True)
f = 2*x

x = Symbol('x', positive=True)
a = 2
f = a*x

リテラルを使っても変数に束縛された値を使っても、プログラム上は同じことです。

Comments

[hamberger]

###
x = Symbol('x', positive = True)
f = a*x**3 -b *x - c
x0 = solve(f, x)

ご回答ありがとうございます。
a, b,c が###の場所で定義されている場合に、f=0を上記のコードで解きたいのですが、（実際に実行すると解が出てこないです。定義場所を変えたくないのですが・・・）

[LouiS0616]

どのように定義されているか書いて下さい。
また、『解が出てこない』という意味が良く分かりません。

[hamberger]

具体的には次のように方程式の解を得たいのですが。

import sympy as sp
a=1.0
b=-1.0
H0 = 0.01
h0 = 1.0

x = sp.Symbol('x', positive = True)
x0 = sp.solve(b * x**3 + a * x - H0 / h0, x)
print(x0)

実行すると[]としか返ってきません。

[LouiS0616]

原因を特定するために、シンボルxの値に制限を付けずに解いてみます。
>>> x = sp.Symbol('x')
>>> sp.solve(b * x**3 + a * x - H0 / h0, x)
[-1.004962991944 + 0.e-23*I, 0.0100010003001201 + 0.e-20*I, 0.994961991643881 + 0.e-23*I]

以上のように、虚数解が導出されています。
しかし虚数部は限りなく0に近いようなので、これは解を求める際の誤差でしょう。

試行錯誤の末、次のようにキーワード引数cubicsをFalseにして実数解を求められました。
>>> sp.solve(b * x**3 + a * x - H0 / h0, x, cubics=False)
[-1.00496299194400, 0.0100010003001201, 0.994961991643881]

三次方程式には解の公式が存在するので、cubics=Trueだとそれを適用するのでしょう。多分。
一方cubicsをFalseにすれば、解析的に解が求められるのだと予想されます。

---
さて、次のようにして正の解を求めることができました。
>>> x = sp.Symbol('x', positive=True)
>>> sp.solve(b * x**3 + a * x - H0 / h0, x, cubics=False)
[0.0100010003001201, 0.994961991643881]

これでいかがでしょう。

[hamberger]

ただ、実際には正の実数解が一つしかないはずなのですが.. 精度が悪いのでしょうか。

[LouiS0616]

-x³+x-0.01 = 0 の解ですよね？

-x³+x = 0 の解が -1, 0, 1 であることを考えると、前述の解もその近傍にあることが予想されます。
問題は0から正負どちらの方向に振れるかです。

先の式を少し変形して、x(-x²+1) = 0.01 とします。
このときxは0に近い数ですから、-x²+1 は正です。x(-x²+1) 全体が正(0.01)になることを考えると、x自体も正の数になるのでは無いでしょうか。

[hamberger]

x^3 - x = 0.01　です！
おっしゃる通りにcubics =falseのオプションを付与することで正の実数解の表示に成功しました。

[0.490673329304963]

大変助かりました。お世話になりました！

[LouiS0616]

2019/09/30 00:37 のコメントに x0 = sp.solve(b * x**3 + a * x - H0 / h0, x) と書いてありますが、これが誤りだったと言うことですね。
確かに質問文だと係数は a, b の順になっていますね。

[hamberger]

はい、なのでLouis様の解釈は正しいです。次からは表記に気を付けます。ありがとうございました！