scipy.interpolate.

BarycentricInterpolator

class scipy.interpolate.BarycentricInterpolator(xi, yi=None, axis=0, *, wi=None, rng=None)

重心（Lagrange 插值，稳定性改进）插值器（C∞ 光滑）。

构建一个通过给定点集的 बहुपद。允许评估多项式及其所有导数，高效更改要插值的 y 值，并通过添加更多 x 和 y 值进行更新。为了数值稳定性，采用重心表示而不是直接计算多项式系数。

参数：

xiarray_like, 形状 (npoints, )

多项式应通过的点的 x 坐标的一维数组

yiarray_like, 形状 (…, npoints, …), 可选

多项式应通过的点的 y 坐标的 N 维数组。如果为 None，y 值将稍后通过 set_y 方法提供。yi 沿插值轴的长度必须等于 xi 的长度。使用 axis 参数选择正确的轴。

axisint, 可选

yi 数组中与 x 坐标值对应的轴。默认为 axis=0。

wiarray_like, 可选

所选插值点 xi 的重心权重。如果缺失或为 None，将从 xi 计算权重（默认）。如果使用相同的节点 xi 计算多个插值，这允许重用权重 wi，而无需重新计算。这也允许对 xi 的某些选择显式计算权重（参见注释）。

rng{None, int, numpy.random.Generator}, 可选

如果 rng 以关键字形式传递，则除了 numpy.random.Generator 之外的类型将传递给 numpy.random.default_rng 以实例化 Generator。如果 rng 已经是 Generator 实例，则使用提供的实例。指定 rng 以实现可重复插值。

如果此参数 random_state 以关键字形式传递，则应用参数 random_state 的传统行为。

• 如果 random_state 为 None（或 numpy.random），则使用 numpy.random.RandomState 单例。

• 如果 random_state 是一个 int，则使用一个新的 RandomState 实例，以 random_state 为种子。

• 如果 random_state 已经是 Generator 或 RandomState 实例，则使用该实例。

版本 1.15.0 中有改动： 作为 SPEC-007 从使用 numpy.random.RandomState 过渡到 numpy.random.Generator 的一部分，此关键字从 random_state 更改为 rng。在过渡期间，两个关键字都将继续工作（只指定其中一个）。过渡期结束后，使用 random_state 关键字将发出警告。random_state 和 rng 关键字的行为如上所述。

属性：

dtype

方法

__call__(x)	在点 x 处评估插值多项式
add_xi(xi[, yi])	向要插值的集合添加更多 x 值
derivative(x[, der])	在点 x 处评估多项式的一个导数。
derivatives(x[, der])	在点 x 处评估多项式的多个导数
set_yi(yi[, axis])	更新要插值的 y 值

注释

此方法是基于 [2] 的 Lagrange 多项式插值 [1] 的变体。该多项式不是使用 Lagrange 或 Newton 公式，而是通过重心公式表示

\[p(x) = \frac{\sum_{i=1}^m\ w_i y_i / (x - x_i)}{\sum_{i=1}^m w_i / (x - x_i)},\]

其中 \(w_i\) 是用一般公式计算的重心权重

\[w_i = \left( \prod_{k \neq i} x_i - x_k \right)^{-1}.\]

这与 AAA 和 FloaterHormannInterpolator 使用的重心形式相同。然而，与此相反，权重 \(w_i\) 的定义使得 \(p(x)\) 是一个多项式而不是一个有理函数。

重心表示避免了许多由浮点运算引起的多项式插值问题。然而，它并不能避免多项式插值固有的问题。即，如果 x 坐标等距，则可以使用 [2] 中的公式显式计算权重

\[w_i = (-1)^i {n \choose i},\]

其中 \(n\) 是 x 坐标的数量。正如 [2] 中指出的，这意味着对于大的 \(n\)，权重会以指数级大的因子变化，导致龙格现象（Runge phenomenon）。

为了避免这种病态，x 坐标应集中在区间的端点。在区间 \([a,b]\) 上，切比雪夫第二类点是绝佳的选择

\[x_i = \frac{a + b}{2} + \frac{a - b}{2}\cos(i\pi/n).\]

在这种情况下，权重可以显式计算为

