Как ускорить символические производные длинных функций с помощью SymPy? - PullRequest
Новая
       101

Как ускорить символические производные длинных функций с помощью SymPy?

2 голосов
/ 16 мая 2019

Я пишу программу на Python для решения уравнения Шредингера с использованием метода Free ICI (хорошо, метод SICI прямо сейчас ... но во что это превратится Free ICI). Если это не звучит знакомо, то это потому, что по этому вопросу очень мало информации и абсолютно нет примеров кода для работы.

Этот процесс включает в себя итеративное достижение решения уравнения в частных производных. При этом существует много символических производных, которые необходимо выполнить. Проблема заключается в том, что при запуске программы функции, которые необходимо дифференцировать, продолжают увеличиваться и увеличиваться, так что на пятой итерации для вычисления символьных производных требуется очень много времени.

Мне нужно ускорить это, потому что я хотел бы достичь как минимум 30 итераций, и я бы хотел, чтобы это было сделано до выхода на пенсию.

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

Вот код, в котором содержится функция, которая вычисляет производные (функция inf_integrate - это просто метод составного Симпсона, поскольку он намного быстрее, чем использование integrate в SymPy, и не вызывает ошибок из-за колебаний функции): 

from sympy import *


def inf_integrate(fun, n, a, b):
    f = lambdify(r, fun)
    h = (b-a)/n
    XI0 = f(a) + f(b)
    XI1 = 0
    XI2 = 0

    for i in range(1, n):
        X = a + i*h

        if i % 2 == 0:
            XI2 = XI2 + f(X)
        else:
            XI1 = XI1 + f(X)

    XI = h*(XI0 + 2*XI2 + 4*XI1)/3

    return XI


r = symbols('r')

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

E1 = symbols('E1')
low = 10**(-5)
high = 40
n = 5000

g = Lambda(r, r)


psi0 = Lambda(r, exp(-1.5*r))

I1 = inf_integrate(4*pi*(r**2)*psi0(r)*H(psi0(r)), n, low, high)
I2 = inf_integrate(4*pi*(r**2)*psi0(r)*psi0(r), n, low, high)

E0 = I1/I2
print(E0)

for x in range(10):

    f1 = Lambda(r, psi0(r))
    f2 = Lambda(r, g(r)*(H(psi0(r)) - E0*psi0(r)))
    Hf1 = Lambda(r, H(f1(r)))
    Hf2 = Lambda(r, H(f2(r)))

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

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

    eqn = Lambda(E1, (H11 - E1*S11)*(H22 - E1*S22) - (H12 - E1*S12)*(H21 - E1*S21))

    roots = solve(eqn(E1), E1)

    E0 = roots[0]

    C = -(H11 - E0*S11)/(H12 - E0*S12)

    psi0 = Lambda(r, f1(r) + C*f2(r))

    print(E0)

Программа работает и точно соответствует ожидаемому результату, но это слишком медленный 1015 *. Любая помощь в ускорении этого процесса очень ценится.

Ответы [ 2 ]

1 голос
/ 16 мая 2019

Здесь вы можете сделать несколько вещей:

  1. Если вы профилируете свой код, вы заметите, что большую часть времени вы проводите в функции интеграции inf_integrate, в основном потому, что выиспользуя ручные петли Python.Это можно исправить, превратив аргумент в векторизованную функцию и используя процедуры интеграции SciPy (которые компилируются и, следовательно, быстро).

  2. Когда вы используете вложенные символьные выражения, это может быть целесообразнопроверка, может ли случайное явное упрощение помочь сдержать взрывную сложность.Похоже, что так и есть.

  3. Все определенные вами функции Lamda не нужны.Вы можете упростить работу с выражениями.Я не проверял, действительно ли это влияет на время выполнения, но, безусловно, помогает в следующем шаге (поскольку у SymEngine еще нет Lambda).

  4. Используйте SymEngine вместо SymPy.SymPy (на данный момент) основан исключительно на Python и, следовательно, медленный.SymEngine является его скомпилированным ядром и может быть значительно быстрее.Он имеет почти все необходимые функции.

  5. На каждом шаге вы решаете уравнение, которое не меняет свою природу: это всегда то же самое квадратное уравнение, меняются только коэффициенты.Решая это один раз в общем, вы экономите много времени, в частности, SymPy, не сталкиваясь со сложными коэффициентами.

Взяв все вместе, я получаю следующее: 

from symengine import *
import sympy
from scipy.integrate import trapz
import numpy as np

r, E1 = symbols('r, E1')
H11, H22, H12, H21 = symbols("H11, H22, H12, H21")
S11, S22, S12, S21 = symbols("S11, S22, S12, S21")
low = 1e-5
high = 40
n = 5000

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)

sampling_points = np.linspace(low,high,n)
def inf_integrate(fun):
    f = lambdify([r],[fun])
    values = f(sampling_points)
    return trapz(values,sampling_points)

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

psi0 = exp(-3*r/2)
I1 = inf_integrate(4*pi*(r**2)*psi0*H(psi0))
I2 = inf_integrate(4*pi*(r**2)*psi0**2)
E0 = I1/I2
print(E0)

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

    H11 = inf_integrate( 4*pi*(r**2)*f1*Hf1 )
    H12 = inf_integrate( 4*pi*(r**2)*f1*Hf2 )
    H21 = inf_integrate( 4*pi*(r**2)*f2*Hf1 )
    H22 = inf_integrate( 4*pi*(r**2)*f2*Hf2 )

    S11 = inf_integrate( 4*pi*(r**2)*f1**2 )
    S12 = inf_integrate( 4*pi*(r**2)*f1*f2 )
    S21 = S12
    S22 = inf_integrate( 4*pi*(r**2)*f2**2 )

    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()

На моей машине это значение достигает значения ½ за несколько секунд.

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

Ответ Wrzlprmft был великолепен.Я пошел дальше и исправил ситуацию, а также поменял неуклюжую интеграционную функцию на интеграцию с SymPy.Это не сработало в моем исходном коде, но прекрасно работает после исправлений / дополнений Wrzlprmft.Программа немного медленнее (все еще на несколько порядков быстрее, чем мой оригинал), но в ней больше нет ошибки, которая ограничивала точность.Вот окончательный код: 

from symengine import *
from sympy import *
import sympy

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

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)
I1 = N(integrate(4*pi*(r**2)*psi0*H(psi0), (r, low, high)))
I2 = N(integrate(4*pi*(r**2)*psi0**2, (r, low, high)))
E0 = I1/I2
print(E0)

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

    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()
...