矢量化双for循环

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 

目前的评价是相当快的，但我不知道是否有向量化双for循环，使之形成方式甚至更快，特别是当使用功能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(:))) 

你也许完全矢量化它，以后会有看看。