拉格朗日插值法的推导
拉格朗日插值法的推导
1、插值的基本定义
设函数 $y = f(x)$ 在区间 $[a, b]$ 上有定义，且已知它在 $n+1$ 个互异点 $a \leq x_0 < x_1 < \ldots < x_n \leq b$ 上的函数值 $y_0, y_1, \ldots, y_n$，若存在一个简单函数 $p(x)$，使得
$$
p(x_i) = y_i, \quad i = 0, 1, 2, \ldots, 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} + \ldots + a_1 x + a_0 \tag{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{cases}
a_0 + a_1 x_0 + a_2 x_0^2 + \ldots + a_n x_0^n = y_0 \
a_0 + a_1 x_1 + a_2 x_1^2 + \ldots + a_n x_1^n = y_1 \
\vdots \
a_0 + a_1 x_n + a_2 x_n^2 + \ldots + a_n x_n^n = y_n \
\end{cases}
\rightarrow
\begin{bmatrix}
1 & x_0 & \ldots & x_0^n \
1 & x_1 & \ldots & x_1^n \
\vdots & \vdots & \ldots & \vdots \
1 & x_n & \ldots & x_n^n \
\end{bmatrix}
\begin{bmatrix}
a_0 \
a_1 \
\vdots \
a^n \
\end{bmatrix}
\begin{bmatrix}
y_0 \
y_1 \
\vdots \
y^n \
\end{bmatrix}
$$
显然，该线性方程组的行列式为
$$
\begin{vmatrix}
1 & x_0 & \ldots & x_0^n \
1 & x_1 & \ldots & x_1^n \
\vdots & \vdots & \ldots & \vdots \
1 & x_n & \ldots & x_n^n \
\end{vmatrix}
= \prod_{i=1}^n \prod_{j=0}^{i-1} (x_i - x_j)
$$
显然，由于 $x_i$ 互不相同，所以上式不为0，所以方程系数 $a_0, \ldots, 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{cases}
1, & i = j \
0, & i \neq j
\end{cases} \tag{2}
$$
图1. 线性插值示意图
2.2、二次插值情况
此时已知 $f(x)$ 上面三个互异点 $(x_i, y_i), i = 0, 1, 2$，要求构造一个不超过2次的代数多项式 $p(x) = ax^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)$，将三个已知点坐标代入，可求得
$$
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)}
$$
此时可得到插值函数为
$$
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
$$
其中，$l_j(x) = \prod_{i=0, i \neq j}^2 \frac{x-x_i}{x_j-x_i}$，显然 $l_j(x)$ 同样具有(2)式所示的性质。该插值称为二次插值或抛物线插值。
2.3、n次插值情况
此时，仿照二次插值的构造方法，令
$$
p(x) = l_0(x)y_0 + \ldots + 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{cases}
1, & i = j \
0, & i \neq j
\end{cases} ，(i, j = 0, 1, \ldots, n) \tag{2}
$$
显然，此时问题归结为构造满足条件的 $n$ 次多项式 $l_j(x)$。事实上，由 $l_j(x) = 0, i \neq j$，知道 $x_0, x_1, \ldots, x_{j-1}, x_{j+1}, \ldots, x_n$ 都是 $l_j(x)$ 的零点，所以可设 $l_j(x) = A(x-x_0)(x-x_1)\ldots(x-x_{j-1})(x-x_{j+1})\ldots(x-x_n)$。其中，$A$ 为待定系数，由条件 $l_j(x_j) = 1$，可确定 $A = \frac{1}{(x_j-x_0)(x_j-x_1)\ldots(x_j-x_{j-1})(x_j-x_{j+1})\ldots(x_j-x_n)}$，所以
$$
l_j(x) = \frac{(x-x_0)\ldots(x-x_{j-1})(x-x_{j+1})\ldots(x-x_n)}{(x_j-x_0)\ldots(x_j-x_{j-1})(x_j-x_{j+1})\ldots(x_j-x_n)} = \prod_{i=0, i \neq 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 \left(\prod_{i=0, i \neq j}^n \frac{x-x_i}{x_j-x_i}\right) \tag{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()