Precision issues when going from Fortran to Matlab

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)

Answer 1

[Thorsten]

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]

Thanks. Which toolbox do you recommend for me and the above example?

[Thorsten]

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.