我正在为学校编写一个程序,我有嵌套的 for 循环来创建一个 4 维数组(坐标为 (x,y) 和 (x',y') 的两点之间的距离),如下所示:
pos_x=1:20;
pos_y=1:20;
Lx = length(pos_x);
Ly = length(pos_y);
Lx2 = Lx/2;
Ly2 = Ly/2;
%Distance function, periodic boundary conditions
d_x=abs(repmat(1:Lx,Lx,1)-repmat((1:Lx)',1,Lx));
d_x(d_x>Lx2)=Lx-d_x(d_x>Lx2);
d_y=abs(repmat(1:Ly,Ly,1)-repmat((1:Ly)',1,Ly));
d_y(d_y>Ly2)=Ly-d_y(d_y>Ly2);
for l=1:Ly
for k=1:Lx
for j=1:Ly
for i=1:Lx
distance(l,k,j,i)=sqrt(d_x(k,i).^2+d_y(l,j).^2);
end
end
end
end
d_x 和 d_y 只是 20x20 矩阵,并且 Lx=Ly 用于试验目的。它非常慢,而且显然不是一种非常优雅的方法。我尝试对嵌套循环进行矢量化,并成功摆脱了两个内部循环:
dx2=zeros(Ly,Lx,Ly,Lx);
dy2=zeros(Ly,Lx,Ly,Lx);
distance=zeros(Ly,Lx,Ly,Lx);
for l=1:Ly
for k=1:Lx
dy2(l,k,:,:)=repmat(d_y(l,:),Ly,1);
dx2(l,k,:,:)=repmat(d_x(k,:)',1,Lx);
end
end
distance=sqrt(dx2.^2+dy2.^2);
它基本上取代了上面的 4 个 for 循环。我现在已经尝试了两天,但找不到一种方法来矢量化所有循环。我想问:
- 是否有可能真正摆脱这两个循环
- 如果是这样,我将不胜感激任何提示和技巧。 到目前为止,我已经尝试在 4 维中再次使用repmat,但是你不能转置 4 维矩阵,所以我尝试在许多不同的组合中一起使用 permute 和 repmat,但无济于事。
任何建议将不胜感激。
<小时/>感谢您的回复。抱歉措辞不好,我基本上想要的是让一组振荡器均匀地位于 x-y 平面上。我想模拟它们的耦合,耦合函数是每个振荡器之间距离的函数。每个振荡器都有一个 x 和 y 坐标,所以我需要找到 osci(1,1)
之间的距离和osci(1,1),..osci(1,N),osci(2,1),..osci(N,N)...
然后 osci(1,2)
也一样和osci(1,1)...osci(N,N)
等等..(所以基本上是所有振荡器和所有其他振荡器之间的距离加上自耦合)如果除了使用 4-D 阵列之外还有更简单的方法来做到这一点,我也绝对想知道它。 .
最佳答案
如果我理解正确的话,振荡器遍布各处,如下所示:
然后你想计算振荡器 1 和振荡器 1 到 100 之间的距离,然后计算振荡器 2 和振荡器 1 到 100 之间的距离等。我相信这可以用 2D 距离矩阵来表示,其中第一个维度是1 到 100,第二个维度从 1 到 100。
例如
%# create 100 evenly spaced oscillators
[xOscillator,yOscillator] = ndgrid(1:10,1:10);
oscillatorXY = [xOscillator(:),yOscillator(:)];
%# calculate the euclidean distance between the oscillators
xDistance = abs(bsxfun(@minus,oscillatorXY(:,1),oscillatorXY(:,1)')); %'# abs distance x
xDistance(xDistance>5) = 10-xDistance; %# add periodic boundary conditions
yDistance = abs(bsxfun(@minus,oscillatorXY(:,2),oscillatorXY(:,2)')); %'# abs distance y
yDistance(yDistance>5) = 10-yDistance; %# add periodic boundary conditions
%# and we get the Euclidean distance
euclideanDistance = sqrt(xDistance.^2 + yDistance.^2);
关于matlab - 在 Matlab 中向量化 4 个嵌套 for 循环,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/6375743/