# How to evaluate this InverseFunction?

Solving an ODE problem

DSolve[{y'[x] == x*(y[x]^2 - 1)^(2/3), y[1] == 3/2}, y[x], x]


{{y[x] -> InverseFunction[( Hypergeometric2F1[1/2, 2/3, 3/2, #1^2] #1 (1 - #1^2)^(2/3))/(-1 + #1^2)^(2/3) &][x^2/2 + 1/2 (-1 + 3 (-1)^(2/3) Hypergeometric2F1[1/2, 2/3, 3/2, 9/4])]}}

, I come to an inverse function and face problems with handling it.

Plot[InverseFunction[(Hypergeometric2F1[1/2, 2/3, 3/2, #1^2] #1 (1 - #1^2)^(2/3))/(-1 + #1^2)^(2/3) &][1/2 (-1 + x^2 +
3 (-1)^(2/3) Hypergeometric2F1[1/2, 2/3, 3/2, 9/4])], {x, 1, 2}]


returns an empty plot. Looking in the documentation/Scope, I try

f = InverseFunction[(Hypergeometric2F1[1/2, 2/3, 3/2, #1^2] #1 (1 - #1^2)^(2/3))/(-1 + #1^2)^(2/3) &]
[ 1/2 (-1 + x^2 +  3 (-1)^(2/3) Hypergeometric2F1[1/2, 2/3, 3/2, 9/4])];
f[1.1]


InverseFunction[( Hypergeometric2F1[1/2, 2/3, 3/2, #1^2] #1 (1 - #1^2)^( 2/3))/(-1 + #1^2)^(2/3) &][ 1/2 (-1 + x^2 + 3 (-1)^(2/3) Hypergeometric2F1[1/2, 2/3, 3/2, 9/4])][1.1]

So how to evaluate this?

• tbl=Table[{x,InverseFunction[...]}, {x,1,3/2,1/16}] followed by ListPlot[tbl,Joined->True] and InputForm[tbl] (and gently fiddling with some of those values) might give you a tiny bit of insight. It looks like it is having a lot of difficulty with some values of x and the large number of different values of x that Plot uses may be trying some of those values.
– Bill
Commented Nov 4, 2022 at 7:06
• @Bill: Thank you. Something strange: y[1]==3/2, but tb1[[1]] is near $-1$. May we trust it? Commented Nov 4, 2022 at 7:19
• Bluntly honest: I do not know. Your DE is complicated enough that after I stared at it for a while I concluded I had no idea what it was or how it behaved. And the solution is much worse. The latest version of MMA even refuses to give an answer, claiming that some parts of this lead to no solution. At first I wondered whether complex values might be part of the reason that Plot was blank, but I'm still not sure. I hope my fumbling might let you see something about this. Good luck with it.
– Bill
Commented Nov 4, 2022 at 7:28
• @Bill: sol = NDSolve[{y'[x] == x*(y[x]^2 - 1)^(2/3), y[1] == 3/2}, y[x], {x, 1, 2}]; seems to be OK. Commented Nov 4, 2022 at 7:35
• When I separate variables and integrate, I obtain a different inverse function expression: $$\frac{y \left(1-y^2\right)^{2/3} \, _2F_1\left(\frac{1}{2},\frac{2}{3};\frac{3}{2};y^2\right)}{\left(y^2-1\right)^{2/3}}=\frac{\text{y0} \left(1-\text{y0}^2\right)^{2/3} \, _2F_1\left(\frac{1}{2},\frac{2}{3};\frac{3}{2};\text{y0}^2\right)}{\left(\text{y0}^2-1\right)^{2/3}}+\frac{x^2}{2}-\frac{1}{2}$$ And this expression, when solved iteratively for y using FindRoot, over x in (1,2) (minus a small residual imaginary part) agrees with the numerical solution via NDSolve.
– josh
Commented Nov 4, 2022 at 15:01

Given:

$$\frac{dy}{dx}=x(y^2-1)^{2/3}, \quad y(x_0)=y_0$$ and separating variables and integrating: $$\int_{3/2}^y \frac{dy}{(y^2-1)^{2/3})}=\int_{1}^x x dx$$ we obtain

$$\left(\frac{y \left(1-y^2\right)^{2/3} \, _2F_1\left(\frac{1}{2},\frac{2}{3};\frac{3}{2};y^2\right)}{\left(y^2-1\right)^{2/3}}\right)\biggr|_{y_0}^y=\frac{x^2}{2}\biggr|_{x_0}^x$$

or:

$$\frac{y \left(1-y^2\right)^{2/3} \, _2F_1\left(\frac{1}{2},\frac{2}{3};\frac{3}{2};y^2\right)}{\left(y^2-1\right)^{2/3}}=\frac{x^2}{2}-\frac{x_0^2}{2}+\frac{y_0 \left(1-y_0^2\right)^{2/3} \, _2F_1\left(\frac{1}{2},\frac{2}{3};\frac{3}{2};y_0^2\right)}{\left(y_0^2-1\right)^{2/3}}=g(x;y_0)$$

which can be written in the Mathematica construct:

$$\text{InverseFunction}\left[\frac{y \left(1-y^2\right)^{2/3} \, _2F_1\left(\frac{1}{2},\frac{2}{3};\frac{3}{2};y^2\right)}{\left(y^2-1\right)^{2/3}}\right]\left[g(x,y_0)\right]$$

In order to evaluate the inverse, first obtain the antiderivative at $$y_0=3/2$$:

intF[y_] = Integrate[1/(y^2 - 1)^(2/3), y];
(* compute integral at y=3/2  *)
intVal = intF[3/2] // N


then use FindRoot over the range $$x\in[1,2]$$

currentY = 3/2;
myPts = Table[
nextY =
Re[y /. FindRoot[(
y (1 - y^2)^(2/3)
Hypergeometric2F1[1/2, 2/3, 3/2, y^2])/(-1 + y^2)^(2/3) ==
x^2/2 - 1/2 + intVal /. x -> xVal, {y, currentY}]];
currentY = nextY;
{xVal, nextY},
{xVal, 1.01, 2, 0.01}
];
p1 = ListPlot[myPts, Joined -> True, PlotStyle -> {Dashed, Red}]


There are likely more efficient ways to do this.

Note I am extracting the real component of the solution as there is a small residual imaginary component. And one way to compare the inverse funtion results with NDSolve is to simply superimpose the inverse function solution over the NDSolve solution:

sol = NDSolveValue[{y'[x] == x*(y[x]^2 - 1)^(2/3), y[1] == 3/2},
y, {x, 1, 2}];
p2 = Plot[sol[x], {x, 1, 2}, PlotStyle -> Black]
Show[{p2, p1},
PlotLabel ->
Style["Superposition of inverse function solution with NDSolve \
solution"]]