如何用Python求解非线性方程组？

Question

我正在尝试用Python解决这个非线性方程组：

2y³ + y² - y⁵ - x⁴ + 2x³ - x² = 0
(-4x³ + 6x² - 2x) / (5y⁴ - 6y² - 2y) = 1

这是我的代码：

from sympy import symbols, Eq, solve
x, y = symbols('x y')
eq1 = Eq(2*y**3 + y**2 - y**5 - x**4 + 2*x**3 - x**2, 0)
eq2 = Eq((-4*x**3 + 6*x**2 - 2*x)/(5*y**4 - 6*y**2 - 2*y), 1)
points = solve([eq1, eq2], [x,y])

我知道这个方程组有5个解，但我用这段代码却没有得到任何结果。有人知道我该怎么做吗？

顺便提一下，第二个方程是第一个方程的导数dy/dx，所以我基本上是在寻找第一个方程曲线上切线斜率为1的所有点。第一个方程的图像是这样的。

方程1的图像

我已经用Geogebra解决了这个问题，并得到了5个答案。我正在尝试用Python，想看看能否用任何Python库来解决这个问题。

使用Geogebra的解法

Answer 1

这些解（5个解中有3个的分母不为0）可以通过这些方程的Groebner基来找到（不过我对如何处理基中存在的两个x-y方程还不是很清楚）：

>>> from sympy import *

# your equations go here, then

>>> n, d = fraction((eq2.lhs-eq2.rhs).normal())
>>> sols = list(groebner((eq1,n),x,y))
>>> yy = real_roots(sols[-1])
>>> xy=list(sorted(set([(s.n(4),yi.n(4)) for e in sols[:2] for yi in yy for s in solve(e.subs(y,yi.n()))])))
>>> [d.subs(dict(zip((x,y),i))) for i in xy]
[6.697, 0, -0.1917, -0.2850, 0, -0.2850, -0.1917, 6.697]
>>> xy = [xy[i] for i in range(len(xy)) if i not in (1,4)]
>>> [i for i in xy if all(abs(v)<1e-4 for v in Tuple(eq1.lhs,n).subs(dict(zip((x,y),i))))]
[(-0.7575, 1.453), (0.2267, -0.5028), (1.106, -0.5731)]

这些只是概念验证用的，精度不高的值。

Answer 2

通过重新写第二个方程，让它没有分母，然后使用 nonlinsolve，我们可以找到所有的解（在这个特定的例子中）：

from sympy import symbols, Eq, solve, fraction, nonlinsolve
x, y = symbols('x y')
eq1 = Eq(2*y**3 + y**2 - y**5 - x**4 + 2*x**3 - x**2, 0)
eq2 = Eq((-4*x**3 + 6*x**2 - 2*x)/(5*y**4 - 6*y**2 - 2*y), 1)
n, d = fraction(eq2.lhs)
eq2 = Eq(n, d)
sol = nonlinsolve([eq1, eq2], [x, y])
print(len(sol))

# 11

几秒钟（或者几分钟，具体看你的电脑性能）后，nonlinsolve 找到了11个解。结果发现其中有6个是复杂解。提取出真实解的一种方法是：

xcoords, ycoords = [], []
for s in sol:
    t1, t2 = s
    if t1.is_real and t2.is_real:
        t1 = t1.n()
        t2 = t2.n()
        xcoords.append(t1)
        ycoords.append(t2.n())
        print("(%s, %s)" % (t1, t2))

# (0, 0)
# (1.00000000000000, 0)
# (1.10625265404821, -0.573060136806016)
# (0.226655046617082, -0.502774766639858)
# (-0.757531019709684, 1.45262440731163)

额外提示。如果你使用 这个绘图模块，你可以快速可视化你正在计算的内容：

from spb import *
graphics(
    implicit_2d(eq1, (x, -2, 2), (y, -3, 3), n=1000),
    implicit_2d(eq2, (x, -2, 2), (y, -3, 3), n=1000),
    list_2d(xcoords, ycoords, scatter=True),
    xlabel="x", ylabel="y"
)