Precision issues when going from Fortran to Matlab
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Thorsten
2015-4-27
编辑: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 个评论
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Fortran with MATLAB 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!