There is `series_reversion` in PR609. On Wednesday, June 10, 2015 at 6:39:02 AM UTC+2, Ondřej Čertík wrote: > > On Tue, Jun 9, 2015 at 7:33 PM, Aaron Meurer <[email protected] > <javascript:>> wrote: > > On Tue, Jun 9, 2015 at 7:09 PM, Ondřej Čertík <[email protected] > <javascript:>> 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] <javascript:>. > >> To post to this group, send email to [email protected] > <javascript:>. > >> 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] <javascript:>. > > To post to this group, send email to [email protected] > <javascript:>. > > 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/9091e3df-f704-47a7-b366-e00c355d7cba%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
