GitHub

Rational function approximation in Julia

Documentation for RationalFunctionApproximation.jl.

This package uses the continuous form of the AAA algorithm to adaptively compute rational approximations of functions on intervals and other domains in the complex plane. See AAA rational approximation on a continuum (or the arXiv version).

Here's a smooth, gentle function on the interval $[-1, 1]$:

using RationalFunctionApproximation, CairoMakie
CairoMakie.update_theme!(size = (400, 250), fontsize=11)
const shg = current_figure

f = x -> exp(cos(4x) - sin(3x))
lines(-1..1, f)
Example block output

To create a rational function that approximates $f$ well on this domain, we use the continuous form of the AAA algorithm:

r = aaa(f)
Barycentric function with 20 nodes and values:
    -1.0=>0.598982,  -0.986111=>0.599042,  -0.944444=>0.605944,  …  1.0=>0.451688

The result is a type (19,19) rational approximant that can be evaluated like a function:

f(0.5) - r(0.5)
3.186340080674199e-14

We see that this approximation has more than 13 accurate digits over most of the interval:

lines(-1..1, x -> f(x)-r(x))
Example block output

The rational approximant interpolates $f$ at greedily selected nodes:

x = nodes(r)
scatter!(x, 0*x, markersize = 8, color=:black)
shg()
Example block output

Here's another smooth example, the hyperbolic secant function:

f = sech
r = aaa(f)
Barycentric function with 7 nodes and values:
    -1.0=>0.648054,  -0.951515=>0.672091,  -0.466667=>0.900188,  …  1.0=>0.648054

We can verify that this is accurate to 14 digits:

x = range(-1, 1, 1000)
extrema(f.(x) - r.(x))
(-3.552713678800501e-15, 6.217248937900877e-15)

Since the sech function has poles in the complex plane, the rational approximant $r$ will have corresponding poles:

using DomainColoring
CairoMakie.update_theme!(size = (360, 360))

domaincolor(r, [-8, 8, -8, 8], abs=true)
Example block output

The poles closest to the interval are found to about 10 digits, while more distant ones are less accurate:

poles(r) / π
6-element Vector{ComplexF64}:
 -2.6605420930314816e-6 - 1.50435605118539im
 -2.6605420930314816e-6 + 1.50435605118539im
  9.914614203220754e-11 - 0.4999999922042446im
  9.914614203220754e-11 + 0.4999999922042446im
  2.6093540880194524e-5 - 2.2138848060268854im
  2.6093540880194524e-5 + 2.2138848060268854im

Here's an example with a more interesting structure of poles and zeros:

f = x -> tanh(10*(x - 0.1)^2)
domaincolor(aaa(f), abs=true)
Example block output