algebraic manipulation - InverseSeries of multiple variables and multiple equations

CONTEXT

Let us consider a bit of the Universe in which we draw spheres

(see a high resolution image here). Astronomers have shown that the density within these spheres could be predicted quite accurately: here is the measured (in red) and predicted (in green) distribution of the density within $50^3$ such spheres at half the age of the Universe: .

Our next purpose is to extend the theory to concentric two or more shells.

In practice, in order to estimate the mildly nonlinear cosmic density of the Universe within concentric shells, I need to find perturbatively the Legendre Transform of a function of two variables in a singular regime (see below).

$$\phi(\lambda_1,\lambda_2)={\rm LT}(\psi(\rho_1,\rho_2))\equiv \sup_{\rho_1,\rho_2}\left[\lambda_1 \rho_1+\lambda_2 \rho_2-\psi(\rho_1,\rho_2)\right]$$

which in turns involves inverting the system $$\partial_\rho \psi =\lambda \quad \quad \quad (1) $$ for $\rho_1,\rho_2$ , and integrating for $\phi(\lambda_1,\lambda_2)$ the system $$ \quad \partial_\lambda \phi =\rho \,. \quad \quad (2) $$ This can in principle be done perturbatively. In my context, it needs to be done in the regime where ${\rm det}| \partial^2_\rho \psi|=0$. I.e. I am interested in Taylor expanding the Legendre transform of $\psi$ near a point (chosen to be zero for simplicity) where one of the eigenvalues of $\partial^2_\rho \psi$ is zero. In physical terms it corresponds to the rare event tail of the density of the Universe in these shells, if you care to know!

ATEMPT

