如果一个变量是一个数组，则集成

Question

我正在尝试做一个双积分，使用 integrate.dblquad 这个想法是传递一个函数，其中一个变量 (q) 是数组：通过数值积分（对于 x 和 y 上的循环，它可以工作，但速度非常慢）。 Scipy 给出以下错误：类型错误：只有大小为 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)

什么是 easiest/most 有效的方法，并将 Int 值的数组分配给一个变量（对于每个 q，但单个数组，我不需要存储 q 值）。我想要这个是因为这是一个非常大的代码的一部分，到目前为止 Int 是积分值的数组。

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

Answer 1

据我所知，没有直接的解决方案可以加速矢量化版本的二重积分。我可以推荐的是通过将 epsabs 增加到 1e-6 或 1e-5

来放宽对 dblquad 的容忍度

另一个有用的选项是减少 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))

如果一个变量是一个数组，则集成

Integration in case one variable is an array

arrays

numpy

integrate

integral

scipy