5.4 Rational function approximation

Rational function basics

Rational function approximation replaces a complicated function with the ratio of two polynomials. Where polynomial approximation struggles with singularities, poles, or functions that change rapidly, rational approximations often succeed with far fewer terms.

Definition of rational functions

A rational function takes the form:

$R(x) = \frac{P(x)}{Q(x)}$

where $P(x)$ is a polynomial of degree $m$ and $Q(x)$ is a polynomial of degree $n$. The notation [m/n] describes the type: $m$ is the numerator degree, $n$ is the denominator degree. So a [3/2] approximation has a cubic on top and a quadratic on the bottom, giving you 6 free coefficients total (after normalizing $Q$).

The denominator is what gives rational functions their real power. Because $Q(x)$ can have roots, the approximation can model poles and asymptotes that polynomials simply cannot represent.

Properties of rational functions

• The domain is all real numbers except where $Q(x) = 0$
• Vertical asymptotes occur at the roots of $Q(x)$
• Horizontal asymptotes exist when $m \leq n$: the asymptote is $a_m / b_n$ if $m = n$, or $0$ if $m < n$
• When $m > n$, the function grows without bound, and you get oblique or higher-order asymptotic behavior
• This built-in ability to model singular behavior is exactly why rational approximations outperform polynomials for many function classes

Advantages over polynomial approximation

• Functions with poles or branch points are natural targets for rational approximation, since polynomials blow up or oscillate wildly near singularities
• Fewer total parameters are often needed. A [4/4] rational approximation (9 parameters) can outperform a degree-15 polynomial (16 parameters) for functions like $\frac{1}{1+25x^2}$
• Rational functions capture asymptotic behavior at $\pm\infty$ correctly, while polynomials always diverge
• Over large intervals, polynomial approximations suffer from Runge's phenomenon; rational approximations are far more resistant to this

Types of rational approximations

Three major approaches dominate the field, each suited to different situations: matching local behavior (Padé), minimizing worst-case error (Chebyshev rational), and interpolating given data (barycentric rational).

Padé approximation

Padé approximation builds a rational function that matches the Taylor series of a target function at a single point (usually $x = 0$). You choose degrees $m$ and $n$, then require that the first $m + n + 1$ Taylor coefficients agree.

This makes Padé approximants especially effective when you have a known series expansion. They're widely used for special functions like $e^x$, $\log(1+x)$, and matrix exponentials. The [2/2] Padé approximant to $e^x$, for example, is far more accurate over a wide range than the degree-4 Taylor polynomial, even though both use the same amount of series information.

Chebyshev rational approximation

Instead of matching at a single point, Chebyshev rational approximation minimizes the maximum error over an entire interval (the minimax principle). The resulting error curve equioscillates, meaning it alternates between its maximum and minimum values at several points.

The Remez algorithm is the standard iterative method for computing these approximations. It provides near-optimal uniform accuracy for continuous functions on finite intervals and is heavily used in library implementations of mathematical functions.

Barycentric rational interpolation

This approach interpolates function values at a set of nodes, but expresses the result in barycentric form:

$R(x) = \frac{\sum_{k=0}^{n} \frac{w_k f_k}{x - x_k}}{\sum_{k=0}^{n} \frac{w_k}{x - x_k}}$

where $w_k$ are barycentric weights, $f_k$ are function values, and $x_k$ are nodes. This form avoids computing polynomial coefficients explicitly, which dramatically improves numerical stability. It works well for both equally spaced and clustered nodes, and remains stable even at high interpolation degrees where other methods break down.

Padé approximation methods

Padé approximants are the most common entry point into rational approximation. They generalize Taylor polynomials by allowing a denominator, which lets them capture behavior that truncated power series miss entirely.

Construction of Padé approximants

The [m/n] Padé approximant to $f(x)$ is written as:

$R_{m,n}(x) = \frac{a_0 + a_1 x + \cdots + a_m x^m}{1 + b_1 x + \cdots + b_n x^n}$

Note the convention: the denominator's leading coefficient is normalized to 1, which removes one degree of freedom and makes the system well-determined.

To construct it:

1. Write out the Taylor series of $f(x)$ to at least order $m + n$: $f(x) = c_0 + c_1 x + c_2 x^2 + \cdots$
2. Multiply both sides of $f(x) \cdot Q(x) = P(x)$ and expand
3. Match coefficients of $x^0, x^1, \ldots, x^{m+n}$ on both sides
4. This yields a system of $m + n + 1$ linear equations. The denominator coefficients $b_1, \ldots, b_n$ are found first (from the higher-order equations), then the numerator coefficients $a_0, \ldots, a_m$ follow directly
5. Solve and assemble $R_{m,n}(x)$

Order conditions

The key property: the Padé approximant agrees with the Taylor series through order $m + n$. That is:

$f(x) - R_{m,n}(x) = O(x^{m+n+1}) \quad \text{as } x \to 0$

This means an [m/n] approximant uses the same information as a degree-$(m+n)$ Taylor polynomial but distributes it between numerator and denominator. Diagonal approximants (where $m = n$) tend to have the best convergence properties and are the most commonly used in practice.

