Avoiding operations on the same component in 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

[Chad Greene]

Is there a reason you can't skip the loops altogether and use gradient or diff? And an aside, just a good habit to get into: avoid using i and j as counters or variables, because sometimes Matlab thinks they're the imaginary unit. Doesn't cause problems very often, but when it does, it can be a real headache.

[Sean de Wolski]

Any reason why you want to get rid of the loops?

[Matt J]

just a good habit to get into: avoid using i and j as counters or variables, because sometimes Matlab thinks they're the imaginary unit.

A better habit is to stop using i and j for the imaginary unit and start using 1i and 1j

    >> for i=1:3, 10*i+1i, end

    ans =

      10.0000 + 1.0000i

    ans =

      20.0000 + 1.0000i

    ans =

      30.0000 + 1.0000i

[Fredrik Carlsson]

Is there a reason you can't skip the loops altogether...

Other than a habit from the Java language where vector operations are limited, no not really, we only had a 45min lecture of the basics in MatLAB in which they stressed the point on trying to avoid for loops as much as possible.

...avoid using i and j as counters or variables, because sometimes Matlab thinks they're the imaginary unit.

This never occured to me, had momentarily forgotten that MatLAB is capable of handling imaginary numbers. =)

Any reason why you want to get rid of the loops?

At first I had three nested loops and managed to get rid of one of them, this increased the speed by a factor of 2, when the plotting was commented out. The loop was only over three components and was between the two loops in the code, I hope that by getting rid of the other two, or at least one more, it will increase the speed even more.

... use gradient or diff?

I will look into this, from a quick glance over the documentation it seems to be what I want. Thanks

Answer 1

[Fredrik Carlsson]

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

Answer 2

[Jan]

This wastes time:

dxdt = zeros(size(W(:,1:3)));

The temporary array "W(:, 1:3)" is created, but not used. Better:

dxdt = zeros(size(W, 1), 3);

norm is not really fast. Look in the FileExchange for "DNorm2". But your vectors are tiny and a simple sqrt(sum(x.*x)) might be more efficient.

Do not calculate "W(i,1:3)-W(j,1:3)" twice. Do this once only and store it in a variable.

[Fredrik Carlsson]

While this optimizes the code to some extent, it does not really answer the question to whether or not there is a way to do the same operations without the loops.