前言

​   最近在学数值分析，把一些算法实现一下

​   具体内容详见 《数值分析》 李庆扬、王能超、易大义（华中科大出版社）

插值的应用与唯一性

• 插值的应用：

1. 对于一些复杂函数，我们往往只能知道某些具体位置的函数值（实际中数据还有误差），当我们需要得到全体定义域上任一点的数值时，便可以采用插值的方法
2. 可以用多项式插值函数来拟合复杂的数据点，简化函数表达 。（多项式历来最被研究者喜欢，可积可导，性质很好！）

• 对于给定的n个点，所有的xi≠xjx_i \neq x_jxi​​=xj​，即所有的采样点均不同，存在且唯一的次数不超过n次的多项式经过所有采样点，称之为插值多项式。（思路：组成的n阶线性方程组的系数矩阵行列式为范德蒙行列式，行列式不为零，所以存在唯一解。）

Lagrange插值法和逐次线性插值

• Lagrange线性插值

  • Lagrange插值的思想是不去求解各x幂次的系数，直接把函数值y作为系数确定下来，反过来确定线性空间的基的形式！

  • 考虑两点的线性插值

    已知f(xk)=yk,f(xk+1)=yk+1f(x_k)=y_k, f(x_{k+1})=y_{k+1}f(xk​)=yk​,f(xk+1​)=yk+1​，设两点的一阶插值多项式为L1(x)L_1(x)L1​(x)，有L1(xk)=yk,L1(xk+1)=yk+1L_1(x_k)=y_k, L_1(x_{k+1})=y_{k+1}L1​(xk​)=yk​,L1​(xk+1​)=yk+1​。

    注意插值多项式一定经过对应的插值节点

    按照Lagrange插值的思想（由直线的点斜式和两点式可以简单推来），我们可以得到：

    L1(x)=x−xk+1xk−xk+1yk+x−xkxk+1−xkyk+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}L1​(x)=xk​−xk+1​x−xk+1​​yk​+xk+1​−xk​x−xk​​yk+1​

    其实可以利用插值多项式的零点和抽样delta的性质来得到多点Lagrange插值多项式的形式：

    • 零点条件：对于lk(x)l_k(x)lk​(x)，在除第k个采样点的点上始终为0
    • 抽样性质：对于lk(x)l_k(x)lk​(x)，在第k个采样点的点上为1

  • 给出lk(x)l_k(x)lk​(x)的一般形式：

    lk(x)=(x−x0)...(x−xk−1)(x−xk+1)...(x−xn)(xk−x0)...(xk−xk−1)(xk−xk+1)...(xk−xn)(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)lk​(x)=(xk​−x0​)...(xk​−xk−1​)(xk​−xk+1​)...(xk​−xn​)(x−x0​)...(x−xk−1​)(x−xk+1​)...(x−xn​)​(k=0,1,...,n)

  • Lagrange插值多项式表达式：

    Ln(x)=∑k=0nyklk(x)=∑i=0nyi∏j=0,j≠inx−xjxi−xj\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}Ln​(x)=∑k=0n​yk​lk​(x)=∑i=0n​yi​∏j=0,j​=in​xi​−xj​x−xj​​

    其中，我们可以看到Lagrange插值多项式的基函数为：

    basen(x)=lk(x)=∏j=0,j≠knx−xjxk−xjbase_n(x)=l_k(x)=\prod_{j=0,j \neq k}^n \frac{x-x_j}{x_k-x_j} basen​(x)=lk​(x)=j=0,j​=k∏n​xk​−xj​x−xj​​

  • 插值余项
    定义截断误差（插值余项）为Rn(x)=f(x)−Ln(x)R_n(x)=f(x)-L_n(x)Rn​(x)=f(x)−Ln​(x)
    根据罗尔定理，可以推得：

    Rn(x)=f(n+1)(ξ)(n+1)!∏i=0n(x−xi)\large R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x-x_i)Rn​(x)=(n+1)!f(n+1)(ξ)​∏i=0n​(x−xi​)

代码实现

