数值分析-多项式插值方法小结
- 前言
- 插值的应用与唯一性
- Lagrange插值法和逐次线性插值
* - 代码实现
- 逐次线性插值
* - 代码实现
- 逼近复杂函数
- Newton插值法
* - 代码实践
- 差分与等距节点差分插值
* - 代码实现
- Hermite插值方法
- 写在最后
前言
最近在学数值分析，把一些算法实现一下
具体内容详见 《数值分析》 李庆扬、王能超、易大义（华中科大出版社）
插值的应用与唯一性
-
插值的应用：
-
对于一些复杂函数，我们往往只能知道某些具体位置的函数值（ 实际中数据还有误差），当我们需要得到全体定义域上任一点的数值时，便可以采用插值的方法
-
可以用多项式插值函数来拟合复杂的数据点，简化函数表达 。（ 多项式历来最被研究者喜欢，可积可导，性质很好！）
-
对于给定的n个点，所有的x i ≠ x j x_i \neq x_j x i =x j ，即所有的采样点均不同， 存在且唯一的次数不超过n次的多项式经过所有采样点，称之为 插值多项式。（ 思路：组成的n阶线性方程组的系数矩阵行列式为范德蒙行列式，行列式不为零，所以存在唯一解。）
Lagrange插值法和逐次线性插值
- Lagrange线性插值
- Lagrange插值的思想是不去求解各x幂次的系数，直接把函数值y作为系数确定下来，反过来确定线性空间的基的形式！
-
考虑两点的线性插值 已知f ( x k ) = y k , f ( x k + 1 ) = y k + 1 f(x_k)=y_k, f(x_{k+1})=y_{k+1}f (x k )=y k ,f (x k +1 )=y k +1 ，设两点的一阶插值多项式为L 1 ( x ) L_1(x)L 1 (x )，有L 1 ( x k ) = y k , L 1 ( x k + 1 ) = y k + 1 L_1(x_k)=y_k, L_1(x_{k+1})=y_{k+1}L 1 (x k )=y k ,L 1 (x k +1 )=y k +1 。 注意插值多项式一定经过对应的插值节点 按照Lagrange插值的思想（由直线的点斜式和两点式可以简单推来），我们可以得到： L 1 ( x ) = x − x k + 1 x k − x k + 1 y k + x − x k x k + 1 − x k y k + 1 \Large L_1(x)=\frac{x-x_{k+1}}{x_k-x_{k+1}}y_k+\frac{x-x_{k}}{x_{k+1}-x_k}y_{k+1}L 1 (x )=x k −x k +1 x −x k +1 y k +x k +1 −x k x −x k y k +1 其实可以利用插值多项式的零点和抽样delta的性质来得到多点Lagrange插值多项式的形式：
- 零点条件：对于l k ( x ) l_k(x)l k (x )，在除第k个采样点的点上始终为0
- 抽样性质：对于l k ( x ) l_k(x)l k (x )，在第k个采样点的点上为1
-
给出l k ( x ) l_k(x)l k (x )的一般形式： l k ( x ) = ( x − x 0 ) . . . ( x − x k − 1 ) ( x − x k + 1 ) . . . ( x − x n ) ( x k − x 0 ) . . . ( x k − x k − 1 ) ( x k − x k + 1 ) . . . ( x k − x n ) ( k = 0 , 1 , . . . , n ) \Large l_k(x)=\frac{(x-x_0)…(x-x_{k-1})(x-x_{k+1})…(x-x_n)}{(x_k-x_0)…(x_k-x_{k-1})(x_k-x_{k+1})…(x_k-x_n)} \,\,\,(k=0,1,…,n)l k (x )=(x k −x 0 )…(x k −x k −1 )(x k −x k +1 )…(x k −x n )(x −x 0 )…(x −x k −1 )(x −x k +1 )…(x −x n )(k =0 ,1 ,…,n )
-
Lagrange插值多项式表达式： L n ( x ) = ∑ k = 0 n y k l k ( x ) = ∑ i = 0 n y i ∏ j = 0 , j ≠ i n x − x j x i − x j \Large L_n(x)=\sum_{k=0}^ny_kl_k(x)=\sum_{i=0}^n y_i \prod_{j=0,j \neq i}^n \frac{x-x_j}{x_i-x_j}L n (x )=∑k =0 n y k l k (x )=∑i =0 n y i ∏j =0 ,j =i n x i −x j x −x j 其中，我们可以看到Lagrange插值多项式的基函数为： b a s e n ( x ) = l k ( x ) = ∏ j = 0 , j ≠ k n x − x j x k − x j base_n(x)=l_k(x)=\prod_{j=0,j \neq k}^n \frac{x-x_j}{x_k-x_j}b a s e n (x )=l k (x )=j =0 ,j =k ∏n x k −x j x −x j
-
插值余项
定义截断误差（插值余项）为R n ( x ) = f ( x ) − L n ( x ) R_n(x)=f(x)-L_n(x)R n (x )=f (x )−L n (x )
根据罗尔定理，可以推得： R n ( x ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ∏ i = 0 n ( x − x i ) \large R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x-x_i)R n (x )=(n +1 )!f (n +1 )(ξ)∏i =0 n (x −x i )
代码实现
import numpy as np
from functools import reduce
class Lagrange_Polynomial_Interpolation(object):
"""Lagrange_Polynomial_Interpolation
base function: Ln(x)=\sum_{k=0}^{n} yk l_k(x)
Author: Junno
Date: 2022-03-01
"""
def __init__(self, x, y):
"""__init__
Args:
x ([np.array]): [x of points]
y ([np.array]): [y of points]
"""
self.N = x.shape[0]
self.x = np.copy(x)
self.y = np.copy(y)
def calculate_lk(self, x, k,):
"""calculate_lk(x)
base function: lk(x)=∏——{j=0,j!=k}^{n} (x-xj)/(xk-xj)
Args:
x ([float]): [x]
k ([int]): [k]
Returns:
[float]: [lk(x)]
"""
ta = reduce(lambda x, y: x*y, [x-self.x[i]
for i in range(self.N) if i != k])
tb = reduce(lambda x, y: x*y, [self.x[k]-self.x[i]
for i in range(self.N) if i != k])
return ta/tb
def predict(self, x):
"""predict
Args:
x ([float]): [x]
Returns:
[float]: [predict y]
"""
return np.sum([self.y[k]*self.calculate_lk(x, k) for k in range(self.N)])
x=np.array([0.32,0.34,0.36])
y=np.array([0.314567,0.333487, 0.352274])
LPI=Lagrange_Polynomial_Interpolation(x,y)
x0=0.3367
LPI.predict(x0)
>>> 0.3303743620375
np.sin(x0)-LPI.predict(x0)
>>> -1.7048187195278786e-07
逐次线性插值
Lagrange插值多项式计算函数近似值时，如果加入了新的采样点，则所有权重表达式都需要重新计算，不是很高效，逐次线性插值法可以在 已有计算的数据上得到高阶的权重表达，更加方便和高效。
一般情况下， 两个k次的插值多项式可以通过线性插值得到一个k+1次的插值多项式 （具体推导详见教材21页）。
记I i 1 , i 2 , . . . , i n ( x ) I_{i_1,i_2,…,i_n}(x)I i 1 ,i 2 ,…,i n (x )表示关于点x i 1 , x i 2 , . . . , x i n x_{i_1},x_{i_2},…,x_{i_n}x i 1 ,x i 2 ,…,x i n 的( n − 1 ) (n-1)(n −1 )次插值多项式 （n个数据点最多组成n-1次），特别地，I i k ( x ) I_{i_k}(x)I i k (x )为零次插值多项式，即函数值f ( x i k ) f(x_{i_k})f (x i k )。
Neville算法：
I 0 , 1 , . . . , k + 1 ( x ) = I 0 , 1 , . . . , k ( x ) + I 1 , . . . , k + 1 ( x ) − I 0 , 1 , . . . , k ( x k + 1 − x 0 ) ( x − x 0 ) \Large I_{0,1,…,k+1}(x)=I_{0,1,…,k}(x)+\frac{I_{1,…,k+1}(x)-I_{0,1,…,k}}{(x_{k+1}-x_0)}(x-x_0)I 0 ,1 ,…,k +1 (x )=I 0 ,1 ,…,k (x )+(x k +1 −x 0 )I 1 ,…,k +1 (x )−I 0 ,1 ,…,k (x −x 0 )
可以理解为要计算0到k+1点的k+1次插值多项式时，可以由0到k点的k次插值多项式和1到k+1点的k次插值多项式通过线性插值得到，这样便组成了一个可以递归求解的公式。
代码实现
class mini_I(object):
"""
class for calculating the interpolation of order k+1 from two component of order k
base funtion: ans=y0+(y1-y0)/(x1-x0)*(x-x0)
Author: Junno
Date: 2022-03-01
"""
def __init__(self, x0=0, x1=0, y0=0, y1=0, first=False):
"""__init__
Args:
x0 (float, optional): [component of the base function]. Defaults to 0.
x1 (float, optional): [component of the base function]. Defaults to 0.
y0 (float, optional): [component of the base function]. Defaults to 0.
y1 (float, optional): [component of the base function]. Defaults to 0.
first (bool, optional): [whether the first element of each order list]. Defaults to False.
"""
self.x0 = x0
self.x1 = x1
self.y0 = y0
self.y1 = y1
self.first = first
self.last_answer = {}
def __call__(self, x=0.):
"""__call__
Args:
x ([float]): [interpolation node]
Returns:
[answer]: [recursive call or the recorded values]
"""
if x not in self.last_answer.keys():
answer = self.y0(x)+(self.y1(x)-self.y0(x))/(self.x1 -
self.x0)*(x-self.x0) if not self.first else self.y0
self.last_answer[x] = answer
return answer
else:
if x in self.last_answer.keys():
return self.last_answer[x]
class Successive_Linear_Interpolation_2D(object):
"""Successive_Linear_Interpolation_2D 逐次线性插值
Author: Junno
Date: 2022-03-01
"""
def __init__(self, x, y, base_order=0):
"""__init__
Args:
x ([np.array]): [x of points]
y ([np.array]): [y of points]
base_order ([int]): [interpolation order]
"""
self.N = x.shape[0]
assert base_order < self.N, 'The base_order should be less than the number of points'
self.x = np.copy(x)
self.y = np.copy(y)
self.base_order = 0
self.calculated = False
self.F = lambda i, j: mini_I(
self.x[i-j], self.x[i], self.I_set[i-1][j-1], self.I_set[i][j-1], first=False)
def calculate_I(self, x, eps=1e-6, order=1):
"""calculate_I calculate y for a given x and order with threshold eps
Args:
x ([float]): [x]
eps ([float], optional): [threshold to stop process]. Defaults to 1e-6.
order ([int], optional): [interpolation order]. Defaults to 1.
"""
assert order < self.N, 'The order should be less than the number of points'
if not self.calculated:
self.base_order = self.N if order is None else order
self.I_set = [[]]
self.I_set[0].append(mini_I(y0=self.y[0], first=True))
last = 0
for i in range(1, self.base_order+1):
self.I_set += [[]]
self.I_set[i].append(mini_I(y0=self.y[i], first=True))
for k in range(1, i+1):
self.I_set[i].append(self.F(i, k))
if i > 1:
new = self.I_set[i][-1](x)
if abs(last-new) eps:
if i < (self.base_order+1)-1:
print(
"Meet with absotion condition of eps={} and the current order of interpolation is {}".format(eps, i))
self.base_order = i
self.calculated = True
return new
last = self.I_set[i][-1](x)
self.calculated = True
return self.I_set[self.base_order][self.base_order](x)
else:
last = self.I_set[0][0](x)
for i in range(1, order+1):
if i > self.base_order:
self.I_set += [[]]
self.I_set[i].append(mini_I(y0=self.y[i], first=True))
for k in range(1, i+1):
self.I_set[i].append(self.F(i, k))
new = self.I_set[i][-1](x)
if abs(last-new) eps:
if check:
print(
"Meet with absotion condition of eps={} and the current order of interpolation is {}".format(eps, i))
return new
last = new
return self.I_set[order][order](x)
def add_point(self, x, y):
"""add_point
Args:
x ([np.array]): [x of points]
y ([np.array]): [y of points]
"""
M = x.shape[0]
self.x = np.concatenate((self.x, x), axis=0)
self.y = np.concatenate((self.y, y), axis=0)
self.N = self.N+M
def print(self, x):
"""print show I_set for a given x
Args:
x ([float]): [x]
"""
for i in range(len(self.I_set)):
print('I_{}:'.format(i), end='\t')
for j in range(len(self.I_set[i])):
print(self.I_set[i][j](x), end='\t')
print('\n')
import numpy as np
x=np.array([0.32,0.34,0.36])
y=np.array([0.314567,0.333487, 0.352274])
SLI=Successive_Linear_Interpolation_2D(x,y)
SLI.calculate_I(0.3345, order=2)
>>> 0.32829725843749996
SLI.calculate_I(0.3345, order=2)-np.sin(0.3345)
>>> 3.3479488331655816e-07
SLI.print(0.3345)
>>>
I_0: 0.314567
I_1: 0.333487 0.32828399999999996
I_2: 0.352274 0.32832057499999995 0.32829725843749996
SLI.add_point(np.array([0.38+i*0.02 for i in range(3)]),np.array([np.sin(0.38+i*0.02) for i in range(3)]))
SLI.x
>>> array([0.32, 0.34, 0.36, 0.38, 0.4 , 0.42])
SLI.calculate_I(0.3678, order=5)
>>>
Meet with stopping condition of eps=1e-06 and the current order of interpolation is 4
0.3595632718504201
SLI.print(0.3678)
>>>
I_0: 0.314567
I_1: 0.333487 0.35978579999999993
I_2: 0.352274 0.35960093000000004 0.35956488035000006
I_3: 0.3709204694129827 0.35954612307106326 0.35956283918438897 0.35956325422139657
I_4: 0.3894183423086505 0.35963676694662533 0.3595637986267979 0.3595632837260384 0.3595632718504201
逼近复杂函数
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
x=np.arange(0,1,0.01)
F=lambda x:np.cos(2*np.pi*x)+np.sin(x)+np.exp(x**2)
xs=np.sort(np.random.choice(x, size=(10), replace=False))
ys=F(xs)
ys_noise=F(xs)+np.random.randn(xs.shape[0],)*0.05
SLI=Successive_Linear_Interpolation(xs,ys)
SLI_n=Successive_Linear_Interpolation(xs,ys_noise)
y_n=np.zeros_like(x)
y=np.zeros_like(x)
for k in range(x.shape[0]):
y[k]=SLI.calculate_I(x[k],order=xs.shape[0]-1,check=False)
y_n[k]=SLI_n.calculate_I(x[k],order=xs.shape[0]-1,check=False)
plt.figure(figsize=(8,5))
S=plt.scatter(xs,ys,s=5, marker='*',c='black')
Predict_without_noise=plt.plot(x,y,c='red')
Predict_with_noise=plt.plot(x,y_n,c='yellow')
plt.xlabel('x')
plt.ylabel('y')
Samples=mpatches.Patch(color='black', label='Samples Points')
Predict_without_noise=mpatches.Patch(color='red', label='Predict_without_noise')
Predict_with_noise=mpatches.Patch(color='yellow', label='Predict_with_noise')
plt.title('Successive_Linear_Interpolation with/withour noise')
plt.legend(handles=[Samples,Predict_without_noise,Predict_with_noise],loc='best')
plt.show()
- 可以看到线性插值对噪声还是非常敏感的， 可以保证的是插值多项式一定过采样点，但是由于噪声两边会剧烈变化。 可见，线性插值方法只有在插值节点数据足够准确时才可以比较好地逼近原函数。
- 当加入的噪声太大时，有可能出现以下的情况：
Newton插值法
- 本质上由于只有n+1个点的情况下，线性插值多项式存在且唯一存在。 不同插值方法，如Lagrange、Neville、Newton等都是从 不同的基函数入手来对插值多项式进行化简总结， 本质上得到的插值多项式以及插值余项在数学上是等价的。 在Lagrange的基础上，Neville和Newton插值法可以利用已计算的阶数数据来求解下一阶的内容，更加高效便捷。
- 一般Newton插值法由 差商概念引入，这里 从最本质的Lagrange插值法推出Newton插值法 首先，Lagrange插值法的n+1点插值多项式表达式如下： L n ( x ) = ∑ k = 0 n y k l k ( x ) = ∑ i = 0 n y i ∏ j = 0 , j ≠ i n x − x j x i − x j \Large L_n(x)=\sum_{k=0}^ny_kl_k(x)=\sum_{i=0}^n y_i \prod_{j=0,j \neq i}^n \frac{x-x_j}{x_i-x_j}L n (x )=∑k =0 n y k l k (x )=∑i =0 n y i ∏j =0 ,j =i n x i −x j x −x j Lagrange插值法由于一旦加入新的插值节点，所有插值基函数都要重新计算，所以我们需要 找到两个相邻阶次的Lagrange插值基函数的差并且可以简洁表达，这样就可以避免重新计算所有数据。 记n阶插值多项式为l n ( x ) l_n(x)l n (x )，n-1阶插值多项式记为l n − 1 ( x ) l_{n-1}(x)l n −1 (x )，记D n ( x ) = l n ( x ) − l n − 1 ( x ) D_n(x)=l_n(x)-l_{n-1}(x)D n (x )=l n (x )−l n −1 (x ) 由前面的内容我们可知，插值基函数在插值节点处的值恒为y i y_i y i ，即在点x i , ( i = 0 , 1 , . . . , n − 1 ) x_i,(i=0,1,…,n-1)x i ,(i =0 ,1 ,…,n −1 )处，l n ( x i ) = l n − 1 ( x i ) l_n(x_i)=l_{n-1}(x_i)l n (x i )=l n −1 (x i )，也就表明D n ( x ) D_n(x)D n (x )有着n个零点，分别是x = x i , ( i = 0 , 1 , . . . , n − 1 ) x=x_i,(i=0,1,…,n-1)x =x i ,(i =0 ,1 ,…,n −1 )，可以确定D n ( x ) D_n(x)D n (x )的表达式可以化简为含有n个1次差项的形式： D n ( x ) = C n ( x ) ⋅ ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) = C n ( x ) ∏ i = 0 n − 1 ( x − x i ) \Large D_n(x)=C_n(x) \cdot (x-x_0)(x-x_1)…(x-x_{n-1})=C_n(x) \prod_{i=0}^{n-1} (x-x_i)D n (x )=C n (x )⋅(x −x 0 )(x −x 1 )…(x −x n −1 )=C n (x )∏i =0 n −1 (x −x i )。 其次，由于D n ( x ) D_n(x)D n (x )是n阶多项式与n-1阶多项式的差，其次数不会超过n阶，所以可以确定C n ( x ) C_n(x)C n (x )为常数项。那么现在的任务就是确定C n C_n C n 的表达式。 在新的插值节点x n x_n x n 处，我们有： C n ∏ i = 0 n − 1 ( x n − x i ) = y n − ∑ i = 0 n − 1 y i ∏ j = 0 , j ≠ i n − 1 x n − x j x i − x j \Large C_n \prod_{i=0}^{n-1} (x_n-x_i)=y_n-\sum_{i=0}^{n-1} y_i \prod_{j=0,j \neq i}^{n-1} \frac{x_n-x_j}{x_i-x_j}C n ∏i =0 n −1 (x n −x i )=y n −∑i =0 n −1 y i ∏j =0 ,j =i n −1 x i −x j x n −x j C n = y n ∏ i = 0 n − 1 ( x n − x i ) − ∑ i = 0 n − 1 y i ∏ j = 0 , j ≠ i n − 1 x n − x j x i − x j ∏ i = 0 n − 1 ( x n − x i ) \Large C_n=\frac{y_n}{\prod_{i=0}^{n-1} (x_n-x_i)}-\frac{\sum_{i=0}^{n-1} y_i \prod_{j=0,j \neq i}^{n-1} \frac{x_n-x_j}{x_i-x_j}}{\prod_{i=0}^{n-1} (x_n-x_i)}C n =∏i =0 n −1 (x n −x i )y n −∏i =0 n −1 (x n −x i )∑i =0 n −1 y i ∏j =0 ,j =i n −1 x i −x j x n −x j 可以发现右边第二项中，分母的连乘部分比分子部分多了一个因式( x n − x i ) (x_n-x_i)(x n −x i )，我们可以进一步化简： C n = y n ∏ i = 0 n − 1 ( x n − x i ) − ∑ i = 0 n − 1 y i ∏ j = 0 , j ≠ i n − 1 x n − x j x i − x j ∏ i = 0 n − 1 ( x n − x i ) ⋅ x n − x i x n − x i \Large \Large C_n=\frac{y_n}{\prod_{i=0}^{n-1} (x_n-x_i)}-\frac{\sum_{i=0}^{n-1} y_i \prod_{j=0,j \neq i}^{n-1} \frac{x_n-x_j}{x_i-x_j}}{\prod_{i=0}^{n-1} (x_n-x_i)} \cdot \frac{x_n-x_i}{x_n-x_i}C n =∏i =0 n −1 (x n −x i )y n −∏i =0 n −1 (x n −x i )∑i =0 n −1 y i ∏j =0 ,j =i n −1 x i −x j x n −x j ⋅x n −x i x n −x i C n = y n ∏ i = 0 n − 1 ( x n − x i ) − ∑ i = 0 n − 1 y i ∏ j = 0 , j ≠ i n ( x i − x j ) \Large \Large C_n=\frac{y_n}{\prod_{i=0}^{n-1} (x_n-x_i)}-\sum_{i=0}^{n-1} \frac{ y_i}{\prod_{j=0,j \neq i}^{n} (x_i-x_j)}C n =∏i =0 n −1 (x n −x i )y n −∑i =0 n −1 ∏j =0 ,j =i n (x i −x j )y i C n = ∑ i = 0 n y i ∏ j = 0 , j ≠ i n ( x i − x j ) \Large C_n=\sum_{i=0}^{n} \frac{ y_i}{\prod_{j=0,j \neq i}^{n} (x_i-x_j)}C n =∑i =0 n ∏j =0 ,j =i n (x i −x j )y i
现在我们可以得到L n ( x ) L_n(x)L n (x )的递推式： L 1 = L 0 + C 1 ( x − x 0 ) L 2 = L 1 + C 2 ( x − x 0 ) ( x − x 1 ) . . . L n ( x ) = C 0 + C 1 ( x − x 0 ) + C 2 ( x − x 0 ) ( x − x 1 ) + . . . + C n ⋅ ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) \begin{aligned} L_1 &= L_0+C_1(x-x_0) \ L_2 &= L_1+C_2(x-x_0)(x-x_1)\ …\ L_n(x) &=C_0+C_1(x-x_0)+C_2(x-x_0)(x-x_1)+…+C_n \cdot (x-x_0)(x-x_1)…(x-x_{n-1}) \ \end{aligned}L 1 L 2 …L n (x )=L 0 +C 1 (x −x 0 )=L 1 +C 2 (x −x 0 )(x −x 1 )=C 0 +C 1 (x −x 0 )+C 2 (x −x 0 )(x −x 1 )+…+C n ⋅(x −x 0 )(x −x 1 )…(x −x n −1 ) 回到Newton插值法，我们可以看到教材中引入的差商概念其实就等价于我们的C n C_n C n ： 一阶差商（类似直线斜率）：f [ x 0 , x k ] = f ( x k ) − f ( x 0 ) x k − x 0 \large f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}f [x 0 ,x k ]=x k −x 0 f (x k )−f (x 0 ) 二阶差商：f [ x 0 , x 1 , x k ] = f [ x 0 , x k ] − f [ x 0 , x 1 ] x k − x 1 \large f[x_0,x_1,x_k]=\frac{f[x_0,x_k]-f[x_0,x_1]}{x_k-x_1}f [x 0 ,x 1 ,x k ]=x k −x 1 f [x 0 ,x k ]−f [x 0 ,x 1 ] k阶差商：f [ x 0 , x 1 , . . . , x k ] = ∑ j = 0 k f ( x j ) ∏ i = 0 , i ≠ j k ( x i − x j ) \large f[x_0,x_1,…,x_k]=\sum_{j=0}^k \frac{f(x_j)}{\prod_{i=0,i \neq j}^k(x_i-x_j)}f [x 0 ,x 1 ,…,x k ]=∑j =0 k ∏i =0 ,i =j k (x i −x j )f (x j ) 差商性质：
1. 差商与节点的排列顺序无关，即有f [ x 0 , x 1 , x 2 ] = f [ x 1 , x 0 , x 2 ] = f [ x 2 , x 1 , x 0 ] f[x_0,x_1,x_2]=f[x_1,x_0,x_2]=f[x_2,x_1,x_0]f [x 0 ,x 1 ,x 2 ]=f [x 1 ,x 0 ,x 2 ]=f [x 2 ,x 1 ,x 0 ]
2. 若f ( x ) f(x)f (x )存在n阶导数，则n阶差商与导数的关系有：f [ x 0 , x 1 , . . . , x n ] = f ( n ) ( ξ ) n ! f[x_0,x_1,…,x_n]=\frac{f^{(n)}(\xi)}{n!}f [x 0 ,x 1 ,…,x n ]=n !f (n )(ξ) Newton插值法的插值多项式： N n ( x ) = C 0 + C 1 ( x − x 0 ) + C 2 ( x − x 0 ) ( x − x 1 ) + . . . + C n ⋅ ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) \large N_n(x)=C_0+C_1(x-x_0)+C_2(x-x_0)(x-x_1)+…+C_n \cdot (x-x_0)(x-x_1)…(x-x_{n-1})N n (x )=C 0 +C 1 (x −x 0 )+C 2 (x −x 0 )(x −x 1 )+…+C n ⋅(x −x 0 )(x −x 1 )…(x −x n −1 ) 刚才我们提到不同方法其实殊途同归，只是选取的线性空间的基函数不同，从上可知，Newton法的基函数为： b a s e n ( x ) = ∏ i = 0 n − 1 ( x − x i ) \large base_n(x)=\prod_{i=0}^{n-1}(x-x_i)b a s e n (x )=i =0 ∏n −1 (x −x i )
* Newton插值法的插值余项 定义Newton插值多项式为N n ( x ) N_n(x)N n (x )，插值余项为R n ( x ) = f ( x ) − N n ( x ) R_n(x)=f(x)-N_n(x)R n (x )=f (x )−N n (x ) 在插值节点x = x i , ( i = 0 , 1 , . . . , n ) x=x_i,(i=0,1,…,n)x =x i ,(i =0 ,1 ,…,n )处，R n ( x ) R_n(x)R n (x )为0，即R n ( x ) R_n(x)R n (x )有n+1个零点，由罗尔定理： R n ( n ) ( ξ ) = f ( n ) ( ξ ) − N n ( n ) ( ξ ) = f ( n ) ( ξ ) − C n ⋅ n ! = 0 \large R_n^{(n)}(\xi)=f^{(n)}(\xi)-N_n^{(n)}(\xi)=f^{(n)}(\xi)-C_n \cdot n!=0 R n (n )(ξ)=f (n )(ξ)−N n (n )(ξ)=f (n )(ξ)−C n ⋅n !=0 即C n = f ( n ) ( ξ ) n ! \large C_n=\frac{f^{(n)}(\xi)}{n!}C n =n !f (n )(ξ) 有前面的内容得到，Lagrange插值法的插值余项为： R n ( x ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ∏ i = 0 n ( x − x i ) \large R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x-x_i)R n (x )=(n +1 )!f (n +1 )(ξ)∏i =0 n (x −x i ) 可以得到Newton的插值余项用C n C_n C n 表达为： R n ( x ) = C n + 1 ∏ i = 0 n ( x − x i ) = f [ x 0 , x 1 , . . . , x n ] ∏ i = 0 n ( x − x i ) R_n(x)=C_{n+1}\prod_{i=0}^n (x-x_i)=f[x_0,x_1,…,x_n]\prod_{i=0}^n (x-x_i)R n (x )=C n +1 ∏i =0 n (x −x i )=f [x 0 ,x 1 ,…,x n ]∏i =0 n (x −x i ) 一般取绝对值： R n ( x ) ≤ ∣ C n + 1 ∏ i = 0 n ( x − x i ) = f [ x 0 , x 1 , . . . , x n ] ∏ i = 0 n ( x − x i ) ∣ R_n(x) \leq |C_{n+1}\prod_{i=0}^n (x-x_i)=f[x_0,x_1,…,x_n]\prod_{i=0}^n (x-x_i)|R n (x )≤∣C n +1 ∏i =0 n (x −x i )=f [x 0 ,x 1 ,…,x n ]∏i =0 n (x −x i )∣ 可以看到，相较于Lagrange插值法，Newton插值法的插值余项也适用于f ( x ) f(x)f (x )导数不存在的情形。
* Newton法的计算复杂度： 由上可知，计算Newton插值多项式需要两个项： C n = ∑ i = 0 n y i ∏ j = 0 , j ≠ i n ( x i − x j ) C_n=\sum_{i=0}^{n} \frac{ y_i}{\prod_{j=0,j \neq i}^{n} (x_i-x_j)}C n =∑i =0 n ∏j =0 ,j =i n (x i −x j )y i N n ( x ) = C 0 + C 1 ( x − x 0 ) + C 2 ( x − x 0 ) ( x − x 1 ) + . . . + C n ⋅ ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) N_n(x)=C_0+C_1(x-x_0)+C_2(x-x_0)(x-x_1)+…+C_n \cdot (x-x_0)(x-x_1)…(x-x_{n-1})N n (x )=C 0 +C 1 (x −x 0 )+C 2 (x −x 0 )(x −x 1 )+…+C n ⋅(x −x 0 )(x −x 1 )…(x −x n −1 ) 即n个点时，C n C_n C n 需要n 2 / 2 n^2/2 n 2 /2次乘除法和n − 1 n-1 n −1次加法，同样地，N n ( x ) N_n(x)N n (x )也需要n 2 / 2 n^2/2 n 2 /2次乘除法和n − 1 n-1 n −1次加法，即O ( N ) = N 2 O(N)=N^2 O (N )=N 2。
代码实践
class Newton_Linear_Interpolation:
"""Newton_Linear_Interpolation 逐次线性插值
Author: Junno
Date: 2022-03-03
"""
def __init__(self,x,y,eps=1e-3,check=False):
"""__init__
Args:
x ([np.array]): [x of points]
y ([np.array]): [y of points]
eps ([float], optional): [threshold to stop process]. Defaults to 1e-3.
check ([bool], optional): [whether to show dc_matrix]. Defaults to False.
"""
self.N=x.shape[0]
self.dc_matrix=np.zeros((self.N,self.N))
self.base_order=0
self.eps=eps
self.x=x
self.y=y
self.Difference_coefficient=lambda a,b,x1,x2:(a-b)/(self.x[x1]-self.x[x2])
for i in range(self.N):
for j in range(i+1):
if j==0:
self.dc_matrix[i][j]=y[i]
else:
self.dc_matrix[i][j]=self.Difference_coefficient(self.dc_matrix[i][j-1],self.dc_matrix[i-1][j-1],i,i-j)
if abs(self.dc_matrix[i][j]-self.dc_matrix[i-1][j])self.eps:
self.base_order=i-1
if self.base_order == 0:
self.base_order=self.N-1
if check:
self.show_dc()
def calulate(self,x,order,show_error=True):
"""__init__
Args:
x ([float]): [x to calculate y]
order ([itn]): [demanding order]
show_error ([bool], optional): [whether to error]. Defaults to True.
Return:
y ([float]): [predict y value]
"""
assert order<self.N,'order should be less than the number of points'
t_order = min(order,self.base_order+1)
error=abs(self.dc_matrix[t_order][t_order]*reduce(lambda a,b: a*b, [x-self.x[k] for k in range(self.dc_matrix.shape[0])]))
if show_error:
print('error:{}'.format(error))
return sum([self.dc_matrix[i][i]*(1 if i==0 else reduce(lambda a,b: a*b, [x-self.x[k] for k in range(i)]) ) for i in range(t_order)])
def add_point(self, x, y):
"""add_point
Args:
x ([np.array]): [x of points]
y ([np.array]): [y of points]
"""
M = x.shape[0]
self.x = np.concatenate((self.x, x), axis=0)
self.y = np.concatenate((self.y, y), axis=0)
self.N = self.N+M
self.dc_matrix=np.pad(self.dc_matrix,((0,M),(0,M)))
for i in range(self.base_order+1,self.N):
for j in range(i+1):
if j==0:
self.dc_matrix[i][j]=self.y[i]
else:
self.dc_matrix[i][j]=self.Difference_coefficient(self.dc_matrix[i][j-1],self.dc_matrix[i-1][j-1],i,i-j)
if abs(self.dc_matrix[i][j]-self.dc_matrix[i-1][j])self.eps:
self.base_order=i-1
def show_dc(self):
"""show difference coefficient matrix
"""
for i in range(self.dc_matrix.shape[0]):
print('{}:'.format(self.y[i]), end='\t')
for j in range(self.dc_matrix.shape[1]):
print(self.dc_matrix[i][j], end='\t')
print('\n')
print('base_order={}'.format(self.base_order))
NSL=Newton_Linear_Interpolation(x,y,check=True)
>>>
0.41075: 0.41075 0.0 0.0 0.0 0.0 0.0
0.57815: 0.57815 1.116 0.0 0.0 0.0 0.0
0.69675: 0.69675 1.1859999999999995 0.2799999999999976 0.0 0.0 0.0
0.88811: 0.88811 1.275733333333333 0.35893333333333377 0.19733333333334047 0.0 0.0
1.02652: 1.02652 1.3841000000000017 0.4334666666666749 0.21295238095240318 0.03123809523812543 0.0
1.25382: 1.25382 1.515333333333332 0.5249333333333217 0.22866666666661706 0.031428571428427754 0.0002930402927728068
base_order=4
NSL.calulate(0.596,5,show_error=True)
>>>
error:4.01693316910874e-09
0.6319175080796159
差分与等距节点差分插值
- 这一节的主要内容是 在各个插值节点间隔相等的情况下，可以利用差分形式对Newton插值法进行简化，实际上，等距节点也是经常会使用的手段。
- 前面我们讲到 差商，其定义为： k阶差商：f [ x 0 , x 1 , . . . , x k ] = ∑ j = 0 k f ( x j ) ∏ i = 0 , i ≠ j k ( x i − x j ) \large f[x_0,x_1,…,x_k]=\sum_{j=0}^k \frac{f(x_j)}{\prod_{i=0,i \neq j}^k(x_i-x_j)}f [x 0 ,x 1 ,…,x k ]=∑j =0 k ∏i =0 ,i =j k (x i −x j )f (x j )
- 在等距节点x k = x 0 + k h ( k = 0 , 1 , . . . , n ) x_k=x_0+kh \,\,\, (k=0,1,…,n)x k =x 0 +k h (k =0 ,1 ,…,n )中，定义以h为步长的 向前差分和 向后差分（中心差分这里不做考虑）： 向前差分：Δ f k = f k + 1 − f k \Delta f_k=f_{k+1}-f_k Δf k =f k +1 −f k 向后差分：∇ f k = f k − f k − 1 \nabla f_k=f_k-f_{k-1}∇f k =f k −f k −1 其中，Δ \Delta Δ和∇ \nabla ∇分别称为 向前差分算子和 向后差分算子。 一般地定义m阶差分形式： Δ m f k = Δ m − 1 f k + 1 − Δ m − 1 f k \Delta ^m f_k=\Delta ^{m-1} f_{k+1}-\Delta ^{m-1}f_k Δm f k =Δm −1 f k +1 −Δm −1 f k
- 为了更好地利用函数值来表示各阶差分形式，引入两个符号算子：
- 不变算子I：I f k = f k If_k=f_k I f k =f k
- 位移算子E：E f k = f k + 1 Ef_k=f_{k+1}E f k =f k +1
则两类差分可以分别用这些算子来表示，m阶可以用二项式展开定理得到，即： - 1阶向前差分：Δ f k = f k + 1 − f k = E f k − I f k = ( E − I ) f k \Delta f_k=f_{k+1}-f_k=Ef_k- If_k=(E-I)f_k Δf k =f k +1 −f k =E f k −I f k =(E −I )f k
- 1阶向后差分：∇ f k = f k − f k − 1 = I f k − E − 1 f k = ( I − E − 1 ) f k \nabla f_k=f_k-f_{k-1}=If_k-E^{-1}f_k=(I-E^{-1})f_k ∇f k =f k −f k −1 =I f k −E −1 f k =(I −E −1 )f k
- n阶向前差分：Δ n f k = ( E − I ) n f k = ∑ j = 0 n C n j ( − 1 ) j E n − j f k = ∑ j = 0 n C n j ( − 1 ) j f k + n − j \Delta ^n f_k=(E-I)^n f_k=\sum_{j=0} ^n C_n ^j (-1)^j E^{n-j} f_k=\sum_{j=0} ^n C_n ^j (-1)^j f_{k+n-j}Δn f k =(E −I )n f k =∑j =0 n C n j (−1 )j E n −j f k =∑j =0 n C n j (−1 )j f k +n −j
-
n阶向后差分：∇ n f k = ( I − E − 1 ) n f k = ∑ j = 0 n C n j I j E − ( n − j ) f k = ∑ j = 0 n C n j ( − 1 ) n − j f k + j − n \nabla ^n f_k=(I-E^{-1}) ^n f_k=\sum_{j=0} ^n C_n^j I^j E^{-(n-j)} f_k=\sum_{j=0} ^n C_n^j (-1)^{n-j} f_{k+j-n}∇n f k =(I −E −1 )n f k =∑j =0 n C n j I j E −(n −j )f k =∑j =0 n C n j (−1 )n −j f k +j −n
-
差分与差商的关系： 可以利用归纳法证明如下：
- 当k = 1 k=1 k =1，有f [ x 0 , x 1 ] = f ( x 1 ) − f ( x 0 ) x 1 − x 0 = Δ 1 f 0 1 ! h 1 \large f[x_0,x_1]=\frac{f(x_1)-f(x_0)} {x_1-x_0}=\frac{\Delta ^1 f_0}{1! h^1}f [x 0 ,x 1 ]=x 1 −x 0 f (x 1 )−f (x 0 )=1 !h 1 Δ1 f 0
- 当k = n − 1 k=n-1 k =n −1，有f [ x 0 , x 1 , . . . , x n − 1 ] = Δ n − 1 f 0 ( n − 1 ) ! h n − 1 \large f[x_0,x_1,…,x_{n-1}]=\frac{\Delta ^{n-1} f_0}{(n-1)! h^{n-1}}f [x 0 ,x 1 ,…,x n −1 ]=(n −1 )!h n −1 Δn −1 f 0 ，f [ x 1 , . . . , x n − 1 , x n ] = Δ n − 1 f 1 ( n − 1 ) ! h n − 1 \large f[x_1,…,x_{n-1},x_n]=\frac{\Delta ^{n-1} f_1}{(n-1)! h^{n-1}}f [x 1 ,…,x n −1 ,x n ]=(n −1 )!h n −1 Δn −1 f 1
- 当k = n k=n k =n，推得f [ x 0 , x 1 , . . . , x n − 1 , x n ] = f [ x 1 , . . . , x n − 1 , x n ] − f [ x 0 , x 1 , . . . , x n − 1 ] x n − x 0 = Δ n − 1 ( f 1 − f 0 ) ( n − 1 ) ! h n − 1 ⋅ 1 n h = Δ n f 0 n ! h n \large f[x_0,x_1,…,x_{n-1},x_n]=\frac{f[x_1,…,x_{n-1},x_n]-f[x_0,x_1,…,x_{n-1}]} {x_n-x_0}=\frac{\Delta ^{n-1} (f_1-f_0)}{(n-1)! h^{n-1}} \cdot \frac{1}{nh}=\frac{\Delta ^n f_0}{n! h^n}f [x 0 ,x 1 ,…,x n −1 ,x n ]=x n −x 0 f [x 1 ,…,x n −1 ,x n ]−f [x 0 ,x 1 ,…,x n −1 ]=(n −1 )!h n −1 Δn −1 (f 1 −f 0 )⋅n h 1 =n !h n Δn f 0
- 差分与n阶导数地关系： 前面我们提到差商与导数的关系为：C n = f ( n ) ( ξ ) n ! \large C_n=\frac{f^{(n)}(\xi)}{n!}C n =n !f (n )(ξ) 由上差分与差商的关系，得到差分与导数的关系：Δ n f 0 = h n f n ( ξ ) \Delta ^n f_0=h^n f^n(\xi)Δn f 0 =h n f n (ξ)
- 差分形式的Newton插值多项式 设节点x k = x 0 + k h x_k=x_0+kh x k =x 0 +k h，待求节点为x = x 0 + t h , 0 ≤ t ≤ 1 x=x_0+th, 0\leq t \leq 1 x =x 0 +t h ,0 ≤t ≤1 上面提到，Newton插值法的插值多项式： N n ( x ) = C 0 + C 1 ( x − x 0 ) + C 2 ( x − x 0 ) ( x − x 1 ) + . . . + C n ⋅ ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) \large N_n(x)=C_0+C_1(x-x_0)+C_2(x-x_0)(x-x_1)+…+C_n \cdot (x-x_0)(x-x_1)…(x-x_{n-1})N n (x )=C 0 +C 1 (x −x 0 )+C 2 (x −x 0 )(x −x 1 )+…+C n ⋅(x −x 0 )(x −x 1 )…(x −x n −1 ) 其中，C n C_n C n 为n阶差商，且C n = Δ n f 0 n ! h n \large C_n=\frac{\Delta ^n f_0}{n! h^n}C n =n !h n Δn f 0 。 代入差分与差商的关系式可以推得： N n ( x 0 + t h ) = ∑ i = 0 n C i ∏ j = 0 i ( x − x j ) = ∑ i = 0 n C i h i ∏ j = 0 i − 1 ( t − j ) = ∑ i = 0 n Δ i f 0 i ! h i h i ∏ j = 0 i − 1 ( t − j ) = ∑ i = 0 n C t i Δ i f 0 \Large N_n(x_0+th)=\sum_{i=0}^n C_i \prod_{j=0}^i (x-x_j)=\sum_{i=0}^n C_i h^i \prod_{j=0}^{i-1} (t-j)= \sum_{i=0}^n \frac{\Delta ^i f_0}{i! h^i} h^i \prod_{j=0}^{i-1} (t-j) = \sum_{i=0}^n C_t^i \Delta ^i f_0 N n (x 0 +t h )=∑i =0 n C i ∏j =0 i (x −x j )=∑i =0 n C i h i ∏j =0 i −1 (t −j )=∑i =0 n i !h i Δi f 0 h i ∏j =0 i −1 (t −j )=∑i =0 n C t i Δi f 0 上式称为Newton前插公式，后插公式同理可推得: 设节点x k = x n + k h x_k=x_n+kh x k =x n +k h，待求节点为x = x n + t h , − 1 ≤ t ≤ 0 x=x_n+th, -1 \leq t \leq 0 x =x n +t h ,−1 ≤t ≤0 N n ( x n + t h ) = ∑ i = 0 n ( − 1 ) i C − t i ∇ i f n \large N_n(x_n+th)= \sum_{i=0}^n (-1)^i C_{-t}^i \nabla ^i f_n N n (x n +t h )=∑i =0 n (−1 )i C −t i ∇i f n 这里可以看到 前插形式的Newton插值公式只跟差分和待求节点距离节点x 0 x_0 x 0 的距离有关，计算机的执行上会更加高效。
代码实现
class Equidistant_Newton_Linear_Interpolation:
"""Equidistant_Newton_Linear_Interpolation 等距节点牛顿插值法
Author: Junno
Date: 2022-03-09
"""
def __init__(self,x,y,h):
"""__init__
Args:
x ([np.array]): [x of points]
y ([np.array]): [y of points]
h ([float]): [gap between of points]
"""
self.N=y.shape[0]
self.base_order=0
self.eps=eps
self.x=x
self.y=y
self.h=h
self.binomial_coefficient=lambda j,k:reduce(lambda x,y:x*y,[j-i for i in range(k)])/reduce(lambda x,y:x*y, [i for i in range(1,k+1)]) if k>0 else 1
self.forward_delta_matrix=self.cal_forward_(self.y,self.N)
self.backward_nabla_matrix=self.cal_backward_(self.y,self.N)
if self.base_order == 0:
self.base_order=self.N-1
def cal_forward_(self,y,N):
"""cal_forward_
Args:
y ([np.array]): [y of points]
N ([int]): [number of points]
Returns:
[np.darray]: [forward_delta_matrix]
"""
temp=np.zeros((N,N))
temp[:,0]=y
for j in range(1,N):
for i in range(N-j):
temp[i][j]=temp[i+1][j-1]-temp[i][j-1]
return temp
def cal_backward_(self,y,N):
"""cal_backward_
Args:
y ([np.array]): [y of points]
N ([int]): [number of points]
Returns:
[np.darray]: [backward_nabla_matrix]
"""
temp=np.zeros((N,N))
temp[:,0]=y
for j in range(1,N):
for i in range(j,N):
temp[i][j]=temp[i-1][j-1]-temp[i][j-1]
return temp
def calulate(self,x,eps=1e-3,direction=None,show_error=True):
"""__init__
Args:
x ([float]): [x to calculate y]
eps ([float], optional): [threshold to stop process]. Defaults to 1e-6.
direction ([str]): [difference direction: 'forward' or 'backward']
show_error ([bool], optional): [whether to error]. Defaults to True.
Return:
y ([float]): [predict y value]
"""
if direction==None:
direction='forward' if x<self.x[0]+self.N/2*self.h else 'backward'
ind=np.argmin(abs(self.x-x))
if direction == 'forward':
t=(x-self.x[ind])/self.h
count=np.sum(np.abs(self.forward_delta_matrix[ind])>eps)
ans= sum([self.binomial_coefficient(t,k)*self.forward_delta_matrix[ind][k] for k in range(count)])
if show_error:
ny=np.sort(np.concatenate((self.y,[ans]),axis=0), axis=0)
nm=self.cal_forward_(ny,self.N+1)
error=np.abs(self.binomial_coefficient(t,count)*np.max(np.abs(nm[:,count])))
print("MAX INTERPOLATION ERROR:{}".format(error))
return ans
else:
t=-(x-self.x[ind])/self.h
count=np.sum(np.abs(self.backward_nabla_matrix[ind])>eps)
ans=sum([(-1)**k*self.binomial_coefficient(-t,k)*self.backward_nabla_matrix[ind][k] for k in range(count)])
if show_error:
ny=np.sort(np.concatenate((self.y,[ans]),axis=0), axis=0)
nm=self.cal_backward_(ny,self.N+1)
error=np.abs(self.binomial_coefficient(-t,count)*np.max(np.abs(nm[:,count])))
print("MAX INTERPOLATION ERROR:{}".format(error))
return ans
def add_point(self, x, y):
"""add_point
Args:
x ([np.array]): [x of points]
y ([np.array]): [y of points]
"""
M = x.shape[0]
self.x = np.concatenate((self.x, x), axis=0)
self.y = np.concatenate((self.y, y), axis=0)
self.N = self.N+M
x = np.arange(0,2, 0.2)
F = lambda x:np.sin(x) + np.cos(x**2)
y = F(x)
EN = Equidistant_Newton_Linear_Interpolation(x, y, 0.2)
EN.calulate(0.79, 'forward', show_error = True)
>>>
MAX INTERPOLATION ERROR:0.01681455917928076
1.5224178420317618
F(v)-DN.calulate(v,show_error=False)
>>> -0.0005751910744202782
EN.calulate(1.59, 'backward', show_error = True)
>>>
MAX INTERPOLATION ERROR:0.05321699743238864
0.18476646350550527
F(v)-DN.calulate(v,show_error=False)
>>> -0.0025930434075897846
Hermite插值方法
- 有些实际插值问题不仅要求在插值节点处插值函数值等于y值，还要求在插值节点处的一阶导数相等。
- 设插值函数为H，则满足H ( x i ) = y i , H ′ ( x i ) = y i ′ ( x ) , ( i = 0 , 1 , . . . , n ) H(x_i)=y_i, H^{‘}(x_i)=y_i^{‘}(x), (i=0,1,…,n)H (x i )=y i ,H ′(x i )=y i ′(x ),(i =0 ,1 ,…,n )。这里一共有2n+2个条件，可以唯一确定一个次数不超过2n+1的多项式，形式如下： H 2 n + 1 ( x ) = a 0 + a 1 x + . . . + a 2 n + 1 x 2 n + 1 H_{2n+1}(x)=a_0+a_1 x+…+a_{2n+1}x^{2n+1}H 2 n +1 (x )=a 0 +a 1 x +…+a 2 n +1 x 2 n +1
- 还是根据Lagrange插值多项式的求解方法，将插值节点处的函数值和一阶导数值作为多项式系数，利用零点和抽样delta性质来求解合适的基函数。记m i = y i ′ ( x ) m_i=y_i^{‘}(x)m i =y i ′(x )，则H的Lagrange多项式形式为： H 2 n + 1 ( x ) = ∑ j = 0 n [ y j α j ( x ) + m j β j ( x ) ] H_{2n+1}(x)=\sum_{j=0}^n[y_j\alpha_j (x)+m_j\beta_j (x)]H 2 n +1 (x )=∑j =0 n [y j αj (x )+m j βj (x )] α j ( x k ) = 1 i f j = k e l s e 0 , α j ′ ( x k ) = 0 0 , ( j , k = 0 , 1 , . . . , n ) \alpha_j(x_k)=1\,\, if \,\,j=k \,\, else \,\, 0,\,\, \alpha^{‘}_j(x_k)=0 \,\, 0, \small (j,k=0,1,…,n)αj (x k )=1 i f j =k e l s e 0 ,αj ′(x k )=0 0 ,(j ,k =0 ,1 ,…,n ) β j ( x k ) = 0 , β j ′ ( x k ) = 1 i f j = k e l s e 0 , ( j , k = 0 , 1 , . . . , n ) \beta_j(x_k)=0, \,\, \beta^{‘}_j(x_k)=1\,\, if \,\,j=k \,\, else \,\, 0, \small (j,k=0,1,…,n)βj (x k )=0 ,βj ′(x k )=1 i f j =k e l s e 0 ,(j ,k =0 ,1 ,…,n ) 即在对应的插值节点x j x_j x j 处，α j ( x ) \alpha_j(x)αj (x )的值为1，其他点出为0，而β j ( x ) \beta_j(x)βj (x )在所有插值节点处都为0，且在对应点出其一阶导数为1。 根据以上性质，可以推理得到，α j ( x ) \alpha_j(x)αj (x ) 至少含有( x − x j ) (x-x_j)(x −x j ) 的平方项（求导后在所有插值节点处都为0），所以我们可以借用Lagrange插值基函数的形式，即l j ( x ) = ∏ k = 0 , k ≠ j n x − x k x j − x k l_j(x)=\prod_{k=0,k \neq j}^n \frac{x-x_k}{x_j-x_k}l j (x )=∏k =0 ,k =j n x j −x k x −x k 但由于α j ( x ) \alpha_j(x)αj (x ) 是不超过2n+1次的多项式，l i ( x ) l_i(x)l i (x ) 是不超过n次的多项式，所以可以将α j ( x ) \alpha_j(x)αj (x )表示为： α j ( x ) = ( a x + b ) l i 2 ( x ) \alpha_j(x)=(ax+b)l_i^2(x)αj (x )=(a x +b )l i 2 (x ) 根据性质有：
- α j ( x j ) = ( a x j + b ) l i 2 ( x j ) = 1 \alpha_j(x_j)=(ax_j+b)l_i^2(x_j)=1 αj (x j )=(a x j +b )l i 2 (x j )=1
- α j ′ ( x j ) = 2 ( a x j + b ) l i ( x j ) l i ′ ( x j ) + a l i 2 ( x j ) = 0 \alpha^{‘}j(x_j)=2(ax_j+b)l_i(x_j)l_i^{‘}(x_j)+al_i^{2}(x_j)=0 αj ′(x j )=2 (a x j +b )l i (x j )l i ′(x j )+a l i 2 (x j )=0 可以解得 α j ( x ) \alpha_j(x)αj (x )的形式为： α j ( x ) = [ 1 − 2 ( x − x j ) ∑ k = 0 , k ≠ j n 1 x j − x k ] l j 2 ( x ) \large \alpha_j(x)=[1-2(x-x_j)\sum{k=0,k \neq j}^{n}\frac{1}{x_j-x_k}]l_j^2(x)αj (x )=[1 −2 (x −x j )∑k =0 ,k =j n x j −x k 1 ]l j 2 (x ) 同理，可以得到β j ( x ) \beta_j(x)βj (x )的形式： β j ( x ) = ( x − x j ) l j 2 ( x ) \large \beta_j(x)=(x-x_j) l_j^2(x)βj (x )=(x −x j )l j 2 (x )
- 插值余项 可以参照Lagrange的插值余项的证明过程，表示为： R ( x ) = f ( 2 n + 2 ) ( ξ ) ( 2 n + 2 ) ! w n + 1 2 ( x ) , w n + 1 ( x ) = ∏ i = 0 n ( x − x i ) 2 \large R(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}w^2_{n+1}(x), \,\, w_{n+1}(x)= \prod_{i=0}^n (x-x_i)^2 R (x )=(2 n +2 )!f (2 n +2 )(ξ)w n +1 2 (x ),w n +1 (x )=∏i =0 n (x −x i )2
写在最后
- 欢迎大家对代码或者内容有问题意见的留言指正，一起进步！当然这也是我闲时的内容，不保证秒回喔！
Original: https://blog.csdn.net/weixin_40519529/article/details/123215263
Author: 小刀来啦
Title: 数值分析-多项式插值方法小结
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/695505/
转载文章受原作者版权保护。转载请注明原作者出处！