Numerical integration 'quad' error

[CAF123]

I have defined a series of functions below with the end function `fA` being inserted for a numerical integration. The integration is with respect with one variable so the other arguments are passed as numeric so that the integration method (quad) may proceed

    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 up to here is runnable and returns a correct function by uncommenting the last print command. I now want to integrate this last function over `ktsq` as follows,

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

    a=0.1
    NN=0.5
    b=-0.2
    w=89
    qbarsq=5
    result = integrated_f(NN,a,b,w,qbarsq)
    print(result)

The problem is where I try to get a number out from this integration by specifying numerical values for each of the other parameters. The error is

`ValueError:
Can't calculate 1st derivative wrt 989.426138911118.`

My only interpretation of this is that the method cannot cope with the complexity of the function (even though I think it is relatively simple in structure) because I did not define any more derivatives and certainly not with respect to this value. If that's correct, is there an easy fix?

 

[Mark44]

What variable or expression has a value of 989.426138911118?
Are all of the values in args necessary to calculate fA() of something? That seems quite complicated to me.
My advice is to see if you can get the integration to work on a simpler function on a specific interval -- a function for which you can calculate the integral by hand so as to compare the results.

 

[CAF123]

What variable or expression has a value of 989.426138911118?

This I am not sure of either, as I mentioned in OP, I did not specify any variable with this value. I interpreted it to mean something underlying in the integration method quad.

Are all of the values in args necessary to calculate fA() of something? That seems quite complicated to me.

I wish to do a numerical integration with numeric NN,a,b,w,qbarsq, and integrated over ktsq. The args specify that the quantities are numeric I think so the numerical integration may proceed.

My advice is to see if you can get the integration to work on a simpler function on a specific interval -- a function for which you can calculate the integral by hand so as to compare the results.

I tried replacing fA with something trivial like ktsq*NN*w and get an error about 'TypeError: 'Mul' object is not callable' , but I do not understand that.

Basically I wish to compute a numerical integration of fA over ktsq over this large interval. Can you suggest a more efficient way or help with the above?
Thanks!

 

[Mark44]

Is the code from the 2nd example you posted in a separate file from that of the 1st example you posted? If so, the 2nd example will need to import the name fA.

Also, the names you have chosen for variables and functions are pretty bad -- e.g., NN, qbarsq, w, MJpsi, FktsqDeriv, xx, Ca, xg, Rg, etc. Names whose purposes can't be easily inferred makes it difficult to comprehend or debug a program.

 

[CAF123]

Is the code from the 2nd example you posted in a separate file from that of the 1st example you posted? If so, the 2nd example will need to import the name fA.

It’s in the same file, I just replaced fA in the integration method with this trivial example.

Also, the names you have chosen for variables and functions are pretty bad -- e.g., NN, qbarsq, w, MJpsi, FktsqDeriv, xx, Ca, xg, Rg, etc. Names whose purposes can't be easily inferred makes it difficult to comprehend or debug a program.

The function names stem from a physics background so they make sense to me, although I can understand why they would be troublesome to an outsider. But regardless, the first block out of two that I posted above before the integration runs so I just need to understand the integration.

 

[Mark44]

I wish to do a numerical integration with numeric NN,a,b,w,qbarsq, and integrated over ktsq. The args specify that the quantities are numeric I think so the numerical integration may proceed.

But you have the following line in the first code example. I'm not familiar with "symbols" which might be part of numpy or scipy. I gather from this line that all of these names (x, NN, etc.) should be treated symbolically, not as numbers. That seems to disagree with your statement above.

x,NN,a,b,ktsq,qbarsq,w = symbols('x NN a b ktsq qbarsq w')

 

[CAF123]

But you have the following line in the first code example. I'm not familiar with "symbols" which might be part of numpy or scipy. I gather from this line that all of these names (x, NN, etc.) should be treated symbolically, not as numbers. That seems to disagree with your statement above.

x,NN,a,b,ktsq,qbarsq,w = symbols('x NN a b ktsq qbarsq w')

oh! I see. I made them symbolic at first so I could do the differentiation but now I wish to pass them as variables unaffected by the numerical integration. But I believe for this integration to be done, all variables have to be defined numeric. Is it possible to make them numeric now? I also thought about using different names, e.g a->aNum but Python might interpret this new variable as still symbolic. Thanks.

Aside: To get a bigger picture, this integration is part of the theoretical model involved in a chi square fit. I will then minimize the chi square for optimal parameters NN,a,b.

 

[Mark44]

oh! I see. I made them symbolic at first so I could do the differentiation but now I wish to pass them as variables unaffected by the numerical integration. But I believe for this integration to be done, all variables have to be defined numeric. Is it possible to make them numeric now? I also thought about using different names, e.g a->aNum but Python might interpret this new variable as still symbolic. Thanks.