import numpy as np
from functools import reduce# Lagrange 线性插值函数
class Lagrange_Polynomial_Interpolation(object):"""Lagrange_Polynomial_Interpolation base function: Ln(x)=\sum_{k=0}^{n} yk l_k(x)Author: JunnoDate: 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/tbdef 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)])

## 测试教材例子
# f(x)=sin(x)
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页）。

记Ii1,i2,...,in(x)I_{i_1,i_2,...,i_n}(x)Ii1​,i2​,...,in​​(x)表示关于点xi1,xi2,...,xinx_{i_1},x_{i_2},...,x_{i_n}xi1​​,xi2​​,...,xin​​的(n−1)(n-1)(n−1)次插值多项式 （n个数据点最多组成n-1次），特别地，Iik(x)I_{i_k}(x)Iik​​(x)为零次插值多项式，即函数值f(xik)f(x_{i_k})f(xik​​)。

Neville算法：

I0,1,...,k+1(x)=I0,1,...,k(x)+I1,...,k+1(x)−I0,1,...,k(xk+1−x0)(x−x0)\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)I0,1,...,k+1​(x)=I0,1,...,k​(x)+(xk+1​−x0​)I1,...,k+1​(x)−I0,1,...,k​​(x−x0​)

可以理解为要计算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 kbase funtion: ans=y0+(y1-y0)/(x1-x0)*(x-x0)Author: JunnoDate: 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 = x0self.x1 = x1self.y0 = y0self.y1 = y1self.first = first# record answer for fast calculationself.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():# here x0,x1,y0,y1 can be a mini_I class to call for input xanswer = self.y0(x)+(self.y1(x)-self.y0(x))/(self.x1 -self.x0)*(x-self.x0) if not self.first else self.y0# recordself.last_answer[x] = answerreturn answerelse:if x in self.last_answer.keys():return self.last_answer[x]class Successive_Linear_Interpolation_2D(object):"""Successive_Linear_Interpolation_2D 逐次线性插值Author: JunnoDate: 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# flag of calculatedself.calculated = False# ref:《数值分析》 李庆扬、王能超、易大义（华中科大出版社）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 epsArgs: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 orderself.I_set = [[]]self.I_set[0].append(mini_I(y0=self.y[0], first=True))last = 0for 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))# early stop when the difference between of two results from k order and k+1 order is less than epsif i > 1:new = self.I_set[i][-1](x)# print(i, abs(last-new) <= eps)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.I_set.pop(-1)self.calculated = Truereturn newlast = self.I_set[i][-1](x)# print(i,)# self.print(x)self.calculated = Truereturn 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 newlast = newreturn 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+Mdef print(self, x):"""print show I_set for a given xArgs: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
# f(x)=sin(x)
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.32829725843749996SLI.calculate_I(0.3345, order=2)-np.sin(0.3345)
>>> 3.3479488331655816e-07# print I_set structure
SLI.print(0.3345)
>>>
I_0:    0.314567    I_1:    0.333487    0.32828399999999996 I_2:    0.352274    0.32832057499999995 0.32829725843749996 # add points
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])# calculate x=0.3678 of order 5-th
SLI.calculate_I(0.3678, order=5)
>>>
Meet with stopping condition of eps=1e-06 and the current order of interpolation is 4
0.3595632718504201SLI.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) # x.shape=200
F=lambda x:np.cos(2*np.pi*x)+np.sin(x)+np.exp(x**2) # original_function# sample x
xs=np.sort(np.random.choice(x, size=(10), replace=False))# with/without measurement error
ys=F(xs)
ys_noise=F(xs)+np.random.randn(xs.shape[0],)*0.05SLI=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)# show results
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')
# Raw=plt.plot(x,F(x), c='g')
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点插值多项式表达式如下：

  Ln(x)=∑k=0nyklk(x)=∑i=0nyi∏j=0,j≠inx−xjxi−xj\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}Ln​(x)=∑k=0n​yk​lk​(x)=∑i=0n​yi​∏j=0,j​=in​xi​−xj​x−xj​​

  Lagrange插值法由于一旦加入新的插值节点，所有插值基函数都要重新计算，所以我们需要找到两个相邻阶次的Lagrange插值基函数的差并且可以简洁表达，这样就可以避免重新计算所有数据。

  记n阶插值多项式为ln(x)l_n(x)ln​(x)，n-1阶插值多项式记为ln−1(x)l_{n-1}(x)ln−1​(x)，记Dn(x)=ln(x)−ln−1(x)D_n(x)=l_n(x)-l_{n-1}(x)Dn​(x)=ln​(x)−ln−1​(x)

  由前面的内容我们可知，插值基函数在插值节点处的值恒为yiy_iyi​，即在点xi,(i=0,1,...,n−1)x_i,(i=0,1,...,n-1)xi​,(i=0,1,...,n−1)处，ln(xi)=ln−1(xi)l_n(x_i)=l_{n-1}(x_i)ln​(xi​)=ln−1​(xi​)，也就表明Dn(x)D_n(x)Dn​(x)有着n个零点，分别是x=xi,(i=0,1,...,n−1)x=x_i,(i=0,1,...,n-1)x=xi​,(i=0,1,...,n−1)，可以确定Dn(x)D_n(x)Dn​(x)的表达式可以化简为含有n个1次差项的形式：

  Dn(x)=Cn(x)⋅(x−x0)(x−x1)...(x−xn−1)=Cn(x)∏i=0n−1(x−xi)\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)Dn​(x)=Cn​(x)⋅(x−x0​)(x−x1​)...(x−xn−1​)=Cn​(x)∏i=0n−1​(x−xi​)。

  其次，由于Dn(x)D_n(x)Dn​(x)是n阶多项式与n-1阶多项式的差，其次数不会超过n阶，所以可以确定Cn(x)C_n(x)Cn​(x)为常数项。那么现在的任务就是确定CnC_nCn​的表达式。

  在新的插值节点xnx_nxn​处，我们有：

  Cn∏i=0n−1(xn−xi)=yn−∑i=0n−1yi∏j=0,j≠in−1xn−xjxi−xj\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}Cn​∏i=0n−1​(xn​−xi​)=yn​−∑i=0n−1​yi​∏j=0,j​=in−1​xi​−xj​xn​−xj​​

  Cn=yn∏i=0n−1(xn−xi)−∑i=0n−1yi∏j=0,j≠in−1xn−xjxi−xj∏i=0n−1(xn−xi)\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)}Cn​=∏i=0n−1​(xn​−xi​)yn​​−∏i=0n−1​(xn​−xi​)∑i=0n−1​yi​∏j=0,j​=in−1​xi​−xj​xn​−xj​​​

  可以发现右边第二项中，分母的连乘部分比分子部分多了一个因式(xn−xi)(x_n-x_i)(xn​−xi​)，我们可以进一步化简：

  Cn=yn∏i=0n−1(xn−xi)−∑i=0n−1yi∏j=0,j≠in−1xn−xjxi−xj∏i=0n−1(xn−xi)⋅xn−xixn−xi\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}Cn​=∏i=0n−1​(xn​−xi​)yn​​−∏i=0n−1​(xn​−xi​)∑i=0n−1​yi​∏j=0,j​=in−1​xi​−xj​xn​−xj​​​⋅xn​−xi​xn​−xi​​

  Cn=yn∏i=0n−1(xn−xi)−∑i=0n−1yi∏j=0,j≠in(xi−xj)\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)}Cn​=∏i=0n−1​(xn​−xi​)yn​​−∑i=0n−1​∏j=0,j​=in​(xi​−xj​)yi​​

  Cn=∑i=0nyi∏j=0,j≠in(xi−xj)\Large C_n=\sum_{i=0}^{n} \frac{ y_i}{\prod_{j=0,j \neq i}^{n} (x_i-x_j)}Cn​=∑i=0n​∏j=0,j​=in​(xi​−xj​)yi​​

  现在我们可以得到Ln(x)L_n(x)Ln​(x)的递推式：

  L1=L0+C1(x−x0)L2=L1+C2(x−x0)(x−x1)...Ln(x)=C0+C1(x−x0)+C2(x−x0)(x−x1)+...+Cn⋅(x−x0)(x−x1)...(x−xn−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} L1​L2​...Ln​(x)​=L0​+C1​(x−x0​)=L1​+C2​(x−x0​)(x−x1​)=C0​+C1​(x−x0​)+C2​(x−x0​)(x−x1​)+...+Cn​⋅(x−x0​)(x−x1​)...(x−xn−1​)​

  回到Newton插值法，我们可以看到教材中引入的差商概念其实就等价于我们的CnC_nCn​：

  一阶差商（类似直线斜率）：f[x0,xk]=f(xk)−f(x0)xk−x0\large f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}f[x0​,xk​]=xk​−x0​f(xk​)−f(x0​)​

  二阶差商：f[x0,x1,xk]=f[x0,xk]−f[x0,x1]xk−x1\large f[x_0,x_1,x_k]=\frac{f[x_0,x_k]-f[x_0,x_1]}{x_k-x_1}f[x0​,x1​,xk​]=xk​−x1​f[x0​,xk​]−f[x0​,x1​]​

  k阶差商：f[x0,x1,...,xk]=∑j=0kf(xj)∏i=0,i≠jk(xi−xj)\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[x0​,x1​,...,xk​]=∑j=0k​∏i=0,i​=jk​(xi​−xj​)f(xj​)​

  差商性质：

  1. 差商与节点的排列顺序无关，即有f[x0,x1,x2]=f[x1,x0,x2]=f[x2,x1,x0]f[x_0,x_1,x_2]=f[x_1,x_0,x_2]=f[x_2,x_1,x_0]f[x0​,x1​,x2​]=f[x1​,x0​,x2​]=f[x2​,x1​,x0​]
  2. 若f(x)f(x)f(x)存在n阶导数，则n阶差商与导数的关系有：f[x0,x1,...,xn]=f(n)(ξ)n!f[x_0,x_1,...,x_n]=\frac{f^{(n)}(\xi)}{n!}f[x0​,x1​,...,xn​]=n!f(n)(ξ)​

  Newton插值法的插值多项式：

  Nn(x)=C0+C1(x−x0)+C2(x−x0)(x−x1)+...+Cn⋅(x−x0)(x−x1)...(x−xn−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})Nn​(x)=C0​+C1​(x−x0​)+C2​(x−x0​)(x−x1​)+...+Cn​⋅(x−x0​)(x−x1​)...(x−xn−1​)

  刚才我们提到不同方法其实殊途同归，只是选取的线性空间的基函数不同，从上可知，Newton法的基函数为：

  basen(x)=∏i=0n−1(x−xi)\large base_n(x)=\prod_{i=0}^{n-1}(x-x_i)basen​(x)=i=0∏n−1​(x−xi​)

