0
$\begingroup$

I consider the following function:

bound[K10_] := xi - (p0*(Q10 - K10)/(p0 - K10 + Sqrt[(K10)^2 - m1^2]*(1 - 2*alpha)))/(u0Max);

This function is basically an implicit function which relates the variable xi to alpha. p0,Q10,m1 and u0Max are constants, which i set to

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

K10P is defined in the following way:

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))];

The implicit function "bound[K10_]" defines a function of the form xi[alpha], but i would like to solve this implicit function for alpha and obtain an expression depending on xi only, so i try:

Solve[bound[K10P] == 0, alpha];

Unfortunately, this does not work in the sense that this command runs forever without solving the equation. Do you have any suggestions how to obtain my desired expression?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
7
$\begingroup$

Could turn it into a polynomial relation using GroebnerBasis.

bound[K10_] := 
  xi - (p0*(Q10 - K10)/(p0 - K10 + 
         Sqrt[(K10)^2 - m1^2]*(1 - 2*alpha)))/(u0Max);
p0 = 1;
m1 = 1/10;
m2 = 1/10;
Q10 = 1/2*(p0 + (m1 - m2)*(m1 + m2)/p0);
Q10Min = (m1^2 + p0^2/4)/p0;
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]

(* Out[43]= -390625 - 37500000 alpha + 37500000 alpha^2 + 1531250 xi + 
 147000000 alpha xi - 147000000 alpha^2 xi - 1875625 xi^2 - 
 217530000 alpha xi^2 + 220410000 alpha^2 xi^2 - 
 5760000 alpha^3 xi^2 + 2880000 alpha^4 xi^2 + 735000 xi^3 + 
 144001200 alpha xi^3 - 149646000 alpha^2 xi^3 + 
 11289600 alpha^3 xi^3 - 5644800 alpha^4 xi^3 - 35985600 alpha xi^4 + 
 38807424 alpha^2 xi^4 - 5698944 alpha^3 xi^4 + 
 2987712 alpha^4 xi^4 - 165888 alpha^5 xi^4 + 55296 alpha^6 xi^4 *)

It turns out that this factors into a quadratic and a quartic in alpha, so if so desired one can get radical solution branches.

Factor[alphaxi]

(* Out[47]= (625 - 1225 xi + 600 xi^2 - 24 alpha xi^2 + 
   24 alpha^2 xi^2) (-625 - 60000 alpha + 60000 alpha^2 + 1225 xi + 
   117600 alpha xi - 117600 alpha^2 xi - 59976 alpha xi^2 + 
   62280 alpha^2 xi^2 - 4608 alpha^3 xi^2 + 2304 alpha^4 xi^2) *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ This works well, thanks for the suggestion. Is it also possible to extend this method such that i can use arbitrary values for p0,m1 and m2? $\endgroup$
    – RealDestructor
    Commented Jul 4, 2019 at 13:42
  • $\begingroup$ You can try that. I think it should work in the same way. $\endgroup$
    – Daniel Lichtblau
    Commented Jul 4, 2019 at 15:32
  • $\begingroup$ It works for sure if i replace the values for p0,m1 or m2 by different values, but i am looking for a final expression which depends on p0,m1 and m2 and can be chosen arbitrarily. If i try this, it does not work. $\endgroup$
    – RealDestructor
    Commented Jul 4, 2019 at 15:45
2
$\begingroup$

First, let's build the function you want to define. We use your equation bound[K10P] == 0 and more general Re[bound[K10P]] == 0

Q10 = 1/2*(p0 + (m1 - m2)*(m1 + m2)/p0);
Q10Min = (m1^2 + p0^2/4)/p0;
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))];
ContourPlot[bound[K10P] == 0, {alpha, -1, 2}, {xi, 0, 1.2}, 
 FrameLabel -> Automatic, PlotPoints -> 150]
ContourPlot[Re[bound[K10P]] == 0, {alpha, -5, 5}, {xi, 0, 1.2}, 
 FrameLabel -> Automatic, PlotPoints -> 150]

Figure 1

We see that the alpha function has several branches for real values of xi. This is the problem. The same can be obtained using

bound[K10P] // FullSimplify

