Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Vectorizing code to calculate radial distance for all points

Asked by Cong Bang Huynh on 12 Jan 2013

Hi all!

I have currently written this code to calculate the radial distance for all points with their x, y and z coordinates.

 for x = xmin:step:xmax
    for y = ymin:step:ymax
        for z = zmin:step:zmax
            r((x-xmin)/step+1,(y-ymin)/step+1,(z-zmin)/step+1) = sqrt(x^2 + y^2 + z^2);

I am wondering if there is a way to vectorize this code so that I don't have to run the nested loops and hence speed up my calculation time, because I have to run the code for quite a number of times to analyze different data sets.

I am new to matlab so this question may seem quite trivial but please do help me :) Thank you in advance.

1 Comment

Roger Stafford on 12 Jan 2013

Aside from the question of vectorizing code, it is unwise to set up for-loops in the way you show here. The problem is that expressions such as (x-xmin)/step+1 used as indices may not turn out to be exact positive integers due to tiny round-off errors and thereby would produce error messages. For example, with xmin = 0, step = 0.1, and xmax = 1, my machine produces a non-integer for (x-xmin)/step+1 at the seventh step and would complain if I used it as an index. With for-loops you should instead do something like the following so that you are guaranteed to have integer values for your indices:

 x = xmin:step:xmax;
 y = ymin:step:ymax;
 z = zmin:step:zmax;
 for kx = 1:length(x)
   for ky = 1:length(y)
     for kz = 1:length(z)
       r(kx,ky,kz) = sqrt(x(kx)^2 + y(ky)^2 + z(kz)^2);
Cong Bang Huynh


No products are associated with this question.

2 Answers

Answer by Matt J on 12 Jan 2013
Accepted answer
 x = (xmin:step:xmax).';
 y = ymin:step:ymax;
 z = reshape(zmin:step:zmax,1,1,[]);
 r=sqrt(bsxfun(@plus, bsxfun(@plus,x.^2,y.^2),z.^2 ));


Matt J
Answer by Roger Stafford on 12 Jan 2013

You can also do:

 [X,Y,Z] = ndgrid(xmin:step:xmax,ymin:step:ymax,zmin:step:zmax);
 r = sqrt(X.^2+Y.^2+Z.^2);


Roger Stafford

Contact us