Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this community
I would like to make the following calculation:
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.)] - 1
Mathematica 8 returns 0. The result is obviously not 0, but my calculation must be losing precision. How can I compute this value with sufficient precision?
For the sake of alternatives, here's another way. You can specify a precision with a backtick. This is called an "arbitrary-precision" number. Such numbers are treated differently than machine precision numbers.
Examples of 3 and 5 digits of precision:
16.`3
16.0
16.`5
16.000
Machine precision depends on the machine, but it's often nearly 16 digits.
$MachinePrecision
15.9546
The results of your calculation at different precisions are given below. Note they are more accurate than machine precision, even though the specified precisions are technically less. That's because Mathematica can use up to $MaxExtraPrecision extra digits of precision in its internal calculations.
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.`3)] - 1
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.`4)] - 1
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.`5)] - 1
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.`15)] - 1
1.3*10^-17 1.25*10^-17 1.250*10^-17 1.2500000000000*10^-17
Which way is better, this one or using N and exact numbers as in the answers of @chuy and @m_goldberg? Well it depends. Normally N[exact] is sufficient, although in this case, it fails -- it gives 0. Arbitrary-precision calculations take longer, but Mathematica does try to keep track of the precision as it goes along. On the other hand, N[exact, prec] (@m_goldberg's answer) tracks the precision and, more importantly, it gives an answer with the target precision, if the internal precision limit $MaxExtraPrecision is not exceeded.
Update: Alternative approach
The difficulty in this calculation is a classic problem in numerical analysis, so please forgive me if this is well-known. The reason can be seen in the numerator of $${1 \over \sqrt{1-a}} -1 = {1 - \sqrt{1 - a} \over \sqrt{1 -a}} \,.$$ If $a$ is very small, then we're subtracting two numbers close to one another, resulting in a loss pf precision. In fact, $\sqrt{1 - a}$ may be equal to $1$ when rounded to machine precision. A solution is to rationalize the numerator to get $${a \over \big(1 + \sqrt{1 - a} \big) \sqrt{1 -a}} = {a \over 1-a+\sqrt{1-a}}\,.$$ In this form, it evaluates with machine precision without a problem.
a/(1 + Sqrt[1 - a] - a) /. a -> (150^2 10^(-4))/(9 10^16.)
1.25*10^-17
Even N has less trouble with this form. Here we limit the extra precision to force the difficulty earlier:
Block[{$MaxExtraPrecision = 4},
N[a/(1 + Sqrt[1 - a] - a) /. a -> (150^2 10^(-4))/(9 10^16), 20]]
1.2500000000000000234*10^-17
Block[{$MaxExtraPrecision = 4},
N[1/Sqrt[1 - a] - 1 /. a -> (150^2 10^(-4))/(9 10^16), 20]]
N::meprec: Internal precision limit $MaxExtraPrecision = 4.` reached while evaluating -1+200000000/Sqrt[39999999999999999]. >> 1.250000*10^-17
This same problem occurs in simple arithmetic:
1. - 1. + 10.^-16
1. + 10.^-16 - 1.
1.*10^-16 0.
In sum, another solution would be to choose a formula that avoids such numerical problems. (I don't know of a way to get Mathematica to figure this out automatically.) While it is rare that N cannot handle a computation, machine precision calculations are fast and to be desired if possible.
$MachineEpsilon is the smallest h such that 1. + h - 1. is nonzero. Smaller than that (as in Michael's example), and you get subtractive cancellation.
$\endgroup$
May 24, 2013 at 13:15
Ditch the decimal point after the 16 (this makes it a machine precision number):
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16)] - 1
yeilds:
-1 + 200000000/Sqrt[39999999999999999]
. I don't get an actual number. I just need better precision.
$\endgroup$
Sqrt[2] is not a number, too? In any event, look up N[], and note well what its second argument does.
$\endgroup$
May 23, 2013 at 14:42
1e-95 or something like that.
$\endgroup$
In V11.0 we try:
N[1/Sqrt[1 - (150^2 10^(-4))/(9 10^16)] - 1]
we get an answer of 0.
As Michael E2 points out we are getting subtractive cancellation causing the loss of significant digits.
So we say
$$a=\frac{22500}{10^4 \left(9 \cdot 10^{16}\right)}$$
In numerical work the Taylor Series is your friend, it is used often. We will try to get a better answer just with a series. We expand a around 0.
Series[1/Sqrt[1 - a] - 1, {a, 0, 3}]//Normal
yields $\frac{5 a^3}{16}+\frac{3 a^2}{8}+\frac{a}{2}$
This is numerically very stable as it has small numerators and no subtraction of nearly equal quantities. We Hornerize it by using
HornerForm[a/2 + (3 a^2)/8 + (5 a^3)/16]
yields $a (1/2 + (3/8 + (5 a)/16) a) $
N[a (1/2 + (3/8 + (5 a)/16) a) /. a -> (22500 10^(-4))/(9 10^16)]
yields $$1.25 \cdot 10^{-17}$$
A better answer.
You can get a high-precision approximation by using N with two arguments, the second argument specifying the precision you desire.
a = 1/Sqrt[1 - (150^2 10^(-4))/(9 10^16)] - 1;
Table[ N[a, p], {p, 17, 22}] // TableForm
1.2500000000000000*10^-17
1.25000000000000002*10^-17
1.250000000000000023*10^-17
1.2500000000000000234*10^-17
1.25000000000000002344*10^-17
1.250000000000000023438*10^-17
.after 16 $\endgroup$