The variables that you are making symbolic won't have numeric values. After the differentiation occurs, just define a new set of variables with the appropriate values.
As far as I know, once you have made a variable symbolic, you can't come back around and say that it should now have a numeric value. I could be wrong in this, though -- I'm more familiar with Python itself, and not at all with numpy, scipy, etc.

Also, I would suggest that you get in the habit of using variable names that are more suggestive of what you're going to do with them. While you understand what your code is supposed to do now, you might find yourself being confused if you look at this code after six months, not to mention that it is hard to understand for a bystander. I will say that as code goes, your code isn't as awful as some code I've seen by physicists, but its readability could be improved with just a little more thought.

 

[CAF123]

The variables that you are making symbolic won't have numeric values. After the differentiation occurs, just define a new set of variables with the appropriate values.

So basically the differentiation will give me a function in terms of symbolic NN,a,b etc. I now want to convert these to numeric variables for use in the integration method. Is there an explicit way to define variables as numeric in python? I know in mathematica in the function argument it is like e.g NN_?NumericQ

I will say that as code goes, your code isn't as awful as some code I've seen by physicists, but its readability could be improved with just a little more thought.

ok thanks for the tip.

 

[Mark44]

So basically the differentiation will give me a function in terms of symbolic NN,a,b etc. I now want to convert these to numeric variables for use in the integration method. Is there an explicit way to define variables as numeric in python?

Ordinary scalar variables by default are numeric in Python, aside from variables that contain characters. By "scalar" I mean non-array, non-tuple, non-dictionary types of variables. You can define a variable simply by giving it a value.

I know in mathematica in the function argument it is like e.g NN_?NumericQ

Python isn't mathematica.

 

[CAF123]

Ordinary scalar variables by default are numeric in Python, aside from variables that contain characters. By "scalar" I mean non-array, non-tuple, non-dictionary types of variables. You can define a variable simply by giving it a value.

I see, so everywhere in my symbolic expression where I have an ‘a’, for example, I need to replace this with ‘aNum’. Would e.g

 fA(ktsq, qbarsq, NN,a,b,w).subs(a,aNum)

work? I don’t have access to python just now so can’t check but it’s my thought. Sorry, still python newbie :)

Python isn't mathematica.

yup I just gave that as an indication of what I was trying to do with replacement of variable name.

 

[Mark44]

I see, so everywhere in my symbolic expression where I have an ‘a’, for example, I need to replace this with ‘aNum’. Would e.g

 fA(ktsq, qbarsq, NN,a,b,w).subs(a,aNum)

work? I don’t have access to python just now so can’t check but it’s my thought. Sorry, still python newbie :)

I don't guarantee that it'll work, but give it a shot!

 

[CAF123]

I don't guarantee that it'll work, but give it a shot!

I tried

def fA(ktsq,qbarsq,NNn,aN,b,w):
    fA(ktsq,qbarsq,NN,a,b,w).subs({NN: NNn, a: aN} )
print(fA(ktsq,qbarsq,NNn,aN,b,w))

but it said that NNn and aN were undefined: 'NameError: name 'aN' is not defined'

To make it work I just write e.g NNn=0.1, aN=0.2 before but I do not want to assign them numbers but just keep them of the numeric type. I looked around and saw something called a 'lamdify' function but this too also seems to rely on a number input to work. I haven't found anything else online yet that will help me.

 

[Mark44]

What does .subs(...) do? It doesn't seem to be part of Python, as I can't find it in the Python 3.5.0 documentation.

I'm having a hard time understanding your code, partly due to its organization. Your integrated_f() function (2nd code snippet in post #1) calls sciquad(), whose first argument is the function to be integrated.

def integrated_f(qbarsq,NN,a,b,w):
        return sciquad(fA,1,(w**2-MJpsi**2)/4, args=(qbarsq,NN, a, b, w))

This seems wrong to me, as integrated_f() doesn't list the function to be integrated as one of its parameters. It always calls sciquad() with the same function -- fA.

Further, to evaluate fA(), Fktsq1() has to be called, which in turn causes FktsqDeriv() to be called. (These last two function names are perfect examples of terrible, uninformative names.) This seems way complicated to me. To calculate an integral, it seems to me that all you need is a function to be evaluated at several points, but your code brings in derivatives to complicate things unnecessarily, as I understand things.

In your last post is the following:

def fA(ktsq,qbarsq,NNn,aN,b,w):
    fA(ktsq,qbarsq,NN,a,b,w).subs({NN: NNn, a: aN} )

