如果一个变量是数组，则进行积分

Question

我正在尝试使用integrate.dblquad 进行二重积分 这个想法是传递一个函数，其中变量 (q) 之一是数组： 通过数值积分（对于 x 和 y 的循环它可以工作，但速度非常慢）。 scipy 给出以下错误： TypeError：只有 size-1 数组可以转换为 Python 标量

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

最简单/最有效的方法是什么， 并将

Int

值数组分配给一个变量（对于每个

q

，但单个数组，我不需要存储

q

值）。我想要这个，因为这是一个非常大的代码的一部分，到目前为止

Int

是积分值的数组。

抱歉问了这个愚蠢的问题，提前谢谢你:)

Answer 1

据我所知，没有直接的解决方案可以加速二重积分的矢量化版本。我可以建议的是放宽

dblquad

的容忍度，将

epsabs

增加到

1e-6

或

1e-5

另一个有用的选项是减少 q 中的样本点数量并使用样条线对它们进行插值：

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