• Newton插值法的插值余项

  定义Newton插值多项式为Nn(x)N_n(x)Nn​(x)，插值余项为Rn(x)=f(x)−Nn(x)R_n(x)=f(x)-N_n(x)Rn​(x)=f(x)−Nn​(x)

  在插值节点x=xi,(i=0,1,...,n)x=x_i,(i=0,1,...,n)x=xi​,(i=0,1,...,n)处，Rn(x)R_n(x)Rn​(x)为0，即Rn(x)R_n(x)Rn​(x)有n+1个零点，由罗尔定理：

  Rn(n)(ξ)=f(n)(ξ)−Nn(n)(ξ)=f(n)(ξ)−Cn⋅n!=0\large R_n^{(n)}(\xi)=f^{(n)}(\xi)-N_n^{(n)}(\xi)=f^{(n)}(\xi)-C_n \cdot n!=0Rn(n)​(ξ)=f(n)(ξ)−Nn(n)​(ξ)=f(n)(ξ)−Cn​⋅n!=0

  即Cn=f(n)(ξ)n!\large C_n=\frac{f^{(n)}(\xi)}{n!}Cn​=n!f(n)(ξ)​

  有前面的内容得到，Lagrange插值法的插值余项为：

  Rn(x)=f(n+1)(ξ)(n+1)!∏i=0n(x−xi)\large R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x-x_i)Rn​(x)=(n+1)!f(n+1)(ξ)​∏i=0n​(x−xi​)

  可以得到Newton的插值余项用CnC_nCn​表达为：

  Rn(x)=Cn+1∏i=0n(x−xi)=f[x0,x1,...,xn]∏i=0n(x−xi)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)Rn​(x)=Cn+1​∏i=0n​(x−xi​)=f[x0​,x1​,...,xn​]∏i=0n​(x−xi​)

  一般取绝对值：

  Rn(x)≤∣Cn+1∏i=0n(x−xi)=f[x0,x1,...,xn]∏i=0n(x−xi)∣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)|Rn​(x)≤∣Cn+1​∏i=0n​(x−xi​)=f[x0​,x1​,...,xn​]∏i=0n​(x−xi​)∣

  可以看到，相较于Lagrange插值法，Newton插值法的插值余项也适用于f(x)f(x)f(x)导数不存在的情形。

