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
add comment

4 Answers

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
add comment

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
add comment

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
add comment

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
add comment

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.