Is there a way to calculate pairwise vector displacement rather than using nested for loops?

Question

I need to calculate the force components on every particle in an N particle system. I can easily find the magnitude of the distance between any two particles using the pdist function, but I also need to compute the xyz displacement between each pair of particles. Pdist only returns magnitudes and does not indicate the direction of displacement. Is there an efficient way to calculate the component wise displacement between vector pairs without having to use nested for loops?

I've used nested for loops to compute the vector displacement between particle pairs, but this is very slow.

function gravity(obj, G)

        obj.rho = squareform(1./((pdist(obj.state(:,1:3))).^3));
        obj.rho = tril(obj.rho) + triu(obj.rho);
        for i = 1:3
           obj.delta_xyz(:,:,i) = squareform(pdist(obj.state(:,i)));
           obj.F_xyz(:,i) = -sum(obj.rho.*obj.delta_xyz(:,:,i), 2);

The above code runs much faster than the nested for loop solution, but pdist computes the magnitudes of the component displacement, and so the direction of the force calculation on each particle is incorrect.

Iban Cereijo Iban Cereijo · Accepted Answer · 2019-05-12T21:00:38

In general, with P having the positions of N particles in your desired shape Nx3:

N = 10;
P = rand(N, 3);

you can generate all the index pairs with meshgrid, keeping the unique pairs i < j:

[i, j] = meshgrid(1:N, 1:N);
mask = i < j;
i = i(mask);
j = j(mask);

and compute the displacements:

D = P(i, :) - P(j, :);

The arrangement should be equivalent to that of pdist, which is equivalent to:

dist = sqrt(sum(D.^2, 2)).'

Is there a way to calculate pairwise vector displacement rather than using nested for loops?

1 Answers