• Newton法的计算复杂度：

  由上可知，计算Newton插值多项式需要两个项：

  Cn=∑i=0nyi∏j=0,j≠in(xi−xj)C_n=\sum_{i=0}^{n} \frac{ y_i}{\prod_{j=0,j \neq i}^{n} (x_i-x_j)}Cn​=∑i=0n​∏j=0,j​=in​(xi​−xj​)yi​​

  Nn(x)=C0+C1(x−x0)+C2(x−x0)(x−x1)+...+Cn⋅(x−x0)(x−x1)...(x−xn−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})Nn​(x)=C0​+C1​(x−x0​)+C2​(x−x0​)(x−x1​)+...+Cn​⋅(x−x0​)(x−x1​)...(x−xn−1​)

  即n个点时，CnC_nCn​需要n2/2n^2/2n2/2次乘除法和n−1n-1n−1次加法，同样地，Nn(x)N_n(x)Nn​(x)也需要n2/2n^2/2n2/2次乘除法和n−1n-1n−1次加法，即O(N)=N2O(N)=N^2O(N)=N2。

代码实践

class Newton_Linear_Interpolation:"""Newton_Linear_Interpolation 逐次线性插值Author: JunnoDate: 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=0self.eps=epsself.x=xself.y=yself.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-1if self.base_order == 0:self.base_order=self.N-1if 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+Mself.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-1def 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=4NSL.calulate(0.596,5,show_error=True)
>>>
error:4.01693316910874e-09
0.6319175080796159

差分与等距节点差分插值

• 这一节的主要内容是在各个插值节点间隔相等的情况下，可以利用差分形式对Newton插值法进行简化，实际上，等距节点也是经常会使用的手段。

• 前面我们讲到差商，其定义为：

  k阶差商：f[x0,x1,...,xk]=∑j=0kf(xj)∏i=0,i≠jk(xi−xj)\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[x0​,x1​,...,xk​]=∑j=0k​∏i=0,i​=jk​(xi​−xj​)f(xj​)​

• 在等距节点xk=x0+kh(k=0,1,...,n)x_k=x_0+kh \,\,\, (k=0,1,...,n)xk​=x0​+kh(k=0,1,...,n)中，定义以h为步长的向前差分和向后差分（中心差分这里不做考虑）：

  向前差分：Δfk=fk+1−fk\Delta f_k=f_{k+1}-f_kΔfk​=fk+1​−fk​

  向后差分：∇fk=fk−fk−1\nabla f_k=f_k-f_{k-1}∇fk​=fk​−fk−1​

  其中，Δ\DeltaΔ和∇\nabla∇分别称为向前差分算子和向后差分算子。

  一般地定义m阶差分形式：

  Δmfk=Δm−1fk+1−Δm−1fk\Delta ^m f_k=\Delta ^{m-1} f_{k+1}-\Delta ^{m-1}f_kΔmfk​=Δm−1fk+1​−Δm−1fk​

• 为了更好地利用函数值来表示各阶差分形式，引入两个符号算子：

  • 不变算子I：Ifk=fkIf_k=f_kIfk​=fk​
  • 位移算子E：Efk=fk+1Ef_k=f_{k+1}Efk​=fk+1​

  则两类差分可以分别用这些算子来表示，m阶可以用二项式展开定理得到，即：

  • 1阶向前差分：Δfk=fk+1−fk=Efk−Ifk=(E−I)fk\Delta f_k=f_{k+1}-f_k=Ef_k- If_k=(E-I)f_kΔfk​=fk+1​−fk​=Efk​−Ifk​=(E−I)fk​
  • 1阶向后差分：∇fk=fk−fk−1=Ifk−E−1fk=(I−E−1)fk\nabla f_k=f_k-f_{k-1}=If_k-E^{-1}f_k=(I-E^{-1})f_k∇fk​=fk​−fk−1​=Ifk​−E−1fk​=(I−E−1)fk​
  • n阶向前差分：Δnfk=(E−I)nfk=∑j=0nCnj(−1)jEn−jfk=∑j=0nCnj(−1)jfk+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}Δnfk​=(E−I)nfk​=∑j=0n​Cnj​(−1)jEn−jfk​=∑j=0n​Cnj​(−1)jfk+n−j​
  • n阶向后差分：∇nfk=(I−E−1)nfk=∑j=0nCnjIjE−(n−j)fk=∑j=0nCnj(−1)n−jfk+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}∇nfk​=(I−E−1)nfk​=∑j=0n​Cnj​IjE−(n−j)fk​=∑j=0n​Cnj​(−1)n−jfk+j−n​

• 差分与差商的关系：

  可以利用归纳法证明如下：

  1. 当k=1k=1k=1，有f[x0,x1]=f(x1)−f(x0)x1−x0=Δ1f01!h1\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[x0​,x1​]=x1​−x0​f(x1​)−f(x0​)​=1!h1Δ1f0​​
  2. 当k=n−1k=n-1k=n−1，有f[x0,x1,...,xn−1]=Δn−1f0(n−1)!hn−1\large f[x_0,x_1,...,x_{n-1}]=\frac{\Delta ^{n-1} f_0}{(n-1)! h^{n-1}}f[x0​,x1​,...,xn−1​]=(n−1)!hn−1Δn−1f0​​，f[x1,...,xn−1,xn]=Δn−1f1(n−1)!hn−1\large f[x_1,...,x_{n-1},x_n]=\frac{\Delta ^{n-1} f_1}{(n-1)! h^{n-1}}f[x1​,...,xn−1​,xn​]=(n−1)!hn−1Δn−1f1​​
  3. 当k=nk=nk=n，推得f[x0,x1,...,xn−1,xn]=f[x1,...,xn−1,xn]−f[x0,x1,...,xn−1]xn−x0=Δn−1(f1−f0)(n−1)!hn−1⋅1nh=Δnf0n!hn\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[x0​,x1​,...,xn−1​,xn​]=xn​−x0​f[x1​,...,xn−1​,xn​]−f[x0​,x1​,...,xn−1​]​=(n−1)!hn−1Δn−1(f1​−f0​)​⋅nh1​=n!hnΔnf0​​

• 差分与n阶导数地关系：

  前面我们提到差商与导数的关系为：Cn=f(n)(ξ)n!\large C_n=\frac{f^{(n)}(\xi)}{n!}Cn​=n!f(n)(ξ)​

  由上差分与差商的关系，得到差分与导数的关系：Δnf0=hnfn(ξ)\Delta ^n f_0=h^n f^n(\xi)Δnf0​=hnfn(ξ)

• 差分形式的Newton插值多项式

  设节点xk=x0+khx_k=x_0+khxk​=x0​+kh，待求节点为x=x0+th,0≤t≤1x=x_0+th, 0\leq t \leq 1x=x0​+th,0≤t≤1

  上面提到，Newton插值法的插值多项式：

  Nn(x)=C0+C1(x−x0)+C2(x−x0)(x−x1)+...+Cn⋅(x−x0)(x−x1)...(x−xn−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})Nn​(x)=C0​+C1​(x−x0​)+C2​(x−x0​)(x−x1​)+...+Cn​⋅(x−x0​)(x−x1​)...(x−xn−1​)

  其中，CnC_nCn​为n阶差商，且Cn=Δnf0n!hn\large C_n=\frac{\Delta ^n f_0}{n! h^n}Cn​=n!hnΔnf0​​。

  代入差分与差商的关系式可以推得：

  Nn(x0+th)=∑i=0nCi∏j=0i(x−xj)=∑i=0nCihi∏j=0i−1(t−j)=∑i=0nΔif0i!hihi∏j=0i−1(t−j)=∑i=0nCtiΔif0\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_0Nn​(x0​+th)=∑i=0n​Ci​∏j=0i​(x−xj​)=∑i=0n​Ci​hi∏j=0i−1​(t−j)=∑i=0n​i!hiΔif0​​hi∏j=0i−1​(t−j)=∑i=0n​Cti​Δif0​

  上式称为Newton前插公式，后插公式同理可推得:

  设节点xk=xn+khx_k=x_n+khxk​=xn​+kh，待求节点为x=xn+th,−1≤t≤0x=x_n+th, -1 \leq t \leq 0x=xn​+th,−1≤t≤0

  Nn(xn+th)=∑i=0n(−1)iC−ti∇ifn\large N_n(x_n+th)= \sum_{i=0}^n (-1)^i C_{-t}^i \nabla ^i f_nNn​(xn​+th)=∑i=0n​(−1)iC−ti​∇ifn​

  这里可以看到前插形式的Newton插值公式只跟差分和待求节点距离节点x0x_0x0​的距离有关，计算机的执行上会更加高效。

代码实现

class Equidistant_Newton_Linear_Interpolation:"""Equidistant_Newton_Linear_Interpolation 等距节点牛顿插值法Author: JunnoDate: 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=0self.eps=epsself.x=xself.y=yself.h=h# function to calculate binomial coefficientsself.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# forward difference matrixself.forward_delta_matrix=self.cal_forward_(self.y,self.N)# backward difference matrixself.backward_nabla_matrix=self.cal_backward_(self.y,self.N)if self.base_order == 0:self.base_order=self.N-1def 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]=yfor 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]=yfor 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:# half to be forward and half to be backwarddirection='forward' if x<self.x[0]+self.N/2*self.h else 'backward'# find nearest points indexind=np.argmin(abs(self.x-x))if direction == 'forward':t=(x-self.x[ind])/self.hcount=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 anselse: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 ansdef 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)# predict and show maximum interpolation error
EN.calulate(0.79, 'forward', show_error = True)
>>>
MAX INTERPOLATION ERROR:0.01681455917928076
1.5224178420317618# delta between true value and predict value
F(v)-DN.calulate(v,show_error=False)
>>> -0.0005751910744202782EN.calulate(1.59, 'backward', show_error = True)
>>>
MAX INTERPOLATION ERROR:0.05321699743238864
0.18476646350550527F(v)-DN.calulate(v,show_error=False)
>>> -0.0025930434075897846

