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

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

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} $$

Я не integral2 функции на моем MATLAB, так что я не проверял, если это повысит скорость много.

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

Answer 2

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

У меня также нет интеграла2 (я использую октаву, integ2 - это новая функция 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)); 

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