scipy интегрировать по массиву с ограничениями переменных

Question

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

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

импорт NumPy в нп

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, где вы должны предоставить значения подынтегрального выражения напрямую и dx-интервал, но это будет слишком неточно для моей сложной функции (функция marcum Q не может найти никакой реализации, это может быть другим способом ее усечения).

Answer 1

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

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

1. Учитывая интеграл (т), которые я буду считать, является интегралом функции Р (х) от 0 до Т [это может быть обобщена на произвольный нижний предел], написать это как дифференциальное уравнение 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. Нет гарантии, что этот подход будет численно стабильным для вашей проблемы.