Following this question I know how to do invert (1) for $\rho(\lambda)$ in 1D (one shell) for the regular

 nn = 3;

 ρofλr = InverseSeries[Series[ψ'[ρ], {ρ, 0,nn}]] /. ρ -> λ /. ψ'[0] -> 0 /. 
   Derivative[n_][ψ][0] :>  Subscript[ψ, n] // Normal

(note the division by ψ''[0])

and the singular case (for which ψ''[0]=0)

 nn = 3;
  ρofλs = InverseSeries[Series[ψ'[ρ], {ρ, 0,nn}]/. ψ''[0]-> 0] /. ρ -> λ /. ψ'[0] -> 0 /. 
    Derivative[n_][ψ][0] :>  Subscript[ψ, n] // Normal

note the square root. It is clear from this expansion that the Legendre transform, ϕ[λ] will have a very different algebraic form in the singular case compared to the regular case.

 ϕ[λ_]=Integrate[ρofλ,λ]

But this is not enough for concentric shells: I need to be able to carry out such Legendre transform when one coordinate is singular.

SOLUTION to first order and for the regular case only

The following works for the regular case only. Let us expand the first system

   nn=1;
   eqn = Thread[{λ1, λ2} == 
   Series[{Derivative[1, 0][ψ][ ρ1,  ρ2], 

   Derivative[0, 1][ψ][ ρ1,  ρ2]}, { ρ1, 0, 
   nn}, { ρ2, 0, nn}]] // Normal;

and reoder it in power of ϵ

  eqn2 = eqn /. 
  Derivative[i_, j_][ψ][0, 
   0] -> ϵ^(i + j - 1) Derivative[i, j][ψ][0, 0] // 
 Series[#, {ϵ, 0, 1}] & // Simplify;

Solve for ρ1,ρ2

 sol = Solve[Normal[eqn2], { ρ1,  ρ2}][[1]] /. 
 Derivative[n_, p_][ψ][__] :>  Subscript[ψ, n, p] // 
 FullSimplify;

We can then integrate the system (2):

 eqn3 = {D[ϕ[λ1, λ2], λ1] == ρ1,D[ϕ[λ1, λ2], λ2] ==  ρ2} /. sol
 ϕ[λ1, λ2]/. DSolve[eqn3,ϕ[λ1, λ2], {λ1, λ2}][[1, 1]] // Apart

But this is not good enough because I am interested in a singular expansion, i.e. near a point where the Jacobian $\partial^2 \psi \partial \rho_i \partial \rho_j$ has zero determinant! (when the above first order solution becomes singular because ψ02 ψ20 -ψ11^2 =0)

QUESTION

I am interested in doing this in 2D (or 3D...) for the singular case, i.e. when the Jacobian $\partial^2_\rho \psi$ has one null eigenvalue.

The main difficulty is that InverseSeries does not work for series of two variables.

The gist of the problem is the following: assuming x1, x2 y1 and y2 are small, how does one invert pertubatively

 y1= x1+ x1^2 + x2^2 + x1 x2 + x2^3...
 y2= x2^2 + x1^2 +x1 x2 + x1^3 +... 

in order to write

 x1= y1+...
 x2= sqrt(y2)+...

This might sound like a technical question, but the core of the problem is fairly general: how does one InverseSeries of multiple variables and multiple equations? in Mathematica. Any suggestions would be very welcome!

Answer

Algebraically, the way to construct the inverse of a series is straightforward. The basic iterative step is to add the next degree terms and solve for the coefficients.

If we have $u = U(x, y)$, $v = V(x, y)$, then we are looking for $x = X(u, v)$, $y = Y(u, v)$ such that $u = U(X(u,v), Y(u,v))$ and $v = V(X(u,v), Y(u,v))$ identically. If the coefficients of the partial series $$ X_n(u,v) = \sum_{k=0}^n\sum_{i=0}^k a_{i,k-i} u^i v^{k-i}, \quad Y_n(u,v) = \sum_{k=0}^n\sum_{i=0}^k b_{i,k-i} u^i v^{k-i} $$ have been determined, then plug $$ x = X_n(u,v) + \sum_{i=0}^{n+1} a_{i,n+1-i} u^i v^{n+1-i}, \quad y = Y_n(u,v) + \sum_{i=0}^{n+1} b_{i,n+1-i} u^i v^{n+1-i} $$ into $$ u = U(X(u,v), Y(u,v)), \quad v = V(X(u,v), Y(u,v)) $$ and set the coefficients of the monomials $u^iv^j$ to zero and solve for the unknown $a$ and $b$.

That's the inductive step, more or less, but how to get started? The original question assumes (justifiably -- by a translation, if necessary) that $(0,0)$ is mapped to $(0,0)$. If the linear terms form an invertible system, then the inductive step can be applied repeatedly up to the desired degree. In the case of the question, the assumptions are that the linear system is nonzero but degenerate ($U_x \ne 0$, $U_y=0$) and the second order system is parabolic, aligned in an independent direction ($U_{yy}(0,0) \ne 0$, $U_{xx}=U_{xy}=0$). In this case, the solutions will be power series in $u$, $\sqrt{v}$.

The function series2DSolveParabolic constructs the series solution up to whatever degree desired. Notes on the code can be found below the example.

terms[poly_, deg_, vars_] := FromCoefficientRules[

   Select[CoefficientRules[poly, vars], Total@First@# <= deg &], vars];

series2D[series2DData[f_, {x_, y_}, dataRules_], n_] := 
  terms[Normal@Series[f[x, y], {x, 0, n}, {y, 0, n}] /. dataRules, 
    n, {x, y}];
series2D[series2DSolData[a_, {u_, v_}], n_] := Sum[u^i*v^(-i + k)*a[i, -i + k],
  {k, 1, n}, {i, 0, k}];

solveCoeff[terms_, pat_] := Solve[Thread[terms == 0], Cases[terms, pat, Infinity]];


series2DSolveParabolic[{u_ == U0_series2DData, v_ == V0_series2DData}, {x_, y_}, deg_] :=
  Module[{a, b, X0, Y0, r},
   a[0, 0] = b[0, 0] = 0;
   X0 = series2DSolData[a, {u, r}];
   Y0 = series2DSolData[b, {u, r}];
   Do[
    solveCoeff[
      DeleteCases[
       Flatten[CoefficientList[#, {u, r}] & /@ {
          u - terms[series2D[U0, n] /.

            {x -> series2D[X0, n], y -> series2D[Y0, n]}, n, {u, r}],
          r^2 - terms[series2D[V0, n + 1] /.
            {x -> series2D[X0, n + 1], y -> series2D[Y0, n + 1]}, n + 1, {u, r}]}],
       _?NumericQ], _a | _b] /. Rule -> Set,
    {n, 1, deg}];

   {x -> Sum[a[i, j - i] u^i v^((j - i)/2), {j, 0, deg}, {i, 0, j}],
    y -> Sum[b[i, j - i] u^i v^((j - i)/2), {j, 0, deg}, {i, 0, j}]}
   ];

Example (a bit messy):

    series2DSolveParabolic[
     {u == series2DData[U, {x, y}, {U[0, 0] -> 0, Derivative[0, 1][U][0, 0] -> 0}],
      v == series2DData[V, {x, y},
        {V[0, 0] -> 0, Derivative[1, 0][V][0, 0] -> 0, 
         Derivative[0, 1][V][0, 0] -> 0, 
         Derivative[1, 1][V][0, 0] -> 0, 
         Derivative[2, 0][V][0, 0] -> 0}]},
        {x, y}, 2]


(* {x -> u/Derivative[1, 0][U][0, 0] -
         (v*Derivative[0, 2][U][0, 0])/
           (Derivative[0, 2][V][0, 0] * Derivative[1, 0][U][0, 0]) -
         (Sqrt[2]*u*Sqrt[v] * Derivative[1, 1][U][0, 0]) /
           (Sqrt[Derivative[0, 2][V][0, 0]] * Derivative[1, 0][U][0, 0]^2) -
         (u^2 * Derivative[2, 0][U][0, 0]) / (2*Derivative[1, 0][U][0, 0]^3), 
   y -> (Sqrt[2]*Sqrt[v]) / Sqrt[Derivative[0, 2][V][0, 0]] -
         (v * Derivative[0, 3][V][0, 0]) / (3*Derivative[0, 2][V][0, 0]^2) -
         (u * Sqrt[v]*Derivative[1, 2][V][0, 0]) / 
           (Sqrt[2]*Derivative[0, 2][V][0, 0]^(3/2) * Derivative[1, 0][U][0, 0]) -

         (u^2*Derivative[2, 1][V][0, 0]) /
           (2 * Derivative[0, 2][V][0, 0]*Derivative[1, 0][U][0, 0]^2)}  *)

$$\left\{x\to -\frac{u^2 U^{(2,0)}(0,0)}{2 U^{(1,0)}(0,0)^3}-\frac{\sqrt{2} u \sqrt{v} U^{(1,1)}(0,0)}{U^{(1,0)}(0,0)^2 \sqrt{V^{(0,2)}(0,0)}}+\frac{u}{U^{(1,0)}(0,0 )}-\frac{v U^{(0,2)}(0,0)}{U^{(1,0)}(0,0) V^{(0,2)}(0,0)}\right.,$$ $$\left.y\to -\frac{u^2 V^{(2,1)}(0,0)}{2 U^{(1,0)}(0,0)^2 V^{(0,2)}(0,0)}-\frac{u \sqrt{v} V^{(1,2)}(0,0)}{\sqrt{2} U^{(1,0)}(0,0) V^{(0,2)}(0,0)^{3/2}}-\frac{v V^{(0,3)}(0,0)}{3 V^{(0,2)}(0,0)^2}+\frac{\sqrt{2} \sqrt{v}}{\sqrt{V^{(0,2)}(0,0)}}\right\}$$

Notes on the code:

terms[poly, deg, vars] -- return the terms of the polynomial whose total degree is at most deg

series2DData[f, {x, y}, dataRules] -- represents a series in the variables x and y, with conditions on the derivatives of f given in dataRules

series2DSolData[c, {x, y},] -- represents a series in the variables x and y with coefficients c[i, j].

series2D[ser, n] -- returns the n-th degree Taylor polynomial of the series ser; ser may be either series2DData or series2DSolData

solveCoeff[terms, pat] -- sets the terms equal to zero and solves for the coefficients represented by the pattern pat.

series2DSolveParabolic[{u == uSeries, v == vSeries}, vars, deg] solves the system of equations for power series of vars in terms of the variables {u, v} up to degree deg. Note that the equation are not really treated as equations. The variables u and v have to be symbols and be the left-hand sides. Note also that internally in the Do loop, r represents Sqrt[v], and that the exponent is off by one (goes up to n+1); this has to do with the recursive procedure starting with degree 1 in x and degree 2 in y.

---

Here's a function that will invert power series if the linear terms are nonsingular.

series2DSolve[{u_ == U0_series2DData, v_ == V0_series2DData}, {x_, y_}, deg_] :=
  Module[{a, b, X0, Y0},
   a[0, 0] = b[0, 0] = 0;
   X0 = series2DSolData[a, {u, v}];
   Y0 = series2DSolData[b, {u, v}];
   Do[
    solveCoeff[

      DeleteCases[
       Flatten[CoefficientList[#, {u, v}] & /@ {
          u -> terms[series2D[U0, n] /. {x -> series2D[X0, n], y -> series2D[Y0, n]},
                 n, {u, v}],
          v -> terms[series2D[V0, n] /. {x -> series2D[X0, n], y -> series2D[Y0, n]},
                 n, {u, v}]}],
        _?NumericQ], _a | _b] /. 
     Rule -> Set,
    {n, 1, deg}];


   {x -> Sum[a[i, j - i] u^i v^(j - i), {j, 0, deg}, {i, 0, j}],
    y -> Sum[b[i, j - i] u^i v^(j - i), {j, 0, deg}, {i, 0, j}]}
   ];

Example

f[x_, y_] := Log[1 + x];
g[x_, y_] := Sin[y];
series2DSolve[
 {u == series2DData[f, {x, y}, {}],
  v == series2DData[g, {x, y}, {}]},

 {x, y}, 5]
  (* {x -> u + u^2/2 + u^3/6 + u^4/24 + u^5/120, y -> v + v^3/6 + (3 v^5)/40} *)

Normal@Series[ArcSin[v], {v, 0, 6}]
  (* v + v^3/6 + (3 v^5)/40 *)

---

The methods used are from basic algebra. There could be, there ought to be, or perhaps there are built-in functions to do some or all of the above.