Matlab and fortran precision

Question

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.1682185947500601;
    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

Does this answer your question? [Precision not respected](https://stackoverflow.com/questions/30827276/precision-not-respected). There's probably a better dupe target out there somewhere, but I can't seem to find it. Basically you're declaring your variables as `real*8`, but then initialising them with default precision constants, e.g. `80.0` rather than the `real*8` precision constant `80.0_8`. Also, `real*8` is not very portable. See [this question](https://stackoverflow.com/questions/22362211/confusing-double-precision-real-in-fortran) for details. — veryreverie, Oct 14 '21 at 19:27
Thanks so much for the suggestion but it doesn't answer my question. — user3307971, Oct 14 '21 at 19:46
I tried what you just suggested @veryreverire but my results didn't change — user3307971, Oct 14 '21 at 19:49
In that case show your updated code and its output. Did you really apply the change to *all* those numbers? — Vladimir F Героям слава, Oct 14 '21 at 20:40
If I make the changes suggested in the duplicates, and fix the broken Fortran ( `real*8`, trailing comma, IMO not using `Implicit None`) I get exactly the same answer as you report for MATLAB. Voting to close as a duplicate until the OP explains why they didn't address their problem. — Ian Bush, Oct 15 '21 at 07:32

user3307971 · Accepted Answer · 2021-10-15T15:05:29.023

Thanks so much guys for your help. I figured out the answer to my question. Below are my matlab and fortran codes to get the same results. The fortran code is the first one and the matlab code is the second one. If you run the two codes you should get the same result ( -82.670010385176070). If you uncomment the double(single) expressions in matlab then you have to take out all of the D0's in fortran to get the same result as matlab ( -82.670007683885615). Fortran Code

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

Matlab code

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)

Matlab and fortran precision

1 Answers1