Why does this vectorized version does return the same as a for loop?
2 views (last 30 days)
Show older comments
I am trying to sum the inter-electrostatic forces of n points at r=[rx,ry,rz]
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/573692/image.png)
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?
7 Comments
Answers (0)
See Also
Categories
Find more on Sparse Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!