cmlabs.interpolate.newtonfd

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

Newton’s forward interpolation formula.

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

where \(\Delta^k f_0\) is the forward 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 forward differences list. If not provided, the function will compute the forward differences table using the forward_differences.

Returns:

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

Return type:

float

See also

forward_differences, newton

Notes

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

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

So the Newton polynomial can be written as:

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

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

\[\begin{split}\begin{gather} P_n(x_0 + th) = \sum_{k=0}^{n} \frac{t(t-1)\ldots(t-k+1)}{k!} \Delta^k f(x_0) \\ C^k_n = \frac{n(n-1)\ldots(n-k+1)}{k!}, \quad k \geq 2, \quad C^0_n = 1, \quad C^1_n = n \\ P_n(x_0 + th) = \sum_{k=0}^{n} C^k_t \Delta^k f(x_0) \\ \end{gather}\end{split}\]

This method is beneficial for determining values of \(x\) at the beginning of a tabulated equally spaced set of values, i.e. \(0 < t < 1\). Also for extrapolating values of \(x\) a bit backward (i.e. to the left) of \(x_0\).

Examples

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