# Newton's root finding method that copes with stationary points

I am using Newton's method to solve an equation of the form $$y=f(x)$$ for increasing values of $$y$$. The value of $$x$$ found for some $$y$$ is used for the next $$y$$ as a starting approximation and this works fine.

$$f$$ is increasing but may have stationary points, and I do have situations where the range of interest of $$y$$ starts or ends on an extremum or very close. This causes the algorithm to experience a division by zero or yield very inaccurate roots.

Is there a simple way to cope ? Is there a known modified Newton's method able to deal with this situation ?

I know that double roots can be handled by the modified iterations $$x\leftarrow x-2\frac{f(x)}{f'(x)}$$ but this still does not work at the root, and one needs to know when to activate the factor $$2$$.

For information, my function $$f$$ is a cubic polynomial, but I want to avoid the explicit computation of the roots by Cardano's formulas.

Update:

If we restate Newton's method at a stationary point, we get to the next term of Taylor's development,

$$f(x+h)-y\approx f(x)-y+\frac{h^2}2f''(x),$$

so that $$h\approx\sqrt{-\frac{f(x)-y}{f''(x)}}.$$

In fact, the approximation

$$f(x+h)-y\approx f(x)-y+h f'(x)+\frac{h^2}2f''(x)$$ always holds, giving the general formula

$$h\approx \frac{f'(x)\pm\sqrt{f'^2(x)-2f''(x)(f(x)-y)}}{f''(x)}.$$

Due to the increased complexity, I don't want to use that formula all the time.

• Assuming $f(x)= ax^3+bx^2+cx+d$, you know $f'(x)$ and $f''(x)$ and the roots of $f'(x)$. When $x$ approaches the roots of $f'(x)$, could you use l'hopital's rule and switch to $f'(x)/f''(x)$? Mar 15, 2021 at 11:58
• @PierreCarre: come on, in the interval of interest.
– user65203
Mar 15, 2021 at 14:03
• @PierreCarre: good point.
– user65203
Mar 15, 2021 at 14:25
• Incidentally, there' s an old but pretty good book, Numerical Recipes, that covers rootfinding in detail including various alternatives to Newton. You might be particularly interested in en.wikipedia.org/wiki/Brent%27s_method . Mar 15, 2021 at 16:10
• @StevenStadnicki: of course. I need very fast convergence (at most two iterations), otherwise the cost of the iterations might outweigh that of Cardano's approach.
– user65203
Mar 15, 2021 at 16:15

There is an easy modification called modified Newton's method. The idea is easy: every time the derivative is zero, one uses a value of the derivative from the last iteration, when it was nonzero (couple of times if needed). This changes the method to the secant kind of method for these couple of iterations. When it becomes nonzero again, one switches back to Newton's method. In more dimensions it works as well.

Additional thing, which you can always try in case of convergence problems, is the damped version of Newton's method. It means, that you reduce each increment by some fixed value smaller than $$1$$, for example by $$0.9$$.

Can you post the function you solve and your starting point?

• Thanks. But what when you start the iterations at a minimum ? The functions are any cubic in a monotonic section, possibly starting or ending on an extremum. If you want a case, $-x^3+3x$ from $(-1,-2)$.
– user65203
Mar 15, 2021 at 11:29
• Interesting. Then you have to start from a different point. (Do you need to find one solution or all solutions?) You can try to generate starting points randomly in the solved domain and run Newton's method from these points unless you are successful. Mar 15, 2021 at 11:37
• Or you can try to examine functional value in more points (equidistantly spaced) and as a starting point choose the one, in which the value is closest to zero. Then you get a starting point of higher quality. Mar 15, 2021 at 11:47
• Yes, this is an option. In fact, for my application it is more natural to process by increasing $y$ and this is why I would prefer to have a solution that starts from the starting point.
– user65203
Mar 15, 2021 at 13:04

Momentum-based algorithms are good at passing through stationary points while retaining fast convergence.

$$m_i=(1-\beta)f'(x_i)+\beta m_{i-1}=m_{i-1}-(1-\beta)(m_{i-1}-f'(x_i))$$ $$x_{i+1}=x_i-(1-\beta^i)\frac{f(x_i)}{m_i}$$

The $$1-\beta^i$$ factor serves as bias correction, supposing one takes $$m_0=0$$, since we get $$m_i=(1-\beta^i)f'(x_i)$$ if $$x_i$$ is fixed.

As $$\beta\to0$$, you get Newton's method and faster convergence, while for larger $$\beta$$ you will get better "resilience" against stationary points.

For your purposes, you know that $$f'(x)\ge0$$, so rather than setting $$m_0=0$$, you may want to initialize it to a small positive value, or perhaps even use the same $$m_i$$ sequence for every $$y$$, since it is likely you will converge to the solution for $$y=f(x)$$ for a given $$y$$ while $$m_i$$ is sufficiently large that starting from where it left off works well for the next $$y$$.

If you go with the latter approach, you'll want to choose $$\beta$$ based on the density of your $$y$$ around the stationary points. If $$\beta$$ is too small and you have many $$y$$ near a stationary point, then $$m_i$$ will decay rapidly, which leaves you with the original issue.

If I understood correctly, the basic step is to solve the equation $$f(x)-y=0$$, for given $$y$$. If we are unsure of the multiplicity of a root, we can apply the modified Newton's method, which corresponds to the application of the usual Newton's method to $$g(x)=\frac{f(x)-y}{f'(x)}$$. This way, the quadratic convergence is recovered, at least in exact arithmetic. The iteration would be $$x_{n+1} = x_n-\dfrac{\frac{f(x_n)-y}{f'(x_n)}}{\frac{f'(x_n)^2 - (f(x_n)-y)f''(x_n)}{f'(x_n)^2} } = x_n - \dfrac{f'(x_n)(f(x_n)-y)}{f'(x_n)^2-f''(x_n)(f(x_n)-y)}$$ If, by any chance, $$f'(x_0)=0$$, just move it a bit.

• $f'(x_0)=0$ (as well as $f'(x_0)\approx0$) is the topic of my question.
– user65203
Mar 15, 2021 at 15:15
• @YvesDaoust I understand. I'm just saying that using this method $f'(x_0)\approx 0$ does not pose a problem. Mar 15, 2021 at 15:22
• With $f'(x_0)=0$, the iterations don't start.
– user65203
Mar 15, 2021 at 15:26
• @YvesDaoust That is why I proposed that "If, by any chance $f'(x_0)=0$, just move it a bit". It would be part of the algorithm to check if $f'(x_0)=0$ and act accordingly. Mar 15, 2021 at 15:29
• You could probably get away with setting $f'(x_0)$ to $ε$ instead of evaluating $f'(x_0+ε)$ to save a few cpu cycles. Mar 17, 2021 at 2:15