Significantly different results between Mathematica and Wolfram Alpha with large radicals and small numbers [closed]

Question

I have some code which involves tiny numbers being put to the power of very large numbers. The function I'm looking at is

$\varphi = \omega(T) \left(1 - (1 - \epsilon)^{n_{e}(T)} \right)$

when $\epsilon $ is very small (~$10^{-16}$) and $n_{e}$ is large ($> 10^{10}$). Both $n_{e}$ and $\omega$ are known functions of $T$. At a known value ($T = 37$) the value of $\varphi$ is known, and we can in theory estimate $\epsilon$ by

$\epsilon = 1 - \left(1 - \frac{\varphi}{\omega} \right)^{\frac{1}{n_{e}}}$

However, when I put this into Mathematica I'm getting very strange behaviour; here's a MWE when $\varphi(37) = 0.2506$;

gammaex = 0.2506; 
omega[t_] := 2.43163218375*10^7*Exp[1700*(1/298.15 - 1/(273.15 + t))];

w[t_] := (3.414105049212413*10^12)/(omega[t]); 
v[t_] := Sqrt[661.6469313477045*(t + 273.15)];
ne[t_] := (v[t]*5.104757516005496*(10^7));

epsilon = SetAccuracy[ 1 - ( 1 - gammaex/w[37])^(1/ne[37]), 30];
test = w[37]*(1 - (1 - epsilon)^(ne[37]))

Now when I evaluate $\varphi(37)$ via the last line I should get the value 0.2506 but I do not; instead I get $\varphi = 0.289104$, over 15% off the true value. I thought maybe this was a precision problem, so I tried some other commands (Surd, Power...) and got the same wrong value. The output value for the epsilon is $\epsilon \approx 1.11 \times 10^{-16} $ from mathematica.

However I've evaluated this with WolframAlpha and got a markedly different value of $\epsilon \approx 9.61 \times 10^{-17}$, and as this link demonstrates the calculation seems to work with WolframAlpha, returning close to the expected value around 0.25*. Any ideas why the calculations are different, and how I can make Mathematica behave with such extreme values?

Incidentally, if I dump the WolframAlpha equation ( (112608)*(1 - (1 - 9.61369*10^-17)^(2.31246*10^10)) ) directly into Mathematica, it still comes up with the wrong value, so I assume it's a radical issue?

Answer 1

As J. M. already suggested, this is a machine precision issue. So let us do your computation with arbitrary precision numbers. For doing so, all machine numbers have to be replaced with arbitrary precision numbers, otherwise the computation falls back to machine numbers. In the following command I have done this by placing `30 after each machine number.

gammaex=0.2506`30;
omega[t_]:=2.43163218375`30*10^7*Exp[1700*(1/298.15`30-1/(273.15`30+t))];

w[t_]:=(3.414105049212413`30*10^12)/(omega[t]);
v[t_]:=Sqrt[661.6469313477045`30*(t+273.15`30)];
ne[t_]:=(v[t]*5.104757516005496`30*(10^7));

epsilon=1-(1-gammaex/w[37])^(1/ne[37]);
test=w[37]*(1-(1-epsilon)^(ne[37]))
FullForm[test]

Now the result is what you expected:

(* 0.25060000000000000000000000000 *)
(* 0.2506`28.52376207789293 *)

Answer 2

You can use Rationalize to convert numbers to exact numbers

gammaex = 0.2506 // Rationalize[#, 0] &;
omega[t_] = 
  2.43163218375*10^7*Exp[1700*(1/298.15 - 1/(273.15 + t))] // 
   Rationalize[#, 0] &;

w[t_] = (3.414105049212413*10^12)/(omega[t]) // Rationalize[#, 0] &;
v[t_] = Sqrt[661.6469313477045*(t + 273.15)] // Rationalize[#, 0] &;
ne[t_] = (v[t]*5.104757516005496*(10^7)) // Rationalize[#, 0] &;

epsilon = 1 - (1 - gammaex/w[37])^(1/ne[37]);
test = w[37]*(1 - (1 - epsilon)^(ne[37]))

1253/5000

% // N

0.2506