How to deal with loss of significance in the case $f(x) = \sqrt{x+2} -\sqrt{x}$?

Question

I would like to evaluate the expression $f(x) = \sqrt{x+2} -\sqrt{x}$ with cases when $x = 3.2 \times 10^{30}$ and $x= 3.2 \times 10^{16}$. I tried using N[Sqrt[x+2] - Sqrt[x], 100] and

ScientificForm[Sqrt[x+2] - Sqrt[x], 100], both yielding 0 as an output. How can I obtain the desired output?

I also tried the method in Making a calculation with high precision by applying .`100 as a suffix but my mathematica 12 doesn't seem to be recognizing it.

score 19 · Answer 1 · answered Oct 17 '19 at 20:14

19

Your expression can be written in a more numerically stable form:

Sqrt[x + 2] - Sqrt[x] == 2/(Sqrt[x + 2] + Sqrt[x]) // FullSimplify
(* True *)

This evaluates without significant loss of precision.

2/(Sqrt[x + 2] + Sqrt[x]) /. x -> N[32*^15]
(* 5.59017*10^-9 *)

2/(Sqrt[x + 2] + Sqrt[x]) /. x -> N[32*^15, 30]
(* 5.59016994374947415367652880074*10^-9 *)

answered Oct 17 '19 at 20:14

mikado

16,741
2
20
54

I assume that it's numerically more stable because you're not subtracting two large numbers from each other? – march Oct 17 '19 at 20:35
@march That's what I meant. – mikado Oct 17 '19 at 20:50

score 12 · Accepted Answer · answered Oct 17 '19 at 20:14

By using arbitrary precision numbers instead of machine numbers:

x = 3.2`100 10^30;
Sqrt[x + 2] - Sqrt[x]

5.59016994374947424102293417182731712454783504367865447721289523170435 *10^-16

Or, even better, by using exact numbers and by doing the conversion afterwards (this forces Mathematica to compute all 100 leading digits):

x = 32/10 10^30;
N[Sqrt[x + 2] - Sqrt[x],100]

5.590169943749474241022934171827317124547835043678654477212895231704350107190085921247950310402211757*10^-16

march · Answer 3 · 2019-10-22T21:10:50.507

10

Use the asymptotic form of the function, i.e.

Normal@Series[Sqrt[x + 2] - Sqrt[x], {x, ∞, 2}]
% /. x -> 3.2*10^30
(* -(1/(2 x^(3/2))) + 1/Sqrt[x] *)
(* 5.59017*10^-16 *)

Note that we don't even need the second term to this level of precision:

Normal@Series[Sqrt[x + 2] - Sqrt[x], {x, ∞, 1}]
% /. x -> 3.2*10^30
(* 1/Sqrt[x] *)
(* 5.59017*10^-16 *)

edited Oct 22 '19 at 21:10

answered Oct 17 '19 at 20:32

march

23,399
2
44
100

1

First I was going to do @mikado's, then I was going to do this but played with it too long. +1 :) – Michael E2 Oct 17 '19 at 20:33
2

@MichaelE2 This is a surprising question in the sense that such a simple thing has invited four useful (and pretty different) answers. – march Oct 17 '19 at 20:37

score 7 · Answer 4 · answered Oct 17 '19 at 20:17

Replace your approximate input with exact rational numbers, reduce the result to a single exact algebraic number, and evaluate that numerically.

f[xx_] := With[{x = Rationalize[xx, 0]}, N[RootReduce[Sqrt[x + 2] - Sqrt[x]]]]
f[3.2 10^30]
(* 5.59017*10^-16 *)

How to deal with loss of significance in the case $f(x) = \sqrt{x+2} -\sqrt{x}$?

4 Answers4