拉格朗日插值法的推导

拉格朗日插值法的推导

引用

CSDN

1.

https://blog.csdn.net/u011572784/article/details/139280319

1、插值的基本定义

设函数y = f ( x ) 在区间[ a , b ] 上有定义，且已知它在n + 1个互异点a ≤ x 0 < x 1 < . . . < x n ≤ 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 (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阶线性方程组
{ 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 → [ 1 x 0 . . . x 0 n 1 x 1 . . . x 1 n ⋮ ⋮ . . . ⋮ 1 x n . . . x n n ] [ a 0 a 1 ⋮ a n ] = [ y 0 y 1 ⋮ y n ]
显然，该线性方程组的行列式为
∣ 1 x 0 . . . x 0 n 1 x 1 . . . x 1 n ⋮ ⋮ . . . ⋮ 1 x n . . . x n n ∣ = ∏ i = 1 n ∏ 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 x − x 1 x 0 − x 1 + y 1 x − x 0 x 1 − x 0 = l 0 ( x ) y 0 + l 1 ( x ) y 1 = ∑ k = 0 1 l k ( x ) y k
其中，l 0 ( x ) = x − x 1 x 0 − x 1称为点x 0的一次插值基函数,l 1 ( x ) = x − x 0 x 1 − x 0称为点x 1的一次插值基函数, 显然p ( x )是l 0 ( x )和l 1 ( x )的线性组合。另外，上述插值基函数满足
l j ( x i ) = δ i j = { 1 , i = j 0 , i ≠ j (2)

图1. 线性插值示意图

2.2、二次插值情况

此时已知f ( x )上面三个互异点( x i , y i ) , i = 0 , 1 , 2，要求构造一个不超过2次的代数多项式p ( x ) = a x 2 + b x + 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 )，将三个已知点坐标代入，可求得
A = y 0 ( x 0 − x 1 ) ( x 0 − x 2 ) B = y 1 ( x 1 − x 0 ) ( x 1 − x 2 ) C = y 2 ( x 2 − x 0 ) ( x 2 − x 1 )
此时可得到插值函数为
p ( x ) = ( x − x 1 ) ( x − x 2 ) ( x 0 − x 1 ) ( x 0 − x 2 ) y 0 + ( x − x 0 ) ( x − x 2 ) ( x 1 − x 0 ) ( x 1 − x 2 ) y 1 + ( 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 = ∑ j = 0 2 l j ( x ) y j
其中，l j ( x ) = ∏ i = 0 , i ≠ j 2 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 = ∑ j = 0 n l j ( x ) y j
其中，l j ( x )是n次多项式，称为插值基函数，它满足条件
l j ( x i ) = δ i j = { 1 , i = j 0 , i ≠ j ， ( i , j = 0 , 1 , . . . , n ) (2)

显然，此时问题归结为构造满足条件的n次多项式l j ( x )。事实上，由l j ( x ) = 0 , i ≠ 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 = 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 ) = ( 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 ) = ∏ i = 0 , i ≠ j n x − x i x j − x i

由此可得
p ( x ) = ∑ j = 0 n l j ( x ) y j = ∑ j = 0 n y j ( ∏ i = 0 , i ≠ j n x − x i x j − x i ) (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()