16
$\begingroup$

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.

$\endgroup$

2 Answers 2

5
$\begingroup$

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

$\endgroup$
1
  • 1
    $\begingroup$ Thanks, I think this will do. I need to learn to use it. At the moment, I am under time press, so I just used MiniMaxApproximation from Mathematica to get the job done, but I want to learn to do this using opensource tools. $\endgroup$ Oct 5, 2012 at 20:11
11
$\begingroup$

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: enter image description here

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.

$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.