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

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]


1. 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]]


1. Keep replacing the above rule until we get rid of all "L"s in the equation:
     L == ReplaceRepeated[Last[%], %]


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)$$

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)$$

• Thanks for your answer. I was really hoping for an analytical expression like the one provided in the "by-hand" computation. This is a perturbative calculation and so theoretically should only require inversion of the polynomial at zeroth-order in c. Commented May 2 at 11:18
• Also note that the Series is to be taken in 1/c; therefore on the third argument of "g", not the first. Commented May 2 at 11:19

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]

• Often, SolveAlways is used in such manipulations. Commented May 2 at 21:51
• That's great insight @yarchik. I did not know about this function. I shall update my answer to use it. Commented May 4 at 13:21