Преобразование символьного выражения в числовое для использования в quad - используйте lambdify?

Question

Я хочу преобразовать выражение, содержащее символические переменные, в числовое, чтобы впоследствии выражение можно было использовать в методе интеграции 'quad'.

import numpy 
import math as m
import scipy
import sympy

#define constants                                                                                                                                                                    
gammaee = 5.55e-6
MJpsi = 3.096916
alphaem = 1/137
lambdasq = 0.09
Ca = 3
qOsq = 2

def qbarsq(qsq):
    return (qsq+MJpsi**2)/4


def xx(qbarsq, w):
    return 4*qbarsq/(4*qbarsq-MJpsi**2+w**2)

from sympy import *
x,NN,a,b,ktsq,qbarsq,w = symbols('x NN a b ktsq qbarsq w')


def xg(a,b,NN,ktsq,x):
    return NN*(x**(-a))*(ktsq**b)*exp(sqrt((16*Ca/9)*log(1/x)*log((log(ktsq/lambdasq))/(log(qOsq/lambdasq)))))


#prints symbolic derivative of xg                                                                                                                                                    
def func(NN,a,b,x,ktsq):
    return (-x*diff(log(xg(a,b,NN,ktsq,x)),x))

#print(func(NN,a,b,x,ktsq))                                                                                                                                                          

#prints symbolic expression for Rg                                                                                                                                                   
def Rg(NN,a,b,ktsq,x):
     return 2**(2*func(NN,a,b,x,ktsq)+3)/sqrt(m.pi)*gamma(func(NN,a,b,x,ktsq)+5/2)/gamma(func(NN,a,b,x,ktsq)+4)
#print(Rg(NN,a,b,ktsq,x))

#prints symbolic expression for Fktsq                                                                                                                                                
def FktsqDeriv(NN,a,b,x,ktsq):
    return diff(Rg(NN,a,b,ktsq,x)*xg(a,b,NN,ktsq,x),ktsq)
#print(FktsqDeriv(NN,a,b,x,ktsq))      

def Fktsq1(qbarsq,ktsq,NN,a,b,w):
    return FktsqDeriv(NN,a,b,x,ktsq).subs(x,4*qbarsq/(4*qbarsq-MJpsi**2+w**2))

#print(Fktsq1(qbarsq,ktsq,NN,a,b,w))


# symbolic expression for fA                                                                                                                                                         
def fA(ktsq,qbarsq,NN,a,b,w):
    return Fktsq1(qbarsq,ktsq,NN,a,b,w)*1/(qbarsq)*1/(qbarsq+ktsq)
print(fA(qbarsq,ktsq,NN,a,b,w)) 

Код работает здесь и возвращает правильную функцию fA, fA является символическим выражением, которое я хочу передать на четырехугольник, чтобы выполнить интеграцию по (более ktsq)

import scipy.integrate.quadrature as sciquad
def integrated_f(qbarsq,NN,a,b,w):
    return sciquad(fA,1,(w**2-MJpsi**2)/4, args=(qbarsq, NN, a, b, w))

Я понимаю, что это не удается, потому что первый аргумент quadфункция имеет символьный тип и не является числовой (= с плавающей запятой), необходимой для квадрата. Как сделать функцию числовой и, таким образом, позволить мне выполнить интеграцию? я пробовал .subs а также lambdify функционировать, но не смог заставить его работать. Первый, кажется, работает, только если указаны числа (т.е. NN=0.1 например, что я не хочу делать), и я попробовал следующее для lambdify

def test(ktsq):
    return fA(ktsq,qbarsq,NN,a,b,w)

f = lambdify(((qbarsq,NN,a,b,w),), test(ktsq))
#print(f(1,2,3,4,5))

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

TypeError: <lambda>() takes 1 positional argument but 5 were given

Answer 1

Да, вы должны использовать lambdify. Первый аргумент lambdify это кортеж символов, а не кортеж, как в вашем коде. Второй аргумент - это выражение SymPy. Пример:

from sympy import *
a, b, c, d, e = symbols('a b c d e')
expr = a*b + 2*c + d/e
f = lambdify((a, b, c, d, e), expr)
print(f(1, 2, 3, 4, 5))   # prints 8.8

В вашем случае это будет выглядеть так

expr = fA(qbarsq, ktsq, NN, a, b, w)
f = lambdify((qbarsq, ktsq, NN, a, b, w), expr, "mpmath")

Здесь в качестве внутреннего интерфейса выбран mpmath, поскольку он может оценивать гамма-функцию, содержащуюся в вашем выражении. В противном случае, вероятно, было бы быстрее использовать опцию "numpy" backend. Смотрите больше на lambdify.

print(f(1, 2, 3, 4, 5, 6))  #  (-4757.21371513605 + 58978.7828908493j)

Как это работает quad... зависит от того, вы получаете реальные или комплексные числа. Когда они реальны, вы можете интегрировать с scipy.integrate.quad:

from scipy.integrate import quad
quad(f, 3, 4, args=(2, 3, 4, 5, 6))[0]

возвращается 30049812.82526324,

Если они сложны, SciPy's quad будет недоволен mpc Тип, который он не понимает. Но mpmath имеет свой quadтак что используйте это вместо:

import mpmath as mp
mp.quad(lambda x: f(2, x, 3, 4, 5, 6), [1, 3])

возвращается mpc(real='7170810.3848631922', imag='-192389955826656.31'), Здесь [1, 3] - интервал интегрирования.