0

I want to convert an expression containing symbolic variables to a numeric one so that the expression may be subsequently used in an integration method '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))

The code runs up to here and returns the correct function fA. fA is a symbolic valued expression which I want to pass onto quad to perform an integration over (over 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))

My understanding is that this fails because the first argument of quad, i.e. the function, is of symbolic type and not numeric (=floating point) required for quad. How to make the function numeric and thus allow me to perform the integration? I've tried .subs and lambdify function but couldn't get it to work. The former seems to work only if numbers are supplied (i.e set NN=0.1 for example which I don't want to do) and I tried the following for 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))

but this gave error about number of positional arguments when I uncommented the print to check if all was working.

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

CAF
  • 329
  • 3
  • 14

1 Answers1

1

Yes, you should use lambdify. The first argument of lambdify is a tuple of symbols, not a tuple of tuples as in your code. The second argument is a SymPy expression. Example:

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

In your case, this would look like 

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

Here mpmath is chosen as the backend because it can evaluate Gamma function that your expression contains. Otherwise, it would probably be faster to use "numpy" backend option. See more on lambdify. 

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

How it works with quad... depends on whether you get real or complex numbers. When they are real, you can integrate with scipy.integrate.quad:

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

returns 30049812.82526324.

If they are complex, SciPy's quad will be unhappy with mpc type which it doesn't understand. But mpmath has its own quad, so use that instead: 

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

returns mpc(real='7170810.3848631922', imag='-192389955826656.31'). Here [1, 3] is the interval of integration.

  • Thanks! Did you manage to get it to work with my example? Removing those extraneous brackets that I had in the first argument I get an error that I could not understand: `AttributeError: 'Mul' object has no attribute 'log'` Any comments about that? – CAF Jan 30 '18 at 14:11
  • I used the def declarations because my expressions are involved in differentiation which generate very long expressions. To reuse them and not instead call them as a function seems inefficient. Isn't the indenting correct? I'm not sure the meaning of 'Mul' - A search on here refers to a clash of numpy and scipy namespaces but I did not write anything that would cause such a clash, as far as I can see. – CAF Jan 30 '18 at 15:17
  • Thanks very much indeed. Given numeric arguments, can I now pass this `f` onto `sciquad`? In my attempt at integration in second part of my OP, I replaced `fA` with `f(ktsq,qbarsq,NN,a,b,w)` where `f = lambdify((ktsq,qbarsq,NN,a,b,w), expr, "mpmath")`. I understand the integration is done over the first argument. But this gave me an mpf error. What is wrong with this? – CAF Jan 30 '18 at 18:32
  • Thanks! Sorry for asking too much but what happens if I wanted to multiply my `f`with,say, an interpolating function depending on the integration variable? The error is that `ValueError: object arrays are not supported`. If it's better I can accept your answer and post this as another question. thanks! – CAF Feb 01 '18 at 14:35
  • I don't have time to answer now, so a separate question makes more sense: someone else may be able to answer. –  Feb 01 '18 at 15:47
  • Ok thanks! here is the question if you are able to help https://stackoverflow.com/questions/48569794/how-to-extend-numerical-integration-method-to-allow-integration-over-an-object-a – CAF Feb 12 '18 at 12:28