拉格朗日插值法的推导

1、插值的基本定义
  设函数$y=f(x)$在区间$[a,b]$上有定义，且已知它在$n+1$个互异点$a\le x_0<x_1<...<x_n\le b$上的函数值$y_0,y_1,...,y_n$，若存在一个简单函数$p(x)$，使得
$$p(x_i)=y_i,i=0,1,2,..,n$$
  成立，则称$p(x)$为$f(x)$插值函数。显然，除上述已知$n+1$个互异点外，在其他位置上，插值函数$p(x)$和原函数$f(x)$之间并没有明确关系，所以插值总是有误差的。不过，若对原函数和插值函数增加一定的约束，则可能使两者保持一致。下面讨论代数插值情况。
  设$p(x)$是次数不超过$n$的代数多项式，即
$$p(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+...+a_{1}x+a_0\text{\hspace{1em}(1)}$$
该函数满足$p(x_i)=y_i$，则称$p(x)$为原函数$f(x)$的$n$次代数插值多项式。该种插值称为代数插值，代数插值的几何意义，其实就是找一条过上述$n+1$个互异点的$n$次代数曲线来近似表示曲线$f(x)$。
  可以证明，上述$n$次插值函数有且只有一个。显然，将$p(x_i)=y_i$的条件代入(1)式，得到下面的$n+1$阶线性方程组
$$\{\begin{array}{c}{a}_0+a_1{x}_0+a_2{x}_0^2+...+a_n{x}_0^n=y_0\\ a_0+a_1{x}_1+a_2{x}_1^2+...+a_n{x}_1^n=y_1\\ ⋮\\ a_0+a_1{x}_n+a_2{x}_n^2+...+a_n{x}_n^n=y_n\end{array}\to \left[\begin{array}{cccc}1& x_0& ...& x_0^n\\ 1& x_1& ...& x_1^n\\ ⋮& ⋮& ...& ⋮\\ 1& x_n& ...& x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a^n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y^n\end{array}\right]$$
  显然，该线性方程组的行列式为
$$\left|\begin{array}{cccc}1& x_0& ...& x_0^n\\ 1& x_1& ...& x_1^n\\ ⋮& ⋮& ...& ⋮\\ 1& x_n& ...& x_n^n\end{array}\right|=\prod _{i=1}^n\prod _{j=0}^{i-1}(x_i-x_j)$$
 显然，由于$x_i$互不相同，所以上式不为0，所以方程系数$a_0,...,a_n$可被唯一确，即该插值多项式有且只有一个。
  插值问题的关键是求解插值多项式，显然利用上述线性方程组，可直接求得多项式系数的最小二乘解。但计算过程涉及矩阵求逆，计算量较大，后面将探究新的计算方法。

2、拉格朗日插值法

2.1、线性插值情况
  我们从一次多项式开始逐步推导。此时有$p(x_i)=y_i,(i=0,1)$，显然，可过这两个点作一条直线，目的是用直线$p(x)$来近似表示原函数$f(x)$，这种插值称为线性插值。该直线的两点式方程可表示为
$$p(x)=y_0\frac{x-x_1}{x_0-x_1}+y_1\frac{x-x_0}{x_1-x_0}=l_0(x)y_0+l_1(x)y_1=\sum _{k=0}^1{l}_k(x)y_k$$
  其中，$l_0(x)=\frac{x-x_1}{x_0-x_1}$称为点$x_0$的一次插值基函数, $l_1(x)=\frac{x-x_0}{x_1-x_0}$称为点$x_1$的一次插值基函数, 显然$p(x)$是$l_0(x)$和$l_1(x)$的线性组合。另外，上述插值基函数满足
$$l_j(x_i)={\delta }_{ij}=\{\begin{array}{c}1,i=j\\ 0,i\ne j\end{array}\text{\hspace{1em}(2)}$$

图1. 线性插值示意图

2.2、二次插值情况
  此时已知$f(x)$上面三个互异点$(x_i,y_i),i=0,1,2$，要求构造一个不超过2次的代数多项式$p(x)=a{x}^2+bx+c$，满足$p(x_i)=y_i,(i=0,1,2)$。为便于求解，我们可将$p(x)$重新整理为$p(x)=A(x-x_1)(x-x_2)+B(x-x_0)(x-x_2)+C(x-x_0)(x-x_1)$，将三个已知点坐标代入，可求得
$$\begin{gathered} A=\frac{y_0}{(x_0-x_1)(x_0-x_2)} \\ B=\frac{y_1}{(x_1-x_0)(x_1-x_2)} \\ C=\frac{y_2}{(x_2-x_0)(x_2-x_1)} \end{gathered}$$
此时可得到插值函数为
$$\begin{gathered} p(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}{y}_0+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}{y}_1+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}{y}_2 \\ =l_0(x)y_0+l_1(x)y_1+l_2(x)y_2=\sum _{j=0}^2{l}_j(x)y_j \end{gathered}$$
其中，$l_j(x)={\prod }_{i=0,i\ne j}^2\frac{x-x_i}{x_j-x_i}$，显然$l_j(x)$同样具有(2)式所示的性质。该插值称为二次插值或抛物线插值。

2.3、n次插值情况
  此时，仿照二次插值的构造方法，令
$$p(x)=l_0(x)y_0+...+l_n(x)y_n=\sum _{j=0}^n{l}_j(x)y_j$$
其中，$l_j(x)$是$n$次多项式，称为插值基函数，它满足条件
$$l_j(x_i)={\delta }_{ij}=\{\begin{array}{c}1,i=j\\ 0,i\ne j\end{array}，(i,j=0,1,...,n)\text{\hspace{1em}(2)}$$
  显然，此时问题归结为构造满足条件的$n$次多项式$l_j(x)$。事实上，由$l_j(x)=0,i\ne j$，知道$x_0,x_1,...,x_{j-1},x_{j+1},...,x_n$都是$l_j(x)$的零点，所以可设$l_j(x)=A(x-x_0)(x-x_1)....(x-x_{j-1})(x-x_{j+1})...(x-x_n)$。其中，$A$为待定系数，由条件$l_j(x_j)=1$，可确定$A=\frac{1}{(x_j-x_0)(x_j-x_1)...(x_j-x_{j-1})(x_j-x_{j+1})...(x_j-x_n)}$，所以
$$l_j(x)=\frac{(x-x_0)...(x-x_{j-1})(x-x_{j+1})...(x-x_n)}{(x_j-x_0)...(x_j-x_{j-1})(x_j-x_{j+1})...(x_j-x_n)}=\prod _{i=0,i\ne j}^n\frac{x-x_i}{x_j-x_i}$$
  由此可得
$$p(x)=\sum _{j=0}^n{l}_j(x)y_j=\sum _{j=0}^n{y}_j(\prod _{i=0,i\ne j}^n\frac{x-x_i}{x_j-x_i})\text{\hspace{1em}(3)}$$
我们称形如(3)式的插值公式为拉格朗日插值。

3、代码实现
  分别仿真拉格朗日插值法用于线性插值、抛物线插值和三次插值的情况，仿真代码见附录，仿真结果如下图所示，图中蓝色点表示已知点，红色点表示插值点。

图2. 线性插值的结果

图3. 抛物线插值的结果

图4. 三次插值的结果

---

【实现代码】

import numpy as np
import matplotlib.pyplot as plt

def lagrange_interpolate(x, xi, yi):
    """
    本函数用于实现拉格朗日函数
    :param x: 待计算的x坐标
    :param xi: 已知的x坐标
    :param yi: 已知的y坐标
    :return y: 对应x得到的y
    """
    if len(xi) != len(yi):
        raise ValueError('x dimension is not equal to y dimension!')

    # #######################计算拉格朗日基函数##################
    lj = np.ones((len(x), len(yi)))
    for j in range(0, len(yi), 1):
        for i in range(0, len(xi), 1):
            if i != j:
                lj[:, j] *= ((x - xi[i]) / (xi[j] - xi[i]))

    # ####################计算拉格朗日插值结果#####################
    y = np.zeros(len(x))
    for j in range(0, len(yi), 1):
        y += yi[j] * lj[:, j]
    return y


if __name__ == '__main__':
    xi = np.array([0, 6, 7, 11])
    yi = np.array([5, 2, 8, 9])
    x = np.arange(0, 12, 0.5)
    y = lagrange_interpolate(x, xi, yi)
    plt.figure()
    plt.plot(xi, yi, 'bo')
    plt.plot(x, y, 'r.-')
    plt.legend(['original points', 'interpolated points'])
    plt.grid()
    plt.show()