Why do I get a big relative error for my function? (Numerical Analisys - floating point)

Question

When evaluating on the computer the following function:

$$f(x)=\frac{x^2}{(\cos(\sin(x)))^2-1}$$

there is a big relative error for values $x\approx0$ (values very close to zero). I used the Taylor Series for the above function and I got: $$ -1 - \frac{2}{3}x^2 - \frac{2}{15}x^4 + \frac{1}{945}x^6 - \dotsb $$ and for $x=0.1\mbox{e-}7$ I got the result $-1$ when the result should be different but the $x^2$, $x^4$, $x^6$ and $x^8$ terms are ignored. I suppose it is because of the floating point representation. Why is this happening? What's the explanation? And how can I find a method to compute $f$ for values $|x|<1$ at machine precision?

Answer 1

This happens because you are trying to divide two very small numbers, and the denominator has a lot more error in it. Let $x \to 0$, so then $\sin x \approx x$, and $\cos x \approx 1-x^2/2$, so $$ f(x) \approx \frac{x^2}{(1-\frac{1}{2}x^2)^2-1} = \frac{x^2}{-x^2+\frac{1}{4}x^4}$$ When $x=10^{-8}$, then $x^4$ and $x^2$ differ by $10^{16}$, which is below machine precision for double floating point numbers, so $-x^2+x^4 = -x^2$, resulting in $-1$.

You can evaluate this more robustly by using the Taylor series below a certain value of $|x|$:

double x2 = x*x;
if(fabs(x) < tol){
    return -1 - (2./3)*x2;
}else{
    double sx = sin(x);
    double csx = cos(sx);
    return x2 / (csx*csx - 1);
}

You can determine the value of tol using the error bound for an alternating series, but it will be more stringent than necessary. Empirically, a value around $10^{-4}$ or $10^{-5}$ seems appropriate, by asking wolframalpha to plot the difference between $f(x)$ and it's truncated Taylor series.

In practice, such kinds of functions are numerically evaluated using Chebyshev polynomial approximations or Pade rational polynomial expansions which can achieve very low uniform error with on the order of 10th order polynomials.

Answer 2

The problem is that the double $-1$ is the one closest to the actual value of $f\left(1\cdot10^{-8}\right)$

We have $$ f\left(1\cdot10^{-8}\right)+1\approx-6.6666666666666668\cdot10^{-17} $$ Machine precision with the datatype double can represent the number above quite fine, however when we subtract that one, we run into the problem of machine precision with the type double has approximately 15 significant digits, and the number is rounded to be exactly $-1$.

Answer 3

One could remark that $\cos^2A-1=-\sin^2A$ and $$ -\frac{x^2}{\sin(\sin x)^2} $$ should evaluate correctly for smaller values of $x$.

---

The floating point evaluation of $\sin x$ remains exact for $|x|<1e{-16}$, essentially by using the formula $\sin x=x·(1-\frac{x^2}{6}+h.o.t)$, which reduces to $\sin x=x$ for $|x|<\sqrt{6}·10^{-8}$. Thus for those small values the modified formula will give the constant $-1$, and because there is no branching involved, the graph is smooth in the transition to larger values of $x$.