I am pleased to announce the initial release of RationalFunctionApproximation.jl, which provides automatic rational function approximation in the complex plane.
It uses an updated version of the AAA algorithm (Driscoll, Nakatsukasa, Trefethen) to adaptively sample the domain of the function to ensure accurate approximations when a singularity is very close to or even on the domain.
The domain can be an interval or any path in the complex plane, and the algorithm iterates on the degree of the approximation until the domain, and optionally its interior or exterior, is free of poles. This usually is enough to overcome the problem of undesired poles in a rational approximation.
Example of a function with three zeros on the unit circle and a pole exterior to that:
(Plot courtesy of Makie.jl and DomainColoring.jl.) The rational approximant shown here is of type (9,9), interpolating the function at 10 points on the unit circle, and is accurate to more than 13 digits on the circle.
Thanks for this nice package. If you will, what is the rationale for publishing the reference article on arXiv, a free service for scientific articles, and then using proprietary Matlab code in it instead of Julia?
The Matlab code we present in the paper is not itself proprietary, and it’s a language that is well understood by most of the readership of the journal we have targeted. A google search for “matlab” on arxiv.org turns up over 12 million hits, so it’s widely recognized by the arXiv readership, as well.
If we publish Julia code alongside Matlab code, it’s going to do more to draw in current Matlab users than pretending Matlab doesn’t exist and offering only code they can’t read or use.
BaryRational works with the original discrete AAA algorithm, which uses an initial discretization of the domain that does not change. This is often fine, unless poles get as close to the domain as the spacing between the sample points.
This package uses a new, continuous form of AAA referenced in my post that adaptively samples the domain, which is possibly more efficient for smooth functions and more reliable for functions with singularities near the domain.
The approximations are based on least squares minimization, as opposed to the max norm. However, there is an iterative reweighting scheme that can approach the minimax approximation. There’s a page in the docs showing this a bit.
