Sure, Ondrej. I will send a PR for it. On Wed, Jun 10, 2015 at 9:23 PM, Shivam Vats <[email protected]> wrote:
> 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-hGF7-gGGfNg90R2oSyE78T%3DMnsr2CGe5SvSMtZ8hwprc_A%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
