2016-02-11 125 views
1

我需要評估的整體,我的代碼是矢量化雙for循環

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時。

回答

2

你可以通過矩陣參數Ytrapz(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(:))) 

你也許完全矢量化它,以後會有看看。

+0

感謝trapz提示,但代碼有一個bug ... – noir

+0

@noir什麼錯誤?數值結果與機牀數值精度相同。我得到'max(abs(T(:) - T2(:)))= 3.4694e-18'其中'T2'是用我的循環計算的。 – Oleg

+0

'Ti'和'r'應該用'repmat'填充,否則尺寸不匹配。 @Oneg – noir