如何在MatLAB中获得Fortran精度 [英] How to obtain Fortran precision in MatLAB

问题描述

我有一段用Fortran和Matlab编写的代码。他们做了完全相同的计算，即

1. 构建一个 tanh -field并找到它的Laplacian


2. 将一些术语相乘

这个乘法的结果产生一个矩阵， （4,4）th和（6,6）th减去。

• 在Fortran中，它们的区别是〜1e-20

这个问题非常关键，因为我测试这个数字是否小于零。
问题：是否有一种方法可以执行计算，以便在Matlab中获得与Fortran中相同的精度？
$ b

我列出以下代码：

---

MatLAB

 清除全部
 
权重= [4./9，1./9,1./9,1。 /9,1./9，1./36,1./36,1./36,1./36]; 
 dir_x = [0,1,0,1,0,1，-1，-1,1]; 
 dir_y = [0，0，1，0，-1，1，-1，-1，-1]; 
 
 
 
 %%%%%%%%%%%%%%%%%%%%%% CONSTANTS 
 length_y = 11; length_x = length_y; 
 y_center = 5; x_center = y_center; 
 
 
 densityHigh = 1.0; 
 densityLow = 0.1; 
 radius = 3.0; 
 c_width = 1.0; 
 
 average_density = 0.5 *（densityHigh + densityLow）; 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 
 
 
 
 for x = 1：length_x 
 for y = 1：length_y 
 for i = 1：9 
 fIn（i，x，y）=权重（ⅰ）* densityHigh; 
 test_radius = sqrt（（x-x_center）*（x-x_center）+（y-y_center）*（y-y_center））;如果（test_radius <=（radius + c_width））
 fIn（i，x，y）=权重（i）*（average_density_0.5 *（densityHigh-densityLow）* tanh（2.0 *（ radius-sqrt（（x-x_center）*（x-x_center）+（y-y_center）*（y-y_center））/ c_width）））; 
 end 
 end 
 end 
 end 
 
 ref_density_2d = ones（length_x）* average_density; 
 for i = 1：length_x 
 ref_density（：，：，i）= abs（ref_density_2d（:, i）'）; 
 end 
 
 
 rho = sum（fIn）; 
 laplacian_rho =（+ 1.0 *（circshift（rho（1，：，:)，[0，-1，-1]）+ circshift（rho（1，：，:)，[0， -1]）+ circshift（rho（1，：，:)，[0，-1，+1]）+ circshift（rho（1，：，:)，[0，+ 1，+ 1]））+ ... 
 + 4.0 *（circshift（rho（1，：，:)，[0，-1，+0]）+ circshift（rho（1，：，:)，[0， + 0]）+ circshift（rho（1，：，:)，[0，+ 0，-1]）+ circshift（rho（1，：，:)，[0，+ 0，+ 1]））+ ... 
 -20.0 * rho（1，：，:)）; （rho-densityLow）。*（rho-densityHigh）。*（rho-ref_density）-laplacian_rho *（1.851851851851852e-04）/6.0;（b-b-b_b）= 4.0 * 0.001828989483310 * 
 
 psi（1,4,4）-psi（1,6,6）
  

---

Fortran

  PROGRAM main 
 
隐式无
 
 INTEGER，PARAMETER :: DBL = KIND（1.D0）
 REAL（KIND = DBL），DIMENSION（1：11,1：11 ）:: psi，rho 
 INTEGER :: i，j，m，即iw，jn，js 
 REAL（KIND = DBL）:: R，rhon，lapRho 
 
 INTEGER，DIMENSION（1：11,1：11,1：4）:: ni 
 
 REAL（KIND = DBL）:: kappa，kappa_6，kappa_12，kappaEf，beta，beta4 
 
 
 
 beta = 12.D0 * 0.0001 /（1.D0 *（（1.0  -  0.1）** 4））
 kappa = 1.5D0 * 0.0001 * 1。 D0 /（（1.0  -  0.1）** 2）
 
 
！--------定义近邻并初始化密度rho --------- ------- 
 DO j = 1，11 
 DO i = 1，11 
 
