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

Question

我正在尝试用 Sympy 和 Python 求解非线性方程组。

结果几乎是正确的，但虚部总是极小，而且这个过程很耗时。

我也在 Matlab 下尝试了同样的计算，结果非常好而且很快。

我知道可以忽略小的虚数部分。但我认为我的代码中一定有问题导致缓慢而虚幻的部分。谁能帮我这个？

Python:3.6

Sympy:1.1.1

import sympy

A1, 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 * F1

eq2 = 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 * F2

s=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.3419

b =

       -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.)]