How to invert a series with two variables, where the series is expanded in the other variable?

Question

I have defined an expression that is a series in powers of $1/c$:

H[L, c] = -(1/(2 L^2)) + 1/c^2 (-(3/(J L^3)) + (15 - m)/(8 L^4)) + 1/c^4 (-(27/(2 J^2 L^4)) + (5 (-7 + 2 m))/(4 J^3 L^3) - (3 (-35 + 4 m))/(4 J L^5) + (-145 + 15 m - m^2)/(16 L^6)) + O[c, Infinity]^6;

defined for any $L > 0$.

I would like to obtain a series expansion for its inverse in the first argument; namely, a $1/c$-series expression: $$L[H, c] = L_0[H] + \frac{1}{c} L_1[H] + \frac{1}{c^2} L_2[H] + \ldots + \mathcal O[\frac{1}{c}]^6 \text. $$

I have tried using the series-manipulating built-in functions without success. Specifically, I tried InverseSeries, but my series is defined in $1/c$, whereas the variable I want to invert is $L$. I also tried

Series[InverseFunction[L |-> Evaluate[H[L, c]]], {c, Infinity, 4}]

but nope; that gives me a series in $\mathcal O\left( c + 1/(2L^2) \right)$, whatever that means.

I have so far been forced to perform the computation by hand. I would like to know, is there a more direct way to do this, using Mathematica's built-in functionality, that I'm missing?

---

The computation "by hand" goes through the following steps:

1.Split H[L,c] into zeroth-order term in c, plus remainder: $H = H_0 + \delta H$

    H0[L] = SeriesCoefficient[H[L, c], 0]

    δH[L, c] = H[L, c] - H0[L]

2. Solve for the zeroth-order term, and take the positive solution

    Last@Flatten@Solve[H == H0[L] + δH, L]

and then replace $\delta H(L, c)$ into the solution, so that we obtain an equation that gives $L$ as a function of $H$, $L$ and $c$.

    % // ReplaceAll[δH -> δH[L, c]]

4. Keep replacing the above rule until we get rid of all "L"s in the equation:

     L == ReplaceRepeated[Last[%], %]

Answer 1

You can use AsymptoticSolve for this purpose:

Lookup[L] @ First @ AsymptoticSolve[h == H[L, c], {L, 1/Sqrt[-2h]}, {c, Infinity, 5}] +O[c,Infinity]^5 //TeXForm

$ \frac{1}{\sqrt{2} \sqrt{-h}}+\frac{-h J m+15 h J+12 \sqrt{2} \sqrt{-h}}{4 \sqrt{2} c^2 \sqrt{-h} J}-\frac{3 \sqrt{-h} h J^3 m^2+30 \sqrt{-h} h J^3 m+35 \sqrt{-h} h J^3+96 \sqrt{2} h J^2 m-240 \sqrt{2} h J^2+80 \sqrt{2} m-280 \sqrt{2}}{32 c^4 \left(\sqrt{2} J^3\right)}+O\left(\left(\frac{1}{c}\right)^5\right) $

Answer 2

By inversion polynomials of order 6 are involved, so computation is working numerical only, if the abstract root object representation has to be avoided.

    H[L_, J_, c_, m_] = 
   -(1/(2 L^2)) + 1/c^2 (-(3/(J L^3)) + (15 - m)/(8 L^4)) + 
     1/c^4 (-(27/(2 J^2 L^4)) + (5 (-7 + 2 m))/(4 J^3 L^3) - 
    (3 (-35 + 4 m))/(4 J L^5) + (-145 + 15 m - m^2)/(16 L^6))

$$\frac{\frac{5 (2 m-7)}{4 J^3 L^3}-\frac{27}{2 J^2 L^4}-\frac{3 (4 m-35)}{4 J L^5}+\frac{-m^2+15 m-145}{16 L^6}}{c^4}+\frac{\frac{15-m}{8 L^4}-\frac{3}{J L^3}}{c^2}-\frac{1}{2 L^2}$$

     g =InverseFunction[
         Function[{L, J, c, m} , H[L, J, c, m]],
         1, 4]

    Function[{L, J, c, m}, 
    Root[145 J^3 - 15 J^3 m + 
      J^3 m^2 + (-420 J^2 + 48 J^2 m) #1 + (216 J - 30 J^3 c^2 + 
      2 J^3 c^2 m) #1^2 + (140 + 48 J^2 c^2 - 40 m) #1^3 + 
      8 J^3 c^4 #1^4 + 16 L J^3 c^4 #1^6 &, 1]]

     Series[g[x, 2, 3, 7], {x, 0, 2}] 

    % // N

$$-(0.558054\, +0.406841 i)+(0.0193974\, +0.224663 i) (x+0.)+(0.23877\, -0.21318 i) (x+0.)^2+O\left((x+0.)^3\right)$$

Answer 3

I have found an answer using Solve. It's not a single function call as I'd hoped, but it could probably be made into a function. Perhaps other users will have improvements to this answer.

Consider H[L_, c] := ... as defined in the question. Define a series representation for the first-argument inverse, L[H,c]:

L[H, c] = L0[H] + 1/c^2 L2[H] + 1/c^4 L4[H]

Here only the even terms are required as can be seen in the definition of H.

Now we can setup the equation H == H[L[H, c], c] and ask Mathematica to Solve for the coefficients at each order in 1/c. For convenience, to avoid imaginary numbers in the solution, I have defined a dummy variable $\varepsilon = - 2 H$.

eqs = Map[# == 0 &, CoefficientList[(-ɛ/2) - H[L[H, c], c] + O[c, Infinity]^6, 1/c]];
Print /@ Select[eqs, ! MatchQ[#, True] &];

sol = Last@Solve[eqs, {L0[H], L2[H], L4[H]}] /. ɛ -> -2 H

Here, Last was used to take the positive solution, as discussed in the question. The final result is:

L[H, c] /. sol

It can be verified that this solution and the t"by-hand" solution using ReplaceRepeated agree.

---

Edit: As pointed out by @yarchik, SolveAlways is designed precisely for setting polynomial equations to zero and solving for coefficients. While I couldn't get SolveAlways to solve for L0, L2, and L4 directly, an equivalent formulation can be achieved using an Eliminate one-liner:

Solve[Eliminate[(-ɛ/2) == H[L[H, c], c] + O[c, Infinity]^6, c], {L0[H], L2[H], L4[H]}] // Simplify // ReplaceAll[ ɛ -> -2 H]