Precision of calculations

Question

I am doing a calculation in Fortran on a double-precision variable, and after the calculation the variable gets the value -7.217301636365630e-24.

However, when I do the same computation in Matlab, the variable just gets the value 0. Is there a way to increase the precision of MatLAB when doing calculations such that I would also be able to get something on the order of 7e-24?

Ideally, it would be something I could apply to all calculations in the script and not just a single variable. Something similar to when using format long.

For me this kind of precision is crucial as I need to determine if a variable is indeed negative or not.

I have added the code. It is rather long, but I couldn't trim it further without throwing away variables and their precision. The last term, Ax(i,:,:), is the one that I would like to have a very high precision on. So the important stuff occurs only in the last line.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONSTANTS
clc
clear all

sym_weight     = [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];



ly = 11; lx = ly;
xC = 5; yC=xC;    

density_high = 1.0;
density_low = 0.1;
radius  = 2;
interface_w = 1;
sigma_st = 0.0001;


beta  = 12*sigma_st/(interface_w*(density_high-density_low)^4);
kappa = 1.5*sigma_st*interface_w/(density_high-density_low)^2;





saturated_density = 0.5*(density_high+density_low);
for x=1:lx
    for y=1:ly
        for i=1:9
            fIn(i, x, y) = sym_weight(i)*density_high;
            gIn(i, x, y) = 3*sym_weight(i);

            test_radius = sqrt((x-xC)^2 + (y-yC)^2);
            if(test_radius <= (radius+interface_w))
                fIn(i, x, y) = sym_weight(i)*( saturated_density - 0.5*(density_high-density_low)*tanh(2*(radius-sqrt((x-xC)^2 + (y-yC)^2))/interface_w) );
            end
        end
    end
end




density_2d = ones(lx)*saturated_density;
for i=1:lx
    density_aux(:,:,i) = abs(density_2d(:, i)');
end



density_local             = sum(fIn);
L_density_local           = (+1.0*(circshift(density_local(1,:,:), [0, +1, +1]) + circshift(density_local(1,:,:), [0, -1, +1]) + circshift(density_local(1,:,:), [0, +1, -1]) + circshift(density_local(1,:,:), [0, -1, -1])) + ...
                             +4.0*(circshift(density_local(1,:,:), [0, +1, +0]) + circshift(density_local(1,:,:), [0, -1, +0]) + circshift(density_local(1,:,:), [0, +0, +1]) + circshift(density_local(1,:,:), [0, +0, -1])) + ...
                             -20.0*density_local(1,:,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




chem_pot   = 4*beta*(density_local-density_low).*(density_local-density_high).*(density_local-density_aux) - kappa*L_density_local/6;






for i=3
    Ax(i,:,:)            =     (+circshift(chem_pot(1,:,:), [0,-2*dir_x(i),-2*dir_y(i)]) -             chem_pot(1,:,:));                                          
end

You'll have to show us your code. Numbers of that range should be no problem at all. — beaker, Apr 25 '15 at 16:45
@beaker It is a very long script. I can try and make a minimal example if it is required. I thought it could be fixed with a simple command such as "format long" — BillyJean, Apr 25 '15 at 16:47
`format long` is only going to affect the display of the number and not its value. I assume you've already done this and it's still displaying `0`. — beaker, Apr 25 '15 at 16:49
@beaker That is true. I will try and make a small example, because the actual code is really long. I also perform the calculations on many variables and all of them need this high precision, so a single command would be nice (if it exists.) — BillyJean, Apr 25 '15 at 16:51
If you're not already doing it, make your variables doubles. Somewhere in the calculations there's got to be somewhere that an intermediate value is getting set to 0. — beaker, Apr 25 '15 at 16:52
Matlab already does calculation by default on double precision values ... the only way you can increase the precision is by using the [variable precision arithmetic](http://mathworks.com/help/symbolic/vpa.html) — Hoki, Apr 25 '15 at 16:53
@beaker I thought they were double by default? If not, then how do I set a variable to double? — BillyJean, Apr 25 '15 at 16:54
`your_variable=double(your_variable)` will make sure they are treated as double (but this is indeed the default type unless you specify something else) — Hoki, Apr 25 '15 at 16:56
@Hoki I added a small example of the code, showing the term where it goes wrong. Thanks — BillyJean, Apr 25 '15 at 17:22
So everything before the last four or five lines produces identical floating point values in FORTRAN and Matlab? Be aware that functions like `abs`, `sqrt`, `tanh`, etc. (and even division) will possibly not produce identical floating point output for the same inputs in different languages. Summing terms in a different order can also result in different output. I assume also that you've confirmed that this issue is not due to a bug in either version. — horchler, Apr 25 '15 at 18:20
@horchler yes, to the 19th decimal. It is when I perform the command to find `Ax` that some of the elements of `Ax` are rounded to 0 in MatLAB, but have values on the order of `1e-24` in Fortran. Also, it is simply impossible for an entry of `Ax` to be 0, when thinking about how it is constructed mathematically. Small yes, but not 0. — BillyJean, Apr 25 '15 at 18:22
@Hoki Thanks for taking a look at it. I am using `i=3`, so it should shift — BillyJean, Apr 25 '15 at 18:41
@Hoki When I output an entry from `density`, it only gives me 9 digits afterwards. If I use, e.g., `vpa(density(1,1,1), 25)` it then adds 25 digits. How can I make 25 digits the default for all variables in the script? — BillyJean, Apr 25 '15 at 18:48
sorry for the confusion about `i`. Still, if you look at the first term of your two members (the one with circshift and the straight one), there are a lot of equal terms (I count 49). So if these terms are meant to be different, the precision is not lost during the last subtraction, but before. I am not familiar with `vpa`, but i would guess you'll have to declare all the variable that matters individually. — Hoki, Apr 25 '15 at 18:57

score 3 · Answer 1 · answered Apr 26 '15 at 13:02

You have not shown the fortran code, but be aware that in Fortran, when you do this:

 density_low = 0.1

The literal 0.1 is single precision, regardless of the type of density_low. All of those literals need to be expressed as 0.1D0 or 0.1_k where k is the appropriate kind integer.

(Sorry if you knew that, but its a common mistake )

Precision of calculations

1 Answers1

Linked