Matlab and Fortran Precision Issue

조회 수: 10 (최근 30일)
SA
SA 2021년 10월 14일
편집: 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

채택된 답변

David Goodmanson
David Goodmanson 2021년 10월 15일
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
SA
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

댓글을 달려면 로그인하십시오.

추가 답변 (2개)

John D'Errico
John D'Errico 2021년 10월 14일
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
John D'Errico
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
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.

댓글을 달려면 로그인하십시오.


Voss
Voss 2021년 10월 14일
The value of F1 is different in MATLAB vs Fortran.
  댓글 수: 1
SA
SA 2021년 10월 14일
Thanks Ben. It's an error and I have edited it but the results are still the same.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Fortran with MATLAB에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by