！初始化密度
 rho（i，j）= 1.D0 
 R = DSQRT（（DBLE（i）-5.0）** 2 +（DBLE（j）-5.0）** 2 b （bLE（3.0） - （b））（b）（b）（b≤1） R）/1.D0）
 END IF 
 
！生成邻居数组
 ni（i，j，1）= i + 1 
 ni（i，j ，2）= j + 1 
 ni（i，j，3）= i  -  1 
 ni（i，j，4）= j  -  1 
 END DO 
 END DO 
 
 
！ （1，3）= 11 
 ni（11，：，1）= 1 
 ni（：，1,4）= 11 
 ni（：11,2）= 1 
 
 
 
！---------界面力的应力形式的微分项---- -  
 DO j = 1，11 
 DO i = 1，11 
 
！识别邻居
 ie = ni（i，j，1）
 jn = ni（i，j，2）
 iw = ni（i，j，3）
 js = ni（i，j，4）
 
！拉普拉斯算子的密度ρbρbρρρρ= 4.D 0 *（rho（即，j）+ rho（iw，j）+ rho（i，jn）+ rho（i，js））& 
 + rho（即jn）+ rho（即js）+ rho（iw，jn）+ rho（iw，js） -  20.D0 * rho（i，j）
 
！定义化学势Psi 
 psi（i，j）= 4.D0 * beta *（rho（i，j）-0.55）*（rho（i，j）-0.1）*（rho（i，j） ） -  1.0）& 
 -kappa * lapRho / 6.D0 
 END DO 
 END DO 
 
 
写（*，*）psi（6,6）-psi （4,4）
 
 
结束计划
  

解决方案

在整个代码中，您仍然没有使用双精度，例如：

  beta = 12。 D0 * 0.0001 /（1.D0 *（（1.0  -  0.1）** 4））
  

和还有很多。如果我强制编译器使用double精度作为float的默认值（对于 gfortran ，编译选项是 -fdefault-real-8
blockquote
0.00000000000000000000000000000000000


所以你需要修复你的代码。例如，引用的行应为：



  beta = 12.D0 * 0.0001D0 /（1.D0 *（ 1.0D0  -  0.1D0）** 4））
  

[虽然我鄙视符号 D0 ，但这是另一回事]

I have a piece of code written in Fortran and in Matlab. They do exactly the same calculation, namely

1. Construct a tanh -field and find its Laplacian
2. Multiply some terms together

The result of this multiplication yields a matrix, whose (4,4)th and (6,6)th I subtract.

• In Fortran their difference is ~1e-20
• In Matlab their difference is identically zero.

This issue is very critical, as I test if this number is less than zero. Question: Is there a way to perform the calculations such that I get the same precision in Matlab as in Fortran?

I list the codes below:

---

MatLAB

clear all

weights = [4./9, 1./9,1./9,1./9,1./9, 1./36,1./36,1./36,1./36];
dir_x   = [  0,   1,  0, -1,  0,    1,  -1,  -1,   1];
dir_y   = [  0,   0,  1,  0, -1,    1,   1,  -1,  -1];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONSTANTS
length_y = 11; length_x = length_y;
y_center = 5; x_center  = y_center;


densityHigh = 1.0;
densityLow  = 0.1;
radius  = 3.0;
c_width = 1.0;

average_density = 0.5*(densityHigh+densityLow);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





for x=1:length_x
    for y=1:length_y
        for i=1:9
            fIn(i, x, y) = weights(i)*densityHigh;
            test_radius = sqrt((x-x_center)*(x-x_center) + (y-y_center)*(y-y_center));
            if(test_radius <= (radius+c_width))
                fIn(i, x, y) = weights(i)*( average_density - 0.5*(densityHigh-densityLow)*tanh(2.0*(radius-sqrt((x-x_center)*(x-x_center) + (y-y_center)*(y-y_center))/c_width)) );
            end
        end
    end
end 

