2
$\begingroup$

I'm trying to use Newton's method to find roots for the function $A \cos(\Theta_2 - \Theta_1) + B \sin(\Theta_1)$. (That is, iterate $x_{i+1} = x_i - f(x_i) / f'(x_i)$). I've got a working implementation and everything's great, except that I'm not getting all the precision I'd expect.

I'm using 64 bit floating point numbers (doubles) throughout, and I'm expecting effective precision to be at least $1e^{-10}$, but I'm getting only maybe $1e^{-9}$. I traced the problem to the fact that $\cos(\frac{\Pi}{2}) \approx 1e^{-9}$ on my machine. Which means when I evaluate my function it returns 0 at a place that isn't exactly a root of my function.

This isn't entirely unexpected. I'm well aware of numeric issues on machines. I think most modern CPUs use CORDIC for trig evaluations, which interpolate between table lookup values. Some numeric noise on the least significant bits is sort of par for the course.

I'm just wondering if there's some known method to get around the limitations of the machine implementation of sine and cosine for root finding with Newton's Method? Is there some way to maybe transform my initial equation in some way that will give me a function with the same roots in some given neighborhood but better numerical behavior? Taylor series approximations come to mind but they don't always have great numerical behavior (although I don't have a lot of experience in that regard so I might be wrong).

Really any general advice about numeric root finding with trig functions would be helpful.

UPDATE

I was slightly wrong. The real problem is that: $\cos(3.141592650897981) = -1$ on my machine. That angle value is roughly $2.69e^{-9}$ off from actual $\Pi$. You can reproduce this with google, actually. Try googling "cos(3.141592650897981)" and "PI - 3.141592650897981" to see what I mean.

  • 0
    @Mich$a$el Hardy: I see your point. I've updated my post.2012-01-03

3 Answers 3

2

With your update at least the problem is clear : $\cos(\pi-\epsilon)=-\cos(\epsilon) \approx -1+\epsilon^2/2$.
$\epsilon^2/2 \approx 3.6\ 10^{-18}$ (that is the 18 digits you may hope using double precision) but $\epsilon \approx 2.69\ 10^{-9}$.

So that your program is probably correct.
Perhaps that some rewriting like x_{i+1} = (x_i f'(x_i) - f(x_i)) / f'(x_i) could help... (quad. precision surely with the adapted mathematical libraries!).
EDIT : This rewriting isn't useful in this case so let's forget that and resume the situation.

The problem is the evaluation of $\cos$ or $\sin$ function at points were the derivative is $0$ so that no first order term will appear after evaluation but only terms of order two or more. You can't 'cure' this problem inside the cos or sin function because of the constant term : $-1 + \frac{\epsilon^2}{2}$ is considered $-1$ for $\epsilon$ of order $10^{-8}$ or less (IEEE 64 bits double).

If you need more precision you'll have to combine two or more terms : in your specific example you may rewrite $1+\cos(u)$ as $\ 2\cos(u/2)^2$ (or directly as $(u-\pi)^2/2!-(u-\pi)^4/4!$) when $|u-\pi|\lt 10^{-6}$ (for example). Note that there should be no loss of significance for $u$ near $0$.

In the more general case $A \cos(\Theta2-\Theta1) + B \sin(\Theta1)$ you could perhaps proceed at some rewriting when : ( this sum is very small ) and ( if $\Theta2-\Theta1$ is 'near' $k \pi$ or if $\Theta1$ is 'near' $(k+1/2) \pi$ ). Another technical solution could be to use higher precision (slower!) arithmetic in the few cases where the sum is small!

The real problem is to use the Newton(–Raphson) method outside of its applicability limits. Quadratic convergence is not warranted when f'(s)=0 for the solution $s$. In your case $f(s)$ with $s=\pi+\epsilon$ is of order $\epsilon^2$ and f'(s) is of order $\epsilon$ but $\epsilon^2 \lt 10^{-16}$ was replaced by $0$ so that instead of adding a term of order $\epsilon=\frac{\epsilon^2}{\epsilon}$ at each iteration we add merely $0$ and the Newton algorithm will not converge. Another solution could be to try an alternative algorithm when $f(s)$ and f'(s) become too small...

  • 0
    Thanks for your explanation. I at least understand the problem pretty well now. I'll play around with various hacks and see if I can get reasonable behavior.2012-01-04
3

Somewhere your program seems to be using single precision rather than double precision. Any reasonable implementation of double-precision arithmetic should give you $\cos(\pi/2)$ with an error of $10^{-16}$ or so. What do you get for $\pi$?

  • 0
    Ah, you were right to be suspicious. The actual problem is that I'm finding the value of sine/cosine from $\Theta = \frac{-\Pi}{2} + \frac{\Pi}{2} * 2.9999999983$. The calculated cosine of that is apparently -1 exactly, but the sine of that is 0.0000000026918122290084216. I'll update my original post2012-01-03
3

Okay, after quite a bit of research and work I think I've found a robust fix for the problem (fixes all the test cases that I had failing before). It's non trivial enough that I wrote up the explanation and hopefully it'll be useful to someone someday.

Basically, cosine and sine lose digits of precision as you get near to $1$ (because $1 - 1e^{-20} = 1$, even in double precision, because the mantissa can't contain all the digits). So you need to separate out the $1$ from the digits of precision and treat them separately. That is, get the trig call in to a form like $1 - \epsilon$, and then calculate $\epsilon$ directly. So to fix, it, I transform the trig calls thusly:

  • $\cos(\Theta) = 1 - 2 \sin^2\frac{\Theta}{2}$ when $\cos(\Theta)$ is near $1$
  • $\cos(\Theta) = -1 + 2 \cos^2\frac{\Theta}{2}$ when $\cos(\Theta)$ is near $-1$

For sine, I combined the two identities above with the identity $\sin(\Theta) = \cos(\frac{\Pi}{2} - \Theta)$. So basically just plugged in $\frac{\Pi}{2} - \Theta$ for $\Theta$ in the above equations. These transformations are exact for any values of $\Theta$, but they lose digits when they're away from $1$ or $-1$ (same as the normal trig calls do when they're near $1$ or $-1$).

The next problem is that, since the aggregate function evaluates to approximately $0$ in the neighborhood of its root, there's a lot of cancellation issues at work when you actually evaluate the functions. So instead of just naively summing up all the terms, I use something like Knuth's summation scheme (algorithm 4.1 (using algorithm 3.1) in this paper: Accurate Sum and Dot Product) (basically keeping a running sum and a running underflow value, and taking the final sum as the sum of the running sum and the underflow at the end).

The final function evaluation seems accurate to something like at least $1e^{-20}$ or so near the root (that is, basically the full double precision), which is sufficient for Newton's method to actually converge to an answer with a good number of significant digits (I stop the iteration at a step size of $1e^{-12}$, but it could probably converge even further).

Anyway, hope that's useful to someone in the future.

  • 0
    I actually came across versine during my research. Never heard of it before, but various versine identities helped me figure some of the above out.2012-01-07