如何求多元多项式方程组的复根 (Python)
How To Find Complex Roots of A Multivariate Polynomial Equation System (Python)
考虑以下方程组:
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 中解决这个问题。顺便说一句,我没有一个很好的初步猜测来使用数字方法。
您可以使用 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)
考虑以下方程组:
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 中解决这个问题。顺便说一句,我没有一个很好的初步猜测来使用数字方法。
您可以使用 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)