How to improve accuracy?

I am trying to evaluate a matlab code. A snippet is shown. I am facing an issue with accuracy of digits. I expect the answer to be 0 but get a non-zero value. My entire code deals with many more similar calculations. Is there a way to improve the accuracy of the code. I am using double precision. Any help would be appreciated

Code:

Tw = 306.13;
Tsat = 302.63;
rho_l = 1440.8;
hlv = 126e3;
Adis = 2e-21;
Pc = (Tw/Tsat - 1)*hlv*rho_l;
del1 = (Adis/Pc)^(1/3);
Pc - Adis/del1^3 %Expected answer = 0

[Aquatris]

Use vpa function if you have the Symbolic Math Toolbox with your calculations. The vpa function also allows you to define the precision by simply adding a second argument to the function;

i.e., "vpa(1/3+5*6-2,32)"

clear 
Tw = 306.13;
Tsat = 302.63;
rho_l = 1440.8;
hlv = 126e3;
Adis = 2e-21;
Pc = vpa((Tw/Tsat - 1)*hlv*rho_l);
del1 = vpa((Adis/Pc)^(1/3));
vpa(Pc - Adis/del1^3) %Expected answer = 0, results -1.5e-32

[Siddharth Iyer]

Thank you Aquatris!

Answer 1

[Walter Roberson]

If you need an exact 0 in those circumstances, then you need to use all of the constants as symbolic rational numbers:

Q = @(x) sym(x, 'r');
Tw = Q(306.13);
Tsat = Q(302.63);
rho_l = Q(1440.8);
hlv = Q(126e3);
Adis = Q(2e-21);
Pc = (Tw/Tsat - 1)*hlv*rho_l;
del1 = (Adis/Pc)^(Q(1)/Q(3));
Pc - Adis/del1^Q(3)

Aquatris suggests using vpa(), but if you vpa() the result of the cube root, then it is rather unlikely that cubing that will give you exactly the original back.

It's like the fact that if you use 2 digits to represent 1/3 in decimal, getting 0.33, and then you multiply that by 3, you will get 0.99 rather than 1 exactly. You could try to fix that by using (say) 20 digits, which would be 0.33333333333333333333, but multiply that by 3 and you get 0.99999999999999999999 and not 1. No matter how many fixed digits you go to, 50, 1000, 10 million, you will still have the problem of not getting back exactly 1. (The argument fails for infinite representations. It turns out that 0.999.... infinite times 9999 .... is the same as exactly 1 according to the theory of infinitesimals.)

[Siddharth Iyer]

Got it. Thank you Walter!