Uniqueness and existence

The [m/n] Padé approximant is unique when it exists. Existence is guaranteed when a certain Hankel determinant built from the Taylor coefficients is nonzero. Specifically, the $n \times n$ Hankel matrix:

$H = \begin{pmatrix} c_{m-n+1} & c_{m-n+2} & \cdots & c_m \\ c_{m-n+2} & c_{m-n+3} & \cdots & c_{m+1} \\ \vdots & & & \vdots \\ c_m & c_{m+1} & \cdots & c_{m+n-1} \end{pmatrix}$

must be nonsingular. When this determinant vanishes, the approximant is said to have a defect, and the actual approximant that exists will be of lower effective order. In such cases, you may need to adjust $m$ or $n$ or use alternative techniques.

Chebyshev rational approximation

Chebyshev rational approximation targets the best possible uniform accuracy over a whole interval, rather than matching behavior at a single point. This makes it the method of choice for building function libraries and high-accuracy numerical routines.

Chebyshev polynomials review

Chebyshev polynomials of the first kind satisfy the recurrence:

$T_n(x) = 2x\,T_{n-1}(x) - T_{n-2}(x), \quad T_0(x) = 1, \quad T_1(x) = x$

They are orthogonal on $[-1, 1]$ with weight $\frac{1}{\sqrt{1 - x^2}}$ and satisfy $|T_n(x)| \leq 1$ on that interval. Their zeros are:

$x_k = \cos\!\left(\frac{(2k-1)\pi}{2n}\right), \quad k = 1, 2, \ldots, n$

The crucial property for approximation: among all monic polynomials of degree $n$, $2^{1-n} T_n(x)$ has the smallest maximum absolute value on $[-1,1]$. This minimax property is why Chebyshev polynomials appear throughout approximation theory.

Minimax approximation principle

The best rational approximation of type [m/n] to $f(x)$ on $[a, b]$ minimizes:

$\|f - R\|_\infty = \max_{x \in [a,b]} |f(x) - R(x)|$

