Clear Filters
Clear Filters

Precision issues when going from Fortran to Matlab

2 views (last 30 days)
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;
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)

Answers (1)

Thorsten
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
  2 Comments
Niles Martinsen
Niles Martinsen on 27 Apr 2015
Thanks. Which toolbox do you recommend for me and the above example?
Thorsten
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.

Sign in to comment.

Categories

Find more on Fortran with MATLAB in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!