Determine inverse function with Gröbner basis

Question

I am trying to generalize a way to invert a function using Gröbner basis. My functions are defined in the following form:

bound[K10_] := 
  xi - (p0*(Q10 - K10)/(p0 - K10 + 
         Sqrt[(K10)^2 - m1^2]*(1 - 2*alpha)))/(u0Max);
p0 = 1;
m1 = 4/10;
(*m2=4/10;*)
(*Q10=1/2*(p0+(m1-m2)*(m1+m2)/p0);*)
Q10 = 1/2;
Q10Min = (m1^2 + p0^2/4)/p0;
u0F[K1_, alpha_] := 
  p0*(Q10 - K1)/(p0 - K1 + Sqrt[K1^2 - m1^2] (1 - 2 alpha));
(*If[m1 \[Equal] m2,
u0Max = p0*(Q10 - Q10Min)/(p0 - Q10Min - Sqrt[Q10Min^2 - m1^2]),
u0Max = u0F[K10P //.alpha\[Rule]1,1]];*)

u0Max = p0*(Q10 - Q10Min)/(p0 - Q10Min - Sqrt[Q10Min^2 - m1^2]);
K10P = -m1^2*
    Q10*(1 - 2*alpha)^2/(p0^2 - 2*p0*Q10 + 
       Q10^2*(1 - (1 - 2*alpha)^2)) + 
   Sqrt[(-m1^2*
        Q10*(1 - 2*alpha)^2/(p0^2 - 2*p0*Q10 + 
           Q10^2*(1 - (1 - 2*alpha)^2)))^2 + (m1^4*(1 - 2*alpha)^2 + 
        m1^2*(p0 - Q10)^2)/(p0^2 - 2*p0*Q10 + 
        Q10^2*(1 - (1 - 2*alpha)^2))];

gb = GroebnerBasis[bound[K10P], {alpha, xi}];
alphaxi = First[gb]

As you can see in u0F[---], i have got a function which depends on alpha and i would like to invert this function such that i obtain an expression alpha[u0F]. The implicit function bound[K10_] is the starting point relating xi and alpha and i would like to solve this for alpha. Additionally, bound[K10_] depends on arbitrary constants p0,m1,Q10, which i chose in my evaluation. Decomposing this with the Gröbner basis, one obtains

-1562500 - 3515625 alpha + 3515625 alpha^2 + 4250000 xi + 
 9562500 alpha xi - 9562500 alpha^2 xi - 3452500 xi^2 - 
 10563750 alpha xi^2 + 12183750 alpha^2 xi^2 - 3240000 alpha^3 xi^2 + 
 1620000 alpha^4 xi^2 + 765000 xi^3 + 5523300 alpha xi^3 - 
 7726500 alpha^2 xi^3 + 4406400 alpha^3 xi^3 - 2203200 alpha^4 xi^3 - 
 1136025 alpha xi^4 + 2154681 alpha^2 xi^4 - 2223936 alpha^3 xi^4 + 
 1578528 alpha^4 xi^4 - 559872 alpha^5 xi^4 + 186624 alpha^6 xi^4

which can be transformed into an expression for alpha:

{{alpha -> (6 xi^2 - Sqrt[-625 xi^2 + 850 xi^3 - 189 xi^4])/(
   12 xi^2)}, {alpha -> (
   6 xi^2 + Sqrt[-625 xi^2 + 850 xi^3 - 189 xi^4])/(
   12 xi^2)}, {alpha -> 
   1/2 (1 - 
      Sqrt[-489 - 625/xi^2 + 850/xi - Sqrt[
       390625 - 1062500 xi + 1583750 xi^2 - 1171300 xi^3 + 
        314721 xi^4]/xi^2]/(6 Sqrt[2]))}, {alpha -> 
   1/2 (1 + 
      Sqrt[-489 - 625/xi^2 + 850/xi - Sqrt[
       390625 - 1062500 xi + 1583750 xi^2 - 1171300 xi^3 + 
        314721 xi^4]/xi^2]/(6 Sqrt[2]))}, {alpha -> 
   1/2 (1 - Sqrt[
      1 + (-625 + 850 xi - 561 xi^2 + Sqrt[
        390625 - 1062500 xi + 1583750 xi^2 - 1171300 xi^3 + 
         314721 xi^4])/(72 xi^2)])}, {alpha -> 
   1/2 (1 + Sqrt[
      1 + (-625 + 850 xi - 561 xi^2 + Sqrt[
        390625 - 1062500 xi + 1583750 xi^2 - 1171300 xi^3 + 
         314721 xi^4])/(72 xi^2)])}}

This is my desired solution for a chosen set of p0, m1 and Q10. I would like to generalize this problem my choosing an arbitrary m1 while keeping p0 = 1 and Q10 = 1/2, so i would like to end up with a function alpha depending on xi and m1. Is it possible to generalize this method?

Answer 1

Sure. If you clear the value for m1 you can treat it as a parameter and redo the GB computation as such:

Timing[gb = GroebnerBasis[bound[K10P], {alpha, xi}, CoefficientDomain -> RationalFunctions];]

This took around 11 minutes on my machine. The first element will give a polynomial in {alpha,xi} with coefficients that are algebraic functions in m1. If you require a polynomial in all variables, that will take more work.

In[216]:= Timing[gbB = GroebnerBasis[gb[[1]], {alpha, xi, m1}];]

(* Out[216]= {1.412, Null} *)

See which elements use all three variables.

In[223]:= Map[Variables, gbB]

(* Out[223]= {{m1}, {m1}, {alpha, m1, xi}, {alpha, m1, xi}, {alpha, m1, xi}} *)

We'll work with the third one. I poked around and saw that it involves a radical of the form sqrt(4 m^2 - 1)^2) so I use PowerExpand to remove the radical.

bigPoly = PowerExpand[Numerator[Together[gbB[[3]]]]];
PolynomialQ[bigPoly, {alpha, xi, m1}]
LeafCount[bigPoly]

(* Out[221]= True

Out[222]= 16804 *)

I'm not convinced the result will be of much use though. It might be better to invert for specific values as needed. This will almost certainly be the case if a goal is subsequent numerical evaluation, since the general formulas will likely not be stable numerically, given the coefficient sizes and length of expressions.