By the equioscillation theorem (a generalization of Chebyshev's theorem to rational functions), the optimal error curve must alternate in sign at least $m + n + 2$ times. This equioscillation property is both the theoretical characterization and the practical stopping criterion for algorithms.

Remez algorithm

The Remez algorithm computes near-minimax rational approximations iteratively:

1. Initialize with $m + n + 2$ reference points in $[a, b]$ (Chebyshev nodes are a good starting choice)

2. Solve the leveled equations: find a rational function $R(x)$ and a leveled error $\epsilon$ such that $f(x_k) - R(x_k) = (-1)^k \epsilon$ at each reference point

3. Locate the actual points of maximum error of $f(x) - R(x)$ over the full interval

4. Update the reference set with these new extrema

5. Repeat until the error equioscillates to within a specified tolerance

Convergence is typically fast for smooth functions, often requiring only 5-10 iterations. For functions with nearby singularities, convergence can be slower and the algorithm may need careful initialization.

Error analysis

Understanding how and why rational approximations fail (or succeed) is essential for choosing the right method and the right degree.

Error bounds for rational approximation

Error bounds depend heavily on the analyticity of the target function. For a function $f$ analytic in a region of the complex plane, the best [n/n] rational approximation on $[-1,1]$ satisfies:

$\|f - R_{n,n}\|_\infty \leq C \cdot e^{-c\sqrt{n}}$

for some constants $C, c > 0$ depending on the singularity structure. This is a root-exponential rate, which is much faster than polynomial approximation when $f$ has singularities near the interval but is still dramatically slower than the geometric rate you get for entire functions.

For functions with poles at known locations, the error bounds tighten considerably because the rational approximation can place its own poles nearby.

Convergence rates

• Entire functions (no singularities anywhere): both polynomial and rational approximations converge geometrically, so rational offers little advantage
• Functions with poles: rational approximation converges geometrically, while polynomial approximation converges only algebraically. This is where rational methods truly shine
• Functions with branch points: rational approximation achieves root-exponential convergence ($O(e^{-c\sqrt{n}})$), much better than the algebraic rate of polynomials but slower than geometric
• Diagonal Padé approximants to Stieltjes functions converge throughout the domain, a result guaranteed by the Montessus de Ballore theorem when the pole locations are known

Stability considerations

Numerical stability can be a real concern with rational approximations:

• Pole-zero cancellation: when a pole and zero of $R(x)$ are close together, evaluation becomes ill-conditioned. Small perturbations in coefficients can move poles significantly
• Coefficient sensitivity: the monomial basis representation of $P(x)$ and $Q(x)$ can have very large coefficients that nearly cancel, leading to catastrophic loss of precision
• Barycentric evaluation largely avoids these issues by never forming the coefficients explicitly
• The condition number of the linear system used to compute Padé coefficients can grow rapidly with degree, so stable solvers (e.g., SVD-based) are preferred for high-order approximants

Applications in numerical analysis

Numerical integration

Rational approximations handle integrands that give standard polynomial-based quadrature trouble. If the integrand has a singularity inside or near the integration interval, replacing it with a rational approximation before integrating can dramatically improve accuracy. Padé-based quadrature rules are also used for oscillatory integrals and integrals with endpoint singularities.

Solution of differential equations

Padé approximation plays a central role in time-stepping for stiff ODEs. The stability function of many implicit Runge-Kutta methods is a Padé approximant to $e^z$. For example, the implicit midpoint rule has stability function equal to the [1/1] Padé approximant of $e^z$, and the two-stage Gauss method corresponds to [2/2]. These rational stability functions allow A-stable and L-stable methods, which polynomial-based explicit methods cannot achieve.

Rational approximations also appear in spectral methods and collocation schemes for boundary value problems.

Approximation of special functions

Most numerical libraries compute functions like $e^x$, $\text{erf}(x)$, Bessel functions, and the gamma function using rational approximations internally. A typical implementation might use argument reduction to map the input to a small interval, then evaluate a precomputed minimax rational approximation on that interval. This approach achieves full machine precision with just a handful of multiplications and one division.

Implementation techniques

Computation of coefficients

• For Padé approximants, solve the linear system using a stable method (QR or SVD decomposition), not naive Gaussian elimination
• Clenshaw's algorithm efficiently evaluates Chebyshev expansions without forming the polynomial explicitly
• Barycentric formulas compute the rational interpolant directly from function values and weights, bypassing coefficient computation entirely
• When choosing $m$ and $n$, start with diagonal ($m = n$) and adjust based on the function's behavior at the boundaries of the interval

Handling of poles and zeros

Spurious poles are a well-known pitfall. A rational approximation may introduce poles that don't correspond to any singularity of the target function. Strategies include:

• Pole detection: after computing the approximation, find the roots of $Q(x)$ and check whether any fall inside or near the approximation interval
• Partial fraction decomposition: rewriting $R(x)$ in partial fractions can improve evaluation near poles and reveal the residue structure
• Froissart doublets: spurious pole-zero pairs that nearly cancel. These are a sign that the chosen [m/n] type has too many degrees of freedom. Reducing $n$ often eliminates them
• Robust Padé methods (e.g., the SVD-based approach of Gonnet, Güttel, and Trefethen) automatically detect and remove spurious poles

Numerical stability issues

• Normalize the denominator polynomial (e.g., set $b_0 = 1$ or $\|Q\| = 1$) to avoid scaling problems
• Use compensated summation (Kahan summation) when evaluating polynomials in the numerator and denominator
• For very high accuracy requirements, consider using extended precision for the coefficient computation step, then round back to double precision for evaluation
• Monitor the condition number of the coefficient matrix during construction. If it exceeds $10^{12}$ or so in double precision, the resulting approximation is likely unreliable

Comparison with other methods

Rational vs polynomial approximation

The rule of thumb: if your function is smooth and entire, use polynomials. If it has singularities, poles, or needs to be accurate over a large interval, use rational approximations.

Rational vs spline approximation

• Splines provide local control (changing data in one region doesn't affect distant regions) while rational approximations are global
• Splines handle discontinuities and preserve shape properties (monotonicity, convexity) more naturally
• Rational functions are better for approximating analytic functions with known singularity structure
• Splines are the standard choice for data fitting and CAD; rational functions dominate in special function evaluation

Trade-offs in accuracy vs complexity

Rational approximations achieve higher accuracy per parameter but at the cost of:

• A nonlinear approximation problem (the denominator coefficients enter nonlinearly)
• The possibility of spurious poles
• More complex error analysis
• Division in evaluation (which can be costly on some hardware and introduces the possibility of division by zero)

For many practical problems, the accuracy gains justify these costs. But when a polynomial or spline does the job adequately, there's no reason to introduce the extra complexity.

Advanced topics

Multivariate rational approximation

Extending rational approximation to functions of several variables is substantially harder. Tensor product constructions (applying one-dimensional rational approximation along each coordinate direction) are the simplest approach but scale poorly with dimension. Sparse grid methods help mitigate the curse of dimensionality. The theory of multivariate Padé approximants is less complete than the univariate case, and uniqueness is not always guaranteed.

Rational approximation on unbounded domains

For functions defined on $[0, \infty)$ or $(-\infty, \infty)$, two main strategies exist:

• Conformal mappings: transform the unbounded domain to a finite interval (e.g., $x \mapsto \frac{x-1}{x+1}$ maps $[0,\infty)$ to $[-1,1]$), then apply standard rational approximation on the finite interval
• Direct construction: build rational approximations with prescribed asymptotic behavior at infinity, ensuring the approximation decays or grows at the correct rate

Weighted approximation, where you minimize $\|w(x)(f(x) - R(x))\|$ for some weight function $w$, is also common for functions with rapid decay.

Adaptive rational approximation methods

Modern adaptive methods automatically select the type [m/n] and node placement:

• The AAA algorithm (Adaptive Antoulas-Anderson) iteratively selects interpolation points from a candidate set, building a barycentric rational approximation that greedily minimizes the approximation error. It's become a standard tool due to its simplicity and robustness
• Hybrid methods combine rational approximation with piecewise polynomial techniques, using rational functions near singularities and polynomials in smooth regions
• Error estimators based on comparing approximations of different degrees guide the adaptive refinement process