Method to solve system of exponential sums of the form $a^x+b^x=c$ given more equations than variables Cross post with mse
For example, let's say I have the following equations.
\begin{gather*}
a^{x-1}+b^{x-1}=337 \\
a^{x}+b^{x}=1267 \\
a^{x+1}+b^{x+1}=4825 \\
a^{x+2}+b^{x+2}=18751.
\end{gather*}
What we can notice is that we can rewrite the left-hand sides of all the equations as
$a^{x+k}(1+(b/a)^{x+k})$
for our values of $k$: $-1$, $0$, $1$, and $2$.
Let's label those formulas as $f(x)=a^x(1+(b/a)^x)$.
If we divide $f(x)^2$ by $f(x+1)f(x-1)$ we will get a new formula that is dependent on $x$ and $b/a$:
$$(1+(b/a)^x)^2/((1+(b/a)^{x-1})(1+(b/a)^{x+1}).$$
We can label $b/a$ as $d$ and $(b/a)^x$ as $f$ (probably should find better variable names).
And then we will have the equation
$$(1+f)^2/((1+f/d)(1+fd))=1267^2/(337 
* 4825).$$
Similarly. Another equation can be created for $f(x+1)^2/(f(x+2)f(x))$:
$$(1+f d)^2/((1+f * d^2)(1+f))=4825^2/(18571 
* 1267).$$
While my first set of 4 (or only 3) equations was unsolvable by Wolfram|Alpha the second set is solvable and 2 pairs of values for $f$ and $d$ are returned.
One of them is
$f=4.21\ldots$ and $d=1.3333 $.
$\log(f)/\log(d)= 5$ which is the true value of $x$.
$d=a/b$ is $4/3$ with numerical error.
Finding $a$ and $b$ is pretty trivial from here and they are 4 and 3 respectively.
My question is can this method be generalized to more than 3 variables (meaning $x$ and more exponentials)?
On paper, it seems that this could be done with a slightly different technique to generate multiple equations based on the principle of removing the dependency on the least significant exponent. In practice, I wasn't quite able to do so yet.
Does anyone know of any research on this topic or if this question was already solved, and in general is there any technique to analyze the sum of exponentials based on multiple values?
To be more clear.
My idea is to take multiple points of a function of the form for example
$a^x+b^x+c^x$ (preferabliy if the number of exponentials is unkown)
Rewrite it as $a^x((b/a)^x+(c/a)^x+1)$
And find the points of the same x value for the function $(b/a)^x+(c/a)^x+1$ And Continue until only one exponential is left
Does anything like that exist or is even possible?
Like this method https://en.m.wikipedia.org/wiki/Newton_polynomial
But for exponentials
 A: I'm not sure about your method, but such equations can easily be turned into polynomial equations by introducing $u:=a^{x-1}$ and $v:=b^{x-1}$:
\begin{cases}
u+v=337 \\
au+bv=1267 \\
a^2u+b^2v=4825 \\
a^3u+b^3v=18751.
\end{cases}
Eliminating variables $a,b,v$, we obtain the equation for $u$:
$$2163001u^2 - 728931337u + 6879707136=0,$$
which determines its values and the values of $x$. The other values are then follow from the above system. Since $u$ and $v$ are not independent variables, we will need only solutions with $\log_a u = \log_b v$.
This approach can be generalized to more variables.
A: We can use linear algebra to make this easier.
Consider the space of sequences of the form $A a^n + B b^n$. We can use several dual bases for this space. One is the obvious basis of maps that send $A a^n + B b^n$ to its coefficients, $A$ and $B$. If neither of $a$ and $b$ are zero and the two are distinct, then another is the basis of maps that send $A a^n + B b^n$ to its values at $n = 0$ and $n = 1$.
Now consider the self-map "increase $n$ by $1$", i.e. taking the sequences $A a^n + B b^n$ to the sequence $A a^{n + 1} + B b^{n + 1}$. In the basis of coefficients, this can be seen as a diagonal map with eigenvalues $a$ and $b$. In other words, it has the corresponding matrix
$A = \left[\begin{matrix}
a & 0 \\
0 & b
\end{matrix}\right]$
In the basis of values at $n = 0, 1$, it clearly moves the value for $n = 1$ into the value for $n = 0$, while doing something to give the new value for $n = 1$, so it has a matrix of the form
$B = \left[\begin{matrix}
0 & 1 \\
z & y
\end{matrix}\right]$. Note that $B$ is conjugate to $A$ - in fact, it's the companion matrix for the characteristic polynomial of $A$.
Now we can approach your specific values. It's not hard to see that we have the equations
$B \left[\begin{matrix} 337 \\ 1267 \end{matrix}\right] = \left[\begin{matrix} 1267 \\ 4825\end{matrix}\right]$ and $B \left[\begin{matrix} 1267 \\ 4825 \end{matrix}\right] = \left[\begin{matrix} 4825 \\ 18751\end{matrix}\right]$
Putting them together, we have that $B \left[\begin{matrix} 337 & 1267 \\ 1267 & 4825 \end{matrix}\right] = \left[\begin{matrix} 1267 & 4825 \\ 4825 & 18751\end{matrix}\right]$
We can then find $B$:
$B = \left[\begin{matrix} 1267 & 4825 \\ 4825 & 18751\end{matrix}\right]\left[\begin{matrix} 337 & 1267 \\ 1267 & 4825 \end{matrix}\right]^{-1} = \frac{1}{576} \left[\begin{matrix} 0 & 576 \\ -13247 & 5717\end{matrix}\right]$
But then we have the characteristic polynomial of $A$ - which the eigenvalues of $A$, the diagonal elements $a$ and $b$, must satisfy. So $a$ and $b$ are the roots of $c^2 - \frac{5717}{576} c + \frac{13247}{576} = 0$. Then $x$ is whichever value satisfies the original equations.

This process generalizes quite readily. If you have $m$ exponentials and $2m$ consecutive equations of the form $\sum_j a_j^i = C_i$ (or even more generally, $\sum_j A_j a_j^i = C_i$), you can find $a_j$ by writing a square matrix of the form
$\left[\begin{matrix}
C_1 & C_2 & \dots & C_m \\
C_2 & C_3 & \dots & C_{m + 1} \\
\vdots & \vdots & \ddots & \vdots \\
C_m & C_{m + 1} & \dots & C_{2m - 1}
\end{matrix}\right]$
and a vector of the form $\left[\begin{matrix}
C_{m + 1} & C_{m + 2} & \dots & C_{2m}
\end{matrix}\right]$, and multiplying the vector by the inverse of the matrix. This works as long as you actually need $m$ such exponentials - if you find that the first matrix isn't invertible, then decrease $m$ by $1$ and try again.

This is one of those funny times where it turns out to be easier to answer a more general question than to answer the specific question, because the more general question has a nicer structure than the specific question.
