Интеграция в случае, если одна переменная является массивом - PullRequest
Новая
       75

Интеграция в случае, если одна переменная является массивом

0 голосов
/ 19 декабря 2018

Я пытаюсь сделать двойное интегрирование, используя integrate.dblquad. Идея состоит в том, чтобы передать функцию, где одна из переменных (q) является массивом: с числовым интегрированием (для циклов по x и y это работает, но этоочень медленно).Scipy выдает следующую ошибку: TypeError: в скаляры Python могут быть преобразованы только массивы размера 1 

#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 был массивом значений для интеграла.

Извините за глупый вопрос и заранее благодарен:)

1 Ответ

0 голосов
/ 19 декабря 2018

Я не знаю простого решения для ускорения векторизованной версии двойного интеграла.Что я могу порекомендовать, так это ослабить допуск 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))

output

...