On Tue, Jun 9, 2015 at 7:09 PM, Ondřej Čertík <[email protected]> wrote: > Hi, > > I would like to evaluate to higher accuracy (quadruple precision is > enough) an inverse Fermi-Dirac integral of order 1/2. The direct > function is the following integral (and it can also be written as a > polylogarithm): > > I_{1/2}(x) = Integral(sqrt(t) / (1+exp(t-x)), (t, 0, oo)) > = -gamma(S(3)/2) * polylog(S(3)/2, -exp(x)) > > and I need its inverse. SymPy can evaluate either representation: > > In [19]: (-gamma(S(3)/2) * polylog(S(3)/2, -exp(x))).subs(x, 3).n(35, > chop=True) > Out[19]: 3.9769853540479774178558497377805661 > > In [20]: Integral(sqrt(t) / (1+exp(t-x)), (t, 0, oo)).subs(x, 3).n(35) > Out[20]: 3.9769853540479774178558497377805661 > > I need to calculate an inverse, i.e. for 3.9769... it would return > 3.00..., to arbitrary accuracy. Once I have that, I want to use > rational function approximation to obtain very fast and accurate > double precision implementation. > > What is the best way to do that? > > There is a recent publication [1], which does that using Mathematica, > and then provides Fortran implementation. I tested it and it works > great. What I actually need is to find a rational approximation to an > expression that contains both the inverse and direct Fermi-Dirac > integrals and I want to just have a rational approximation for the > final expression. The inverse Fermi-Dirac is the hardest part, so > that's why I am asking. > > In [1], they calculate a series expansion of the Fermi-Dirac integral > and then they reverse (invert) the series, see eq. (9), (10), (14), > (15) and the Mathematica code afterwards. Unfortunately, SymPy raises > an exception for a series of the polylog function: > > In [23]: polylog(S(3)/2, -z).series(z, 0, 5) > > https://github.com/sympy/sympy/issues/9497 > > But let's say we fix it. Here i s a slightly modified expression that works: > > In [6]: polylog(S(3)/2, 1+z).series(z, 0, 5) > zeta(3/2) + z*zeta(1/2) + z**2*(zeta(-1/2)/2 - zeta(1/2)/2) + > z**3*(zeta(1/2)/3 + zeta(-3/2)/6 - zeta(-1/2)/2) + > z**4*(11*zeta(-1/2)/24 + zeta(-5/2)/24 - zeta(-3/2)/4 - zeta(1/2)/4) + > O(z**5) > > How can we invert the series? Mathematica has a function InverseSeries: > > http://reference.wolfram.com/language/ref/InverseSeries.html > > It seems SymPy can't do it yet. > > Shivam, I think that would be a very useful addition to the series > module that you are implementing. We have to figure out an algorithm > for it.
I think there's a closed form for the inversion of a formal power series. See https://en.wikipedia.org/wiki/Formal_power_series#The_Lagrange_inversion_formula. I'm not sure what the algorithm is for series like x*sin(x), which that Mathematica example shows is a series in sqrt(x). I suppose that's because the linear term is 0 (the Wikipedia article claims that the constant needs to be 0 and the linear term needs to be nonzero). Strictly speaking x*sin(x) is not invertible near x=0 (it's an even function). I guess this picks one branch and inverts that. Aaron Meurer > > Ondrej > > > [1] Fukushima, T. (2015). Precise and fast computation of inverse > Fermi–Dirac integral of order 1/2 by minimax rational function > approximation. Applied Mathematics and Computation, 259, 698–707. > doi:10.1016/j.amc.2015.03.015 > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. > To post to this group, send email to [email protected]. > Visit this group at http://groups.google.com/group/sympy. > To view this discussion on the web visit > https://groups.google.com/d/msgid/sympy/CADDwiVCp7rfTUk0TrooNaWqpVrADFbGPb_Ptup2Jn49twTapEg%40mail.gmail.com. > For more options, visit https://groups.google.com/d/optout. -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6K%2BXXF6-Bb7tjxXHMbt74wEG_ufEu_%3D48-2UnQ4YbJjFw%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
