如何找到多元多项式方程系统的复根 (Python)答案

【问题标题】：How To Find Complex Roots of A Multivariate Polynomial Equation System (Python)如何找到多元多项式方程系统的复根 (Python)
【发布时间】：2015-08-27 10:27:44
【问题描述】：

考虑以下方程组：

1623.66790917 * x ** 2 + 468.829686367 * x * y + 252.762128419 * y ** 2 + -1027209.42116 * x + -301192.975791 * y + 188804356.212 = 0

11154.1759415 * x ** 2 + 31741.0229155 * x * y + 32933.5622632 * y ** 2 + -16226174.4037 * x + -26323622.7497 * y + 6038609721.67 = 0

如您所见，系统有两对复杂的解决方案。我试过sympy，但没有成功。我想知道如何在 Python 中解决这个问题。顺便说一句，我没有一个很好的初始猜测来使用数字方法。

【问题讨论】：

这些方程可能实际上是一个四次方程，虽然通常可以用符号求解，但它是通常可以用符号求解的程度的极限。也许您正在寻找一种数值方法？
@skyking 一个揭示复杂根源的象征性解决方案，真是太好了。我尝试了 sympy solve_poly_system 但花了很长时间没有解决方案。
您想找出那个特定的方程组的解吗？但你已经知道了。大概你的意思是某种方程系统。什么课？
@strubbly 不寻找特定的方程组。我想找到一个方程组的所有四个解，该方程组由两个在射影单应性下共面圆的圆锥曲线组成。一对解表示单应性下的圆点（ICP）图像。

标签： python sympy complex-numbers polynomial-math equation-solving

【解决方案1】：

您可以使用 mpmath 的 findroot() 对这些方程进行数值求解。据我所知，没有办法告诉findroot() 找到多个根，但我们可以绕过这个限制：首先，找到一个解对 (xa, ya)，然后将方程除以 (x - xa )*(y - ya)。您确实必须提供一个初步的近似解决方案，但我设法找到了一些只在几次尝试中就可以工作的东西。

from mpmath import mp

mp.dps = 30
prec = 20

f_1 = lambda x, y: 1623.66790917 * x ** 2 + 468.829686367 * x * y + 252.762128419 * y ** 2 + -1027209.42116 * x + -301192.975791 * y + 188804356.212

f_2 = lambda x, y: 11154.1759415 * x ** 2 + 31741.0229155 * x * y + 32933.5622632 * y ** 2 + -16226174.4037 * x + -26323622.7497 * y + 6038609721.67

def solve(f1, f2, initxy):
    return mp.findroot([f1, f2], initxy, solver='muller')

def show(x, y):
    print 'x=', mp.nstr(x, prec)
    print 'y=', mp.nstr(y, prec)
    print mp.nstr(f_1(x, y), prec)
    print mp.nstr(f_2(x, y), prec)
    print

f1a = f_1
f2a = f_2
xa, ya = solve(f1a, f2a, (240+40j, 265-85j))
show(xa, ya)

f1b  = lambda x, y: f1a(x, y) / ((x - xa) * (y - ya))
f2b  = lambda x, y: f2a(x, y) / ((x - xa) * (y - ya))
xb, yb = solve(f1b, f2b, (290+20j, 270+30j))
show(xb, yb)

输出

x= (246.82064795986653023 + 42.076841530787279711j)
y= (261.83565021239842638 - 81.555049135736951496j)
(0.0 + 3.3087224502121106995e-24j)
(0.0 + 0.0j)

x= (289.31873055121622967 + 20.548128321524345062j)
y= (272.23440694481666637 + 29.381152413744722108j)
(0.0 + 3.3087224502121106995e-24j)
(0.0 + 0.0j)

【讨论】：

感谢您的解决方案和示例代码。正如我在 OP 上评论的那样，我想找到一个方程组的所有四个解，包括两个二次曲线的两个方程。这些圆锥曲线是射影单应性下的共面圆。一对解表示单应性下的圆点（ICP）图像。正如 OP 中提到的，我一开始没有很好的估计。我正在寻找一种更通用的方法来自动求解方程组。
@Bahribayli：嗯，是的，提出良好的初步估计并不容易。 :) 如果您可以为这些圆锥曲线创建成对的参数方程，那么您可以轻松地生成曲线上的点列表，然后在列表中搜索近似匹配；本质上是在 CxC 中绘制曲线并寻找它们相交的点。 OTOH，只需使用在圆锥曲线覆盖的 CxC“平面”区域内随机生成的点，您就可以获得可接受的结果。