Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I'm implementing a root-solver for finding x coordinates of a function f(x), after I have an y-coordinate. The function is periodic, roughly sinusoidal with constant amplitude but non-linearly varying frequency; for an inverse I don't have a closed-form (it is an infinite series), so I use the Newton iteration to find the x-value at a given y beginning the iteration at $x_0$ which is rather near the true value by something like $x=newton(x=x_0,f(x)-y)$.

In most cases this works fine, however if the y is in the near of a maximum (or minimum) of f, where the shape is very similar to the maximum of a sinus-curve, the newton-iteration does not converge. The wikipedia gives a bit of information about this, but not a workaround. The last way out would be to resort to binary search, but which I'd like to avoid since the computation of f(x) is (relatively) costly.

Does someone know an improvement in the spirit of the Newton-iteration (which has often quadratic convergence) for this region of y-values?

[update] hmmm... perhaps it needs only a tiny twist? It just occurs to me, that it might be possible to go along the way how I find the easily approximated x for the maximum y: here I use the Newton on the derivative of the function and search for the zero:
$x_{max}=newton(x=x_0,f(x)') \qquad $ and this has the usual quadratic convergence. But how to apply this for some y in the near of the maximum?

share|improve this question
    
You could approximate your function at $x$ by a parabola, using $f(x)$, $f'(x)$ and $f''(x)$, instead of a line using just the first two... –  Jaime Dec 7 '12 at 23:33
    
Do you mean that the starting value is very close to a local extremum or that the actual zero of the function is also an extremum? In the latter case, newton for $\sqrt f$ might help ... –  Hagen von Eitzen Dec 7 '12 at 23:34
    
@Hagen : the $y$ for which it is difficult to find the $x$ are in the near of $\sin(\pm \pi / 2)\cdot \alpha$ where $\alpha$ is the amplitude of my function $f(x)$ –  Gottfried Helms Dec 7 '12 at 23:42
    
@Jaime : looks like a sort of idea which I'm looking for. Would you mind to elaborate this a bit more to help me step in? (Surely I'll also need time to test&adapt...) –  Gottfried Helms Dec 7 '12 at 23:46
    
Maybe you could use quadratic interpolation along with bisection to (hopefully) converge faster? –  copper.hat Dec 8 '12 at 0:05

5 Answers 5

Have you considered implementing a hybrid method? For example, at each step:

IF a Newton step would result in an iteration that is outside the bounds where you have determined the root must lie, then take a bisection step (slower than Newton, but bisection always converges to a root and is not affected by extrema), or a step using a method other than Newton that is not prone to failing near extrema.

ELSE proceed with a Newton step (since it converges quadratically, as you pointed out).

share|improve this answer
    
Yes, I'm considering bisection, but I'd like to find something else with better rate-of-convergence –  Gottfried Helms Dec 7 '12 at 23:38
    
Then replace "bisection" with another method of your choosing? I chose bisection because it is reliable, and it will only be used at steps near a maximum/minimum. –  Eric Angle Dec 7 '12 at 23:56

By popular demand from the OP...

In Newton's method you are replacing your function $y=f(x)$ by a linear approximation around the point $x_0$, $y = f(x_0) + f'(x_0) (x-x_0)$, which intersects the x axis ($y=0$) at $x=x_0-f(x_0)/f'(x_0)$.

