Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
2  
Just remove the . after 16 –  rm -rf May 23 '13 at 14:36

3 Answers 3

up vote 4 down vote accepted

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.

share|improve this answer
    
Thank you very much for your answer! There is really usefull information in your answer! Thank's! –  Thanos May 24 '13 at 6:07
1  
@Thanos Thanks. I added another solution that in some ways is better but more work for the user, and I hope I made it clear that I prefer m_goldberg's way of dealing with the given calculation –  Michael E2 May 24 '13 at 13:14
    
In particular, $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. –  J. M. May 24 '13 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]

share|improve this answer
    
Thank you very much for your answer! If I remove the . I don't get an actual number. I just need better precision. –  Thanos May 23 '13 at 14:40
1  
What do you mean you don't get an actual number, @Thanos? This is a number, symbolically presented. I suppose you will say that Sqrt[2] is not a number, too? In any event, look up N[], and note well what its second argument does. –  J. M. May 23 '13 at 14:42
    
@J.M.: he he he! It is a true number indeed! all I need is something in the format 1e-95 or something like that. –  Thanos May 23 '13 at 14:43
    
Hence my third sentence, @Thanos. Did you read the docs already? –  J. M. May 23 '13 at 14:45
1  
@Thanos, then it seems you don't want the maximum precision (which your question explicitly asks for). –  chuy May 23 '13 at 15:05

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

share|improve this answer
    
Thank you very much for your answer! –  Thanos May 24 '13 at 6:10

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.