ref_density_2d = ones(length_x)*average_density;
for i=1:length_x
    ref_density(:,:,i) = abs(ref_density_2d(:, i)');
end


          rho = sum(fIn);
laplacian_rho = (+1.0*(circshift(rho(1,:,:), [0, -1, -1]) + circshift(rho(1,:,:), [0, +1, -1]) + circshift(rho(1,:,:), [0, -1, +1]) + circshift(rho(1,:,:), [0, +1, +1])) + ...
                 +4.0*(circshift(rho(1,:,:), [0, -1, +0]) + circshift(rho(1,:,:), [0, +1, +0]) + circshift(rho(1,:,:), [0, +0, -1]) + circshift(rho(1,:,:), [0, +0, +1])) + ...
                -20.0*rho(1,:,:));

psi   = 4.0*0.001828989483310*(rho-densityLow).*(rho-densityHigh).*(rho-ref_density) - laplacian_rho*(1.851851851851852e-04)/6.0;

psi(1,4,4)-psi(1,6,6) 

---

Fortran

 PROGRAM main

 IMPLICIT NONE

 INTEGER, PARAMETER :: DBL = KIND(1.D0)
 REAL(KIND = DBL), DIMENSION(1:11,1:11) :: psi, rho
 INTEGER :: i, j, m, ie, iw, jn, js
 REAL(KIND = DBL) :: R, rhon, lapRho

 INTEGER, DIMENSION(1:11,1:11,1:4) :: ni

 REAL(KIND = DBL) :: kappa, kappa_6, kappa_12, kappaEf, beta, beta4



 beta     = 12.D0*0.0001/(1.D0*( (1.0 - 0.1)**4 ))
 kappa    = 1.5D0*0.0001*1.D0/( (1.0 - 0.1)**2 ) 


!-------- Define near neighbours and initialize the density rho ----------------
 DO j = 1, 11
   DO i = 1, 11

! Initialize density
      rho(i,j) = 1.D0
        R =  DSQRT( ( DBLE(i)-5.0 )**2 + ( DBLE(j)-5.0 )**2 )
        IF (R <= (DBLE(3.0) + 1.D0)) THEN
          rho(i,j) = 0.55D0 - 0.5*0.9*TANH(2.D0*(DBLE(3.0) - R)/1.D0)
        END IF

 !Generate neighbors array
      ni(i,j,1) = i + 1
      ni(i,j,2) = j + 1
      ni(i,j,3) = i - 1
      ni(i,j,4) = j - 1
   END DO
 END DO


! Fix neighbours at edges
 ni(1,:,3) = 11
 ni(11,:,1) = 1
 ni(:,1,4) = 11
 ni(:,11,2) = 1



!--------- Differential terms for the stress form of the interfacial force -----
 DO j = 1, 11
   DO i = 1, 11

! Identify neighbors
     ie = ni(i,j,1)
     jn = ni(i,j,2)
     iw = ni(i,j,3)
     js = ni(i,j,4)

! Laplacian of the density rho
     lapRho = 4.D0*( rho(ie,j ) + rho(iw,j ) + rho(i ,jn) + rho(i ,js) )       &
             + rho(ie,jn) + rho(ie,js) + rho(iw,jn) + rho(iw,js) - 20.D0*rho(i,j)

! Define the chemical potential Psi
     psi(i,j) = 4.D0*beta*( rho(i,j) - 0.55 )*( rho(i,j) - 0.1 )*( rho(i,j) - 1.0 ) &
              - kappa*lapRho/6.D0
   END DO
 END DO


write(*,*) psi(6,6)-psi(4,4)


 END PROGRAM

解决方案

You are still not consequently using double precision throughout your code, e.g.:

beta     = 12.D0*0.0001/(1.D0*( (1.0 - 0.1)**4 ))

and many more. If I force the compiler to use double precision as default for floats (for gfortran the compile option is -fdefault-real-8), the result from your code is:

0.00000000000000000000000000000000000

So you need to fix your code. The cited line, for instance, should read:

beta     = 12.D0*0.0001D0/(1.D0*( (1.0D0 - 0.1D0)**4 ))

[Although I despise the notation D0, but that's a different story]

这篇关于如何在MatLAB中获得Fortran精度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！