As long as the series to be inverted has finite precision `n`, has no constant term and finite linear term, one can revert the series iteratively; I do not think there are convergence problems, since these series are formal, or anyway for sufficienly small region around 0, the function is monotonic, so it can be reversed.
On Wednesday, June 10, 2015 at 5:53:16 PM UTC+2, Shivam Vats 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] <javascript:>> > 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/e699d112-7697-4b93-8088-591398ef3ae5%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
