Why does this vectorized version does return the same as a for loop?

2 次查看（过去 30 天）

显示更早的评论

Yi-xiao Liu 2021-4-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/793552-why-does-this-vectorized-version-does-return-the-same-as-a-for-loop

评论： Mohammad Sami 2021-4-6

在 MATLAB Online 中打开

I am trying to sum the inter-electrostatic forces of n points at r=[rx,ry,rz]

However, the vectorized version does not return the same as the loop

rng("shuffle")
n=5;
r=10*rand([n,3])-5;
e=rand([n,1]);
m=rand([n,1]);
Q_pre=rand;
InterElectric(e,m,r,Q_pre)-InterEBencht(e,m,r,Q_pre)
>>ans=
 2.16840434497101e-19	  0.00235284634075158	  0.00929194817544228
 2.08166817117217e-17     0.0990903618666820     -0.102663948053765
-3.46944695195361e-18    -0.0270906896509036	 -0.0250703363558977
-1.38777878078145e-17	 -0.0505941945906215	  0.142937845957793
 1.08420217248550e-19	 -0.000650770705875512	 -0.000346897691602804
function f=InterEBencht(e,m,r,Q_pre)
for i=1:length(e)
    f(i,:)=[0,0,0];
    for j=1:length(e)
        if i~=j
            Qij=e(i)*e(j)/m(i)*Q_pre;
            rij=r(i,:)-r(j,:);
            Rij=vecnorm(rij);
            f(i,:)=f(i)+(Qij/Rij^2)*(rij/Rij);
        end
    end
end
end
function f=InterElectric(e,m,r,Q_pre)
x=r(:,1);
y=r(:,2);
z=r(:,3);
xij=x-x';
yij=y-y';
zij=z-z';
Qij=e.*e'.*Q_pre./m;
Rij=cat(3,xij,yij,zij);
Rij=vecnorm(Rij,2,3);
fMag=Qij./(Rij.^3);
fx=sum(fMag.*xij,2,"omitnan");
fy=sum(fMag.*yij,2,"omitnan");
fz=sum(fMag.*zij,2,"omitnan");
%When i=j, Rij=0, xij=0, xij/Rij=nan,so omitnan ensure that the force a
%bead has on itself(i=j) is not summed, same for y and z
f=[fx,fy,fz];
f(isnan(f))=0;
end

What scratches my head is that the x component of f is the same between 2 versions, but the y and z component is not despite they are computed using the same formula as x. Can anyone help me?