provide robust two (geometric entity) numerical solver

[smichr]

As raised in the issue in which a light weight solution is given, a good numerical solver for equations that one might encounter when dealing with 2D geometric entities would be nice.

Here is such an equation set:

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 = 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)

And here is the start of a routine that I put (cobbled) together.

def n2roots(eq1, eq2, x, y, n=3):
  """Return numerical approximations of solution for coupled equations
  that are polynomial of degree 2 or less in variables x and y."""
  from sympy.core.containers import Tuple
  from sympy.solvers.solvers import unrad, solve
  from sympy import nfloat
  eqs = Tuple(*[e.expand().n(n) for e in (eq1,eq2)])
  if any(e.is_number for e in eqs) or not set([x,y]) & eqs.free_symbols:
    raise ValueError('one or more equations is missing symbols')
  for e in eqs:
    if any(e.diff(v,2).free_symbols for v in (x,y)):
      raise ValueError('only quadratics (or less) in %s and %s supported' % (x, y))
  # check for independent case: one are both are independent of the other
  free = [e.free_symbols for e in eqs]
  if any(len(i)==1 for i in free):
    if len(free[0]|free[1]) == 1:
      raise ValueError('uncoupled system')
    if len(free[0])==1:
      i = 0 if free[0].pop()==x else 1
      j = 0 if i else 1
    else:
      i = 1 if free[1].pop()==x else 0
      j = 0 if i else 1
    rv = []
    for xi in roots(eqs[i]):
      for yi in roots(eqs[j].subs(x, xi)):
        rv.append({x:xi, y:yi})
    print('case0')
    return nfloat(rv, n)
  # check for linear case: one is linear in x or y
  for i, e in enumerate(eqs):
    for j in (x, y):
      if j not in free[i]:
         continue
      if j not in e.diff(j).free_symbols:
        ji = solve(e, j)[0]
        s = x if j == y else y
        o = eqs[0 if i else 1]
        print('case1')
        rv = [{j: ji.subs(s, si).expand(), s: si} for si in roots(o.subs(j, ji))]
        return nfloat(rv, n)
  # biquadratic
  anx = solve(eqs[0], x)[0]
  yeq = eqs[1].subs(x, anx)
  z = unrad(yeq)
  if not z:
    z = yeq
  else:
    z = z[0]
  yy = roots(z)
  def norm(x,y):
    return abs((x**2+y**2).n(2))
  got=[]
  for yi in yy:
    ty = eqs.subs(y, yi)
    for xi in roots(ty[0],x):
      got.append((norm(*ty.subs(x,xi)),(xi,yi)))
  print('case2')
  rv = sorted(got, key=lambda x:x[0])
  return sorted(nfloat([t for i, t in rv[:len(rv)//2]], n))

But sorting out the use of roots, real_roots or nroots is not done. The above was tested with the following:

for e in [
  (x-2,y-3),
  (x-2,x+y-1),
  (x+y-2,x-y-3),
  (x**2+y**2-2,x**2+y-1),
  (x**2+y**2-1, (x-2)**2+(y-.1)**2-4.1),
  (x**2+2*x*y-y**2-1,x**2+y**2-5),
  (eq1.subs({k:x,b:y}),eq2.subs({k:x,b:y}))]:
  args = e + (x, y)
  print(n2roots(*args))

giving output

case0
[{x: 2.00, y: 3.00}]
case0
[{x: 2.00, y: -1.00}]
case1
[{x: 2.50, y: -0.500}]
case1
[{y: -0.618, x: -1.27}, {y: -0.618, x: 1.27}, {y: 1.62, x: -0.786*I}, {y: 1.62, x: 0.786*I}]
case2
[(0.178, 0.984), (0.276, -0.961)]
case2
[(-2.12, 0.707), (-1.00, -2.00), (1.00, 2.00), (2.12, -0.707)]
case2
[(1.07, -27.3), (1.79, -76.9), (2.34, -19.2), (5.21, -106.)]

You can see that imaginary numbers slipped by in one of the case1 cases. So there is still more work to be done.

Anyway, this is here if anyone --even me, perhaps -- would like to take it up.

cf #3571

[smichr]

Note also the following existing solution but it includes imaginary components:

>>> def g2(e, x, y, n=3):
...    sol = solve_poly_system([Poly(i, x, y) for i in e], x, y, domain=QQ)
...    if n is not None:
...        sol = [tuple([nfloat(i, n) for i in j]) for j in sol]
...    return sol
...
>>> g2((eq1,eq2),k,b)
[
 (1.07 + 2.29e-17*I, -27.3 + 0.e-15*I),
 (5.21 + 3.29e-18*I, -106.0 + 8.82e-18*I),
 (2.34 - 1.8e-17*I, -19.2 - 0.e-14*I),
 (1.79 - 8.15e-18*I, -76.9 - 0.e-15*I)]

[smichr]

There is a method given for solving two quadratic/conic expressions here.