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]