Matlab and Fortran Precision Issue
7 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
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.
More Answers (2)
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
John D'Errico
on 14 Oct 2021
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.
See Also
Categories
Find more on Fortran with MATLAB in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!