Vectorize двойной для петель

Question

мне нужно оценить интеграл, и мой код

r=0:25; 
t=0:250; 
Ti=exp(-r.^2); 
T=zeros(length(r),length(t)); 
for n=1:length(t) 
    w=1/2/t(n); 
    for m=1:length(r) 
    T(m,n)=w*trapz(r,Ti.*exp(-(r(m).^2+r.^2)*w/2).*r.*besseli(0,r(m)*r*w)); 
    end 
end 

В настоящее время оценка довольно быстро, но мне интересно, если есть способ векторизации двойной для цикла и сделать его даже быстрее, особенно когда используется функция trapz.

Answer 1

Вы можете оптимизировать его, передавая матрицу аргумент Y в trapz(A,Y), и с помощью dim = 2, то есть цикл становится:

r = 0:25; 
t = 0:250; 
Ti = exp(-r.^2); 

tic 
T = zeros(length(r),length(t)); 
for n = 1:length(t) 
    w = 1/2/t(n); 
    for m = 1:length(r) 
     T(m,n) = w*trapz(r,Ti.*exp(-(r(m).^2+r.^2)*w/2).*r.*besseli(0,r(m)*r*w)); 
    end 
end 
toc 

tic 
T1 = zeros(length(r),length(t)); 
for n = 1:length(t) 
    w = 1/2/t(n); 
    Y = bsxfun(@times,Ti.*r, exp(-bsxfun(@plus,r'.^2,r.^2)*w/2).*besseli(0,bsxfun(@times,r',r*w))); 
    T1(:,n) = w* trapz(r,Y,2); 
end 
toc 

max(abs(T(:)-T1(:))) 

Вы могли бы, вероятно, векторизации его полностью, будет иметь вид позже.