Precision issues when going from Fortran to Matlab

Question

Niles Martinsen 2015-4-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/213791-precision-issues-when-going-from-fortran-to-matlab

评论： Thorsten 2015-4-27

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

Construct a tanh -field and find its Laplacian
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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Thorsten 2015-4-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/213791-precision-issues-when-going-from-fortran-to-matlab#answer_176859

编辑：Thorsten 2015-4-27

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 个评论
显示无隐藏无

Niles Martinsen 2015-4-27

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

Thorsten 2015-4-27

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.

请先登录，再进行评论。