Hermite插值方法

• 有些实际插值问题不仅要求在插值节点处插值函数值等于y值，还要求在插值节点处的一阶导数相等。

• 设插值函数为H，则满足H(xi)=yi,H′(xi)=yi′(x),(i=0,1,...,n)H(x_i)=y_i, H^{'}(x_i)=y_i^{'}(x), (i=0,1,...,n)H(xi​)=yi​,H′(xi​)=yi′​(x),(i=0,1,...,n)。这里一共有2n+2个条件，可以唯一确定一个次数不超过2n+1的多项式，形式如下：

  H2n+1(x)=a0+a1x+...+a2n+1x2n+1H_{2n+1}(x)=a_0+a_1 x+...+a_{2n+1}x^{2n+1}H2n+1​(x)=a0​+a1​x+...+a2n+1​x2n+1

• 还是根据Lagrange插值多项式的求解方法，将插值节点处的函数值和一阶导数值作为多项式系数，利用零点和抽样delta性质来求解合适的基函数。记mi=yi′(x)m_i=y_i^{'}(x)mi​=yi′​(x)，则H的Lagrange多项式形式为：

  H2n+1(x)=∑j=0n[yjαj(x)+mjβj(x)]H_{2n+1}(x)=\sum_{j=0}^n[y_j\alpha_j (x)+m_j\beta_j (x)]H2n+1​(x)=∑j=0n​[yj​αj​(x)+mj​βj​(x)]

  αj(xk)=1ifj=kelse0,αj′(xk)=00,(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​(xk​)=1ifj=kelse0,αj′​(xk​)=00,(j,k=0,1,...,n)

  βj(xk)=0,βj′(xk)=1ifj=kelse0,(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​(xk​)=0,βj′​(xk​)=1ifj=kelse0,(j,k=0,1,...,n)

  即在对应的插值节点xjx_jxj​处，α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−xj)(x-x_j)(x−xj​)的平方项（求导后在所有插值节点处都为0），所以我们可以借用Lagrange插值基函数的形式，即lj(x)=∏k=0,k≠jnx−xkxj−xkl_j(x)=\prod_{k=0,k \neq j}^n \frac{x-x_k}{x_j-x_k}lj​(x)=∏k=0,k​=jn​xj​−xk​x−xk​​但由于αj(x)\alpha_j(x)αj​(x)是不超过2n+1次的多项式，li(x)l_i(x)li​(x)是不超过n次的多项式，所以可以将αj(x)\alpha_j(x)αj​(x)表示为：

  αj(x)=(ax+b)li2(x)\alpha_j(x)=(ax+b)l_i^2(x)αj​(x)=(ax+b)li2​(x)

  根据性质有：

  1. αj(xj)=(axj+b)li2(xj)=1\alpha_j(x_j)=(ax_j+b)l_i^2(x_j)=1αj​(xj​)=(axj​+b)li2​(xj​)=1
  2. αj′(xj)=2(axj+b)li(xj)li′(xj)+ali2(xj)=0\alpha^{'}_j(x_j)=2(ax_j+b)l_i(x_j)l_i^{'}(x_j)+al_i^{2}(x_j)=0αj′​(xj​)=2(axj​+b)li​(xj​)li′​(xj​)+ali2​(xj​)=0

  可以解得 αj(x)\alpha_j(x)αj​(x)的形式为：

  αj(x)=[1−2(x−xj)∑k=0,k≠jn1xj−xk]lj2(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−xj​)∑k=0,k​=jn​xj​−xk​1​]lj2​(x)

  同理，可以得到βj(x)\beta_j(x)βj​(x)的形式：

  βj(x)=(x−xj)lj2(x)\large \beta_j(x)=(x-x_j) l_j^2(x)βj​(x)=(x−xj​)lj2​(x)

• 插值余项

  可以参照Lagrange的插值余项的证明过程，表示为：

  R(x)=f(2n+2)(ξ)(2n+2)!wn+12(x),wn+1(x)=∏i=0n(x−xi)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)^2R(x)=(2n+2)!f(2n+2)(ξ)​wn+12​(x),wn+1​(x)=∏i=0n​(x−xi​)2