\[\begin{split}w_i = \begin{cases} (-1)^i/2 & i = 0,n \\ (-1)^i & \text{otherwise} \end{cases}.\end{split}\]

更多信息请参见 [2]。请注意，对于大的 \(n\)，显式计算权重（参见示例）将比通用公式更快。

参考文献

[1]

https://en.wikipedia.org/wiki/Lagrange_polynomial

[2] (1,2,3,4)

Jean-Paul Berrut 和 Lloyd N. Trefethen，《Barycentric Lagrange Interpolation》，SIAM Review 2004 46:3, 501-517 DOI:10.1137/S0036144502417715

示例

在浏览器中试用！

生成一个五次重心插值函数，使用 \((0, \pi/2)\) 中六个随机间隔的节点来近似函数 \(\sin x\) 及其前四个导数。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import BarycentricInterpolator
>>> rng = np.random.default_rng()
>>> xi = rng.random(6) * np.pi/2
>>> f, f_d1, f_d2, f_d3, f_d4 = np.sin, np.cos, lambda x: -np.sin(x), lambda x: -np.cos(x), np.sin
>>> P = BarycentricInterpolator(xi, f(xi), random_state=rng)
>>> fig, axs = plt.subplots(5, 1, sharex=True, layout='constrained', figsize=(7,10))
>>> x = np.linspace(0, np.pi, 100)
>>> axs[0].plot(x, P(x), 'r:', x, f(x), 'k--', xi, f(xi), 'xk')
>>> axs[1].plot(x, P.derivative(x), 'r:', x, f_d1(x), 'k--', xi, f_d1(xi), 'xk')
>>> axs[2].plot(x, P.derivative(x, 2), 'r:', x, f_d2(x), 'k--', xi, f_d2(xi), 'xk')
>>> axs[3].plot(x, P.derivative(x, 3), 'r:', x, f_d3(x), 'k--', xi, f_d3(xi), 'xk')
>>> axs[4].plot(x, P.derivative(x, 4), 'r:', x, f_d4(x), 'k--', xi, f_d4(xi), 'xk')
>>> axs[0].set_xlim(0, np.pi)
>>> axs[4].set_xlabel(r"$x$")
>>> axs[4].set_xticks([i * np.pi / 4 for i in range(5)],
...                   ["0", r"$\frac{\pi}{4}$", r"$\frac{\pi}{2}$", r"$\frac{3\pi}{4}$", r"$\pi$"])
>>> for ax, label in zip(axs, ("$f(x)$", "$f'(x)$", "$f''(x)$", "$f^{(3)}(x)$", "$f^{(4)}(x)$")):
...     ax.set_ylabel(label)
>>> labels = ['Interpolation nodes', 'True function $f$', 'Barycentric interpolation']
>>> axs[0].legend(axs[0].get_lines()[::-1], labels, bbox_to_anchor=(0., 1.02, 1., .102),
...               loc='lower left', ncols=3, mode="expand", borderaxespad=0., frameon=False)
>>> plt.show()

接下来，我们将展示如何使用切比雪夫第二类点来避免龙格现象。在此示例中，我们还显式计算了权重。

>>> n = 20
>>> def f(x): return np.abs(x) + 0.5*x - x**2
>>> i = np.arange(n)
>>> x_cheb = np.cos(i*np.pi/(n - 1))  # Chebyshev points on [-1, 1]
>>> w_i_cheb = (-1.)**i  # Explicit formula for weights of Chebyshev points
>>> w_i_cheb[[0, -1]] /= 2
>>> p_cheb = BarycentricInterpolator(x_cheb, f(x_cheb), wi=w_i_cheb)
>>> x_equi = np.linspace(-1, 1, n)
>>> p_equi = BarycentricInterpolator(x_equi, f(x_equi))
>>> xx = np.linspace(-1, 1, 1000)
>>> fig, ax = plt.subplots()
>>> ax.plot(xx, f(xx), label="Original Function")
>>> ax.plot(xx, p_cheb(xx), "--", label="Chebshev Points")
>>> ax.plot(xx, p_equi(xx), "--", label="Equally Spaced Points")
>>> ax.set(xlabel="$x$", ylabel="$f(x)$", xlim=[-1, 1])
>>> ax.legend()
>>> plt.show()

返回在新标签页中打开