# Interpolation

Implementation: `eko.interpolation`

In order to obtain the operators in an PDF independent way we use approximation theory. Therefore, we define the basis grid

from which we define our interpolation

Thus each grid point \(x_j\) has an associated interpolation polynomial \(p_j(x)\)
(represented by `eko.interpolation.BasisFunction`

).
We interpolate in \(\ln(x)\) using Lagrange interpolation among the nearest
\(N_{degree}+1\) points, which renders the \(p_j(x)\): polynomials of order
\(O(\ln^{N_{degree}}(x))\).

## Algorithm

First, we split the interpolation region into several areas (represented by
`eko.interpolation.Area`

), which are bound by the grid points:

Note, that we include the right border point into the definition, but not the left which keeps all areas disjoint. This assumption is based on the physical fact, that PDFs do have a fixed upper bound (\(x=1\)), but no fixed lower bound.

Second, we define the interpolation blocks, which will build the interpolation polynomials and contain the needed amount of points:

Now, we construct the interpolation in a bottom-up approach for a given point
\(\bar x\) for the interpolating function \(\bar f(\bar x)\).
(The actual implementation is done in `eko.interpolation.InterpolatorDispatcher`

)

Determine the relevant active area \(\bar A(\bar x)\) :

Determine the associated block \(\bar B(\bar A)\) : Choose the block in which \(\bar A(\bar x)\) is located most central. For an odd number of \(N_{degree}\) choose the block which is located higher, i.e. closer to \(x=1\) (following the choice of borders for the areas). At the border of the grid, choose the “most central”.

Only the points and their associated polynomials which are in \(\bar B\) participate in the interpolation of \(\bar f(\bar x)\). All other polynomials vanish:

The active polynomials build a Lagrange interpolation over their associated points

The generated polynomials are orthonormal in the following sense

and thus continuous, however, they are not necessarily differentiable. They are also complete in the following sense

## Example

To outline the interpolation algorithm, let’s consider the following example configuration:

The mapping between areas and blocks are then given by

```
A0 -> B0, A1 -> B0, A2 -> B1, A3 -> B2, A4 -> B3, A5-> B4, A6 -> B5, A7 -> B5
```

The generated polynomials then look like:

## Change Interpolation Basis

It is always possible to change interpolation basis, that corresponds to change interpolation grid.

The way it is done is the following:

So the change of basis to apply to coefficients/polynomials is:

that corresponds to evaluate the old polynomials \(p_j\) (i.e. took the place of the continuous function in \(\bar f\)) on the new points \(\tilde x_k\).

Take care that in the above sections it has been explained that the target
function is interpolated by approximating it with a **piecewise polynomial**.

Piecewise Polynomials Interpolation

It is very relevant to notice that it is **not a polynomial**, indeed being
defined piecewise only the *function is* **continuous**, but **not its
derivatives**.

Then even interpolating on a bigger but different grid might cause a further approximation.

Indeed no polynomial is able to produce a discontinuity in the derivatives, thus if an old grid point falls in the middle of a new grid area an approximation is expected (that intuitevely it would happen only if the grid were smaller, but it is also happening for larger poorly chosen ones).

Since further interpolating will cause a loss of accuracy it is *recommended* to
have the new grid being a subset or a superset of the previous one, or in
general *to share as many points as possible*.