写在最后

• 欢迎大家对代码或者内容有问题意见的留言指正，一起进步！当然这也是我闲时的内容，不保证秒回喔！

数值分析-多项式插值方法小结相关推荐

1. 数值分析——多项式插值

   多项式插值 这段时间关注了一个数值分析的课程,是华东师范大学潘建瑜老师的课,看了一遍课件,将内容梳理一下,做个笔记. 什么是插值 已知一个函数 f ( x ) f(x) f(x)在 [ a , b ] ...

2. 数值分析 | 多项式求值（Horner法则）

   Horner法则 又名秦九韶算法或嵌套乘法,是一种高效的多项式求值算法. 对于 n n n次多项式 f ( x ) = a 0 + a 1 x + . . . + a n − 1 x n − 1 + ...

3. 分段二次插值函数表达式_数值分析（拟合、插值和逼近）之数据插值方法（线性插值、二次插值、Cubic插值、埃米尔特zz...

   插值.拟合和逼近的区别 据维基百科,科学和工程问题可以通过诸如采样.实验等方法获得若干离散的数据,根据这些数据,我们往往希望得到一个连续的函数(也就是曲线)或者更加密集的离散方程与已知数据相吻合,这过 ...

4. 数值分析：数据插值方法

   http://blog.csdn.net/pipisorry/article/details/62227459 插值.拟合和逼近的区别 据维基百科,科学和工程问题可以通过诸如采样.实验等方法获得若干离 ...

