I have a program to calculate the evolution of a system of N particles with equal mass attracted by gravity. When I execute a profile of the program, I get the first 7 functions that take the longest are for the calculation of the acceleration of the system (91% of the time execution), which is okay because it is the hardest part of the problem. However, I want to know if there is a better way to do it.
My acceleration program is as follows:
function [A] = aceleracion(R)
Rij=R-permute(R,[2,1,3]);
A = sum(((vecnorm(Rij,2,3) + 0.05).^-3).*Rij , 2);
endfunction
Here, R is the variable where the data of the position is saved and it has 3 dimensions: rows separate particles, columns separate timesteps, and the third dimension is for xyz coordinates.