如何加快 nsolve 或二分法？
How to speed up nsolve or the bisection method?
我正在编写一个需要某种根查找器的程序，但我使用的每个根查找器都非常慢。我正在寻找加快速度的方法。
我使用了 SymPy 的 nsolve，虽然它产生了非常精确的结果，但它非常慢（如果我对我的程序进行 12 次迭代，则需要 12 多个小时才能 运行）。我编写了自己的二分法，效果更好，但仍然很慢（12 次迭代需要大约 1 小时到 运行）。我一直找不到符号引擎求解器，或者这就是我要使用的。我将 post 我的两个程序（使用二分法和 nsolve）。非常感谢任何有关如何加快速度的建议。
这是使用 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(10):
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)
这是使用我的二分法的代码：
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 bisection(fun, a, b, tol):
NMax = 100000
f = Lambdify(E, fun)
FA = f(a)
for n in range(NMax):
p = (b+a)/2
FP = f(p)
if FP == 0 or abs(b-a)/2 < tol:
return p
if FA*FP > 0:
a = p
FA = FP
else:
b = p
print("Failed to converge to desired tolerance")
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(11):
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 = bisection(eqn, E0 - 1, E0, 10**(-15))
print(roots)
正如我所说，它们都按预期工作，但速度非常慢。
以下是对您的代码的一些优化，
- 使用
Lambdify(E, fun, cse=True)使用公共子表达式消除
- 添加
pi = sympify(sympy.N(pi))以使用数值pi。由于大型表达式，将 pi 保持为符号会造成伤害。
- 将
.simplify 调用更改为 .expand 调用。
- 您的积分表达式具有特殊形式。它们具有特殊的形式，
integrate(r**n * exp(-p*r), (r, 0, inf)，可以很容易地集成。
In [21]: var("n, r, p", positive=True)
Out[21]: (n, r, p)
In [22]: integrate(q*r**n*exp(-p*r), (r, 0, oo))
Out[22]: p**(-n)*q*gamma(n + 1)/p
您可以使用下面这样的 hack 来利用它。 （理想情况下 sympy 应该能够更快地完成此操作，但 sympy 在这方面做得不好。去年夏天，我 运行 在尝试象征性地求解狄拉克和薛定谔方程以调试我的数字代码时遇到了同样的问题。我假设你正在尝试做类似的事情）
def integrate(*args):
args = list(args)
expr = args[0].expand()
r = sympy.S(args[1][0])
limits = args[1][1:]
p = sympy.Wild("p")
n = sympy.Wild("n")
q = sympy.Wild("q")
pattern = q * r**n * sympy.exp(p*r)
terms = expr.args
if not expr.is_Add:
terms = [expr]
result = 0
for arg in terms:
d = sympy.S(arg).match(pattern)
if d is None:
result += sympy.N(sympy.integrate(arg, args[1]))
continue
if d[p].is_number and d[q].is_number and d[n].is_number:
ex = d[q]*(-d[p])**(-d[n])/d[p]*sympy.lowergamma(d[n]+1, -d[p]*r)
result += sympify(sympy.factorial(d[n])*d[q]/(-d[p])**(d[n]+1))
else:
result += sympy.N(sympy.integrate(arg, args[1]))
return result
这 4 个更改对我来说将时间减少到 16 秒。
我正在编写一个需要某种根查找器的程序，但我使用的每个根查找器都非常慢。我正在寻找加快速度的方法。
我使用了 SymPy 的 nsolve，虽然它产生了非常精确的结果，但它非常慢（如果我对我的程序进行 12 次迭代，则需要 12 多个小时才能 运行）。我编写了自己的二分法，效果更好，但仍然很慢（12 次迭代需要大约 1 小时到 运行）。我一直找不到符号引擎求解器，或者这就是我要使用的。我将 post 我的两个程序（使用二分法和 nsolve）。非常感谢任何有关如何加快速度的建议。
这是使用 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(10):
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)
这是使用我的二分法的代码：
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 bisection(fun, a, b, tol):
NMax = 100000
f = Lambdify(E, fun)
FA = f(a)
for n in range(NMax):
p = (b+a)/2
FP = f(p)
if FP == 0 or abs(b-a)/2 < tol:
return p
if FA*FP > 0:
a = p
FA = FP
else:
b = p
print("Failed to converge to desired tolerance")
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(11):
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 = bisection(eqn, E0 - 1, E0, 10**(-15))
print(roots)
正如我所说，它们都按预期工作，但速度非常慢。
以下是对您的代码的一些优化，
- 使用
Lambdify(E, fun, cse=True)使用公共子表达式消除 - 添加
pi = sympify(sympy.N(pi))以使用数值pi。由于大型表达式，将pi保持为符号会造成伤害。 - 将
.simplify调用更改为.expand调用。 - 您的积分表达式具有特殊形式。它们具有特殊的形式，
integrate(r**n * exp(-p*r), (r, 0, inf)，可以很容易地集成。
In [21]: var("n, r, p", positive=True)
Out[21]: (n, r, p)
In [22]: integrate(q*r**n*exp(-p*r), (r, 0, oo))
Out[22]: p**(-n)*q*gamma(n + 1)/p
您可以使用下面这样的 hack 来利用它。 （理想情况下 sympy 应该能够更快地完成此操作，但 sympy 在这方面做得不好。去年夏天，我 运行 在尝试象征性地求解狄拉克和薛定谔方程以调试我的数字代码时遇到了同样的问题。我假设你正在尝试做类似的事情）
def integrate(*args):
args = list(args)
expr = args[0].expand()
r = sympy.S(args[1][0])
limits = args[1][1:]
p = sympy.Wild("p")
n = sympy.Wild("n")
q = sympy.Wild("q")
pattern = q * r**n * sympy.exp(p*r)
terms = expr.args
if not expr.is_Add:
terms = [expr]
result = 0
for arg in terms:
d = sympy.S(arg).match(pattern)
if d is None:
result += sympy.N(sympy.integrate(arg, args[1]))
continue
if d[p].is_number and d[q].is_number and d[n].is_number:
ex = d[q]*(-d[p])**(-d[n])/d[p]*sympy.lowergamma(d[n]+1, -d[p]*r)
result += sympify(sympy.factorial(d[n])*d[q]/(-d[p])**(d[n]+1))
else:
result += sympy.N(sympy.integrate(arg, args[1]))
return result
这 4 个更改对我来说将时间减少到 16 秒。