5. matlab 平滑曲线连接_平滑轨迹插值方法之多项式插值（附代码）

   前言 今天我们来聊聊轨迹插值,在机器人的运动规划和控制领域,参考轨迹的生成是一个历史悠久的问题,已经发展出了一系列的方法.今天我们就来聊一聊轨迹插值领域中最常见的轨迹插值方法:多项式插值. 说明:本文 ...

6. 数学建模准备 插值（拉格朗日多项式插值，牛顿多项式插值，分段线性插值，分段三次样条插值，分段三次Hermite插值）

   文章目录 摘要(必看) 0 基础概念 什么是插值 插值用途 什么是拟合 插值和拟合的相同点 插值和拟合的不同点 1 常用的基本插值方法 1.1 多项式插值法 1.1.1 拉格朗日多项式插值法 多项式插 ...

7. 多项式插值之Lagrange、PCHIP与Spline以及BD-Rate和BD-PSNR的计算

   做视频编解码算法的话一般都会接触到一个目前比较通用的客观评价指标,即 Bjøntegaard Delta(BD),如果关注码率变化的话就是 BD-Rate,如果关注 PSNR 变化的话就是BD-PSN ...

8. 【数值分析】拉格朗日插值法

   在数值分析中,拉格朗日插值法是以法国十八世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法.许多实际问题中都用函数来表示某种内在联系或规律,而不少函数都只能通过实验和观测来了解.如对实践中的某个物理量 ...

