Cómo obtener la precisión Fortran en MatLAB

double-precision

Question

Tengo un código escrito en Fortran y en Matlab. Ellos hacen exactamente el mismo cálculo, es decir

1. Construya un tanh -field y encuentre su Laplacian
2. Multiplicar algunos términos juntos

El resultado de esta multiplicación produce una matriz, cuya (4,4) th y (6,6) th restamos.

• En Fortran su diferencia es ~ 1e-20
• En Matlab, su diferencia es idénticamente cero.

Este problema es muy crítico, ya que pruebo si este número es menor que cero. Pregunta : ¿Hay alguna manera de realizar los cálculos de manera que obtenga la misma precisión en Matlab que en Fortran?

Enumero los códigos a continuación:

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

Answer 1

Por lo tanto, todavía no está utilizando una precisión doble en todo su código, por ejemplo:

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

y muchos más. Si -fdefault-real-8 al compilador a usar doble precisión como predeterminado para flotantes (para gfortran la opción de compilación es -fdefault-real-8 ), el resultado de su código es:

0.00000000000000000000000000000000000

Entonces necesitas arreglar tu código. La línea citada, por ejemplo, debe leer:

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

[Aunque desprecio la notación D0 , pero esa es una historia diferente]