-1

我有一个多项式，我试图找到所有使用 nsolve 进行数值计算的根。当我尝试使用 nsolve 查找最低根（这是我真正需要的唯一根，但我不介意找到所有根）时，我收到一条错误消息，提示“无法从 -0.499923944877944 创建 mpf”。

我尝试过使用多个不同的求解器。当我使用 SymPy 的求解时，它只找到了 5 个根（应该有 6 个）。使用solve也花了很长时间，因为我相信它最初是在尝试象征性地解决它。我尝试了求解集，但没有给出正确的答案。

以下是我的所有代码。一切都按预期进行，直到最底部的 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，但它甚至没有达到它给我一个根的程度。

4

2 回答 2

0

需要考虑两件事：nsolve如果您知道根所在的区间，将接受 solver='bisect' 关键字：

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

此外，real_roots在这种情况下，这可能是获取根的有用方法。在您发布的代码中使用 x in range(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]
于 2019-05-24T18:28:14.317 回答
0

使用 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 花费的时间太长（当我在 for 循环中进行 11 次迭代时，我一直在等待它解决它，这太荒谬了） .

于 2019-05-21T04:36:29.847 回答