9. 数值分析c语言编程辛普森公式,数值计算方法：矩形法、梯形法与辛普森公式...

   在数值分析中,数值积分是计算定积分数值的方法和理论.在数学分析中,给定函数的定积分的计算不总是可行的.许多定积分不能用已知的积分公式得到精确值.数值积分是利用黎曼积分和积分中值等数学定义和定理,用数值 ...

最新文章

1. 开源用于寻找系外行星的代码
2. 利用mysql做信息管理_利用MySql实现学生信息管理系统的后台数据管理
3. 使用NuGet发布自己的类库包（Library Package）
4. [转]关于sizeof()的一些思考
5. 51单片机雾化片自动扫频程序_单片机简介
6. 用递归解决冒泡排序问题
7. python语言接收信息的内置函数_python接收信息的内置函数是
8. PE学习.动手写PE.见缝插针
9. 一位36岁程序员的困惑(转)
10. 基于jquery横向手风琴效果
11. 一个老程序员“伯伯”的独白
12. 第十节 直流变直流电路（DCDC）芯片选型
13. angular6添加子路由_如何将Ionicons添加到Angular 6应用
14. 通用代码：发送短信并显示倒计时
15. python计算圆周率_用python程序求圆周率到任意位
16. shopnc route.php,shopnc自动结算的问题
17. 【bzoj3065】: 带插入区间K小值 详解——替罪羊套函数式线段树
18. 西南大学计算机学院推免,2019年西南大学计算机与信息科学学院硕士研究生拟录取名单的公示（不含推免生）...
19. 论文投稿指南——中文核心期刊推荐（力学）
20. 最新Java资源整理，大多数人的选择

热门文章

1. matlab编程刀尖频响,用半理论法预测主轴系统刀尖点频响函数
2. 天猫店铺半自动商品详情数据
3. UnityShader初级篇——渐变纹理
4. 硬盘坏了数据可以恢复吗?
5. 宁夏-银川地区geojson数据
6. 飞秋（FeiQ）关闭好友上下线提示功能
7. 北大才女笔记：这样学习线性回归和梯度下降(上篇)
8. python项目对接钉钉SDK
9. 从零开始的python学习Day4
10. A connection attempt failed because the connected party did not properly respond after a period of……