Matlab and Fortran Precision Issue

Question

SA 2021-10-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1563821-matlab-and-fortran-precision-issue

编辑： SA 2021-10-15

I am running a code in Matlab and fortran 90 but I get different results althought the codes are the same. Is this due to different precisions in the languages?

My code is posted below

XSTART = 2.0;
EPSA = 1.0;
EPSW = 80.0;
BULK_STRENGTH = 9.42629*1.0;
KAPPA = 8.486902807*BULK_STRENGTH/EPSW; 
VK = sqrt(KAPPA);
EC2_KBT = (332.06364/0.5921830)*1.0;
F1 = 1.1682217268107287;
UNC1 = F1 - ((EC2_KBT*1.0)/(EPSA*XSTART));
FREE_ENERGY = (0.50)*(1.0)*(UNC1)*(0.5921830*1.0);
ORIGINAL_FE = (0.50)*(1.0^2)*(332.06364)*(0.50)* ...
              (1.0/((VK*XSTART+1.0)*EPSW) - 1.0/EPSA)
          abs(FREE_ENERGY-ORIGINAL_FE);

for ORIGINAL_FE, I get -82.670010385176070 in matlab but -82.670007683885615 in fortran 90 and I am not sure why. My fortran code is posted below (you still get the results I had using implicit double precision (A-H,O-Z)

          PROGRAM MIB_HDM
          IMPLICIT real*8 (A-H,O-Z)
          
          REAL*8 :: EPSW,VK,XSTART
          REAL*8 :: EC2_KBT,KAPPA
          REAL*8 :: UNC1,BULK_STRENGTH
          REAL*8 :: ORIGINAL_FE
          REAL*8 :: EPSA
          
          XSTART = 2.0
          EPSA = 1.0
          EPSW = 80.0
          BULK_STRENGTH = 9.42629*1.0
          KAPPA = 8.486902807*BULK_STRENGTH/EPSW
          VK = sqrt(KAPPA)
          EC2_KBT = (332.06364/0.5921830)*1.0
          F1 = 1.1682217268107287
          
          
          UNC1 = F1 - ((EC2_KBT*1.0)/(EPSA*XSTART))
          FREE_ENERGY = (0.50)*(1.0)*(UNC1)*(0.5921830)
          ORIGINAL_FE = (0.50)*(1.0**2)*(332.06364)*(0.50)* &
          (1.0/((VK*XSTART+1.0)*EPSW) - 1.0/EPSA)
          print *, original_fe
          
          end program

Any explantion will be greatly appreciated. Thanks

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2021-10-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1563821-matlab-and-fortran-precision-issue#answer_809026

Hello SA

It's the constants. First of all for simplicity I commented out everythng not concerned with calculation of ORIGINAL_FE. I don't have a fortran complier but I tried to simulate what is going on by replacing three constants with double(single(constant)). If you uncomment the three double(single)) lines in the code below, you get exactly the fortran result. I was surprised to get exactly that result to all decimal places, but that's what happens.

XSTART = 2.0;
EPSA = 1.0;
EPSW = 80.0;
BULK_STRENGTH = 9.42629*1.0;
% BULK_STRENGTH = double(single(BULK_STRENGTH));
c1 = 8.486902807; 
% c1 = double(single(c1));
KAPPA = c1*BULK_STRENGTH/EPSW; 
VK = sqrt(KAPPA);
             % EC2_KBT = (332.06364/0.5921830)*1.0;
             % F1 = 1.1682217268107287;
             % F1 = double(single(F1))
             % UNC1 = F1 - ((EC2_KBT*1.0)/(EPSA*XSTART));
             % FREE_ENERGY = (0.50)*(1.0)*(UNC1)*(0.5921830*1.0);
c2 = 332.06364;
% c2 = double(single(c2));
ORIGINAL_FE = (0.50)*(1.0^2)*c2*(0.50)* ...
    (1.0/((VK*XSTART+1.0)*EPSW) - 1.0/EPSA)

