矢量的雙for循環包括兩個變量

Question

我需要數值評估一些積分;這些都是在該圖像中所示的形式的正弦：

這些積分是N x N矩陣的矩陣元素，所以我需要以評估它們的所有可能的組合n和m的範圍爲1到N。積分是n和m對稱的，我在我目前的嵌套for環方案已經實施：

function [V] = coulomb3(N, l, R, R0, c, x) 
r1 = 0.01:x:R; 
r2 = R:x:R0; 
r = [r1 r2]; 
rl1 = r1.^(2*l); 
rl2 = r2.^(2*l); 
sines = zeros(N, length(r)); 
V = zeros(N, N); 
for i = 1:N; 
    sines(i, :) = sin(i*pi*r/R0); 
end 

x1 = length(r1); 
x2 = length(r); 
for nn = 1:N 
    for mm = 1:nn 
     f1 = (1/6)*rl1.*r1.^2.*sines(nn, 1:x1).*sines(mm, 1:x1); 
     f2 = ((R^2/2)*rl2 - (R^3/3)*rl2.*r2.^(-1)).*sines(nn, x1+1:x2).*sines(mm, x1+1:x2); 
     value = 4*pi*c*x*trapz([f1 f2]); 
     V(nn, mm) = value; 
     V(mm, nn) = value; 
    end 
end 

我想，在循環中調用sin(x)是個壞主意，所以我計算所有需要的值，並將其儲存。爲了評估積分我使用trapz，但由於第一和第二/第三積分有不同的範圍，功能值需要分開計算，然後組合。

我試過了幾種不同的矢量化方法，但唯一能給出正確結果的方法比上面的循環花費的時間要長得多（使用gmultiply但創建的數組很有效）。我也做了一個解析解（假設m和n是整數和R0 > R > 0），但這些解決方案包含一個餘弦積分（MATLAB中的cosint）函數，對於大的N來說，這個函數非常慢。

我不確定整個事情可以矢量化而不創建非常大的數組，但內部循環至少應該是可能的。任何想法將不勝感激！

我目前使用的輸入是：

​​

使用這些值

V(1, 1:3) = 873,379900963549 -5,80688363271849 -3,38139152472590 

雖然對角線值永遠不會隨R0因此它們更有趣的收斂。

Answer 1

您將失去從我的方法的問題的對稱性收益，但這意味着2損失的因素。可能性是你最終還是會受益的。

這個想法是使用多維數組，使用trapz支持這些輸入。我會證明你的身材的第一項，因爲兩個人應該同樣做，並且這個點的技術：

r1 = 0.01:x:R; 
r2 = R:x:R0; 
r = [r1 r2].'; 
rl1 = r1.'.^(2*l); 
rl2 = r2.'.^(2*l); 
sines = zeros(length(r),N);  %// CHANGED!! 
%// V = zeros(N, N); not needed now, see later 
%// you can define sines in a vectorized way as well: 
sines = sin(r*(1:N)*pi/R0);  %//' now size [Nr, N] ! 

%// note that implicitly r is of size [Nr, 1, 1] 
%// and sines is of size [Nr, N, 1] 
sines2mat = permute(sines,[1, 3, 2]); %// size [Nr, 1, N] 

%// the first term in V: perform integral along first dimension 
%//V1 = 1/6*squeeze(trapz(bsxfun(@times,bsxfun(@times,r.^(2*l+2),sines),sines2mat),1))*x; %// 4*pi*c prefactor might be physics, not math 
V1 = 1/6*permute(trapz(bsxfun(@times,bsxfun(@times,r.^(2*l+2),sines),sines2mat),1),[2,3,1])*x; %// 4*pi*c prefactor might be physics, not math 

關鍵的一點是，bsxfun(@times,r.^(2*l+2),sines)是大小[Nr,N,1]的矩陣，這又是乘以sines2mat，使用bsxfun，結果的大小爲[Nr,N,N]，並且元素(k1,k2,k3)對應於在徑向點​​,n=k2和m=k3處的被積函數。使用trapz()明確地說，第一維（這將是默認值）將其減少到大小爲[1,N,N]的數組，這正好是您在squeeze()後所需要的數組。更新：按照@Dev-iL's comment你應該使用permute而不是squeeze擺脫領先的單身維度，因爲這可能更有效。

其他兩個術語可以用同樣的方式處理，當然如果你基於重疊和不重疊的部分重構積分，它當然可能有幫助。