Выполнение двойного интеграла по матрице пределов

Question

Я недавно изучал, как выполнять двойные интегралы в python. Это то, что я использую:

myint = dblquad(lambda xa,xb: np.exp(-(xa-xb)**2),-np.inf,x1,lambda x: -np.inf, lambda x: x2) 

и для целей тестирования я выбрал x1 и x2 будет сказать, 5 и 10. Это похоже на работу.

Но в действительности мои x1 = [1,2,3,4,5] и x2 = [5,6,7,8,9], и я хочу, чтобы двойной интеграл выполнялся над каждой комбинацией x1 и х2, т. е. матрица. Я мог бы сделать это с 2 для петель, я думаю, но я думал, что может быть лучший способ.

Итак, мой вопрос - как мне выполнить двойную интеграцию по матрице пределов.

спасибо.

---

редактировать:

я получил следующее предупреждение:

UserWarning: The maximum number of subdivisions (50) has been achieved. 
If increasing the limit yields no improvement it is advised to analyze 
the integrand in order to determine the difficulties. If the position of a 
local difficulty can be determined (singularity, discontinuity) one will 
probably gain from splitting up the interval and calling the integrator 
on the subranges. Perhaps a special-purpose integrator should be used. 

Означает ли это, что он не сходится? Я не очень понимаю сообщение.

Когда я сюжет:

y = exp(-(x-5)^2) 

, например, это выглядит как гауссовой кривой, так что нет никаких проблем интегрирования по правильно? Является ли проблема из-за двойного интеграла?

спасибо.

---

редактировать:

Ах, я вижу. Спасибо Raman Shah, я сейчас понимаю проблему.

Answer 1

Используя itertools, вы можете создать итератор ограничений, чтобы пройти. Это по существу представляет собой двойную петлю, но для более расширяемой, как вы можете иметь произвольное число входов с itertools.product и не хранить все границы сразу:

import numpy as np 
from scipy.integrate import dblquad 
import itertools 

f = lambda xa,xb: np.exp(-(xa-xb)**2) 
intg = lambda (x1,x2): dblquad(f,-np.inf,x1, 
           lambda x:-np.inf, 
           lambda x:x2) 


X1 = np.arange(1,6) 
X2 = np.arange(5,10) 
for limit in itertools.product(X1,X2): 
    print limit, intg(limit) 

Если вам нужно больше скорости, вы можете посмотреть в модуль multiprocessing для параллельных вычислений, поскольку каждый процесс независим.

Answer 2

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