You could instead approximate by a parabola as $y=f(x_0) + f'(x_0)(x-x_0) +\frac{1}{2}f''(x_0)(x-x_0)^2$, which intercepts the x-axis at $x = x_0 -\frac{f'(x_0)\mp\sqrt{f'(x_0)-2f(x_0)f''(x_0)}}{f''(x_0)}$. You will of course have the issue of having two, not one, possible next iteration point, but there are multiple ways to get around these: choose the closest one, always move up (or down), choose the one with a smallest $f(x_0)$...

share|improve this answer
    
Ahh, very nice - it looks, that with my update I was even near already... So I'll just try it. That shall need some time, then I'll come back to this. Thanks so far! –  Gottfried Helms Dec 8 '12 at 0:14
    
Hmm, I couldn't manage to make this working with my specific application, maybe simply programming errors. Perhaps if I get the solver working by more detailed study of the process in the extreme cases I'll come back to this and try to locate my errors. Thanks anyway for that answer! –  Gottfried Helms Jan 3 '13 at 15:37

I couldn't tell you the cost of this idea, but maybe you could work it out:

You are trying to solve for a root of $g$, where $g(x)=f(x)-y$ and you want your solution in a neighborhood of $x_0$. Let's call that solution $x_s$. The problem is that $g'(x)$ is too small near $x_s$, so in the Newton algorithm you divide by very small things, yielding big changes from $x_i$ to $x_{i+1}$; possibly so big that the algorithm converges to a different solution or not at all.

So what if you engineered a substitute function $\tilde{g}$, who still satisfied the demand that $\tilde{g}(x_s)=0$, but has $\tilde{g}'(x_s)$ not so small. For example, $\tilde{g}(x)$ could equal $g(x)\cdot\ln|g(x)|$. This has $\lim_{x\to x_s}\tilde{g}(x)=0$ and has $\tilde{g}'(x)=g'(x)\cdot(1+\ln|g(x)|)$. The absolute value of $\tilde{g}'(x_s)$ will be quite larger than that of $g'(x_s)$.

share|improve this answer

In most cases this works fine, however if the y is in the near of a maximum (or minimum) of f, where the shape is very similar to the maximum of a sinus-curve, the newton-iteration does not converge.

It doesn't converge in what way? Knowing the failure mode is important in fixing them. There are three main possibilities: the iteration is converging slowly/asymptotically, the iteration is converging to the wrong value (ie diverging from the answer you want), or the iteration is approaching a limit cycle (bouncing around among 2, 3, or more values repeatedly).

If your problem is that convergence near a minimum is too slow, and the second derivative doesn't vanish there, too, you could try Halley's method. The trouble with Halley's method is that when your state is far from the root it can cause you to overshoot, reducing your basin of convergence. So, I take a cue from "Numerical Recipes in C," when using it I clamp the effect so that it isn't too small or too large:

\begin{align} $x_{n+1} = x_n - \left(\frac{f_n}{f_n'}\right) \cdot \mathrm{CLAMP}\left(1 - \frac{f_nf_n''}{2[f'_n]^2}, \frac{1}{2}, 2\right)^{-1}. \end{align}

The ability to speed up convergence by using higher order methods in this class (Househoder's methods) is limited, however, by the number of derivatives that are zero. A simple example of a function that no Newton like method will find the root of quickly: $e^{-1/x^2}$, (all of its derivatives vanish at $x=0$).

If the problem is that the iteration is diverging from the root you want to find then there's a really easy answer: iterate in the other direction. See, in any iterative process like this you will generically have attractors (points the iteration will converged to, if the starting point is in its basin), repulsors (points the iteration will diverge from), and limit cycles (attractors of the double step, triple step, etc). I just recently had a problem with Newton's method always trying to converge to the solution I didn't want in a two solution problem. Flipping the sign to change the direction of iteration fixed that up nicely, because reversed iteration flips the roles of attractors and repulsors.

If the problem is that you're hitting a limit cycle, then that's more tricky. See, if $x_{n+1} = g(x_n)$ is your update rule, then you can define a fixed point as $g(x) = x$. Where the trouble begins is $g$ implies a whole family of update rules, call them $g^{[n]}$, that are the same as applying $g$ $n$ times. The fixed points of $g$ will all be fixed points of $g^{[n]}$, but $g^{[n]}$ will generically have more. Because iterating in $g$ is the same as iterating in $g^{[n]}$, the fixed points of $g^{[n]}$ can become the attractors that dominate the behavior of the sequence, pulling it into a limit cycle of length $n$. If we assume, as is often the case, that the fixed point you want is a reuplsor within the basin of attraction of a limit cycle, then you might be able to find that point if you: 1, detect the cycle; 2, synthesize a guess that is not to close to the limit cycle points (eg average them, or use the $x_0$ that converged to the cycle in the first place); and 3, iterate backwards.

I can't offer any guarantees without further information, though, because those basins of attraction are generically fractals in shape.

share|improve this answer

If the minimum/maximum is negative, then an $x_0$ such that $f(x_0)$ is positive is preferred or if the minimum/maximum is positive then an $x_0$ such that $f(x_0)$ negative is more reasonable but if the root is at the minimum/maximum then I don't think there should be a problem.

Something else with a better rate of convergence is the Secant Method. This has a convergence rate of $\cfrac{1+\sqrt 5}{2}=1.618...$ mostly due to the fact that it starts with two initial values.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.