Python中的拉格朗日插值

Question

我想使用拉格朗日插值法来插值一个多项式，但是这段代码无法正常工作：

def interpolate(x_values, y_values):
    def _basis(j):
        p = [(x - x_values[m])/(x_values[j] - x_values[m]) for m in xrange(k + 1) if m != j]
        return reduce(operator.mul, p)

    assert len(x_values) != 0 and (len(x_values) == len(y_values)), 'x and y cannot be empty and must have the same length'

    k = len(x_values)
    return sum(_basis(j) for j in xrange(k))

我跟随维基百科的方法，但当我运行它时，在第三行收到了一个IndexError错误!

谢谢

- rubik

---

如果那些点踩的人能解释一下为什么要点踩，我会非常感激。 - rubik

Answer 1

尝试

def interpolate(x, x_values, y_values):
    def _basis(j):
        p = [(x - x_values[m])/(x_values[j] - x_values[m]) for m in xrange(k) if m != j]
        return reduce(operator.mul, p)
    assert len(x_values) != 0 and (len(x_values) == len(y_values)), 'x and y cannot be empty and must have the same length'
    k = len(x_values)
    return sum(_basis(j)*y_values[j] for j in xrange(k))

你可以按照以下步骤进行确认:

>>> interpolate(1,[1,2,4],[1,0,2])
1.0
>>> interpolate(2,[1,2,4],[1,0,2])
0.0
>>> interpolate(4,[1,2,4],[1,0,2])
2.0
>>> interpolate(3,[1,2,4],[1,0,2])
0.33333333333333331

因此，结果是根据给定点的多项式插值值。在这种情况下，3个点定义了一个抛物线，前3个测试显示给定x_value时返回了指定的y_value。

- smichr

---

要在Python3中使其工作：pip3 install future，然后您可以添加以下代码：from past.builtins import reduce, xrange; import operator - SurpriseDog

Answer 2

虽然我很晚才发现这个，但是我在寻找Lagrange插值的简单实现时找到了这个。@smichr的回答非常好，但是Python有点过时，我也想要一些能够良好使用np.ndarrays并且易于绘图的东西。也许其他人会发现这个有用：

import numpy as np
import matplotlib.pyplot as plt


class LagrangePoly:

    def __init__(self, X, Y):
        self.n = len(X)
        self.X = np.array(X)
        self.Y = np.array(Y)

    def basis(self, x, j):
        b = [(x - self.X[m]) / (self.X[j] - self.X[m])
             for m in range(self.n) if m != j]
        return np.prod(b, axis=0) * self.Y[j]

    def interpolate(self, x):
        b = [self.basis(x, j) for j in range(self.n)]
        return np.sum(b, axis=0)


X  = [-9, -4, -1, 7]
Y  = [5, 2, -2, 9]

plt.scatter(X, Y, c='k')

lp = LagrangePoly(X, Y)

xx = np.arange(-100, 100) / 10

plt.plot(xx, lp.basis(xx, 0))
plt.plot(xx, lp.basis(xx, 1))
plt.plot(xx, lp.basis(xx, 2))
plt.plot(xx, lp.basis(xx, 3))
plt.plot(xx, lp.interpolate(xx), linestyle=':')
plt.show()

- jds

---

有人能解释一下 xx vs lp.basis(xx,0),(xx,1).... (xx,3) 的图形代表什么吗？ - Devansh Mishra

它正在绘制拉格朗日多项式的基函数，重建一个类似于这样的图形：https://en.wikipedia.org/wiki/File:Lagrange_polynomial.svg - jds

Answer 3

检查一下索引，维基百科上说“k+1个数据点”，但是你设置的是k = len(x_values)，如果你完全按照公式进行操作，那么应该是k = len(x_values) - 1。

- AndiDog

---

好的，如果我这样做：interpolate（[1,2,3]，[1,4,9]）是为什么会返回 -0.5x^2 + 1.5x？看看这个链接：http://i.imgur.com/MkATz.gif - rubik

@rubik：抱歉，如果不知道插值算法，我无法帮助您解决这样具体的问题（而且我也不会去查阅相关资料）。请再次检查您的逻辑或搜索现有的实现。如果您在如何应用插值方面发布更多代码（例如，在您的问题中缺少“x”的定义/初始值），那么可能会有人能够为您提供进一步的帮助。 - AndiDog

我正在使用pypol（http://pypol.altervista.org/），x是单项式（x=1）（http://pypol.altervista.org/functions.html#pypol.monomial）。 - rubik

Answer 4

这段代码与Python 3兼容：

def Lagrange (Lx, Ly):
    x=sympy.symbols('x')
    if  len(Lx)!= len(Ly):
        return 1
    y=0
    for k in range ( len(Lx) ):
        t=1
        for j in range ( len(Lx) ):
            if j != k:
                t=t* ( (x-Lx[j]) /(Lx[k]-Lx[j]) )
        y+= t*Ly[k]
    return y

- nima ansari