用 sympy 求解非线性方程，但我得到的结果虚部很小

Question

我正在尝试用 Sympy 和 Python 求解非线性方程组。结果几乎是正确的，但虚部总是极小，而且这个过程很耗时。我也在 Matlab 下尝试了同样的计算，结果非常好而且很快。我知道可以忽略小的虚数部分。但我认为我的代码中一定有问题导致缓慢而虚幻的部分。谁能帮我这个？Python:3.6Sympy:1.1.1import sympyA1, B1, C1, D1, E1, F1 = (0.0019047619047619048,                          -1.7494954273533616e-19,                          0.0004761904761904762,                          -8.747477136766808e-18,                          0.047619047619047616,                          1.0)A2, B2, C2, D2, E2, F2 = (8.264462809917356e-05,                          -0.0,                          0.00033057851239669424,                          -0.008264462809917356,                          -0.03305785123966942,                          1.0)k, b = sympy.symbols('k b')eq1 = B1 ** 2 * b ** 2 + 2 * B1 * D1 * b - 2 * B1 * E1 * b * k - 4 * F1 * B1 * k + D1 ** 2 + 2 * D1 * E1 * k + \      4 * C1 * D1 * b * k + E1 ** 2 * k ** 2 - 4 * A1 * E1 * b - 4 * A1 * C1 * b ** 2 - 4 * C1 * F1 * k ** 2 - 4 * A1 * F1eq2 = B2 ** 2 * b ** 2 + 2 * B2 * D2 * b - 2 * B2 * E2 * b * k - 4 * F2 * B2 * k + D2 ** 2 + 2 * D2 * E2 * k + \      4 * C2 * D2 * b * k + E2 ** 2 * k ** 2 - 4 * A2 * E2 * b - 4 * A2 * C2 * b ** 2 - 4 * C2 * F2 * k ** 2 - 4 * A2 * F2s=sympy.solve([eq1,eq2],[k,b])print(s)这就是我在 Python 和 Sympy 下得到的，虚部非常小。它几乎需要10秒。这对我的整个项目来说是不可接受的。[(1.07269682322063 + 2.8315655624133e-28*I, -27.3048937553762 + 0.e-27*I), (1.79271658724978 - 2.83156477591471e-28*I, -76.8585791921325 - 0.e-27*I), (2.34194482854222 + 2.83156702952074e-28*I, -19.2027508047623 - 0.e-26*I), (5.20930842765403 - 2.83156580622397e-28*I, -105.800442914396 - 7.59430998293648e-28*I)]这就是我在 MATLAB 下使用“solve”得到的结果。它非常快。这就是我想要的。k =       5.2093       1.7927       1.0727       2.3419b =       -105.8      -76.859      -27.305      -19.203

Answer 1

SymPy 打算成为一个符号包，而不是数字，所以结果应该是预期的。但是，有一个函数nsolve可以用来找到数值解。在 SymPy/mpmath 中，没有专门的方法（我知道）可以计算多项式的所有根——在这种情况下，是一对二次方。您将必须一一找到根：

>>> list(nsolve((eq1, eq2), (k,b), (1, 1)))

[1.07269682322063, -27.3048937553762]

但是可以使用现有工具为此类方程构建求解器。以下是一个示例（在极端情况下可能会出现大量数值问题）：

def n2solve(eq1, eq2, x, y, n):

  """Return numerical solutions for 2 equations in

  x and y where each is polynomial of order 2 or less

  as would be true for equations describing geometrical

  objects.

  Examples

  ========

  >>> n2solve(x**2 + y, y**2 - 3*x*y + 4, x, y, 3)

  (-2.82, -7.96)

  (-1.34, -1.80)

  """

  from sympy.core.containers import Tuple

  from sympy.solvers.solvers import unrad, solve

  eqs = Tuple(eq1, eq2)

  sym = set([x, y])

  assert all(i.free_symbols == sym for i in eqs)

  anx = solve(eq1, x)[0]

  yeq = eq2.subs(x, anx)

  z = unrad(yeq)

  z = z[0] if z else yeq

  yy = real_roots(z)

  def norm(x,y):

    return abs((x**2+y**2).n(2))

  got=[]

  for yi in yy:

    yi = yi.n(n)

    ty = eqs.subs(y, yi)

    for xi in real_roots(ty[0]):

      xi = xi.n(n)

      got.append((norm(*ty.subs(x, xi)), xi, yi))

  return sorted([(x,y) for e,x,y in sorted(got)[:len(got)//2]])

这为问题中提出的方程提供了以下解决方案：

[(1.07, -27.3),

(1.79, -76.9),

(2.34, -19.2),

(5.21, -106.)]