0 个投票

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

2 个评论
显示 无 隐藏 无

Aquatris
Aquatris 2018-6-8
编辑：Aquatris 2018-6-8
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
Siddharth Iyer 2018-6-12
Thank you Aquatris!

请先登录，再进行评论。

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

 采纳的回答

Walter Roberson
Walter Roberson 2018-6-8
编辑：Walter Roberson 2018-6-12

0 个投票

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.)

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Number Theory 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by