Order of elements inside sqrt/exponentiation changes ode output

I am running a code that runs an ode in which the equations of motion are exactly the same for both cases (except the line where I am making the little change). With the same input, I obtain different output.

Here's a simplified working version of the main code. The only difference between solver_ex1 and solver_ex2 is that in the calculation of r23, I use (x(1)+mu-1)^2 versus (x(1)-1+mu)^2. mu is a constant.

mu = 1.21506683e-2;
x0 = [0.8 0.3 0.1 0 0.2 0];
tspan = linspace(0,200,1000);

options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[~,y1] = ode113(@(t,y) solver_ex1(t,y,mu),tspan,x0,options);
[~,y2] = ode113(@(t,y) solver_ex2(t,y,mu),tspan,x0,options);

figure; hold on;
plot(y1(:,1),y1(:,2),'g-','LineWidth',0.5);
plot(y2(:,1),y2(:,2),'m--','LineWidth',0.5);

function xout = solver_ex1(~,x,mu)

% Distances (in r^(3/2) format)
r13 = ( (x(1)+mu-0)^2 + x(2)^2 + x(3)^2 )^1.5;
r23 = ( (x(1)+mu-1)^2 + x(2)^2 + x(3)^2 )^1.5;

% New state (nondimensional equations of motion)
dx = [  x(4)
        x(5)
        x(6)
        x(1) + 2*x(5) - (1-mu)*(x(1)+mu)/r13 - mu*(x(1)-1+mu)/r23
        x(2) - 2*x(4) - (1-mu)*x(2)/r13      - mu*x(2)/r23
        0             - (1-mu)*x(3)/r13      - mu*x(3)/r23];

xout(1:6) = dx(1:6);
xout = xout';

end

function xout = solver_ex2(~,x,mu)

% Distances (in r^(3/2) format)
r13 = ( (x(1)+mu-0)^2 + x(2)^2 + x(3)^2 )^1.5;
r23 = ( (x(1)-1+mu)^2 + x(2)^2 + x(3)^2 )^1.5;

% New state
dx = [  x(4)
        x(5)
        x(6)
        x(1) + 2*x(5) - (1-mu)*(x(1)+mu)/r13 - mu*(x(1)-1+mu)/r23
        x(2) - 2*x(4) - (1-mu)*x(2)/r13      - mu*x(2)/r23
        0             - (1-mu)*x(3)/r13      - mu*x(3)/r23];

xout(1:6) = dx(1:6);
xout = xout';

end

The output differs in each case. Plotting them as in the code also gives different visual trajectories after a while. That is, they are similar in the beginning, but start to diverge later on.

I get that it might be the tolerance, and the fact that it goes near the singularity (for r23), but shouldn't MATLAB deal with both cases the same way?

I am a bit confused at the moment. Could you give me a hint on what's going on?

Answer 1

[dpb]

Addition is mathematically associative so that x-1+mu == x+mu-1 identically but that isn't necessarily an identity in floating point; it is not guaranteed that operations are completely associative to the last bit (as you've discovered another case thereof).

An optimizing compiler might choose to rearrange the two expressions into the same identical code for execution, but MATLAB for the most part as an interpreting language generally operates on expressions as it finds them. Hence, it's a case where the order of the two operations does, in fact, cause a difference in rounding at some point in the calculation and then that difference eventually propogates to large enough difference to be observable.

[Mitsu]

Thanks you, dpb! I will keep this in mind.