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?
7 个评论
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
2021-4-6
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
2021-4-6
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
2021-4-6
编辑：Yi-xiao Liu
2021-4-6
DGM
2021-4-6
Awesome
Mohammad Sami
2021-4-6
Great.
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
2021-4-6
2021-4-6
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
