# Precision issues when going from Fortran to Matlab

2 views (last 30 days)
Niles Martinsen on 27 Apr 2015
Commented: Thorsten on 27 Apr 2015
Hi
I have a piece of code written in Fortran that I have now finished writing in Matlab. They do exactly the same calculation, namely
1. Construct a tanh -field and find its Laplacian
2. Multiphy some terms together
The result of this multiplication yields a matrix, whose (4,4)th and (6,6)th element are identical according to Matlab. However, in Fortran their difference is ~1e-20. This is very critical for me, as I test if this number is less than zero or not.
Is there a way to perform the calculations such that I get the same precision as in Fortran?
A minimal version of my code is the following, where the last line, psi(1,4,4)-psi(1,6,6), is the one that incorrectly gives 0:
clear all
format long
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;
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;
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)

Thorsten on 27 Apr 2015
Edited: Thorsten on 27 Apr 2015
eps = 2^(-52) = 2.2204e-16 is Matlab's precision for double. e-20 is below this precision. In fact, it is below the precision of double-precision binary floating-point according to IEEE 754. So you need higher precision, which is not supported by native Matlab. You have to use multiprecision extensions (aka toolboxes) for Matlab, like http://www.mathworks.com/matlabcentral/fileexchange/36534-hpf-a-big-decimal-class or http://www.mathworks.com/matlabcentral/fileexchange/6446-multiple-precision-toolbox-for-matlab
Niles Martinsen on 27 Apr 2015
Thanks. Which toolbox do you recommend for me and the above example?
Thorsten on 27 Apr 2015
I would recommend the first, but you have to try yourself which fits your needs best. There are also commercial solutions around. Please do not forget to accept the answer.