首页 > 解决方案 > Matlab 和 fortran 精度

问题描述

我在 Matlab 和 fortran 90 中运行代码,但尽管代码相同,但我得到不同的结果。这是由于语言的精度不同吗?我的代码发布在下面

    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);

对于 ORIGINAL_FE,我在 matlab 中得到 -82.670010385176070,但在 fortran 90 中得到 -82.670007683885615,我不知道为什么。我的 fortran 代码发布在下面(您仍然可以得到我使用隐式双精度 (AH,OZ) 的结果

          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

标签: matlabfortran

解决方案


非常感谢你们的帮助。我想出了我的问题的答案。下面是我的 matlab 和 fortran 代码以获得相同的结果。第一个是fortran代码,第二个是matlab代码。如果您运行这两个代码,您应该得到相同的结果(-82.670010385176070)。如果您在 matlab 中取消注释双(单)表达式,那么您必须在 fortran 中取出所有D0以获得与 matlab 相同的结果(-82.670007683885615)。Fortran 代码

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代码

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)

推荐阅读