MATLAB: Avoiding operations on the same component in nested for loops

nested for loops

I'm trying to optimize an implementation of Newtons gravity law using a Runge-Kutta of the fourth order for an assignment in a course I'm taking, so far it works as it should but is rather slow when it is to plot/make video the orbits of the objects. Aside from using a different numerical integrator, as I am aware of the problems of Runge-Kutta, but the assignment was specific about implementing the method, I have tried to take away as many loops as I can during the calculations of the derivatives.
So far I have two nested for loops in a function that will calculate the derivatives as shown below. The W argument is a matrix containg the objects location and velocity in x y z direction on a row as,
object1 (x,y,z,vx,vy,vz);
object2 (x,y,z,vx,vy,vz); …
i.e. each row is a planet or star in this case.
The Gm parameter is an array containing the masses of the objects multiplied by G.
The t should be a scalar but is not used.
As shown in the code the dvdt operation is only to be done when i~=j since it would otherwise mean division with 0. I'm wondering if there is a way to get rid of the second loop, I have tried it by using vector operations as MatLAB is constructed to do, but I cannot find a way to make sure that if statement is included either explicitly or implicitly.
 function dwdt = g2(W,Gm,t)
 dxdt = zeros(size(W(:,1:3)));
 dvdt = zeros(size(W(:,1:3)));
 for i=1:size(W,1) %Iterate over number of objects

     dxdt(i,:) = W(i,4:6);
     for j=1:size(W,1) %Iterate over number of objects
         if j~=i
             r3 = norm(W(i,1:3)-W(j,1:3)).^3;
             dvdt(i,:) = dvdt(i,:) - Gm(j).*(W(i,1:3)-W(j,1:3))./r3;
         end
     end
 end
 dwdt = [dxdt dvdt];
 end

Best Answer

I mananged to get rid of the loop using bsxfun function and made it faster, though not by a ridiculus amount, I thought I should post my solution here
 function dwdt2 = g2(W,Gm,t)
 dvdt = zeros(size(W(:,1:3)));
 for i=1:size(W,1) %Iterate over number of objects
    diff = bsxfun(@minus,W(i,1:3),W(:,1:3));
    r3 = sqrt(sum(diff.*diff,2)).^3;
    b = bsxfun(@times,Gm,diff);
    c = bsxfun(@rdivide,b,r3);
    c(isnan(c)) = 0;
    dvdt(i,:) = - sum(c,1);
 end
 dwdt2 = [W(:,4:6) dvdt];
 end
I suspect this is not an optimal use of bsxfun but it got rid of the inner loop and I solved the dividing by 0 problem which the if statement originally solved
Related Question