Mario, the algorithm looks quite similar to the Newton method you have used for functional inverses for expanding tan and tanh. We could also expand `sin/cos`. Which one do think should be faster?
Also, can the convergence of the method be proven? I could not find any reference to it (except the algorithm in Computer Arithmetic by Paul Zimmerman). Shivam Vats On Wed, Jun 10, 2015 at 9:00 PM, mario <[email protected]> wrote: > > 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]> 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/CAGd-hGGYuApMHGcp1QXi%2B0gE1CJEQwdd%2B0qe_%2BkM7iX03tCP3A%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
