Sympy 的solve() 函数不会产生所有解决方案

Question

我需要求解以下非线性方程组：

我正在使用sympy的solve()函数来获得方程组的解。通常，solve() 会产生所有解决方案。我已经指定了常数 a、b、c、d、g、h 的值并期望获得所有解决方案。

但是solve()只产生一种解决方案，而且这也是错误的。我的问题是：

如何使用单个 python 实用程序获得所有符号解决方案？我不需要数值解。

这里是实现solve()的代码部分：

from sympy import  solve

from sympy import symbols
x,y,a,b,c,d,g,h = symbols('x y a b c d g h')
a = 0.12
b = -2.031529100521498e-30
c = b
d = 1
g = 11
h = 0
F = (y + a/2)*x*(2*x**2 - 2*y*a + 2*d**2 + b + c) - (x*y*g**2)
G = (x**2 - y*a + b + d**2)*(x**2 - y*a + c+d**2) - 4*((y+a/2)**2)*x**2 - (h**2 + (x**2)*(g**2))
sol = solve([F,G],(x,y),dict=True,simplify=False)
print(sol)

Answer 1

这些是你的方程：

In [45]: from sympy import  solve
    ...: 
    ...: from sympy import symbols
    ...: x,y,a,b,c,d,g,h = symbols('x y a b c d g h')
    ...: a = 0.12
    ...: b = -2.031529100521498e-30
    ...: c = b
    ...: d = 1
    ...: g = 11
    ...: h = 0
    ...: F = (y + a/2)*x*(2*x**2 - 2*y*a + 2*d**2 + b + c) - (x*y*g**2)
    ...: G = (x**2 - y*a + b + d**2)*(x**2 - y*a + c+d**2) - 4*((y+a/2)**2)*x**2 - (h**2 + (x**2)*(g*
    ...: *2))

我将把浮点数转换为有理数，因为浮点数和多项式一般不能很好地混合。我只是在这里指出，您的

b

值非常小，实际上为零（在标准 64 位浮点中）。如果您希望

b

非零，那么您需要从一开始就避免浮动，因为在执行上面所示的代码后它已经消失了。

In [46]: F, G = nsimplify(F), nsimplify(G) # don't use floats with polynomials

In [47]: F
Out[47]: 
                        ⎛   2   6⋅y    ⎞
-121⋅x⋅y + x⋅(y + 3/50)⋅⎜2⋅x  - ─── + 2⎟
                        ⎝        25    ⎠

In [48]: G
Out[48]: 
                                            2
     2           2        2   ⎛ 2   3⋅y    ⎞ 
- 4⋅x ⋅(y + 3/50)  - 121⋅x  + ⎜x  - ─── + 1⎟ 
                              ⎝      25    ⎠ 

方程 F 因式分解，所以让我们依次计算两个因式：

In [56]: F.factor()
Out[56]: 
  ⎛      2         2        2               ⎞
x⋅⎝1250⋅x ⋅y + 75⋅x  - 150⋅y  - 74384⋅y + 75⎠
─────────────────────────────────────────────
                     625                     

In [57]: _, F1, F2 = F.factor().args

In [58]: F1
Out[58]: x

In [59]: solve([F1, G], [x, y])
Out[59]: [(0, 25/3)]

我相信这就是您提到的解决方案。您将其描述为不正确，但对于代码中所示的方程来说，它绝对是正确的。

这很简单。现在对于

F2

，我们计算 Groebner 基：

In [61]: gb = groebner([F2, G], [x, y], method='f5b')

In [62]: print(gb[0])
x**2 + 16*y**4/121 + 595216*y**3/9075 + 892608*y**2/75625 + 1844017554*y/1890625 - 75643/75625

In [63]: print(gb[1])
y**5 + 74411*y**4/150 + 148777*y**3/1250 + 1845583341*y**2/250000 + 691424307*y/781250 - 113451/125000

我们可以看到

gb[1]

是

y

的五次方程。这个五次方程是不可约的，并且可能没有其根的激进公式（Abel-Ruffini），但是当系数是显式有理数时，SymPy 可以将根表示为 RootOf：

In [66]: r1, r2, r3, r4, r5 = Poly(gb[-1]).all_roots()

In [67]: r1
Out[67]: 
       ⎛          5               4               3                 2                              
CRootOf⎝18750000⋅x  + 9301375000⋅x  + 2231655000⋅x  + 138418750575⋅x  + 16594183368⋅x - 17017650, 0

⎞
⎠

每个根都可以用在第一个方程中来求解

x

，例如：

In [69]: print(solve([gb[0],y-r1], [x, y]))
[(-sqrt(-2250000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**4 - 200836800*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**2 + 17019675 - 16596157986*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0) - 1116030000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**3)/4125, CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)), (sqrt(-2250000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**4 - 200836800*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**2 + 17019675 - 16596157986*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0) - 1116030000*CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0)**3)/4125, CRootOf(18750000*x**5 + 9301375000*x**4 + 2231655000*x**3 + 138418750575*x**2 + 16594183368*x - 17017650, 0))]

还有其他方法可以做到这一点，但这给出了 RootOf 方面的所有解决方案。我不知道这是否是您想要的，但您可能希望 SymPy 返回一些不可能返回的东西。这是全套解决方案：

In [71]: sx1, sx2 = solve(gb[0], x)

In [72]: sols = solve([F1, G], [x, y], dict=True)

In [73]: sols.extend([{x:sx1.subs(y, r), y:r} for r in [r1,r2,r3,r4,r5]])

In [74]: sols.extend([{x:sx2.subs(y, r), y:r} for r in [r1,r2,r3,r4,r5]])

这就是它们的数字：

In [76]: for s in sols:
    ...:     print(s[x].n(3), s[y].n(3))
    ...: 
0 8.33
-0.0610 -496.
-10.9 -0.121
-0.0917 0.00102
-7.71 + 0.091*I -0.045 - 3.86*I
-7.71 - 0.091*I -0.045 + 3.86*I
0.0610 -496.
10.9 -0.121
0.0917 0.00102
7.71 - 0.091*I -0.045 - 3.86*I
7.71 + 0.091*I -0.045 + 3.86*I

关于您是否有兴趣以 a、b 等作为符号而不是明确的数字来获得解决方案，这个问题是不明确的。事实上

gb[-1]

是不可约五次方程，因此对于符号系数的情况不太可能存在任何有意义的显式解：

https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem

一般来说，面对此类问题时，最好考虑采用不同的方法。显式解决方案可能是不可能的（如此处）或者太复杂而没有多大用处（例如，当使用卡尔达诺的四次公式时）。当我提出这个建议时，似乎很少有人会注意到，但实际上，在这种情况下，最好是“不”寻求“明确”的解决方案。您的原始多项式方程已经比任何根式公式可以更好地表示解：对于符号工作，通常最好使用这些多项式作为系统解的“隐式”表示。我无法建议如何做到这一点，因为我不知道您下一步打算做什么（如果您想得到答案，请提出一个新问题）。