Out[]= -((25 (96 - Sqrt[(
        1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 1/(
        alpha - alpha^2)))/(12 (196 + 1/(1 - alpha) + 1/alpha - Sqrt[(
        1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 
        Sqrt[2] (1 - 
           2 alpha) Sqrt[((1 - 
            2 alpha)^2 (1 + (-1 + alpha) alpha (-48 + Sqrt[(
               1 - 96 (-1 + alpha) alpha)/((-1 + 
                  alpha)^2 alpha^2)])))/((-1 + alpha)^2 alpha^2)]))) +
  xi

From here we immediately find (xif=xi)

xif[alpha_] := (
 25 (96 - Sqrt[(
    1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 1/(
    alpha - alpha^2)))/(
 12 (196 + 1/(1 - alpha) + 1/alpha - Sqrt[(
    1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 
    Sqrt[2] (1 - 
       2 alpha) Sqrt[((1 - 
        2 alpha)^2 (1 + (-1 + alpha) alpha (-48 + Sqrt[(
           1 - 96 (-1 + alpha) alpha)/((-1 + 
              alpha)^2 alpha^2)])))/((-1 + alpha)^2 alpha^2)]))

Get the same curves

{Plot[xif[x], {x, 0, 1.2}, PlotPoints -> 200, 
  AxesLabel -> {"alpha", "xi"}], 
 Plot[Re[xif[x]], {x, -5, 5}, PlotPoints -> 200, 
  AxesLabel -> {"alpha", "xi"}]}

Figure2

Define the inverse function

\[Alpha] = 
 InverseFunction[
  Function[{alpha}, (
   25 (96 - Sqrt[(
      1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 1/(
      alpha - alpha^2)))/(
   12 (196 + 1/(1 - alpha) + 1/alpha - Sqrt[(
      1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 
      Sqrt[2] (1 - 
         2 alpha) Sqrt[((1 - 
          2 alpha)^2 (1 + (-1 + alpha) alpha (-48 + Sqrt[(
             1 - 96 (-1 + alpha) alpha)/((-1 + 
                alpha)^2 alpha^2)])))/((-1 + alpha)^2 alpha^2)]))]]

If you evaluate the alpha function, then there will be no answer, since the system will not be able to select a branch. But you can choose the branch of the solution of the equation. Suppose we chose branch 0<alpha<1. Define the function

xif[alpha_ /; alpha > 0 && alpha < 1] := (
 25 (96 - Sqrt[(
    1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 1/(
    alpha - alpha^2)))/(
 12 (196 + 1/(1 - alpha) + 1/alpha - Sqrt[(
    1 - 96 (-1 + alpha) alpha)/((-1 + alpha)^2 alpha^2)] + 
    Sqrt[2] (1 - 
       2 alpha) Sqrt[((1 - 
        2 alpha)^2 (1 + (-1 + alpha) alpha (-48 + Sqrt[(
           1 - 96 (-1 + alpha) alpha)/((-1 + 
              alpha)^2 alpha^2)])))/((-1 + alpha)^2 alpha^2)]))

Make sure that we select the desired branch:

{Plot[xif[x], {x, 0, 1.2}, PlotPoints -> 200, 
  AxesLabel -> {"alpha", "xi"}, PlotRange -> All], 
 Plot[Re[xif[x]], {x, -5, 5}, PlotPoints -> 200, 
  AxesLabel -> {"alpha", "xi"}, PlotRange -> {0, 1}]}

Figure 3

Define and plot the inverse function

alph = InverseFunction[xif]

Plot[alph[x], {x, .3, 0.99}, PlotRange -> All, 
 AxesLabel -> {"xi", "alpha"}]

Figure 4

Improve this answer
$\endgroup$
3
  • $\begingroup$ How can i choose the branch in particular? Lets suppose that alpha is between 0 and 1 as you already plotted above. $\endgroup$
    – RealDestructor
    Commented Jul 3, 2019 at 20:23
  • $\begingroup$ @RealDestructor See update to my answer . $\endgroup$
    – Alex Trounev
    Commented Jul 4, 2019 at 3:40
  • $\begingroup$ This works perfectly fine, thanks for that. Is it also possible to obtain an explicit analytical expression for the inverse function? $\endgroup$
    – RealDestructor
    Commented Jul 4, 2019 at 8:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.