qmat.lagrange

Base module for Barycentric Lagrange Approximation, based on [Berrut & Trefethen, 2004]. Allows to easily build integration / interpolation / derivative matrices, from any list of node points.

Examples

>>> # Base usage to generate a quadrature matrix
>>> from qmat.lagrange import LagrangeApproximation, np
>>>
>>> grid = np.linspace(0, 1, num=5)
>>> approx = LagrangeApproximation(grid)
>>> Q = approx.getIntegrationMatrix([(0, tau) for tau in grid])
>>>
>>> # Interpolation
>>> fGrid = np.linspace(0, 1, num=200)
>>> u = np.exp(grid)
>>> P = approx.getInterpolationMatrix(fGrid)
>>> uFine = P @ u
>>>
>>> # Alternative interpolation using the object as a function
>>> uFine = approx(fGrid)
>>>
>>> # Derivative
>>> D = approx.getDerivativeMatrix()
>>> du = D @ u

Classes

LagrangeApproximation

Class approximating any function on a given set of points using barycentric

Functions

`computeFejerRule`(n)	Compute a Fejer rule of the first kind, using DFT [Waldvogel, 2006].
`getSparseInterpolationMatrix`(inPoints, outPoints, order)	Get a sparse interpolation matrix from inPoints to outPoints of order

Module Contents

computeFejerRule(n)[source]

