Dear Dima, I'm trying to offer a tool for "engineering" problems, similar to (and along the lines of) our taylor() and series() functions for symbolic expressions. Therefore, few assumptions can be made. Accordingly, we can accept some failures (as we accept failures of taylor() or series() as in (for example):
sage: sqrt(x).series(x, 5) 1 + Order(x^5) sage: sqrt(x).taylor(x, 0, 5) sqrt(x) (Note : the second one is more informative than the first...). Since anyway rational fractions approximations need some serious thinking before use (e. g. introducing poles when the original function is smooth...), this tool is more of an help than a black-box solution... Le mercredi 13 novembre 2019 13:57:57 UTC+1, Dima Pasechnik a écrit : > > In the multivariate case a lot depends on input. > E.g., do you know something about zeros of your function? > Not in general > E.g. do you have derivatives easily available? > I'm postulating that either taylor() or series() can be used to get them. > If derivates are hard, you probably would like to avoid them all together, > using something known as Newton-Pade > approximation: https://link.springer.com/article/10.1007/BF01390129 > Thanks for this reference. I think that this is an useful tool for getting numerical approximations in caes where symbolic solutions (whch I'm aiming for) aren't available... I'll keep them in mind for further work. > > By the way, have you looked at the survey > https://www.sciencedirect.com/science/article/pii/S037704270000337X > <https://www.google.com/url?q=https%3A%2F%2Fwww.sciencedirect.com%2Fscience%2Farticle%2Fpii%2FS037704270000337X&sa=D&sntz=1&usg=AFQjCNH19JTDvc5Wg3RMaf0k3DHR1rf-hw> > > ? > That's the paper I quoted in my first mail... ;-). Thanks a lot ! > > Dima > > On Wed, 13 Nov 2019, 11:16 Emmanuel Charpentier, <emanuel.c...@gmail.com > <javascript:>> wrote: > >> Le mardi 12 novembre 2019 20:42:16 UTC+1, rjf a écrit : >>> >>> Since Maxima is free and open source and gpl, why not just read the >>> algorithm implemented there >>> and rewrite it in Python? >>> >> >> That can be done. But I had other interest in mind: >> >> - multivariate case (my solution of iterative univariate >> approximations may not be optimal...) >> - other (newer) algorithms, possibly more efficient. >> >> The hint given by Vincent is interesting, and I'm following it. I'll keep >> this thread posted. >> Anyway, wrapping the relevant Maxima calls in a Python function is >> trivial, and may serve as a first implementation. >> >> >>> >>> RJF >>> \ >>> >>> On Monday, November 11, 2019 at 1:29:56 AM UTC-8, Emmanuel Charpentier >>> wrote: >>>> >>>> Dear Vincent, >>>> >>>> a very quick answer (limbic system level. :-). Thank you for this hint. >>>> It seems really interesting. But I'll need time to explore it and find my >>>> way. >>>> >>>> Of course, I will have to work on PolynomialRing(SR) in order to be >>>> able to work on one variable while ignoring the rest... >>>> >>>> I'll keep you posted. >>>> >>>> Le lundi 11 novembre 2019 09:28:58 UTC+1, Vincent Neiger a écrit : >>>>> >>>>> Dear Emmanuel, >>>>> >>>>> You may be interested in taking a look at the following function: >>>>> Matrix_polynomial_dense.minimal_approximant_basis >>>>> >>>>> This only supports the univariate case. This solves a problem which >>>>> generalizes Padé approximation (the documentation gives a precise >>>>> description of what it computes; taking notation from there, instead of >>>>> thinking in terms of polynomials you may view the matrix F as having >>>>> power >>>>> series entries, with the column j truncated at order d_j ). >>>>> >>>>> Coming back to our problem: if you have a power series S(x) at some >>>>> order d and you want to find a Padé approximation f(x) / g(x) of it, you >>>>> can find it by calling the above function on the 2 x 1 polynomial matrix >>>>> [[S(x).polynomial()], [-1]] with the input order being d. Specifically on >>>>> this input the above function will return a 2 x 2 matrix of univariate >>>>> polynomials, and its first row will be [[g(x) , f(x)]], your sought >>>>> approximant. >>>>> >>>>> By default you will get the approximant with f and g of degree about >>>>> d/2 and deg(g) > deg(f) (well, at least on typical instances), but by >>>>> manipulating the optional argument "shift" you can get any other type >>>>> that >>>>> you want. This shift is basically equivalent to the notion of "defects" >>>>> often found in the literature on approximants. >>>>> >>>>> This is based on an algorithm described by Van Barel and Bultheel and >>>>> by Beckermann and Labahn in the early 1990s. A much faster algorithm >>>>> exists, using a divide and conquer approach and fast polynomial >>>>> multiplication, but for it to be interesting we would need Sage to have a >>>>> faster polynomial matrix multiplication (for the moment this faster >>>>> algorithm is not really interesting except for small matrix dimensions... >>>>> so it could actually bring substantial speedups for this 2 x 1 case). >>>>> >>>>> As you can guess from the documentation of the above function, it has >>>>> been mainly designed with an "exact arithmetic" context in mind, >>>>> typically >>>>> working over a finite field or the rationals. So it is my turn to write >>>>> that I would be very interested in reading your comments on how this >>>>> existing method behaves in your context. >>>>> >>>>> >>>>> Le dimanche 10 novembre 2019 14:32:52 UTC+1, Emmanuel Charpentier a >>>>> écrit : >>>>>> >>>>>> Dear list, >>>>>> >>>>>> IMHO, Sage may use an implementation of Padé approximants (similar t >>>>>> its implementation of Taylor series), if only for numerical use reasons. >>>>>> Currently, this can be done: >>>>>> >>>>>> - by wrapping the Maxima functions taylor and pade (Maxima's pade >>>>>> needs a Maxima taylor development object, which does not cleanly >>>>>> translate to Sage); >>>>>> - by using the PowerSeriesRing.pade method. >>>>>> >>>>>> Some trials have shown that the latter method, as advertised in its >>>>>> documentation, is indeed unsuitable even for moderately large degrees of >>>>>> the numerator and denominator: the expression thus obtained are >>>>>> extremely >>>>>> unwieldy large and are slow to evaluate numerically. >>>>>> In contrast, the algorithm used by Maxima, gives easily usable >>>>>> results (even if they can be enhanced by expansion and possibly >>>>>> factorization and/or simplification). But using it worsens our >>>>>> dependence >>>>>> on Maxima. >>>>>> Hence, a couple of questions: >>>>>> >>>>>> *Algorithms:* >>>>>> >>>>>> Do you have pointers to better algorithms for Padé approximants >>>>>> computation (especially for the multivariate case ? (These might also be >>>>>> helpful in the implementation of PowerSeriesRing.pade ...) >>>>>> >>>>>> *User interface:* >>>>>> >>>>>> We can follow our current taylor() convention, which is a bit of a >>>>>> straightjacket in the multivariate case,imposing the same degree for all >>>>>> the developments wrt the different variables. >>>>>> We can also allow so specify different degrees for the development >>>>>> wrt the different variables (this can make sense for very asymetric >>>>>> functions). >>>>>> >>>>>> Suggestions ? >>>>>> >>>>>> *Multivariate case:* >>>>>> >>>>>> An elementary implementation (see (Huard and Guillaume, 2000) >>>>>> <https://www.sciencedirect.com/science/article/pii/S037704270000337X> >>>>>> for example) is to iteratively create a Padé approximant for the >>>>>> successive >>>>>> variables. i. e. if p(f, x) denotes the Pade approximant wrt the >>>>>> variable x, you end up computing (p(p(p(f,x),y),z) (the >>>>>> implementation is trivial). The paper I quoted hints that there are >>>>>> better >>>>>> solutions, but is a bit above my pay grade (my initial formation is in >>>>>> dentistry and surgery...). >>>>>> >>>>>> Do you have suggestions on this point ? >>>>>> >>>>>> More generally, any comments are welcome ! >>>>>> >>>>>> -- >> You received this message because you are subscribed to the Google Groups >> "sage-devel" group. >> To unsubscribe from this group and stop receiving emails from it, send an >> email to sage-...@googlegroups.com <javascript:>. >> To view this discussion on the web visit >> https://groups.google.com/d/msgid/sage-devel/5b020b3f-fee8-4503-b190-4dbcfb897675%40googlegroups.com >> >> <https://groups.google.com/d/msgid/sage-devel/5b020b3f-fee8-4503-b190-4dbcfb897675%40googlegroups.com?utm_medium=email&utm_source=footer> >> . >> > -- You received this message because you are subscribed to the Google Groups "sage-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/5556869a-7ed9-42f5-8f18-d403d446c2d9%40googlegroups.com.