Open source implementation of rational approximation to a function

Question

I am looking for some open source implementation (any of Python, C, C++, Fortran is fine) of rational approximation to a function. Something along the article [1]. I give it a function and it gives me back two polynomials, whose ratio is the approximation on the given interval and the error is oscillating with the same amplitude and it is an optimal approximation, or close to it.

Here is what I found:

• It looks like chebfun can do that, but I don't have access to Matlab *.

• There is a simple program in the section 5-13 "Rational Chebyshev Approximation" in Numerical Recipes (NR).

• Mathematica has EconomizedRationalApproximation and MiniMaxApproximation

I was wondering whether there is something newer (possibly better tested) than the NR code.

My application is that I have a set of special functions, about 10, which are given either as hypergeometric series, or some formula which has numerical cancellations, and I want to have a robust, fast and accurate evaluation function, that is being called in the inner most loop of calculating two particle matrix elements in Hartree Fock calculations. I put a simple example of a function that works for me at [2]. As you can see, it's either using a direct formula, or a series around x=0, that I calculated using SymPy. It sort of works, but the accuracy is not great, around x=1 about half of the significant digits are lost (but for x=0.1 as well as x=1e5 it is accurate to almost all significant digits). I am looking for a better approximation.

[1] Deun, J., & Trefethen, L. N. (2011). A robust implementation of the Carathéodory-Fejér method for rational approximation. BIT Numerical Mathematics, 51(4), 1039–1050. doi:10.1007/s10543-011-0331-7 (June 2010 e-print)

[2] https://gist.github.com/3831580

(*) Nor to the chebfun website, which gives me 404, but Pedro suggested that must be my provider issue.

Answer 1

Might this be of any help? http://www.alglib.net/interpolation/rational.php

Answer 2

Doing one-off best rational approximations can often be accomplished by "manual" iterations of the Remez algorithm: interpolate a rational approximation with (relative or absolute) alternating sign errors at an initial guess for interpolation points, locate one (or more) points where the actual error exceeds that of the interpolation points and pivot (swapping out one or more of the guessed points for these, in a way that preserves alternating sign errors). It's a good way to understand the theory from the ground up. When the required accuracy is modest, I find a spreadsheet implementation helps me feel my way through the calculations and adapt problems through cut-and-paste.

To give a simple illustration, consider the non-polynomial function $y=1/x$ on interval $1 \le x \le 2$. The secant line through the endpoints looks like:

If we lower the straight line by half the distance at the point of "maximum error", then we will have three "equi-oscillatory" points of absolute error counting the two endpoints and this critical point in the interior. Thus a best linear polynomial approximation to a smooth convex function can be found almost by inspection.

The Chebfun software itself is open source (BSD-style license) since v4.0. Not sure why links to their site may have given a 404 error back at that time, but they are now at v5.2.1. (They also have a GitHub repository.) The Chebfun project aims at more than best rational approximations, with the central concept being a Matlab vector-like syntax for (univariate real) functions on an interval.

To complete the circle would require getting it working under Octave rather than Matlab. There was an Octave-maintainers thread about this started in 2010. A 2012 paper by Chebfun maintainers/authors suggests some effort has gone into Octave compatibility from their side.

Looking at the Chebfun-related projects page suggests that because the Chebfun software is open source, it might be possible to do this port, or to make use of one of the open source packages listed there, such as Olivier Verdier's pychebfun, a Python Chebyshev Functions project hosted at GitHub.