Compute a Fejer rule of the first kind, using DFT [Waldvogel, 2006]. Inspired from quadpy (https://github.com/nschloe/quadpy @Nico_Schlömer)

Parameters:

n (int) – Number of points for the quadrature rule.

Returns:

nodes (np.1darray(n)) – The nodes of the quadrature rule
weights (np.1darray(n)) – The weights of the quadrature rule.

class LagrangeApproximation(points, weightComputation='AUTO', scaleWeights=False, scaleRef='MAX', duplicates='USE_LEFT', fValues=None)[source]

Class approximating any function on a given set of points using barycentric Lagrange interpolation.

Let note \((t_j)_{0\leq j<n}\) the set of points, then any scalar function \(f\) can be approximated by the barycentric formula :

\[p(x) = \frac{\displaystyle \sum_{j=0}^{n-1}\frac{w_j}{x-x_j}f_j} {\displaystyle \sum_{j=0}^{n-1}\frac{w_j}{x-x_j}},\]

where \(f_j=f(t_j)\) and

\[w_j = \frac{1}{\prod_{k\neq j}(x_j-x_k)}\]

are the barycentric weights. The theory and implementation is inspired from [1].

Parameters:

points (list, tuple or np.1darray) – The given interpolation points, no specific scaling, but must be ordered in increasing order.
weightComputation (str, optional) – Algorithm used to compute the barycentric weights. Can be :
- ‘FAST’ : uses the analytic formula (unstable for large number of points)
- ‘STABLE’ : uses logarithmic difference and scaling of the weights
- ‘CHEBFUN’ : uses the same approach as in the chebfun package
The default is ‘AUTO’ : it tries the ‘FAST’ algorithm, and if an overflow is detected, it switches to the ‘STABLE’ algorithm.
scaleRef (str, optional) – Scaling used in the ‘STABLE’ algorithm for weight computation. Can be :
- ‘ZERO’ : scaling based on the weight for the value closest to \(t=0\).
- ‘MAX’ : scaling based on the maximum weight value.
The default is ‘MAX’.
duplicates (str) – Which strategy to use in case of duplicated values within the interpolation points. Can be :
- ‘USE_LEFT’ : uses the first value from the left in the values vector
- ‘USE_RIGHT’ : uses the first value from the right in the values vector
The default is ‘USE_LEFT’.
fValues (list, tuple or np.1darray) – Function values to be used when evaluating the LagrangeApproximation as a function

Variables:

points (np.1darray) – The interpolating points
weights (np.1darray) – The associated barycentric weights
nPoints (int (property)) – The number of points, can also be retrieve with n (legacy alias)
uniquePoints (np.1darray) – The unique interpolating points. When there is no duplicates, points == uniquePoints.
nUniquePoints (int (property)) – The number of unique points
duplicates (str) – The strategy used when there is duplicated interpolation points

References

weights

weightComputation = 'AUTO'

fValues = None

property nPoints: int: The number of points

n: Legacy alias for nPoints

property nUniquePoints: int: The number of unique points

property hasDuplicates: bool: Wether the points have duplicates or not

getInterpolationMatrix(times, duplicates=True) → numpy.ndarray[source]

Compute the interpolation matrix for a given set of discrete “time” points.

For instance, if we note \(\vec{f}\) the vector containing the \(f_j=f(t_j)\) values, and \((\tau_m)_{0\leq m<M}\) the “time” points where to interpolate. Then \(I[\vec{f}]\), the vector containing the interpolated \(f(\tau_m)\) values, can be obtained using :

\[I[\vec{f}] = P_{Inter} \vec{f},\]

where \(P_{Inter}\) is the interpolation matrix returned by this method.

Parameters:

times (list-like or np.1darray) – The discrete “time” points where to interpolate the function.
duplicates (bool) – Wether or not take into account duplicates in the points. This has no impact if all interpolating points are distinct. Default is True.

Returns:

P – The interpolation matrix, with \(M\) rows (size of the times parameter) and \(n\) columns.

Return type:

np.2darray(M, n)

getIntegrationMatrix(intervals, numQuad='FEJER', duplicates=True) → numpy.ndarray[source]

Compute the integration matrix for a given set of intervals.

For instance, if we note \(\vec{f}\) the vector containing the \(f_j=f(t_j)\) values, and \((\tau_{m,left}, \tau_{m,right})_{0\leq m<M}\) the different interval where the function should be integrated using the barycentric interpolant polynomial. Then \(\Delta[\vec{f}]\), the vector containing the approximations of

\[\int_{\tau_{m,left}}^{\tau_{m,right}} f(t)dt,\]

can be obtained using :

\[\Delta[\vec{f}] = P_{Integ} \vec{f},\]

where \(P_{Integ}\) is the interpolation matrix returned by this method.

Parameters:

intervals (list of pairs) – A list of all integration intervals.
numQuad (str, optional) – Quadrature rule used to integrate the interpolant barycentric polynomial. Can be :
- ‘LEGENDRE_NUMPY’ : Gauss-Legendre rule from Numpy
- ‘LEGENDRE_SCIPY’ : Gauss-Legendre rule from Scipy
- ‘FEJER’ : internally implemented Fejer-I rule
The default is ‘FEJER’.
duplicates (bool) – Wether or not take into account duplicates in the points. This has no impact if all interpolating points are distinct. Default is True.

Returns:

Q – The integration matrix, with \(M\) rows (number of intervals) and \(n\) columns.

Return type:

np.2darray(M, n)

getDerivativeMatrix(order=1, duplicates=True) → numpy.ndarray[source]

Generate derivative matrix of first or second order (or both) based on the Lagrange interpolant. The first order differentiation matrix \(D^{(1)}\) approximates

\[D^{(1)} u \simeq \frac{du}{dx}\]

on the interpolation points. The formula is :

\[D^{(1)}_{ij} = \frac{w_j/w_i}{x_i-x_j}\]

for \(i \neq j\) and

\[D^{(1)}_{jj} = -\sum_{i \neq j} D^{(1)}_{ij}`\]

The second order differentiation matrix \(D^{(2)}\) approximates

\[D^{(2)} u \simeq \frac{d^2u}{dx^2}\]

on the interpolation points. The formula is :

\[D^{(1)}_{ij} = -2\frac{w_j/w_i}{x_i-x_j}\left[ \frac{1}{x_i-x_j} + \sum_{k \neq i}\frac{w_k/w_i}{x_i-x_k} \right]\]

for \(i \neq j\) and

\[D^{(2)}_{jj} = -\sum_{i \neq j} D^{(2)}_{ij}\]

⚠️ If you want a derivative matrix with many points (~1000 or more), favor the use of weightComputation=”STABLE” when initializing the LagrangeApproximation object. If not, some (very small) weights could be approximated by zeros, which would make the computation of the derivative matrices fail …

Note

There is a typo in the formula for \(D^{(2)}\) given in the paper of Berrut and Trefethen. The formula above is the correct one.

Parameters:

order (int or str, optional) – The order of the derivative matrix, use “ALL” to retrieve both. The default is 1.
duplicates (bool) – Wether or not take into account duplicates in the points. This has no impact if all interpolating points are distinct. Default is True.

Returns:

D – Derivative matrix. If order=”ALL”, return a tuple containing all derivative matrices in increasing derivative order.

Return type:

np.2darray or tuple of np.2darray

getDerivationMatrix(*args, **kwargs) → numpy.ndarray[source]

getSparseInterpolationMatrix(inPoints, outPoints, order, gridPeriod=-1)[source]

Get a sparse interpolation matrix from inPoints to outPoints of order order using barycentric Lagrange interpolation.

The matrix will have order entries per line, and tends to be banded when both inPoints and outPoints are equispaced and cover the same interval.

Parameters:

inPoints (np.1darray) – The points you want to interpolate from
outPoints (np.1darray) – The points you want to interpolate to
order (int) – Order of the interpolation
grid_period (float) – Period of the grid. Negative values indicate non-periodic grids

Returns:

A – Sparse interpolation matrix

Return type:

scipy.sparse.csc_matrix(len(outPoints), len(inPoints))