Integration in case one variable is an array

Question

I am trying to do a double integrale, using integrate.dblquad The idea is to pass a function where one of the variables (q) is array: With the numerical integration (for loops over x and y it works, but it is very slow). Scipy gives the following error: TypeError: only size-1 arrays can be converted to Python scalars

#set of values for the variables:
q=np.linspace(0.0001, 0.6, num=200)
rho1=0.2
rho2=0.5
rho_s=0.340
a = 20.1
b = 11.12
c = 6.18
ta=6.0
tb=5.5
tc=2.2

import numpy as np
from scipy import integrate

#equation simplifier:
def Bessel_like(z):
    Bes = 3 * (np.sin(z) - z * np.cos(z)) / (z**3.)
    return Bes

def Intensity(rho1, rho2, rho_s, a, b, c, ta, tb, tc, q):


    V1        =  a * b* c 
    V1pV2     = (a+ta) * (b+tb) * tc
    factorV1    = V1    * (rho1-rho2)
    factorV1pV2 = V1pV2 * (rho2-rho_s)

    def f(x,y):

        t1_1  = np.square(a * np.cos(np.pi * x/3))
        t1_2  = np.square(b * np.sin(np.pi * x/3)) * (1 - np.square(y))
        t1_3  = np.square(c*y)
        t1    = q * np.sqrt(t1_1 + t1_2 + t1_3)

        t2_1  = np.square( (a+ta) * np.cos(np.pi * x/3) )
        t2_2  = np.square( (b+tb) * np.sin(np.pi * x/3) ) * (1 - np.square(y))
        t2_3  = np.square( (c+tc)*y )
        t2    = q * np.sqrt(t2_1 + t2_2 + t2_3)


        return np.square(factorV1 * Bessel_like(t1)  + factorV1pV2 * Bessel_like(t2) )

    Int = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1)

    return Int[0]

# latter on, calling integral
Icalc = Intensity(rho1, rho2, rho_s, a, b, c, ta, tb, tc, q)

What would be the easiest/most efficient way to do this, and assign the array of Int values to one variable (for each q, but single array, I don't need a q value stored). I want this because this is a part of a really big code, and so far Int was array of values for the integral.

Sorry for the stupid question and thank you in advance :)

Answer 1

There is no straightforward solution that I know of to speed up a vectorized version of double integral. What I can recommend is to loosen the tolerance of dblquad, by increasing epsabs to, say, 1e-6 or 1e-5

An additional useful option is to reduce the number of sample points in q and interpolate them using a Spline:

from scipy.interpolate import InterpolatedUnivariateSpline as IUS
def Intensity(q, rho1, rho2, rho_s, a, b, c, ta, tb, tc):
   # I reversed your variable order putting q first, so you can vectorize on q
    V1        =  a * b* c 
    V1pV2     = (a+ta) * (b+tb) * tc
    factorV1    = V1    * (rho1-rho2)
    factorV1pV2 = V1pV2 * (rho2-rho_s)

    def f(x,y):

        t1_1  = np.square(a * np.cos(np.pi * x/3))
        t1_2  = np.square(b * np.sin(np.pi * x/3)) * (1 - np.square(y))
        t1_3  = np.square(c*y)
        t1    = q * np.sqrt(t1_1 + t1_2 + t1_3)

        t2_1  = np.square( (a+ta) * np.cos(np.pi * x/3) )
        t2_2  = np.square( (b+tb) * np.sin(np.pi * x/3) ) * (1 - np.square(y))
        t2_3  = np.square( (c+tc)*y )
        t2    = q * np.sqrt(t2_1 + t2_2 + t2_3)


        return np.square(factorV1 * Bessel_like(t1)  + factorV1pV2 * Bessel_like(t2) )
    Int = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1)
    return Int[0]

# latter on, calling integral
args = [rho1, rho2, rho_s, a, b, c, ta, tb, tc]
Icalc = Intensity(q[0], *args)
print(Icalc)
# construct spline
ius = IUS(q[::10], np.vectorize(Intensity)(q[::10])
plt.plot(q, np.vectorize(Intensity)(q), 'go')
plt.plot(q, ius(q))