Sqrt returns complex values

Question

I am trying to find sine from a cosine relation using the following commands:

CosEta[x_,y_] := 1/2 Cosh[Log[3]/2] (-Sqrt[y^2 + (x - Sech[Log[3]/2])^2] + Sqrt[y^2 + (x + Sech[Log[3]/2])^2]);    

sinsquare[x_,y_] := 1 - (CosEta[x,y])^2;    

SinEta[x_,y_] := Sqrt[sinsquare[x,y]];    

At some point, i.e. x-> 1.1, y-> -0.000000001 , CosEta[x,y] becomes -4.44089 * 10^-16, and as a consequence, SinEta returns complex number as 0. + 2.10734 * 10^-8 I . What can be done to get rid of this complex number return?

I have to implement this commands in my full code which includes multiple 1st and 2nd derivatives of SinEta and CosEta for which I can not use Re[] in defining the SinEta since Re' and Re'' shows up in the final outcome.

Answer 1

it is better to pass exact numbers to eliminate the chance of round off errors and at the very end apply N to the result. This way any accumulated round off are eliminated.

Added FullSimplify and now the complex value is gone

CosEta[x_, y_] := 
  1/2   Cosh[
    Log[3]/2]   (-Sqrt[y^2 + (x - Sech[Log[3]/2])^2] + 
     Sqrt[y^2 + (x + Sech[Log[3]/2])^2]);
sinsquare[x_, y_] := (1 - (CosEta[x, y])^2)
SinEta[x_, y_] := FullSimplify[Sqrt[sinsquare[x, y]]]
{CosEta[x, y], sinsquare[x, y], SinEta[x, y]} /. {x -> 1.1, 
   y -> -0.000000001} // Chop

To elaborate, the problem was this

CosEta[x_,y_]:=FullSimplify[1/2  Cosh[Log[3]/2]  (-Sqrt[y^2+(x-Sech[Log[3]/2])^2]+Sqrt[y^2+(x+Sech[Log[3]/2])^2])];
sinsquare[x_,y_]:=(1-(CosEta[x,y])^2)
SinEta[x_,y_]:=Sqrt[sinsquare[x,y]]

And now

SinEta[x, y]

Plugging {x -> 1.1, y -> -0.000000001} into the above gives complex

But if we first simplify the expression

 SinEta[x, y] // FullSimplify

And if we now plug the same x and y values it now gives small real value and no complex.

 (SinEta[x, y] // FullSimplify) /. {x -> 1.1, y -> -0.000000001}

Answer 2

The OP's code for CosEta suffers from subtractive/catastrophic cancellation when y is a small floating point number. This may be remedied by rationalization. Simplification will change the round-off error, change the edge cases, but is not designed to eliminate the problem.

Here is a version of CosEta, obtained by rationalizing the OP's code for it.

(* Num. stable *)
CosEta[x_, y_] := (2 x)/(
  Sqrt[y^2 + (x - Sech[Log[3]/2])^2] + 
   Sqrt[y^2 + (x + Sech[Log[3]/2])^2]);

sinsquare[x_, y_] := (1 - CosEta[x, y]^2)(*(1+CosEta[x,y])*);

SinEta[x_, y_] := Sqrt[sinsquare[x, y]];

OP's example:

SinEta[x, y] /. {x -> 1.1, y -> -0.000000001}

(* 0. *)

SinEta[1.1, -0.000000001]

(* 0. *)

More edge cases. The above code:

Table[
   SeedRandom[n, Method -> "MersenneTwister"];(* for reproducibility *)
   With[{x = RandomReal[{-6, 20}]}, 
    SinEta[x, x*$MachineEpsilon*RandomReal[2^28]]],
   {n, 200}] // Im // Counts

(* <|0 -> 200|> *)

With OP's code for CosEta:

Table[
   SeedRandom[n, Method -> "MersenneTwister"];(* for reproducibility *)
   With[{x = 10^RandomReal[{-6, 20}]}, 
    SinEta[x, x*$MachineEpsilon*RandomReal[2^24]]],
   {n, 200}] // Im // Counts // KeySort
(*
<|0 -> 140, 2.10734*10^-8 -> 2, 3.65002*10^-8 -> 4, 
 1.47514*10^-7 -> 2, 4.67906*10^-7 -> 5, 1.83557*10^-6 -> 2, 
 0.0000113751 -> 2, 0.0000258281 -> 1, 0.000181464 -> 3, 
 0.00031899 -> 8, 0.00236808 -> 10, 0.0207007 -> 10, 0.0518223 -> 4, 
 0.144338 -> 4, 0.57735 -> 3|>
*)

A numerically stable refactoring of SinEta:

SinEta[x_, y_] := With[{yy = y^2},
   With[{u = 4 yy + 3, w = 4 yy - 3, xx = x^2},
    With[{z = (u + Sqrt[u^2 + (8 w + 16  xx) xx])},
     Sqrt[(z - 4 xx)/(z + 4 xx)]
     ]]];

SinEta[x, y] - Sqrt[sinsquare[x, y]] // TrigToExp // 
 FullSimplify[#, {x, y} \[Element] Reals] &

(* 0 *)