Incidentally it's not needed for this calculation but your fortran version there is no declaration of F1 as REAL*8.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

SA 2021-10-15

编辑：SA 2021-10-15

Thanks so much David. You are right and with this eplanation I was able to make some changes to my fortran code and had the same value (-82.670010385176070) as the matlab code (without uncommenting the double(single) expressions). The D0's attatched to the numbers in fortran indicate double precision. If you take out the D0's in fortran then you have to uncomment the double(single) expressions in matlab to get the same result (-82.670007683885615). Below is the new fortran code to get the same result as matlab. Thanks so much for the help and explantion. Greatly appreciated.

PROGRAM MIB_HDM
IMPLICIT real*8 (A-H,O-Z)
REAL*8 :: EPSW,VK,XSTART
REAL*8 :: EC2_KBT,KAPPA
REAL*8 :: UNC1,BULK_STRENGTH
REAL*8 :: ORIGINAL_FE
REAL*8 :: EPSA,F1
XSTART = 2.D0
EPSA = 1.D0
EPSW = 80.D0
BULK_STRENGTH = 9.42629D0
KAPPA = 8.486902807D0*BULK_STRENGTH/EPSW
VK = sqrt(KAPPA)*1.D0
EC2_KBT = (332.06364D0/0.5921830D0)
F1 = 1.1682217268107287D0
UNC1 = F1 - ((EC2_KBT*1.D0)/(EPSA*XSTART))
FREE_ENERGY = (0.5D0)*(1.D0)*(UNC1)*(0.5921830)
ORIGINAL_FE = (0.5D0)*(1.D0**2)*(332.06364D0)*(0.5D0)* &
(1.D0/((VK*XSTART+1.D0)*EPSW) - 1.D0/EPSA)
print *, original_fe
end program

请先登录，再进行评论。

Answer 2

John D'Errico 2021-10-14

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1563821-matlab-and-fortran-precision-issue#answer_808761

No.

In the most common case, such a disagreement is because you copied over some numbers from Fortran to MATLAB (or the other way around) but only used a limited number of significant digits.

For example, we see computations like this: (332.06364/0.5921830)*1.0. Are those the EXACT values used in both cases? Do you know those to be EXACTLY the numbers that were used in both cases?

The reason i suggest this, is because the results that you show as different are the same, down to about the same number of significant digits. It is a very common mistake that we see.

The precision of MATLAB and FORTRAN is the same. They will both be using the same IEEE standard form for double precision variables. So assuming that those variables were declared in Fortran to be doubles, the precision will be the same. In MATLAB, the default is a double. And this brings up a second possibility. Could you have you declared something to be a single in the Fortran code? We don't know.

3 个评论
显示 1更早的评论隐藏 1更早的评论

John D'Errico 2021-10-14

I don't have your code, as written in BOTH languages to see if you are telling the truth, or just what you THINK to be true. However, when I see things written like:

BULK_STRENGTH = 9.42629*1.0;

or

ORIGINAL_FE = (0.50)*(1.0^2)*(332.06364)*(0.50)* ...

it very much implies you are a novice, and therefore easily subject to the type of mistakes I have alluded to. The presence of many mutiplications by 1.0 or even 1 squared, suggests a novice to programming at work.

My expectation is that one (or more) of those numbers is not known to be exactly what you represent it as.

The fact is, this is NOT a question of precision carried, IF both languages are using double precision. Nothing you have done is even remotely close to creating a precision problem.

SA 2021-10-14

The 1.0's in the code are from fortran which was 1.0d0 and it's not because I am now learnign to code in matlab. Should you try out the code yourslef by removing the 1.0's and trying it out in a fortran compiler with double precision you will still get the same result. You can try it out and if your reults are the same then I will agree I am wrong but until then the results I posted are true from my computation.

请先登录，再进行评论。