On Tue, Jun 9, 2015 at 7:33 PM, Aaron Meurer <[email protected]> wrote: > 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.
I don't quite get it with the constant term, but it seems to be true, e.g. the inversion of exp(x) = 1 + x + .... doesn't exist, because that would be log(x), which doesn't have an expansion around x=0. One can only expand e.g. log(1+x), as the inverse is exp(x)-1. I found a pretty good intro here: http://mathworld.wolfram.com/SeriesReversion.html with further references to literature. > > 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. -- 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/CADDwiVBkeVs3TmYxH5MVXYXurkLxUnwBqraX%3D41oA-tBKTBq2Q%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
