Rational functions

(updated )

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

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 $$

Adaptive improvement

Increase $p$

Increase $q$

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

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

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