2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Are you familiar with Newton's Method? $\endgroup$ – mweiss Jul 20 '14 at 14:12
  • $\begingroup$ As something in conjunction with @mweiss's answer, note that $x=1$ gives a decent seed for Newton's method especially as $a$ gets bigger. $\endgroup$ – Semiclassical Jul 20 '14 at 14:15
  • 1
    $\begingroup$ The question being is how to find all the solutions (or roots), provided a=16 there are 3. Standard Newton's method gives you one solution at a time and you have to provide starting point (or seed if you like). I do not mind Newton's method (I implemented this a couple of times) but I would to make sure if it is worth further investigation. How to modify it or combine with a different method, that is the question. $\endgroup$ – algorytmus Jul 20 '14 at 14:30
  • $\begingroup$ Your equation can be rewritten as $x=a^{-a^{-x}}$, which gives us a simple recurrence formula, easily implemented in many programming languages. I've just tested it in C, and it works like a charm. But, of course, it only returns a single value of all possible three. $\endgroup$ – Lucian Jul 20 '14 at 18:55
  • 1
    $\begingroup$ I think the fact that this problem has a pitchfork bifurcation of one root into three at $a\approx 15$ is the principally interesting part of this. So you may be able to attract more interest if you emphasize that in your question text. $\endgroup$ – Semiclassical Jul 21 '14 at 0:57
4
$\begingroup$

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: enter image description here

#!/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
$\endgroup$
  • $\begingroup$ I'm a bit curious as to what would happen in that latter method if you started with $a$ above the bifurcation and then continued it lower. Presumably one would see them all collide. $\endgroup$ – Semiclassical Jul 21 '14 at 16:24
  • $\begingroup$ @Semiclassical - they should indeed collide, we'll see :) $\endgroup$ – nbubis Jul 21 '14 at 16:25
  • $\begingroup$ Since the solutions presumably don't disappear but rather become complex below the bifurcation, probably one won't see a true 'collision' (but eh, when are numerical methods ever exact?) $\endgroup$ – Semiclassical Jul 21 '14 at 16:27
  • $\begingroup$ @Semiclassical - they indeed collide, see above. $\endgroup$ – nbubis Jul 21 '14 at 17:29
  • 2
    $\begingroup$ The bifurcation occurs at $a = e^e$. $\endgroup$ – Antonio Vargas Jul 21 '14 at 18:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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