Solving equation $a^{-x} + \log x/\log a = 0$

Question

Please can you instruct me how should I start writing an algorithm (pseudo-code, to be implemented) for finding all solutions for the following equation: $a^{-x} + \log x/\log a = 0$

where $a$ ($a$ will be between $(2, 20)$ ) is a constant and the equation should be solved for $x$. I already know a couple of algorithms for solving simple polynomial equations but I would like to write code in any programing language (using basic operations, no special functions available like find_root) which finds the solutions for the equation with no special assumptions. First step I want to try to assume that my $a=16$ and hence implement algorithm for finding solutions for this particular form.

Answer 1

There are two problems involved here: One, is a general local root finding problem, for which many method's (such as Newton's) can be used. The second, is the problem of finding all (at most three) solutions.

One method I could suggest is the "Homotopy Continuation Method". To use it, first find all three solutions for some $a> e^e\sim15.15$. Now, vary $a$ slightly $a\to a+\Delta a$, and use the previous roots $x_i$ as initial values for your solver. Continue this process literately, until you find the solution curves for all $a$.

You can see how the solutions converge around the bifurcation $a= e^e$, as calculated using the attached python script.

Addendum: The bifurcation is at $a= e^e$, since three solution can only exist when $df/dx=0$. This happens when $a^x = \log^2 (a) x$, and since the Lambert W function is not real valued below $-1/e$, we get that $a\ge e^e$.

Code:

#!/usr/bin/python
from math import log

def f(x,a):
    return a**(-x) + log(x) / log(a)
def f_over_fprime(x,a):
    return x * (log(a) + (a**x) * log(x) ) / ( (a**x) - x * log(a) * log(a))

# find initial roots for a = 16

a = 16.0

roots = []

for x in [x * 0.1 for x in range(1, 10)]:
    y = x
    while (y > 0 and abs(f(y,a)) > 1.0e-14):
        y = y - f_over_fprime(y,a)
        if (y < 0): y = 0
    found = False
    for x in roots:
        if (x == round(y,14)): found = True
    if (found == False): roots.extend([round(y,12)])

print "Roots for a = 16: ", roots

for a in [a * 0.1 for a in range(160, 20, -1)]:
    y = roots
    for i in [0,1,2]:
        while (y[i] > 0 and abs(f(y[i],a)) > 1.0e-14):
            y[i] = y[i] - f_over_fprime(y[i],a)
            if (y[i] < 0): y[i] = 0
    roots = y
    print "Roots for a = ", a, "are: ", roots