Как получить точность Fortran в MatLAB

Question

У меня есть часть кода, написанная в Фортране и в Matlab. Они делают точно такой же расчет, а именно

1. Построить tanh -поле и найти ее лапласианом
2. Multiply некоторые термины вместе

В результате этого умножения дает матрицу, которой (4,4) th и (6,6) th я вычитаю.

• В Fortran их разность составляет ~ 1e-20
• В Matlab их разность тождественно равна нулю.

Эта проблема очень важна, поскольку я тестирую, если это число меньше нуля. Вопрос: Есть ли способ выполнить вычисления, чтобы получить такую ​​же точность в Matlab, как в Fortran?

I список кодов ниже:

---

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

Вы все еще не следовательно, при использовании двойной точности на протяжении всего кода, например .:

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

и многое другое. Если я заставить компилятор использовать двойную точность по умолчанию для поплавков (для gfortran опция компиляции является -fdefault-real-8), результат от вашего кода:

0,00000000000000000000000000000000000

Так что вам нужно исправить код , Приведенные строки, например, следует читать:

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

[Хотя я презираю обозначение D0, но это другая история]