Series approximation(s) of a difficult recursive equation

Question

New user here. I'm working on trying to get asymptotic solutions to the following recursive function:

$f(r)=\frac{1}{r-k}\lgroup\sqrt{\frac{2}{k^2-1}}+\sqrt{\frac{1}{2k^2-1}}\rgroup$ (Eqn. 1)

where $k$ is such that

$r=\frac{1}{k}\lgroup\frac{\sqrt{2}(k^2-1)^{-1/2}+(2k^2-1)^{-1/2}}{\sqrt{2}(k^2-1)^{-3/2}+(2k^2-1)^{-3/2}}\rgroup+k$. (Eqn. 2)

Now I have solved this function numerically (in matlab), basically by taking the value of $r$ and solving Eqn 2 numerically for $k$, and then using that value of $k$ in Eqn 1 to find $f(r)$. The function goes to positive infinity at $r=1$, and goes to 0 as $r\rightarrow\infty$ (I would post an image of it here, but it's not letting me post an image as a new user).

The function becomes complex for $r<1$, but that isn't relevant because the physical system this is based on is not defined for $r<1$.

I want to get series approximations for $f(r)$ as $r\rightarrow1$ and $r\rightarrow\infty$. I plotted the same function on a log-log plot, while shifting the x-axis over 1 so that I could see the asymptotic behavior near $r=1$ (again I would post an image of the plot, but I cannot). Basically, on the log-log plot it goes to a straight line with a negative slope as $r\rightarrow1$, and then a slightly different slope as it goes to $r\rightarrow\infty$.

Based on this, it certainly should have a well-behaved series expansion at both of these extremes, if I can just figure out how to evaluate it.

My idea so far is to do the following:

1. Get series expansions for Eqn 2 for $r(k)$ at the two asymptotes.
2. Invert series expansions for Eqn 2 to get series solutions for $k(r)$.
3. Substitute these series solutions into Eqn 1 and then evaluate series expansions for $f(r)$.

I've done #1 (using the series expansion function in the symbolic math toolbox in matlab), and gotten the two following series:

As $r\rightarrow1$, $r(k)=1+3(k-1)+2(k-1)^{3/2}-(k-1)^2+O(k^{5/2})$

Similarly, as $r\rightarrow\infty$, $r(k)=2k-\frac{5}{6k}-\frac{1}{12k^3}-\frac{5}{96k^5}+O(k^{-7})$.

After this though, I don't know what to do. I've looked up how to invert standard Taylor series, but these are Laurent and Puiseux series, and I can't figure out how to invert them. When inverting Taylor series you already know what powers will be in the inverted function, so then you just solve for the appropriate coefficients. In these though, I don't know a priori what the powers of the inverted series would be so I can't just solve for them using the same method.

I'm pretty sure that this could be easily done in Mathematica or Magma (since they both have a series inversion capability) but I don't have access to such software, and the symbolic math toolbox in matlab, which I do have, doesn't have such capability.

Answer 1

I think Groebner bases could be helpful by proceeding as follows. First transform your equations into algebraic equations. For instance, your first equation implies $$((r-k)^2f^2-\frac{2}{k^2-1}-\frac{1}{2k^2-1})^2-2\frac{2}{(k^2-1)(2k^2-1)}=0$$ (without error of my part) which becomes a polynomial equation after multiplication by $((k^2-1)(2k^2-1))^2$.

Similar operations yield a polynomial from the second equation. Use Groebner base techniques for eliminating $k$ in your polynomial system. This yields a polynomial $\mathbb C[f,r]$ and choosing the correct branches you get series approximations (or Puiseux-series approximations at branch points).