Newton–Raphson method impossible when a fuction is almost "flat"

Question

I'm trying to compute the value of x^x.
I have calculated it's derivative and I'm
using the Newton–Raphson method to find the root of f(x)=x^x-a for a given a.

This is the function I use to calculate the next approximation:

double xkp1(double xk, double a){
    double lnxp1 = log(xk) +1;
    return xk -( (fx(xk, a))  /  (exp(xk*log(xk)) * (log(xk) + 1)  ));
}

where function fx is defined this way:

double fx(double x, double a){
    return exp(x*log(x))-a;
}

Now, this approach works correctly only if the starting value of xk is already close to the root.
If it differs by +-0.5, xk just explodes to a really high value and f(x) becomes infinity.
Now i think what might be wrong here - the derivative of x^x is very small in comparison to the actual value of x^x, so the whole (fx(xk, a)) / (exp(xk*log(xk)) * (log(xk) + 1) ) just becomes +infinity - the tangent line overshoots the root.
This means I'd need higher double precision, but is it possible in any way? If not, can something be modified in the method to make it applicable in this situation?

What are your test values of a? Also, you could instead take the logarithm of both sides and solve x*log(x) = log(a) - far less likely to overshoot and end up with astronomical values. — meowgoesthedog

Yngve Moe Yngve Moe · Accepted Answer · 2019-10-12T19:30:19

I'll assume that you are only interested in positive x values, not negative integers (which also have real values).

I have two solutions that may help you with your problem. Firstly, the numerical scheme might be more stable if you first use a log transform. Your equation would then be to find a such that x log(x) = log(a). The derivative of x log(x) is log(x) + 1, which may be more stable than your exponential function.

Otherwise, you can use you knowledge in that f(x) = x^x is monotonically increasing on x>1/e to run a bisecting search. Start with x_l = 1/e and x_r = 2/e. If f(x_l) > a, then there is no solution. Then while f(x_r) < a, set x_l = x_r and x_r = c * x_r with c > 0(or any other scheme to increase x_r). Once you have a pair, x_l, x_r that follow f(x_l) < a < f(x_r), then you can begin a bisecting search to find a small interval that contain the solution to f(x) = a. Once this interval is small enough, you can start with iterations of Newton's method.

You may even combine the two methods above, and perform a bisecting search to solve x log(x) = log(a).

Newton–Raphson method impossible when a fuction is almost "flat"

1 Answers