cmlabs.interpolate.newtonbd

cmlabs.interpolate.newtonbd(xvals, x, yvals=None, coef=None)[source]

Newton’s backward interpolation formula.

\[\begin{split}\begin{gather} P_n(x) = \sum_{k=0}^{n} \frac{\nabla^k f(x_n)}{k!h^k} \omega_k(x) \\ P_n(x_n + th) = \sum_{k=0}^{n} C^k_t \nabla^k f(x_n) \end{gather}\end{split}\]

where \(\nabla^k f_n\) is the backward difference of the function.

Parameters:
  • xvals (array_like, 1-D) – The sorted x-coordinates of the data points.

  • x (float) – The x-coordinate at which to evaluate the polynomial.

  • yvals (array_like, 1-D, optional, default: None) – The y-coordinates of the data points, i.e., f(\(x\)). Only used if coef is not provided.

  • coef (array_like, 1-D, optional, default: None) – The backward differences list. If not provided, the function will compute the backward differences table using the backward_differences.

Returns:

res – The value of the polynomial at \(x\).

Return type:

float

Notes

The relation between divided difference and backward difference is proved by induction.

\[\begin{gather} f(x_0, x_1, \ldots, x_k) = \frac{\nabla^k f(x_n)}{k! h^k}, \quad h = x_i - x_{i-1} = const \end{gather}\]

So the Newton polynomial can be written as:

\[\begin{split}\begin{gather} P_n(x) = \sum_{k=0}^{n} \frac{\nabla^k f(x_n)}{k!h^k} \omega_k(x) \\ \end{gather}\end{split}\]

Using a variable \(s = \frac{x - x_n}{h}\) we can rewrite the polynomial as:

\[\begin{split}\begin{gather} P_n(x_n + th) = \sum_{k=0}^{n} \frac{t(t+1)\ldots(t+k-1)}{k!} \nabla^k f(x_n) \\ C^k_{-n} = \frac{-n(-n-1)\ldots(-n-k+1)}{k!} = (-1)^k \frac{n(n+1)\ldots(n+k-1)}{k!} \\ P_n(x) = \sum_{k=0}^{n} (-1)^k C^k_{-t} \nabla^k f(x_n) \\ \end{gather}\end{split}\]

This method is beneficial for determining values of \(x\) at the end of the tabulated equally spaced set of values, i.e. \(-1 < t < 0\). Also for extrapolating values of \(x\) a little ahead (to the right) of \(x_n\).

Examples

>>> import numpy as np
>>> from cmlabs.interpolate import newtonbd
>>> xvals = np.array([0, 1, 2, 3])
>>> yvals = np.array([1, 3, 2, 5])
>>> x = np.float32(2.5)
>>> newtonbd(xvals, x, yvals=yvals)
2.5625