Different results by changing the order of operations

Hi all,

Replicating the cardiac nodal cell model (set of ODEs) from MATLAB (2018a, 64bit) to Simulink I found different results in solution of a particulal equation (without derivatives). Later, I also found the same difference in Matlab results either by changing the order of mathematical operations (e.g. (x./y).*z ~= x.*z./y), or calculating parts of the equation separately first, and then combining them. Example of the latter case is below (see also attached file). The problem appears in equation in the line 451. Temporal variables al57 and al571 show different results, though the difference is not too big. For the first time step, 1.1324 and 1.0209, respectively, and the second value seems correct being checked with calculator, with the same (2nd) result in Simulink.

A(:,57) = 1.0+ (S(:,6)./C(:,19)).*(1.0+exp(( - C(:,16).*S(:,1))./C(:,49))+C(:,14)./C(:,23))+ C(:,14)./C(:,20)+power(C(:,14), 2.0)./ (C(:,20).*C(:,21))+power(C(:,14), 3.0)./( C(:,20).*C(:,21).*C(:,22));

al57 = A(:,57)

al570 = (1.0+exp(( - C(:,16).*S(:,1))./C(:,49))+C(:,14)./C(:,23))+C(:,14)./C(:,20)+power(C(:,14), 2.0)./ (C(:,20).*C(:,21))+ power(C(:,14), 3.0)./( C(:,20).*C(:,21).*C(:,22));

al571 = 1.0+ (S(:,6)./C(:,19)).*al570

-----

al57 =

1.1324

al571 =

1.0209

Can anyone comment or give some advice? Thanks in advance.

[Stephen23]

"..difference in Matlab results either by changing the order of mathematical operations"

As operations on binary floating point numbers are neither associative nor commutative, why do you particularly expect the outputs to be exactly the same? Floating point mathematics is definitely not the same as symbolic or algebraic mathematics.

Search this forum for "floating point" for hundreds of threads disscussing this topic.

Most likely your algorithm is numerically unstable, e.g. incudes differences of very similar numbers. If you explained the algorithm them someone might be able to help you with that.

[Maxim]

Hi Stephen,

Obviously, I am familiar with binary floating point operations and possible difference in results with changing order of operations. Usually, the difference is in 5th-6th decimal position. But in my example shown above, “not too big” meant 5-10% error (sorry for unclear description), sometimes acceptable for general computer simulation results, but for my case – 10% difference in oscillation period.

The result in al571 seems to be correct, as it is confirmed by a calculator and Simulink (with any order of operstions), which, I think, use the same floating point operations, except working with real numbers, not matrixes. Moreover, below there is a result of changing the order of operations I mentioned above - (x./y).*z -> (x.*z)./y (moved C(:,19) to the end of the expression)

al572 = 1.0+ (S(:,6).*(1.0+exp(( - C(:,16).*S(:,1))./C(:,49))+C(:,14)./C(:,23))+ C(:,14)./C(:,20)+ power(C(:,14), 2.0)./( C(:,20).*C(:,21))+power(C(:,14), 3.0)./( C(:,20).*C(:,21).*C(:,22)))./C(:,19)

-------

al57 =

1.1324

al571 =

1.0209

al572 =

6.4486

Here the difference is large. Upgraded to version 2019b, but results are exactly the same.

Answer 1

[James Tursa]

You haven't broken up the calculations properly. I.e., you are comparing different calculations. Your code is essentially:

    A(:,57) = 1.0+ (S(:,6)./C(:,19)).*(1.0+exp((  - C(:,16).*S(:,1))./C(:,49))+C(:,14)./C(:,23))+ C(:,14)./C(:,20)+power(C(:,14), 2.0)./ (C(:,20).*C(:,21))+ power(C(:,14), 3.0)./( C(:,20).*C(:,21).*C(:,22));
    al570 =                           (1.0+exp((  - C(:,16).*S(:,1))./C(:,49))+C(:,14)./C(:,23))+ C(:,14)./C(:,20)+power(C(:,14), 2.0)./ (C(:,20).*C(:,21))+ power(C(:,14), 3.0)./( C(:,20).*C(:,21).*C(:,22));
    al57  = A(:,57);
    al571 =   1.0+ (S(:,6)./C(:,19)).*al570;

In the A(:,57) line and thus the al57 calculation, the factor (S(:,6)./C(:,19)) is only applied to the following sub-expression:

    (1.0+exp((  - C(:,16).*S(:,1))./C(:,49))+C(:,14)./C(:,23))

But in the al570 and al571 lines, the factor (S(:,6)./C(:,19)) gets applied to this entire expression:

    (1.0+exp((  - C(:,16).*S(:,1))./C(:,49))+C(:,14)./C(:,23))+ C(:,14)./C(:,20)+power(C(:,14), 2.0)./ (C(:,20).*C(:,21))+ power(C(:,14), 3.0)./( C(:,20).*C(:,21).*C(:,22))

You can see this by going into the editor and clicking on the opening parenthesis in the (1.0+exp... and seeing that the associated closing parenthesis is at the ...C(:,23)) part of the expression and not at the end of the line.

You need to correct your code and make sure it is calculating your source equations properly. Your differences are not due to changing the order of calculations, but are caused by doing different calculations altogether.

[Maxim]

OMG... You are right. Thank you!

Answer 2

[Etsuo Maeda]

"Validated numerics" will help you to understand your problem.

https://en.wikipedia.org/wiki/Validated_numerics

My best guess is that "vpa" in Symbolic Math Toolbox will show you a proper result.

HTH