# Rational functions

hpq-FEM: approximate functions using a mesh (h) of rational functions of numerator degree p and denominator degree q.

$$p(x) = \frac{a_0 + a_1 x + a_2 x^2 \cdots a_p x^p}{1 + b_1 x + b_2 x^2 \cdots b_q x^q}$$

The main goal is to reduce the number of degrees-of-freedom in the problem while still approximating the solution within the error tolerance.

## Evaluating Polynomials

Two methods are optimal.

### Langrange

$$p(x) = \left(\prod_i (x - x_i)\right) \sum_i \frac{w_j }{(x - x_i)} f(x_i)$$

with

$$w_j = \frac{1}{\prod_{1\neq j} x_j - x_i}$$

### Barycentric

$$p(x) = \frac {\sum_i \frac{w_j }{x - x_i} f(x_i)} {\sum_i \frac{w_j }{x - x_i}}$$

with the same weigths $w_i$.

The barycentric one becomes unstable when $x$ deviates to much from $x_i$. The Langrange version does not suffer from this instability.

Rational functions can be evaluated using two Langrange evaluations, dividing both results in the end.

## Solving

We want to approximate $f(x)$ with $p(x)$. We do this by equating the values on a set of points $x_i$. Take $y_i = f(x_i)$ and set $p(x_i) = y_i$.

$$\frac{a_0 + a_1 x_i + a_2 x_i^2 \cdots a_p x_i^p}{ 1 + b_1 x_i + b_2 x_i^2 \cdots b_q x_i^q} = y_i$$

$$a_0 + a_1 x_i + a_2 x_i^2 \cdots a_p x_i^p = y_i (1 + b_1 x_i + b_2 x_i^2 \cdots b_q x_i^q)$$

$$a_0 + a_1 x_i + a_2 x_i^2 \cdots a_p x_i^p = y_i + b_1 y_ix_i + b_2 y_ix_i^2 \cdots b_q y_i x_i^q$$

$$a_0 + a_1 x_i + a_2 x_i^2 \cdots a_p x_i^p • b_1 y_i x_i - b_2 y_ix_i^2 \cdots b_q y_i x_i^q = y_i$$

or equivalently

$$a_0 + a_1 x_i + a_2 x_i^2 \cdots a_p x_i^p - b_1 y_i x_i - b_2 y_i x_i^2 \cdots b_q y_i x_i^q) = y_i$$

This can be solved using the system

$$\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots &x_0^p & -y_0 x_0 & -y_0 x_0^2 & \cdots & -y_0 x_0^q \\ \vdots & \vdots &\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & x_n^2 & \cdots &x_n^p & -y_n x_n & -y_n x_n^2 & \cdots & -y_n x_n^q \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_p \\ b_1 \\ \vdots \\ b_q \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \\ \end{bmatrix}$$

where $n = p + q + 1$.

Method from Numerical Recipes

$$y_i = y_i + e \frac{e_i}{|e_i|}$$

$$\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^p & -y_0' x_0 & -y_0' x_0^2 & \cdots & -y_0' x_0^q \\ \vdots & \vdots &\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & x_n^2 & \cdots &x_n^p & -y_n x_n & -y_n x_n^2 & \cdots & -y_n x_n^q \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_p \\ b_1 \\ \vdots \\ b_q \end{bmatrix} = \begin{bmatrix} y_0' \\ y_1 \\ \vdots \\ y_n \\ \end{bmatrix}$$

## Degrees of freedom

The degrees of freedom of an element is

$$p + q + 1$$

### Subdivide

An element $(p,q)$ can subdivided in two elements of degree $(p_a, q_a)$ and $(p_b, q_b)$ if we allow the total degrees of freedom to increase by one, we are constraint by

$$p_a + q_a + p_b + q_b = p + q$$

The number of subdivision options is equal to the number of ways $p+q$ can be distributed into four bins. This is the stars-and-bars problem with solution:

$$\binom{p + q + 3}{p + q}$$

## Benchmark problem

\begin{aligned} u(0) &= γ_0 \\ u(1) &= γ_1 \\ -u'' &= f \end{aligned}

with $γ_0$, $γ_1$ and $f$ such that exact solution is

$$u(x) = \mathrm{atan} \left(60\left(x - \frac{\pi}{3}\right) \right)$$

## References

• Rational approximation. https://arxiv.org/pdf/math/0101042.pdf

• Numerical Recipes 3rd edition page 248.

• Comentary on the numerical recipes algorithm. https://chasethedevil.github.io/post/rational_fit/

• Alternative design matrix. https://math.stackexchange.com/questions/924482/least-squares-regression-matrix-for-rational-functions

• The remez / minimax implementation in Chebfun appears to be the state of the art for rational approximation. https://people.maths.ox.ac.uk/trefethen/remezR2_arxiv.pdf http://www.chebfun.org/ http://www.chebfun.org/docs/guide/guide04.html http://www.chebfun.org/examples/approx/RationalAbsx.html https://github.com/chebfun/chebfun/blob/master/minimax.m https://jncf2018.lip6.fr/files/slides/filip.pdf

• https://people.maths.ox.ac.uk/trefethen/barycentric.pdf

• https://people.maths.ox.ac.uk/trefethen/publication/PDF/2011_138.pdf

• AAA algorithm. https://arxiv.org/abs/1612.00337 http://guettel.com/rktoolbox/examples/html/example_aaa.html http://www.chebfun.org/examples/approx/AAAApprox.html https://jncf2018.lip6.fr/files/slides/filip.pdf https://people.maths.ox.ac.uk/trefethen/AAAfinal.pdf

• AAA-Lawson. https://arxiv.org/pdf/1705.10132.pdf

• History of rational approximation and Floater-Hormann. http://www.alglib.net/interpolation/rational.php https://cs.fit.edu/~dmitra/SciComp/13Spr/MikeChirstian-NRPresentation.pdf

https://github.com/chebfun/chebfun/blob/a26ae3cab130646c05fca428c1cbdfdd95208a49/minimax.m

• Optimiality of Barycentric evaluation. https://www.researchgate.net/profile/Richard_Baltensperger/publication/226684962_Recent_Developments_in_Barycentric_Rational_Interpolation/links/00b4951518a41531de000000.pdf http://eprints.maths.ox.ac.uk/1400/1/NA-11-11.pdf http://faculty.cs.tamu.edu/schaefer/research/pyramid.pdf

Remco Bloemen
Math & Engineering
https://2π.com