Matlab and Fortran Precision Issue

6 views (last 30 days)
SA
SA on 14 Oct 2021
Edited: SA on 15 Oct 2021
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

Accepted Answer

David Goodmanson
David Goodmanson on 15 Oct 2021
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 Comment
SA
SA on 15 Oct 2021
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

Sign in to comment.

More Answers (2)

John D'Errico
John D'Errico on 14 Oct 2021
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 Comments
SA
SA on 14 Oct 2021
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.

Sign in to comment.


Benjamin
Benjamin on 14 Oct 2021
The value of F1 is different in MATLAB vs Fortran.
  1 Comment
SA
SA on 14 Oct 2021
Thanks Ben. It's an error and I have edited it but the results are still the same.

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!