pythonfunctionscipynumerical-integrationquad

Defining a integral of multi variable function as a second function python


Defining a integral of multi variable function as a second function python

I am using python to integrate a multivariable function over only one of the variables (it is a function of x and theta and I am integrating over theta from 0 to 2*pi, hence the result is a function of x). I have attempted the following:

import numpy as np
import scipy.integrate as inte
d=10.0
xvals=np.linspace(-d,d,1000)
def aIntegrand(theta,x):
    return 1/(2*np.pi)*np.sin(2*np.pi*x*np.sin(theta)/d)**2
def A(x):
    return (inte.quad(aIntegrand,0,2*np.pi,args=(x,))[0])**(1/2)

plt.plot(xvals,A(xvals))
plt.xlabel("x")
plt.ylabel("A(x)")
plt.show()

and I get the following error:

TypeError: only size-1 arrays can be converted to Python scalars

I assume this is because the result of the quad integrator is an array with two elements and python doesn't like defining the function based on an indexed array? That is a complete guess of the problem though. If anyone knows how I can fix this and could let me know, that'd be great :)


A second attempt

I have managed to get the plot of the integral successfully using the following code:

import numpy as np
import scipy.integrate as inte
import matplotlib.pyplot as plt
d=10.0
xvals=np.linspace(-d,d,1000)
thetavals=np.linspace(0.0,2*np.pi,1000)
def aIntegrand(theta,x):
    return 1/(2*np.pi)*np.sin(2*np.pi*x*np.sin(theta)/d)**2
def A(x):
    result=np.zeros(len(x))
    for i in range(len(x)):
        result[i]=(inte.quad(aIntegrand,0,2*np.pi,args=(x[i],))[0])**(1/2)
    return result
def f(x,theta):
    return x**2* np.sin(theta)

plt.plot(xvals,A(xvals))
plt.xlabel("x")
plt.ylabel("A(x)")
plt.show()

Integral Plot

But this does not give A(x) as a function, due to the way I have defined it, it requires an array form input. I need the function to be of the same form as aIntegrand, where when given parameters returns a single value & so the function can be integrated repeatedly.


Solution

  • I do not think that what you seek exists within Scipy. However, you have at least two alternatives.

    1. First, you can create an interpolator using interp1d. This means that the range of x values you initially give determines the range for which you will be able to call the interpolator. This interpolator can then be called for any value of x in this interval.
    2. You can perform the second integral even though you do not have a callable. The function simps only requires values and their location to provide an estimate of the integration process. Otherwise, you can perform the double integration in one call with dblquad.