How to use Newton's method to find the roots of an oscillating polynomial?

Question

Use Newton’s method to find the roots of $32x^6 − 48x^4 + 18x^2 − 1 = 0$ accurate to within $10^{-5}$.

Newton's method requires the derivative of this function, which is easy to find. Problem is, there are several roots near each other:

These roots are so close that they seem to prevent Newton's method from working correctly, but I don't know why. How to modify the method so that it still works in this case? And how to code this modified method (in Matlab) so that the algorithm find these roots efficiently?

Answer 1

Here you go:

alpha = 0.05;
NM = @(x,f)x-alpha*prod(([1,0]*f([x,1;0,x])).^[1,-1]);
f = @(x)32*x^6-48*x^4+18*x^2-x^0;
x = 0.5;
eps = 1;
while eps > 1e-5
    xn = NM(x,f);
    eps = abs(xn-x);
    x = xn;
end
[x f(x)]

I caution you to not turn in this code to your professor. He will probably ask you how it works.

I provide this so that you may see how the method works. However, the underlying mechanics of how this particular code works are somewhat complicated and likely beyond your level. You will have to write your own function. If you turn this in, and the professor reviews your code, I guarantee you will get penalized for plagiarism.

---

In an earnest attempt at an answer, the above code includes one trick.

An iteration of Newton's method looks like this:

$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.$$

However, if we plot our polynomial, we find that it oscillates somewhat like a sine wave centered at $x=0$. This causes a unique behavior of Newton's method: each projection of the slope throws the value to another part of the wave, and you get a "back and forth" action.

In order to accommodate this, we add a scaling parameter, $\alpha$:

$$x_{n+1} = x_n - \alpha\frac{f(x_n)}{f'(x_n)}.$$

Setting this scaling parameter such that $0 < \alpha < 1$ attenuates the effect of the slope. This slows down the convergence, but it ensures that we don't project "too far" ahead and get caught in an infinite loop.

Examine the difference between alpha = 0.05 and alpha = 1.