0
J'ai un code écrit en Fortran et en Matlab. Ils font exactement le même calcul, à savoirComment obtenir la précision Fortran dans MatLAB
- construire un
tanh-field et trouver sa Laplacien - Multiplier certains termes ensemble
Le résultat de cette multiplication donne une matrice, dont (4,4) th et (6,6) th je soustrais.
- En Fortran leur différence est ~ 1e-20
- En Matlab leur différence est identique à zéro.
Ce problème est très critique, car je teste si ce nombre est inférieur à zéro. Question: Existe-t-il un moyen d'effectuer les calculs de manière à obtenir la même précision dans Matlab que dans Fortran?
I Liste des codes ci-dessous:
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
En dépit d'avoir été informé par ailleurs dans la réponse à une autre version de cette question (http : //stackoverflow.com/questions/29867870/precision-of-calculations) vous utilisez toujours des nombres à simple précision dans votre code Fortran. C'est '0.1' est, par défaut, un nombre réel à simple précision dans Fortran (sauf si vous avez utilisé les options du compilateur pour ajuster le comportement habituel et je parie que vous n'avez pas, avez-vous vous petit coquin?). Cela pourrait expliquer les différences que vous observez. À ce stade de votre quête, vous devriez au moins avoir étudié cette possibilité. –
Je l'ai changé donc il utilise une double précision, mais la conclusion est la même - les résultats diffèrent encore. J'ai passé beaucoup de temps à enquêter, croyez-moi ... – BillyJean
Je suis désolé, mais personne n'est susceptible de passer par beaucoup de code (en deux langues), l'exécuter, et le déboguer.Affinez votre problème. Découvrez le premier endroit où vos résultats divergent. En plus de ce que @HighPerformanceMark dit ci-dessus, voir aussi mon commentaire à votre question précédente - vous devez vérifier que les fonctions comme 'abs',' sqrt', 'tanh', etc. (et même la division) retournent le même résultat dans tous cas. Et par "vérifier" je ne veux pas simplement imprimer à quelques décimales. – horchler