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

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?

[DGM]

I think it might be valuable to know which of the two results is the correct one. Is there a simple test case that can be used? Consider:

n=3;
r=[1 1 1; 0 0 0; -1 -1 -1]
e=rand([n,1]);
m=rand([n,1]);
Q_pre=rand;
A=InterElectric(e,m,r,Q_pre)
B=InterEBencht(e,m,r,Q_pre)
A-B

gives essentially zero difference between the two. I'm wondering if there's a symmetry problem, though for n=2, the results match for any r.

[Yi-xiao Liu]

Try this condition, frankly I would think they both should be correct, probably I should first do it by hand:

n=3;
r=50*[1,1,1;0,1,0;1,-1,-1];
e=100*ones([n,1]);
m=ones([n,1]);
Q_pre=1;

[DGM]

I'm not so sure about that. Consider

r=[1 0 0; 0 0 0; 0 1 0]
e=ones([n,1]);
m=ones([n,1]);
Q_pre=1;

gives this for the loop method

B =
   1.353553390593274   0.646446609406726   1.000000000000000
  -1.000000000000000  -2.000000000000000  -1.000000000000000
  -0.353553390593274   0.646446609406726  -0.353553390593274

and gives this for the vectorized method

A =
   1.353553390593274  -0.353553390593274                   0
  -1.000000000000000  -1.000000000000000                   0
  -0.353553390593274   1.353553390593274                   0

Since the location vectors share the same z-component, I would expect that fz would be zero. Also, notice the nice symmetry of the x and y components of A. If I understand the math, that would suggest that the loop method is giving incorrect results.

Is your goal to have both versions working, or do you just want the vectorized version?

[Mohammad Sami]

Can you check if this is this a mistake in your for loop

f(i,:)=f(i)+(Qij/Rij^2)*(rij/Rij); % used f(i) instead of f(i,:)
%        ^

f(i,:)=f(i,:)+(Qij/Rij^2)*(rij/Rij); % change f(i) to f(i,:)

[Yi-xiao Liu]

Mohammad, you are right!

I missed the index again, so annoyed at myself

After correcting the index in the for loop, now both version returns the same and correct (I believe) results. Thank you guys!

[DGM]

Awesome

[Mohammad Sami]

Great.