Scipy интегрировать по массиву с переменными границами

Question

Я пытаюсь интегрировать функцию по списку точек и передать весь массив функции интегрирования, чтобы векторизовать объект. Для начала, вызов scipy.integrate.quad слишком медленный, так как у меня есть что-то вроде 10 000 000 точек для интеграции. Использование scipy.integrate.romberg делает трюк намного быстрее, почти мгновенно, в то время как квад работает медленно, так как вы должны зациклить его или векторизовать его.

Моя функция довольно сложна, но для демонстрации, скажем, я хочу интегрировать x^2 из a в b, но x - это скалярный массив для оценки x. Например

импортировать numpy как np

from scipy.integrate import quad, romberg

def integrand(x, y):   

    return x**2 + y**2


quad(integrand, 0, 10, args=(10) # this fails since y is not a scalar

romberg(integrand, 0, 10)  # y works here, giving the integral over
                           # the entire range 

Но это работает только для фиксированных границ. Есть ли способ сделать что-то вроде

z = np.arange(20,30)
romberg(integrand, 0, z)  # Fails since the function doesn't seem to
                          # support variable bounds

Единственный способ, которым я вижу это, состоит в том, чтобы повторно реализовать сам алгоритм в numpy и использовать его вместо этого, чтобы я мог иметь переменные границы. Любая функция, которая поддерживает что-то вроде этого? Существует также romb, где вы должны указать значения integrand напрямую и интервал dx, но это будет слишком неточно для моей сложной функции (функция marcum Q не смогла найти какую-либо реализацию, которая могла бы быть другим способом расставить точки над ней)).

Answer 1

Наилучший подход при попытке оценить специальную функцию - написать функцию, которая использует свойства функции для быстрой и точной оценки ее во всех режимах параметров. Маловероятно, что один подход даст точные (или даже стабильные) результаты для всех диапазонов параметров. Прямая оценка интеграла, как в этом случае, почти наверняка сломается во многих случаях.

При этом общую проблему оценки интеграла по многим диапазонам можно решить, превратив интеграл в дифференциальное уравнение и решив его. Грубо говоря, шаги будут

1. Учитывая, что интеграл I(t), который я буду считать интегралом от функции f(x) от 0 до t [это можно обобщить до произвольного нижнего предела], запишите его как дифференциальное уравнение dI/dt = f(x).).
2. Решите это дифференциальное уравнение, используя scipy.integrate.odeint() для некоторых начальных условий (здесь I(0)) в некотором диапазоне времен от 0 до t. Этот диапазон должен содержать все пределы интереса. То, насколько точно это сделано, зависит от функции и точности ее оценки.
3. Результатом будет значение интеграла от 0 до t для набора t, который мы вводим. Мы можем превратить это в "непрерывную" функцию, используя интерполяцию. Например, используя сплайн, мы можем определить i = scipy.interpolate.InterpolatedUnivariateSpline(t,I),
4. Задан набор верхних и нижних пределов в массивах b а также aсоответственно, тогда мы можем оценить их все сразу как res=i(b)-i(a),

Для того, чтобы этот подход работал в вашем случае, вам потребуется тщательно изучить его по всем параметрам. Также отметим, что функция Маркума Q содержит полубесконечный интеграл. В принципе, это не проблема, просто преобразовать интеграл в единицу в конечном диапазоне. Например, рассмотрим преобразование x->1/x. Нет гарантии, что этот подход будет численно устойчивым для вашей проблемы.