Решить нелинейные уравнения с симпой, но я получил результаты с небольшой мнимой частью

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 с помощью «решения». Это довольно быстро. Это то, что я хотел.

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