[Math] Newton’s method and trig functions on a computer

numerical methodsrootstrigonometry

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.

Best Answer

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$?