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

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