Sympy 解决 returns 一个空集

Question

我正在尝试使用 solve 求解以下多项式方程但是它 returns 是一个空集。

    `3*d*x**2 + 2*w*x = 0,` 
    `3*d*y**2 + 2*w*y = 0,`
    `-d + 2*w*z + 1 =0,`
    `x**3 + y**3 - z = 0,`
    `x**2 + y**2 + z**2 - 1 = 0`

wolfram-alpha returns 遵循真正的解决方案。

    `d≈0.417238, w≈-0.516977, x = 0, y≈0.826031, z≈0.563624`
    `d≈0.417238, w≈-0.516977, x≈0.826031, y = 0, z≈0.563624`
    `d≈0.417238, w≈0.516977, x = 0, y≈-0.826031, z≈-0.563624`
    `d≈0.417238, w≈0.516977, x≈-0.826031, y = 0, z≈-0.563624`
    `d≈0.528689, w≈-0.492358, x≈0.620853, y≈0.620853, z≈0.478626`

如果存在解决方案，我如何确保获得解决方案？有没有其他方法可以更稳健地处理 0 和其他条件？

Answer 1

嗯，sympy 的 solve 寻找完全符号化的解决方案。当前版本 returns 如果方程太难以符号方式求解（不输出错误或警告），则为空列表。

要获得数值解，sympy 有 nsolve 需要初始猜测，并且（当前）只有 returns 一个解。

from sympy import nsolve, symbols

d, w, x, y, z = symbols('d w x y z', real=True)

nsolve([3 * d * x ** 2 + 2 * w * x,
        3 * d * y ** 2 + 2 * w * y,
        -d + 2 * w * z + 1,
        x ** 3 + y ** 3 - z,
        x ** 2 + y ** 2 + z ** 2 - 1], (d, w, x, y, z), (1, 1, 1, 1, 1))

报告的解决方案：

Matrix([[ 0.528689459190352],
        [-0.492357687731282],
        [ 0.620853041008598],
        [ 0.620853041008598],
        [ 0.478626161989451]])

Answer 2

如果您将 SymPy 用作纸笔增强功能并像手动求解方程那样处理方程，您会发现整个系统可以简化为一个单一的多项式表达式（在这种情况）有一个明确的解决方案或可以通过数值求解；然后可以将该单个变量的解反代到其他方程中。有一些 hints SymPy 有一天会更自动地解决这样的系统，但我希望自己看到这些关系会有一些快乐，就像我对以下内容所做的那样：

定义方程

>>> from sympy.abc import *
>>> eqs = (3*d*x**2 + 2*w*x,
... 3*d*y**2 + 2*w*y ,
... -d + 2*w*z + 1 ,
... x**3 + y**3 - z ,
... x**2 + y**2 + z**2 - 1 )

在某个变量中寻找一个线性方程并求解它，存储解决方案

>>> reps=[(w,solve(eqs[0],w)[0])]

更新方程式

>>> eqs=[i.subs(reps) for i in eqs]; eqs
[0, -3*d*x*y + 3*d*y**2, -3*d*x*z - d + 1, x**3 + y**3 - z, x**2 + y**2 + z**2 - 1]

重复

>>> reps.append((d, solve(eqs[2],d)[0]))
>>> eqs=[i.subs(reps) for i in eqs]; eqs
[0, -3*x*y/(3*x*z + 1) + 3*y**2/(3*x*z + 1), -3*x*z/(3*x*z + 1) + 1 - 1/(3*x*z + 1), x**3 + y**3 - z, x**2 + y**2 + z**2 - 1]
>>> [cancel(i) for i in eqs]
[0, -(3*x*y - 3*y**2)/(3*x*z + 1), 0, x**3 + y**3 - z, x**2 + y**2 + z**2 - 1]

>>> reps.append((z, solve(eqs[-2],z)[0]))
>>> [cancel(i.subs(reps)) for i in eqs]
[0, -(3*x*y - 3*y**2)/(3*x**4 + 3*x*y**3 + 1), 0, 0, x**6 + 2*x**3*y**3 + x**2 + y**6 + y**2 - 1]

>>> reps.append((x, solve(eqs[1],x)[0]))
>>> [cancel(i.subs(reps)) for i in eqs]
[0, 0, 0, 0, 4*y**6 + 2*y**2 - 1]

直到只有 1 个等式。它是 y**2 的立方体，所以会有六个解

>>> yy=solve(eqs[-1].subs(reps),y)
>>> [i.n(2) for i in yy]
[-0.62, 0.62, 0.55 - 0.71*I, 0.55 + 0.71*I, -0.55 - 0.71*I, -0.55 + 0.71*I]

现在为每个 y 建立完整的解决方案

>>> for i in yy:
...     sol=[(y,i)]
...     for v,e in reversed(reps):
...         sol.append((v,e.subs(sol).n(2)))
...     sol[0] = sol[0][0],sol[0][1].n(2)
...     print(sol)
[(y, -0.62), (x, -0.62), (z, -0.48), (d, 0.53), (w, 0.49)]
[(y, 0.62), (x, 0.62), (z, 0.48), (d, 0.53), (w, -0.49)]
[(y, 0.55 - 0.71*I), (x, 0.55 - 0.71*I), (z, -1.3 - 0.59*I), (d, -0.26 - 0.2*I), (w, 0.43 - 0.12*I)]
[(y, 0.55 + 0.71*I), (x, 0.55 + 0.71*I), (z, -1.3 + 0.59*I), (d, -0.26 + 0.2*I), (w, 0.43 + 0.12*I)]
[(y, -0.55 - 0.71*I), (x, -0.55 - 0.71*I), (z, 1.3 - 0.59*I), (d, -0.26 + 0.2*I), (w, -0.43 - 0.12*I)]
[(y, -0.55 + 0.71*I), (x, -0.55 + 0.71*I), (z, 1.3 + 0.59*I), (d, -0.26 - 0.2*I), (w, -0.43 + 0.12*I)]

这会让您顺利地分析这个方程组。我构建的解决方案仅用于演示目的；如果您不在每次添加时进行评估而是等到结束然后在打印前进行sol = {k:v.n(2) for k,v in sol}，将会获得更准确的解决方案。

通过求解方程式，您发现基本上有 4 个线性关系和一个变量的多项式。

Sympy 解决 returns 一个空集

Sympy Solve returns an empty set

python

solver

sympy