Why nsolve is saying "cannot create mpf from -0.499923944877944."?

Question

I have a polynomial that I am trying to find all of the roots for numerically using nsolve. When I try to use nsolve to find the lowest root (that is the only one I really need, but I would not mind finding all of them), i get an error that says "cannot create mpf from -0.499923944877944."

I've tried using multiple different solvers. When I used SymPy's solve, it only found 5 of the roots(when there should be 6). Using solve also took a very long time, as I believe it is trying to solve it symbolically at first. I tried solveset, and that did not give the correct answer.

Below is all of my code. Everything works as it is supposed to up until the nsolve at the very bottom.

from symengine import *
import sympy
from sympy import Matrix
from sympy import nsolve

trial = Matrix()

r, E1, E = symbols('r, E1, E')
H11, H22, H12, H21 = symbols("H11, H22, H12, H21")
S11, S22, S12, S21 = symbols("S11, S22, S12, S21")
low = 0
high = oo

integrate = lambda *args: sympy.N(sympy.integrate(*args))

quadratic_expression = (H11-E1*S11)*(H22-E1*S22)-(H12-E1*S12)*(H21-E1*S21)
general_solution = sympify(sympy.solve(quadratic_expression, E1)[0])


def solve_quadratic(**kwargs):
    return general_solution.subs(kwargs)


def H(fun):
    return -fun.diff(r, 2)/2 - fun.diff(r)/r - fun/r


psi0 = exp(-3*r/2)
trial = trial.row_insert(0, Matrix([psi0]))
I1 = integrate(4*pi*(r**2)*psi0*H(psi0), (r, low, high))
I2 = integrate(4*pi*(r**2)*psi0**2, (r, low, high))
E0 = I1/I2
print(E0)

for x in range(5):
    f1 = psi0
    f2 = r * (H(psi0)-E0*psi0)
    Hf1 = H(f1).simplify()
    Hf2 = H(f2).simplify()

    H11 = integrate(4*pi*(r**2)*f1*Hf1, (r, low, high))
    H12 = integrate(4*pi*(r**2)*f1*Hf2, (r, low, high))
    H21 = integrate(4*pi*(r**2)*f2*Hf1, (r, low, high))
    H22 = integrate(4*pi*(r**2)*f2*Hf2, (r, low, high))

    S11 = integrate(4*pi*(r**2)*f1**2, (r, low, high))
    S12 = integrate(4*pi*(r**2)*f1*f2, (r, low, high))
    S21 = S12
    S22 = integrate(4*pi*(r**2)*f2**2, (r, low, high))

    E0 = solve_quadratic(
            H11=H11, H22=H22, H12=H12, H21=H21,
            S11=S11, S22=S22, S12=S12, S21=S21,
        )
    print(E0)

    C = -(H11 - E0*S11)/(H12 - E0*S12)
    psi0 = (f1 + C*f2).simplify()
    trial = trial.row_insert(x+1, Matrix([[psi0]]))

# Free ICI Part

h = zeros(x+2, x+2)
HS = zeros(x+2, 1)
S = zeros(x+2, x+2)

for s in range(x+2):
    HS[s] = H(trial[s]).simplify()

for i in range(x+2):
    for j in range(x+2):
        h[i, j] = integrate(4*pi*(r**2)*trial[i]*HS[j], (r, low, high))

for i in range(x+2):
    for j in range(x+2):
        S[i, j] = integrate(4*pi*(r**2)*trial[i]*trial[j], (r, low, high))

m = h - E*S
eqn = m.det()

roots = nsolve(eqn, E0)

print(roots)

The smallest root should be greater than or equal to -0.5, but it's not even getting to the point where it gives me a root.

Answer 1

Two things to consider: nsolve will accept a solver='bisect' keyword if you know the interval in which your root lies:

>>> nsolve(x**2-3,x,(0.,2.),solver='bisect')
1.73205080756888

Also, real_roots may be a useful method of obtaining roots in cases like this. With x in range(6) in the code you posted, a polynomial of order 88 is obtained and solved in a short order with:

>>> eq=nsimplify(eqn, rational=True).as_numer_denom()[0]
>>> [i.n(3) for i in real_roots(Poly(eq,E))][:7]
[-0.836, -0.298, -0.117, -0.0821, 0.0854, 0.181, 0.399]

Answer 2

The initial guess when using nsolve needs to be a float, but I was putting in a symbolic number. Below is the corrected code.

from symengine import *
import sympy
from sympy import Matrix
from sympy import nsolve

trial = Matrix()

r, E1, E = symbols('r, E1, E')
H11, H22, H12, H21 = symbols("H11, H22, H12, H21")
S11, S22, S12, S21 = symbols("S11, S22, S12, S21")
low = 0
high = oo

integrate = lambda *args: sympy.N(sympy.integrate(*args))

quadratic_expression = (H11-E1*S11)*(H22-E1*S22)-(H12-E1*S12)*(H21-E1*S21)
general_solution = sympify(sympy.solve(quadratic_expression, E1)[0])


def solve_quadratic(**kwargs):
    return general_solution.subs(kwargs)


def H(fun):
    return -fun.diff(r, 2)/2 - fun.diff(r)/r - fun/r


psi0 = exp(-3*r/2)
trial = trial.row_insert(0, Matrix([psi0]))
I1 = integrate(4*pi*(r**2)*psi0*H(psi0), (r, low, high))
I2 = integrate(4*pi*(r**2)*psi0**2, (r, low, high))
E0 = I1/I2
print(E0)

for x in range(6):
    f1 = psi0
    f2 = r * (H(psi0)-E0*psi0)
    Hf1 = H(f1).simplify()
    Hf2 = H(f2).simplify()

    H11 = integrate(4*pi*(r**2)*f1*Hf1, (r, low, high))
    H12 = integrate(4*pi*(r**2)*f1*Hf2, (r, low, high))
    H21 = integrate(4*pi*(r**2)*f2*Hf1, (r, low, high))
    H22 = integrate(4*pi*(r**2)*f2*Hf2, (r, low, high))

    S11 = integrate(4*pi*(r**2)*f1**2, (r, low, high))
    S12 = integrate(4*pi*(r**2)*f1*f2, (r, low, high))
    S21 = S12
    S22 = integrate(4*pi*(r**2)*f2**2, (r, low, high))

    E0 = solve_quadratic(
            H11=H11, H22=H22, H12=H12, H21=H21,
            S11=S11, S22=S22, S12=S12, S21=S21,
        )
    print(E0)

    C = -(H11 - E0*S11)/(H12 - E0*S12)
    psi0 = (f1 + C*f2).simplify()
    trial = trial.row_insert(x+1, Matrix([[psi0]]))

# Free ICI Part

h = zeros(x+2, x+2)
HS = zeros(x+2, 1)
S = zeros(x+2, x+2)

for s in range(x+2):
    HS[s] = H(trial[s]).simplify()

for i in range(x+2):
    for j in range(x+2):
        h[i, j] = integrate(4*pi*(r**2)*trial[i]*HS[j], (r, low, high))

for i in range(x+2):
    for j in range(x+2):
        S[i, j] = integrate(4*pi*(r**2)*trial[i]*trial[j], (r, low, high))

m = h - E*S
eqn = m.det()

roots = nsolve(eqn, float(E0))

print(roots)

If anyone knows how to speed up nsolve, I would love to hear it. I am probably just going to write my own script for the bisection method, as nsolve is taking way too long (I've been waiting hours for it to solve it when I do 11 iterations in that for loop, and that is quite ridiculous).