Как векторизовать функцию для интеграла2?

Question

Я хочу оценить двойной интеграл вида $$\int_{-\infty}^a \int_{-\infty}^b \sum_{i,j}^K a_ia_jx^iy^j\exp(-x^2 - y^2 + xy)dx dy $$

где $ a_i $ и $ a_j $ являются константами. Поскольку интеграл линейный, я могу поменять местами суммирование и интегрирование, но в этом случае мне нужно вычислить интегралы $K^2$, и это займет слишком много времени. В этом случае я делаю следующее:

for i = 1:K
   for j = 1:K
       fun = @(x,y) x.^i.*y.^j.*exp(-2.*(x.^2 + y.^2 - 2.*x.*y))
       part(i,j) = alpha(i)*alpha(j)*integral2(fun,-inf,a,-inf,b)
   end
end

Это занимает слишком много времени, поэтому я хочу оценить только один интеграл, но я не знаю, как векторизовать $\sum_{i,j}^K a_ia_jx^iy^j\exp(-x^2 - y^2 + xy)$, а именно, как поставить его на интеграл2. Буду очень признателен за любую помощь.

Answer 1

Похоже, вам нужно, чтобы i и j были третьим и четвертым измерениями, чтобы был шанс, что код будет работать.

У меня также нет интеграла 2 (я использую октаву, интеграл 2 - это новая функция matlab, которой еще нет в октаве), поэтому я не могу проверить ее, но я думаю, что-то вроде этого может сработать:

alphaset=zeros(1,1,K,K);
alphaset(1,1,1:K,1:K)=alpha(1:K)'*alpha(1:K);
i_set=zeros(1,1,K,1);
j_set=zeros(1,1,1,K);
i_set(:)=1:K;
j_set(:)=1:K;
fun=@(x,y)  x.^i_set.*y.^j_set.*exp(-2.*(x.^2 + y.^2 - 2.*x.*y));
part = squeeze(alphaset.*integral2(fun,-inf,a,-inf,b));

Как я уже сказал, я не могу обещать, что это сработает, потому что я не знаю, как работает integra2. Но если вы замените интеграл2 просто "сумма (сумма (вес ([1,2,4],[3,-1,2])))", то он будет работать так, как задумано для этой операции (то есть он суммирует по значениям x и y, и результатом является матрица над набором индексов).

Answer 2

1. Если вы просто хотите улучшить скорость, попробуйте parfor.

2. Пусть $X=(x,x^2,\cdots,x^K)$, $Y=(y,y^2,\cdots,y^K)$, $A=(a_{ij})$ $ матрица с $ a_ {ij} = a_ {i} a_ {j} $, затем $$\sum_{i,j}^K a_{i}a_{j}x^iy^j=XAY^{T}$$

У меня нет функции integra2 на моем matlab, поэтому я не проверял, сильно ли она улучшит скорость.

Кроме того, я думаю, что вам нужно использовать символы x и y, после того как вы вычислите $$XAY^{T}$$, а затем используйте matlabFunction для преобразования символического выражения в дескриптор функции. Вот мой тестовый код: syms x y; Х =[х, х ^2]; Y=[у, у ^2]; Z=X*Y'; fun =matlabFunction(Z); ff=@(x,y) x^2+y^2;. = Гг весело (х, у) * и далее (х, у);