Function in parfor loop produces incorrect output, in normal for-loop output is correct

I believe I am being hampered by the knowledge of how a parfor loop work, hopefully someone knowns in which direction to point me.

The problem: I have a script which uses a function to calculate the crash location of an object with {x,y,z} and {vx,vy,vz} coordinates and speeds. I use the simple script below to evaluate the effect of a decreasing time-step dt on the performance and precission. Put in a normal for-loop, the output is as expected. (I understand that I write over my output in the main loop, but ignore that for now).

for dt = [0.1 0.01 0.001 0.0001 0.00001 0.000001]
    for l = 1:10  
        [x(l), y(l)] = f_crash_location(0,0,100,10,0,0,0.0038,dt);
    end
end

The output is then x = 41.244... and y = 0, which is physically correct :). However, in the main simulation where I am using the function, l is not 10, but more like a couple million runs, so I thought, perfect for parallelization. But when I change the for to a parfor-loop like this:

for dt = [0.1 0.01 0.001 0.0001 0.00001 0.000001]
    parfor l = 1:10  
        [x(l), y(l)] = f_crash_location(0,0,100,10,0,0,0.0038,dt);
    end
end

The code is executed insanely fast, which seems nice, but the output is useless. x = 10^-5 is plainly wrong. (I figured out the parfor loop only does dt*x_dot once.)

The code of the function f_crash_location.m is:

 function [x_n,y_n] = f_crash_location(x_n,y_n,z_n,x_dot_n,y_dot_n,z_dot_n,C_D,dt)
    %%Declare variables
    g = 9.807;

    %%Numerical model
    while z_n > 0
        % Vector velocity
        V = sqrt(x_dot_n^2 + y_dot_n^2 + z_dot_n^2);

        % Accelerations
        x_2dot_n = -V*x_dot_n*C_D;
        y_2dot_n = -V*y_dot_n*C_D;    
        z_2dot_n = -g + -V*z_dot_n*C_D;

        % Velocities
        x_dot_n = x_dot_n + x_2dot_n*dt;
        y_dot_n = y_dot_n + y_2dot_n*dt;
        z_dot_n = z_dot_n + z_2dot_n*dt;

        % Positions
        x_n = x_n + x_dot_n*dt;
        y_n = y_n + y_dot_n*dt;
        z_n = z_n + z_dot_n*dt;
    end
end

After playing around a bit I believe it has something to do with how variables are used differently in parfor then in for-loops, but I am stuck on how to proceed. What I am doing wrong that the parfor loop spits out, well, that brown stuff?

[Stephen23]

I cannot replicate the behavior on R2012b using default settings:

for dt = [0.1 0.01 0.001 0.0001 0.00001 0.000001]
    parfor l = 1:10
        [x(l), y(l)] = f_crash_location(0,0,100,10,0,0,0.0038,dt);
    end
    x,y % display output
end

Where parfor prints this:

x =
   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446
y =
     0     0     0     0     0     0     0     0     0     0

And for prints this:

x =
   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873   40.8873
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553   41.2553
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401   41.2401
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444   41.2444
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446
y =
     0     0     0     0     0     0     0     0     0     0
x =
   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446   41.2446
y =
     0     0     0     0     0     0     0     0     0     0

They look the same to me. What MATLAB version are you using?

[Wouter Grimme]

Thank you for trying to help me btw. I think I have pinpointed the problem:

In my main script I use:

global g;
g = 9.807;

And then in my function I use:

global g;

to call the variable again. If I do it like this, I get the behavior I was describing, when I rewrite the function file to:

g=9.807;

Then it works like a charm.

Probably I use global variables in a wrong way together with a par-for loop?

[Stephen23]

"Probably I use global variables in a wrong way together with a par-for loop?"

Global variables are the wrong way to do anything.

• If you want to know what experienced programmers think about globals, then search the 'net.
• If you want to know what experienced MATLAB users think of globals, then search this forum.
• If you want to know what the writers of MATLAB think about globals, then read their blogs.
• If you want to know why they all think that, then look at your own question: it is a classic example of how globals cause bugs that are hard to track down.

Simple rule to make your life easier and waste less time chasing down hard-to-find bugs: do NOT use global variables. You might find this useful, it describes much better alternatives to global variables:

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Wouter Grimme]

Really learned a few things today, thank you for your help and explanation Stephen!

Answer 1

[Wouter Grimme]

Short answer:

Do not use global variables!