Is there any chance of having a formula or approximation to inverse the Modified Bessel function of the first kind? I mean to solve $I_M(x)$ with respect to x: $I^{-1}_M(x)$?
2 Answers
I've come across a more precise method, in Mathematica you can use
InverseSeries[Series[BesselI[0,x],{x,0,20}]]
which gives \begin{equation} 2\sqrt{x-1}-\frac{1}{4}(x-1)^{3/2}+\frac{47}{576}(x-1)^{5/2}-\frac{161}{4608}(x-1)^{7/2}+\frac{565571}{33177600}(x-1)^{9/2}-\cdots \end{equation}

- 4,321
The behaviour of $I_0(x)$ closely follows $\cosh(x)$, \begin{equation} I_0(x)=1+\frac{x^2}{4}+\frac{x^4}{64}+\frac{x^6}{2304}+\cdots,\\ \cosh(x)=1+\frac{x^2}{2}+\frac{x^4}{24}+\frac{x^6}{720}+\cdots, \end{equation}
so you could take inspiration from an expansion of $\mathrm{acosh}(x)$ about $x=1$, \begin{equation} \mathrm{acosh}(x)=\sqrt{2} \sqrt{x-1}-\frac{(x-1)^{3/2}}{6 \sqrt{2}}+\frac{3 (x-1)^{5/2}}{80\sqrt{2}}-\frac{5(x-1)^{7/2}}{448 \sqrt{2}}+\cdots, \end{equation} and find some coefficients $a_n$ \begin{equation} I^{-1}_0(x)=\sum_{n=1}^\infty a_n(x-1)^{(2n-1)/2} \end{equation} By fitting some values from an appropriate root finder (Wolfram|Alpha will do) with a truncated series the first few terms appear to be \begin{equation} a_1=2\\ a_2=-\frac{1}{4} \end{equation} Below are the results of a parameter fitting from $\mathrm{gnuplot}$. \begin{equation} \begin{matrix} a_1 & 2 & \pm 6.13\cdot10^{-8} \\ a_2 & -0.249991 & \pm 1.3\cdot10^{-6} \\ a_3 & 0.0815024 & \pm 8.409\cdot10^{-6} \\ a_4 & -0.0345023 & \pm 2.578\cdot10^{-5} \\ a_5 & 0.0159056 & \pm 4.474\cdot10^{-5} \\ a_6 & -0.00706415 & \pm 4.802\cdot10^{-5} \\ a_7 & 0.00270988 & \pm3.337\cdot10^{-5} \\ a_8 & -0.000821109 & \pm 1.528\cdot10^{-5} \\ a_9 & 0.000182402 & \pm 4.578\cdot10^{-6} \\ a_{10} & -2.74386\cdot10^{-5} & \pm 8.628\cdot10^{-7} \\ a_{11} & 2.47422\cdot10^{-6} & \pm 9.278\cdot10^{-8} \\ a_{12} & -1.00428\cdot10^{-7} & \pm 4.339\cdot10^{-9} \end{matrix} \end{equation} we can set these first two values and then fit again \begin{equation} \begin{matrix} a_3 & 0.0814389 & \pm 7.448\cdot10^{-6} \\ a_4 & -0.0340809 & \pm 2.745\cdot10^{-5} \\ a_5 & 0.0150122 & \pm 4.322\cdot10^{-5} \\ a_6 & -0.00607734 & \pm 3.866\cdot10^{-5} \\ a_7 & 0.00205892 & \pm 2.205\cdot10^{-5} \\ a_8 & -0.000552547 & \pm 8.5\cdot10^{-6} \\ a_9 & 0.000113841 & \pm 2.285\cdot10^{-6} \\ a_{10} & -1.76097\cdot10^{-5} & \pm 4.334\cdot10^{-7} \\ a_{11} & 1.99748\cdot10^{-6} & \pm 5.785\cdot10^{-8} \\ a_{12} & -1.6062\cdot10^{-7} & \pm 5.316\cdot10^{-9} \\ a_{13} & 8.64948\cdot10^{-9} & \pm 3.203\cdot10^{-10} \\ a_{14} & -2.79347\cdot10^{-10} & \pm 1.139\cdot10^{-11} \\ a_{15} & 4.08662\cdot10^{-12} & \pm 1.812\cdot10^{-13} \end{matrix} \end{equation}
This expansion of $I^{-1}_0(x)$ seems worst case accurate up to $\varepsilon=\pm 1.2\cdot10^{-6}$ in the interval $[0,10]$, (see residues plot).
Residues of the series fitting against sampled values of $I^{-1}_0(x)$

- 4,321