Question

Я пишу программу, которая требует какого-либо корневого поиска, но каждый используемый мной корневой поиск неудовлетворительно медленный.Я ищу способ ускорить это.

Я использовал nsolve SymPy, и, хотя он дает очень точные результаты, он очень медленный (если я делаю 12 итераций моей программы, это занимает 12+часов бегать).Я написал свой собственный метод деления пополам, и он работает намного лучше, но все еще очень медленно (12 итераций занимает ~ 1 час).Мне не удалось найти решатель сименгина, или я бы этим воспользовался.Я опубликую обе мои программы (с помощью метода деления пополам и с помощью 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)

Как я уже сказал, они оба работают так, как и должны, но делают это очень медленно.

isuruf · Answer 1 · 22 мая 2019

Вот некоторые оптимизации для вашего кода,

Используйте Lambdify(E, fun, cse=True), чтобы использовать Common Subexpression Elmination
Добавьте 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

Вы можете воспользоваться этим преимуществом, используя хак, подобный приведенному ниже.(В идеале 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 секунд для меня.

Как ускорить nsolve или метод деления пополам?

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Как ускорить nsolve или метод деления пополам?

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Похожие темы