Почему nsolve говорит "не может создать mpf из -0.499923944877944."? - PullRequest
Новая
       29

Почему nsolve говорит "не может создать mpf из -0.499923944877944."?

0 голосов
/ 20 мая 2019

У меня есть многочлен, который я пытаюсь найти все корни для численного использования nsolve.Когда я пытаюсь использовать nsolve, чтобы найти самый низкий корень (это единственный, который мне действительно нужен, но я не против найти их все), я получаю сообщение об ошибке «не могу создать mpf из -0.499923944877944».

Я пробовал использовать несколько разных решателей.Когда я использовал решение SymPy, он нашел только 5 корней (когда их должно быть 6).Использование решения также заняло очень много времени, так как я считаю, что сначала оно пытается решить символически.Я попытался решить, и это не дало правильного ответа.

Ниже приведен весь мой код.Все работает так, как и положено до nsolve в самом низу. 

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)

Наименьший корень должен быть больше или равен -0,5, но он даже не доходит до того, что дает мнекорень.

Ответы [ 2 ]

0 голосов
/ 24 мая 2019

Следует учесть две вещи: nsolve примет ключевое слово solver = 'bisect', если вы знаете интервал, в котором находится ваш корень: 

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

Кроме того, real_roots может быть полезным методом получения корней в подобных случаях. С х в ​​диапазоне (6) в коде, который вы разместили, полином порядка 88 получается и решается в коротком порядке с: 

>>> 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]
0 голосов
/ 21 мая 2019

Первоначальное предположение при использовании nsolve должно быть числом с плавающей точкой, но я вводил символическое число. Ниже приведен исправленный код. 

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)

Если кто-нибудь знает, как ускорить nsolve, я бы хотел это услышать. Я, вероятно, просто собираюсь написать свой собственный скрипт для метода деления пополам, так как nsolve занимает слишком много времени (я часами жду, когда он решит эту проблему, когда я выполню 11 итераций в цикле for, и это довольно нелепо) .

...