拉格朗日插值原理及实现(Python)
一. 前言
Lagrange插值是利用n次多项式来拟合(n+1)个数据点从而得到插值函数的方法。(注意n次多项式的定义是未知数最高次幂为n,但是多项式系数有n+1个,因为还有个常数项)
Lagrange插值和Newton插值本质上相同,都是用(n-1)次多项式来拟合n个数据点。所以这两种插值方法得到的插值函数相同,因为多项式拟合的基本定理:同时通过n个数据点,且最高次幂小于(n-1)的多项式函数唯一。下面顺手证明一下这个重要的定理。
如果已知n+1个数据点\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\),假设\(L_1=k_0+k_1x+k_2x^2+\cdots+k_nx^n\)和\(L_2=k_0’+k_1’x+k_2'x^2+\cdots+k_n’x^n\)都是通过这n个数据点的插值函数。那么应该有\(L_1 - L_2\)通过所有\((x_1,0),(x_2,0),\cdots,(x_n,0)\)
代入这些点得到齐次线性方程组:
\[\begin{bmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^n\\
1 & x_1 & x_1^2 & \cdots & x_1^n\\
\vdots & \vdots & \vdots & & \vdots \\
1 & x_n & x_n^2 & \cdots & x_n^n\\
\end{bmatrix}
\begin{bmatrix}
k_0-k_0'\\
k_1-k_1'\\
k_2-k_2'\\
\vdots \\
k_n-k_n'\\
\end{bmatrix}
=
\begin{bmatrix}
0 \\
0 \\
0 \\
\vdots \\
0
\end{bmatrix}
\]它的系数行列式是Vandermonde行列式,所以:
\[\begin{vmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^n\\
1 & x_1 & x_1^2 & \cdots & x_1^n\\
\vdots & \vdots & \vdots & & \vdots \\
1 & x_n & x_n^2 & \cdots & x_n^n\\
\end{vmatrix}
=
\prod_{n\ge i > j \ge 0 }(x_i-x_j)
\]因为每个点都是不同的,所以\(x_i \ne x_j,i\ne j\),所以齐次线性方程组的系数行列式不等于0,故方程解唯一且为0解:
\[k_i - k_i' = 0\\
k_i = k_i' \\
L_1 = L_2
\]证毕。
二. 3种形式的Lagrange插值函数推导
1. 原始形态的Lagrange插值
为了用n次多项式拟合n+1个数据点:\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)
Lagrange插值函数采用的方法是构建一个这样的函数:
\]
也就是用一组基函数\(\{l_0(x),l_1{x}, \cdots ,l_n(x)\}\)去构建插值函数\(L(x)\),那么不难想到这样的基函数需要满足这样的条件:
0 , & j\ne i \\
1 , & j = i
\end{cases}
\]
这样对于\(L(x)\)就会有:
L(x_i) &= l_0(x_i)y_0 + \cdots l_i(x_i)y_i + \cdots +l_n(x)y_n \\
&=y_i
\end{split}
\]
这样就实现了\(L(x)\)通过所有的数据点\((x_i,y_i)\)。接下来就构建这样的基函数\(l_i(x)\)。
首先实现\(l_i(x_j) = 0,i\ne j\) :
\cdots(x-x_n)
\]
这个函数实现了\(l_i(x)\)在所有非\((x_i,y_i)\)的点处为0,但是在\(x=x_i\)处:
\]
为了让\(l_i(x_i)=1\),我们可以将\(l_i(x)\)除以这个值进行修正:
l_i(x_i)&= \frac{(x-x_0)(x -x_1)\cdots (x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}\\
&=\prod_{n\ge j \ge0,j\ne i} \frac{x-x_j}{x_i -x_j}
\end{split}
\tag{2}
\]
将(2)代入(1)就得到了原始形态的Lagrange插值函数.
L(x) &= l_0(x)y_0 + l_1(x)y_1 + \cdots +l_n(x)y_n \\
&= \prod_{n\ge j \ge0,j\ne 0} \frac{x-x_j}{x_0 -x_j}y_0 +\prod_{n\ge j \ge0,j\ne 1} \frac{x-x_j}{x_1 -x_j}y_1 +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{x-x_j}{x_n -x_j}y_n
\end{split}
\tag{3}
\]
例:已知点(1,1),(2,2),(3,3),用原始Lagrange插值计算插值函数。
首先计算三个插值基函数:
\[l_0(x) = \frac{(x-2)(x-3)}{(1-2)(1-3)}=\frac{(x-2)(x-3)}{2}\\
l_1(x) = \frac{(x-1)(x-3)}{(2-1)(2-3)}=-{(x-1)(x-3)}\\
l_2(x) = \frac{(x-1)(x-2)}{(3-1)(3-2)}=\frac{(x-1)(x-2)}{2}\\
\]从而得到插值函数:
\[L(X) = l_0(x)+ 2l_1(x)+3l_2(x)=\frac{(x-2)(x-3)}{2}-2{(x-1)(x-3)}+ \frac{3(x-1)(x-2)}{2}
\]插入点x=4试一下:\(L(4)=4\) 。
原始模式的Lagrange插值函数,每次计算插值点时需要进行n(n-1)次乘法,时间复杂度为\(O(n^2)\)
2. 第一形式Lagrange插值
为了降低计算的时间复杂度,我们对原始的Lagrange插值函数进行改进。
为了书写方便,我们令\(w_i\)为
w_i &= \prod_{n\ge j \ge0,j\ne i} (x_i -x_j) \\
&= (x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)
\end{split}
\tag{4}
\]
观察(3)式,我们可以提取一个公因数\((x-x_0)(x-x_1)\cdots (x-x_n) = g(x)\):
L(x)
&= (x-x_0)(x-x_1)\cdots (x-x_n)[\prod_{n\ge j \ge0,j\ne 0} \frac{y_0}{(x_0 -x_j)(x-x_0)} +\prod_{n\ge j \ge0,j\ne 1} \frac{y_1}{(x_1 -x_j)(x-x_1)} +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{y_n}{(x_n -x_j)(x-x_n)}]\\
&=g(x)[\frac{y_0}{w_0(x-x_0)}+\frac{y_1}{w_1(x-x_1)}+\cdots+\frac{y_n}{w_n(x-x_n)}]\\
&= g(x)\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)}
\end{split}
\tag{5}
\]
这样就得到了第一形式Lagrange插值。这一形式的Lagrange插值计算的流程如下:
- 预处理:根据已知数据点,利用公式(4)计算\(w_0,w_1,\cdots,w_n\);
- 插值:计算\(g(x)=(x-x_0)(x-x_1)\cdots(x-x_n)\),然后根据公式(5),先计算中括号内的加式,然后乘以\(g(x)\)得到插值点的值,这样计算一个新的点的时间复杂度就是\(O(n)\);
- 补充数据点:如果数据集有更新,只需要更新\(w_i\)即可,加入一个新的点到数据集,只需要将每个\(w_i\)乘以\((x_i-x_{n+1})\),此外再增加一个\(w_{n+1}\)。
例:已知点(1,1),(2,2),(3,3),使用Lagrange第一形式计算插值函数。
首先计算\(w_i\):
\[w_0=(1-2)(1-3)=2\\
w_1=(2-1)(2-3)=-1\\
w_2=(3-1)(3-2)=2
\]插值函数就是:
\[L(x)=(x-1)(x-2)(x-3)[\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}]
\]插入点x=4试一下:\(L(4)=4\) 。
3. 第二形式的Lagrange插值(重心插值公式)
第一形式的Lagrange插值还要计算\(g(x)\),可以再进行优化。
第一形式:
\]
为了消掉\(g(x)\),我们取\(y=1\)这条直线上的点进行插值,即取n+1个点:\((x_0,1),(x_1,1),\cdots,(x_n,1)\)那么这n+1个点的插值函数就是:
\]
而\(L'(x)=1,x=x_0,x_1,\cdots,x_n\) ,所以我们可以用(6)除以(7)消去g(x):
\]
这样就得到了第二形式的Lagrange插值,也称为重心插值,通常Lagrange插值采用这种形式。
它的计算过程如下:
- 预处理:计算\(w_i\);
- 插值:对需要插入的点(x,y),计算(x-x_i)(可以存到一个列表中,计算时直接取用);
- 补充数据点:补充新的点到数据集只需要更新\(w_i\) 。
例:已知点(1,1),(2,2),(3,3),利用重心插值公式计算插值函数。
首先计算\(w_i\):
\[w_0=(1-2)(1-3)=2\\
w_1=(2-1)(2-3)=-1\\
w_2=(3-1)(3-2)=2
\]得到\(L(x)\):
\[L(x)=\frac{\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}}{\frac{1}{2(x-1)}-\frac{1}{(x-2)}+\frac{1}{2(x-3)}}
\]插入点x=4试一下:\(L(4)=4\) 。
三. 利用Python编程实现这三种Lagrange插值
import numpy as np
class LagrangeInterpolation:
"""
There are three modes of the LagrangeInterpolation.The input data is supposed
to be two lists.
---------------------------------------------------------------------------------------
self.data : Contain the points known before.
self.dataLength : Indicate the number of points in the data list.
self.weight : The self.weight[i] is (xi - x1)(xi - x2)...(xi - x(i-1))(xi - x(i+1))
...(xi - xn)
self.items : The self.items[i] is (x - x1)(x - x2)...(x - x(i-1))(x - x(i+1))...
(x - xn)
---------------------------------------------------------------------------------------
"""
def __init__(self, x, y):
self.data = {'x': list(x), 'y': list(y)}
self.dataLength = len(self.data['x'])
self.weight = []
self.items = []
# control is a flag indicating if there is anything wrong with the data.
self.control = True
if len(self.data['x']) != len(self.data['y']):
print("The length of x isn't equal to the length of y!")
self.control = False
else:
self.__preprocess(order=0)
# Appending function is used to add points to the data list.
def data_append(self, ap):
if ap[0] in self.data['x']:
print("The point already exist.")
self.control = False
else:
self.control = True
self.data['x'].append(ap[0])
self.data['y'].append(ap[1])
self.dataLength = self.dataLength + 1
self.__preprocess(order=1)
"""
Preprocessing is used to update the self.weight and self.items.
order = 0 : Initialize the self.weight.
order = 1 : Update the self.weight when new point is added to the data list.
order = 2 : Calculate the self.items for each point waiting for interpolation.
"""
def __preprocess(self, order=0, x=0):
if order == 0:
self.weight = list(np.zeros(self.dataLength))
for i in range(self.dataLength):
weight_temp = 1
for j in range(self.dataLength):
if i == j:
pass
else:
weight_temp = weight_temp * (self.data['x'][i] - self.data['x'][j])
self.weight[i] = weight_temp
elif order == 1:
self.weight.append(1)
for i in range(self.dataLength - 1):
self.weight[i] = self.weight[i] * (self.data['x'][i] - self.data['x'][-1])
self.weight[-1] = self.weight[-1] * (self.data['x'][-1] - self.data['x'][i])
elif order == 2:
self.items = list(np.zeros(self.dataLength))
for i in range(self.dataLength):
self.items[i] = x - self.data['x'][i]
# The mode1 is the initial mode of Lagrange interpolation.
def mode1(self, px):
self.__preprocess(order=2, x=px)
if self.control:
dataCheck = False
for w in self.weight:
if w != 0:
dataCheck = True
else:
dataCheck = False
if dataCheck:
py = 0.0
for i in range(self.dataLength):
py_temp = 1
for j in range(self.dataLength):
if i == j:
pass
else:
py_temp = py_temp * self.items[j]
py = py + py_temp * self.data['y'][i] / self.weight[i]
return py
else:
print("There is a same x!")
return None
else:
return None
def mode2(self, px):
self.__preprocess(order=2, x=px)
itemsProd = np.prod(self.items)
itemsSum = 0
for i in range(self.dataLength):
itemsSum = itemsSum + self.data['y'][i] / (self.weight[i] * self.items[i])
py = itemsProd * itemsSum
return py
def mode3(self, px):
self.__preprocess(order=2, x=px)
denomSum = 0
numeSum = 0
for i in range(self.dataLength):
dtemp = self.weight[i] * self.items[i]
numeSum = numeSum + self.data['y'][i] / dtemp
denomSum = denomSum + 1 / dtemp
py = numeSum / denomSum
return py
demo.py :
from lagrange_interpolation import LagrangeInterpolation as lag
x = [1, 2, 3, 4]
y = [1, 2, 3, 4]
inter1 = lag(x, y)
inter1.data_append((5, 5)) #往数据集中追加一个点
z1 = inter1.mode1(6)
z2 = inter1.mode2(6)
z3 = inter1.mode3(6)
print(z1)
print(z2)
print(z3)
"""
Result:
6.0
6.000000000000002
6.0000000000000036
"""
本站文章如无特殊说明,均为本站原创,如若转载,请注明出处:拉格朗日插值原理及实现(Python) - Python技术站