I don't see how this could possibly work. You are defining fA() by calling itself, with some fiddling around with the parameters via subs() (whose purpose I don't understand, as it doesn't seem to be part of Python).

A more reasonable partitioning of chores in regard to the function is to have two separate function definitions -- one to be treated symbolically for finding its derivative, and another to be used to find function values. The two functions would be different names for the same operation. As a simple example, if the function is supposed to square its single input value, one function could be named fSymbolic() and the other could be named just f().

 

[ChrisVer]

That number is pretty random, I copy-pasted your code in my python IDE and that's the error I got...

Traceback (most recent call last):
  File "C:/Users/PainSama00017/PycharmProjects/learningTutorial/learningTutorial.py", line 83, in <module>
    result = integrated_f(NN,a,b,w,qbarsq)
  File "C:/Users/PainSama00017/PycharmProjects/learningTutorial/learningTutorial.py", line 76, in integrated_f
    return sciquad(fA,1,(w**2-MJpsi**2)/4, args=(qbarsq,NN, a, b, w))
  File "C:\Users\PainSama00017\venv\learningTutorial\lib\site-packages\scipy\integrate\quadrature.py", line 190, in quadrature
    newval = fixed_quad(vfunc, a, b, (), n)[0]
  File "C:\Users\PainSama00017\venv\learningTutorial\lib\site-packages\scipy\integrate\quadrature.py", line 85, in fixed_quad
    return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None
  File "C:\Users\PainSama00017\venv\learningTutorial\lib\site-packages\scipy\integrate\quadrature.py", line 115, in vfunc
    return func(x, *args)
  File "C:/Users/PainSama00017/PycharmProjects/learningTutorial/learningTutorial.py", line 69, in fA
    return Fktsq1(qbarsq, ktsq, NN, a, b, w) * 1 / (qbarsq) * 1 / (qbarsq + ktsq)
  File "C:/Users/PainSama00017/PycharmProjects/learningTutorial/learningTutorial.py", line 61, in Fktsq1
    return FktsqDeriv(NN, a, b, x, ktsq).subs(x, 4 * qbarsq / (4 * qbarsq - MJpsi ** 2 + w ** 2))
  File "C:/Users/PainSama00017/PycharmProjects/learningTutorial/learningTutorial.py", line 54, in FktsqDeriv
    return diff(Rg(NN, a, b, ktsq, x) * xg(a, b, NN, ktsq, x), ktsq)
  File "C:\Users\PainSama00017\venv\learningTutorial\lib\site-packages\sympy\core\function.py", line 1837, in diff
    return Derivative(f, *symbols, **kwargs)
  File "C:\Users\PainSama00017\venv\learningTutorial\lib\site-packages\sympy\core\function.py", line 1113, in __new__
    Can\'t calculate %s%s derivative wrt %s.''' % (count, ordinal, v)))
ValueError:
Can't calculate 1st derivative wrt 2.42613891111800.

Also checking the expression that is used in the sympy function that fails (__new__) I got this:
"1.02074930014737e+33*2**(-1.10544201885086e-33*x**0.8*(3.61846205571064e+32*x**(-0.8)*exp(0.567675746179508*sqrt(log(1/x))) - 5.13528286874444e+32*x**(-0.8)*exp(0.567675746179508*sqrt(log(1/x)))/sqrt(log(1/x)))*exp(-0.567675746179508*sqrt(log(1/x))) + 3)*x**0.2*exp(0.567675746179508*sqrt(log(1/x)))*gamma(-5.52721009425429e-34*x**0.8*(3.61846205571064e+32*x**(-0.8)*exp(0.567675746179508*sqrt(log(1/x))) - 5.13528286874444e+32*x**(-0.8)*exp(0.567675746179508*sqrt(log(1/x)))/sqrt(log(1/x)))*exp(-0.567675746179508*sqrt(log(1/x))) + 2.5)/gamma(-5.52721009425429e-34*x**0.8*(3.61846205571064e+32*x**(-0.8)*exp(0.567675746179508*sqrt(log(1/x))) - 5.13528286874444e+32*x**(-0.8)*exp(0.567675746179508*sqrt(log(1/x)))/sqrt(log(1/x)))*exp(-0.567675746179508*sqrt(log(1/x))) + 4)"

This looks pretty nasty...

 

[Tom.G]

@CAF123 Perhaps I'm way off-base because I've never used Python, but here goes.

It seems your call to integrated_f() is trying to pass values by keyword reference. This appearance is because the definition and the call to it have their arguements in a different order... however the syntax does not seem to match the example in this reference:
https://www.learnpython.org/en/Multiple_Function_Arguments

That reference, and many others, were found with: https://www.google.com/search?&q=python+function+arguments

Tom