拉格朗日插值法与牛顿插值多项式

多项式插值

先有一个函数 f ( x ) f(x) f(x)，如果给定在区间 [ a , b ] [a,b] [a,b]上的 n + 1 n+1 n+1个点 a < = x 0 < x 1 < ⋅ ⋅ ⋅ < x n < = b a<=x_0<x_1<···<x_n<=b a<=x0<x1<⋅⋅⋅<xn<=b，已知 f ( x ) f(x) f(x)在这些点的函数值 y 0 , y 1 , ⋅ ⋅ ⋅ , y n y_0,y_1,···,y_n y0,y1,⋅⋅⋅,yn，现在要求一个次数不超过n次的多项式函数P(x)，使得 P ( x i ) = y i P(x_i)=y_i P(xi)=yi，这样的一个过程叫做多项式插值。现在问题是，是否一定存在这样的多项式，如果存在，这样的多项式是否唯一？
求解这个问题就等价于求解线性方程组： { 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 \begin{cases} a_0+a_1x_0+a_2x_0^2+···+a_nx_0^n=y_0\\ a_0+a_1x_1+a_2x_1^2+···+a_nx_1^n=y_1\\·····················\\ a_0+a_1x_n+a_2x_n^2+···+a_nx_n^n=y_n \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧a0+a1x0+a2x02+⋅⋅⋅+anx0n=y0a0+a1x1+a2x12+⋅⋅⋅+anx1n=y1⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅a0+a1xn+a2xn2+⋅⋅⋅+anxnn=yn
这是一个关于系数 ( a 0 , a 1 , ⋅ ⋅ ⋅ , a n ) (a_0,a_1,···,a_n) (a0,a1,⋅⋅⋅,an)的线性方程组，考虑它的系数矩阵行列式： ∣ 1 x 0 x 0 2 x 0 3 ⋅ ⋅ ⋅ x 0 n 1 x 1 x 1 2 x 1 3 ⋅ ⋅ ⋅ x 1 n 1 x 2 x 2 2 x 2 3 ⋅ ⋅ ⋅ x 2 n ⋅ ⋅ ⋅ ⋅ ⋅ 1 x n x n 2 x n 3 ⋅ ⋅ ⋅ x n n ∣ \begin{vmatrix} 1 & x_0 &x_0^2 & x_0^3 ··· &x_0^n \\ 1 & x_1 &x_1^2 & x_1^3 ··· &x_1^n \\ 1 & x_2 &x_2^2 & x_2^3 ··· &x_2^n\\ ·&·&·&·&·\\ 1 & x_n &x_n^2 & x_n^3 ··· &x_n^n\\ \end{vmatrix} ∣∣∣∣∣∣∣∣∣∣111⋅1x0x1x2⋅xnx02x12x22⋅xn2x03⋅⋅⋅x13⋅⋅⋅x23⋅⋅⋅⋅xn3⋅⋅⋅x0nx1nx2n⋅xnn∣∣∣∣∣∣∣∣∣∣
这是一个范德蒙行列式，它的值为： ∏ 0 < = i < j < = n ( x j − x i ) \prod_{0<=i<j<=n}(x_j-x_i) ∏0<=i<j<=n(xj−xi)，由于 x j ≠ x i x_j ≠ x_i xj̸=xi，这个行列式的值不等于零，因此系数矩阵满秩，该线性方程组有且只有一个解。
因此，给定一个点值序列 ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋅ ⋅ ⋅ , ( x n , y n ) (x_0,y_0),(x_1,y_1),···,(x_n,y_n) (x0,y0),(x1,y1),⋅⋅⋅,(xn,yn)，如果这些点的横坐标互不相同，那么这个序列就可以被映射到一个唯一的多项式，反过来，如果给定互不相同的n+1个x，对于任意的两个不同的n次多项式，它们在这些x下的取值一定不完全相同，否则上述方程组就不应该只有一个解。因此点值序列与多项式之间存在一个一一对应的关系。
直接求解上述方程组可以得到结果，但是这个方法效率太低。

拉格朗日插值法

先从最简单的情况出发，即n=1时，这个时候对应的多项式应该是一条直线，由直线的两点式表示法有： L ( x ) = y k x k + 1 − x x k + 1 − x k + y k + 1 x − x k x k + 1 − x k L(x) = y_k\frac{x_{k+1}-x}{x_{k+1}-x_{k}} + y_{k+1}\frac{x-x_k}{x_{k+1}-x_k} L(x)=ykxk+1−xkxk+1−x+yk+1xk+1−xkx−xk
观察这个式子，可以看到 L ( x ) L(x) L(x)实际上是两个线性函数 l k ( x ) = x k + 1 − x x k + 1 − x k , l k + 1 ( x ) = x − x k x k + 1 − x k l_k(x)=\frac{x_{k+1}-x}{x_{k+1}-x_k},l_{k+1}(x)=\frac{x-x_k}{x_{k+1}-x_k} lk(x)=xk+1−xkxk+1−x,lk+1(x)=xk+1−xkx−xk的线性组合
可以看到这两个线性函数的特点： l k ( x k ) = 1 , l k ( x k + 1 ) = 0 l_k(x_k) = 1, l_k(x_{k+1}) =0 lk(xk)=1,lk(xk+1)=0 l k + 1 ( x k ) = 0 , l k + 1 ( x k + 1 ) = 1 l_{k+1}(x_k)=0,l_{k+1}(x_k+1)=1 lk+1(xk)=0,lk+1(xk+1)=1
由此可以试想，如果第k个函数 l k l_k lk只在 x = x k x=x_k x=xk的时候取值为1，其他的时候取值都为0，那么如果代入 x = x k x=x_k x=xk进去，就一定可以得到函数值为 l k l_k lk的系数 y k y_k yk，这样就可以构造出想要的多项式。
称类似 l k l_k lk的函数为线性插值基函数，求n+1个点值的插值多项式，实际上就是求n+1个线性插值基函数，然后分别乘上对应的系数 y k y_k yk即可，这样就是拉格朗日插值法的思想。引入记号 w n + 1 ( x ) = ( x − x 0 ) ( x − x 1 ) ⋅ ⋅ ⋅ ( x − x n ) w_{n+1}(x)=(x-x_0)(x-x_1)···(x-x_n) wn+1(x)=(x−x0)(x−x1)⋅⋅⋅(x−xn) 那么就可以用它来表达 l k l_k lk： l k = y k w n + 1 ( x ) ( x − x k ) w n + 1 ′ ( x k ) l_k=y_k\frac{w_{n+1}(x)}{(x-x_k)w'_{n+1}(x_k)} lk=yk(x−xk)wn+1′(xk)wn+1(x)
于是最终插值的多项式为： L ( x ) = ∑ k = 0 n y k w n + 1 ( x ) ( x − x k ) w n + 1 ′ ( x k ) L(x)=\sum_{k=0}^{n}y_k\frac{w_{n+1}(x)}{(x-x_k)w'_{n+1}(x_k)} L(x)=k=0∑nyk(x−xk)wn+1′(xk)wn+1(x)
计算 w n + 1 ( x ) w_{n+1}(x) wn+1(x)复杂度为 o ( n 2 ) o(n^2) o(n2)，一共需要累加n次，每一次中计算 y k w n + 1 ′ ( x k ) \frac{y_k}{w'_{n+1}(x_k)} wn+1′(xk)yk复杂度为 o ( n ) o(n) o(n)，除以 ( x − x k ) (x-x_k) (x−xk)需要 o ( n ) o(n) o(n)，因此总体复杂度为 o ( n 2 ) o(n^2) o(n2)

洛谷P4781是拉格朗日插值法的模板题，这里给出自己的代码（写的不好常数较大仅供参考）

#include <stdio.h>const int maxN = 2010;typedef long long ll;
ll w_coeffcient[maxN];
ll tmp[maxN];
ll l_coeffcient[maxN];
const ll P = 998244353;
int xx[maxN];
int yy[maxN];ll fast_pow( ll x , ll y)
{ll res = 1;while (y){if (y & 1) res = res * x % P;y >>= 1;x = x * x % P;}return res;
}int main()
{int n, k;scanf("%d%d", &n, &k);int x, y;scanf("%d%d", &x, &y);xx[0] = x;yy[0] = y;w_coeffcient[1] = 1;w_coeffcient[0] = (-x + P ) % P;for (int i = 1; i < n; i++){scanf("%d%d", &x, &y);xx[i] = x;yy[i] = y;for (int j = 0 ; j < i + 1 ; j ++ ){//计算 多项式 *（x - xk ）//多项式 * x 整体系数往后偏移1//多项式 * xk 直接所有系数对应相乘tmp[1 + j] = w_coeffcient[j];w_coeffcient[j] = w_coeffcient[j] * ((-x + P) % P) % P;}for (int j = 0; j <= i + 1; j++){w_coeffcient[j] = (w_coeffcient[j] + tmp[j] + P) % P;}}for (int i = 0; i < n; i ++ ){ll con = yy[i];ll s = 1;//求w'(x_k) for (int j = 0 ; j < n ; j ++ ){if (j == i) continue;s = (s * ((xx[i] - xx[j] + P) % P )) % P;}//因为要除以s 模P意义下求s的逆元ll inv = fast_pow(s, P - 2);//得到yi / w'(x_i)con = con * inv % P;ll new_c = w_coeffcient[n];//让w(x) / (x - x_i) 然后所得系数加至最终得到的结果多项式上for (int j = n - 1 ; j >= 0 ; j -- ){l_coeffcient[j] = (l_coeffcient[j] + (new_c * con % P)) % P  ;new_c = ( (w_coeffcient[j] + (new_c * xx[i] % P))  % P);}}ll res = 0;for (int i = n - 1 ; i >= 0 ; i -- ){res = (res + fast_pow(k, i) * l_coeffcient[i] % P) % P;}printf("%d\n", res);return 0;
}

拉格朗日插值法的误差估计
如果说给定的n个点值本来就位于一个多项式函数曲线上时，那么拉格朗日插值法得到的就是精确的结果，而如果原本的函数并不是一个多项式函数的话，得到的结果只是原函数的一个多项式逼近，这之间是存在一定的误差的，那么下面来估计其误差，设 R n ( x ) = f ( x ) − L ( x ) R_n(x) = f(x)-L(x) Rn(x)=f(x)−L(x)，现有 R n ( x i ) = 0 ( i = 0 , 1 , 2 , ⋅ ⋅ ⋅ , n ) R_n(x_i) = 0 ( i = 0 , 1 , 2 , ··· , n ) Rn(xi)=0(i=0,1,2,⋅⋅⋅,n)，因此可以设 R n ( x ) = K ( x ) ( x − x 0 ) ( x − x 1 ) ⋅ ⋅ ⋅ ( x − x n ) = K ( x ) w n + 1 ( x ) R_n(x) = K(x)(x-x_0)(x-x_1)···(x-x_n)=K(x)w_{n+1}(x) Rn(x)=K(x)(x−x0)(x−x1)⋅⋅⋅(x−xn)=K(x)wn+1(x) K ( x ) K(x) K(x)是一个与x有关的函数。下面令 ϕ ( t ) = f ( t ) − L n ( t ) − K ( x ) ( t − x 0 ) ( t − x 1 ) ⋅ ⋅ ⋅ ( t − x n ) \phi(t) = f(t) - L_n(t) - K(x)(t-x_0)(t-x_1)···(t-x_n) ϕ(t)=f(t)−Ln(t)−K(x)(t−x0)(t−x1)⋅⋅⋅(t−xn)
可知 ϕ ( t ) \phi(t) ϕ(t)至少存在n+2个零点 x , x 0 , x 1 , ⋅ ⋅ ⋅ , x n x,x_0,x_1,···,x_n x,x0,x1,⋅⋅⋅,xn 由罗尔中值定理可知， ϕ ( n + 1 ) \phi^{(n+1)} ϕ(n+1)至少存在一个零点，即 ϕ ( n + 1 ) ( c ) = f ( n + 1 ) ( c ) − ( n + 1 ) ! K ( x ) , c ∈ ( a , b ) 且 c 依赖于 x \phi^{(n+1)}(c)=f^{(n+1)}(c)-(n+1)!K(x) , c∈(a,b)且c依赖于x ϕ(n+1)(c)=f(n+1)(c)−(n+1)!K(x),c∈(a,b)且c依赖于x因此有 K ( x ) = f ( n + 1 ) ( c ) ( n + 1 ) ! K_(x) = \frac{f^{(n+1)}(c)}{(n+1)!} K(x)=(n+1)!f(n+1)(c)
如果说能够算出 f ( n + 1 ) ( c ) f^{(n+1)}(c) f(n+1)(c)的话，误差的最大值就可以得到。

均差与牛顿插值多项式

拉格朗日插值法的一个缺点是当插值的节点有删减或增添时，计算要全部重新进行，因此可以考虑一种逐次生成多项式的方法。先来看n=1时的情形，这个时候所得到的仍然是一条直线，考虑用直线的点斜式表示： P 1 ( x ) = f ( x 0 ) + f ( x 1 ) − f ( x 0 ) x 1 − x 0 ( x − x 0 ) P_1(x) = f(x_0) + \frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0) P1(x)=f(x0)+x1−x0f(x1)−f(x0)(x−x0)它可以看做由 P 0 ( x ) = f ( x 0 ) P_0(x) = f(x_0) P0(x)=f(x0)经过 P 1 ( x ) = P 0 ( x ) + a 1 ( x − x 0 ) P_1(x) = P_0(x) + a_1(x-x_0) P1(x)=P0(x)+a1(x−x0)得来，其中 a 1 = f ( x 1 ) − f ( x 0 ) x 1 − x 0 a_1 = \frac{f(x_1)-f(x_0)}{x1-x0} a1=x1−x0f(x1)−f(x0) 称之为差商。
如果再往上考虑二次曲线，即 P 2 ( x ) = P 1 ( x ) + a 2 ( x − x 0 ) ( x − x 1 ) P_2(x) = P_1(x) + a_2(x-x_0)(x-x_1) P2(x)=P1(x)+a2(x−x0)(x−x1)它已经满足 P 2 ( x 1 ) = 0 ， P 2 ( x 0 ) = 0 P_2(x_1) = 0，P_2(x_0) = 0 P2(x1)=0，P2(x0)=0 还需要有 P 2 ( x 2 ) = f ( x 2 ) P_2(x_2) = f(x_2) P2(x2)=f(x2) 因此 a 2 = P 2 ( x 2 ) − P 1 ( x 2 ) ( x 2 − x 0 ) ( x 2 − x 1 ) = f ( x 2 ) − f ( x 0 ) − f ( x 1 ) − f ( x 0 ) x 1 − x 0 ( x 2 − x 0 ) ( x 2 − x 0 ) ( x 2 − x 1 ) a_2 = \frac{P_2(x_2)-P_1(x_2)}{(x_2-x_0)(x_2-x_1)}=\frac{f(x_2) - f(x_0) - \frac{f(x_1)-f(x_0)}{x1-x0}(x_2-x_0)}{(x_2-x_0)(x_2-x_1)} a2=(x2−x0)(x2−x1)P2(x2)−P1(x2)=(x2−x0)(x2−x1)f(x2)−f(x0)−x1−x0f(x1)−f(x0)(x2−x0) = f ( x 2 ) − f ( x 0 ) x 2 − x 0 − f ( x 1 ) − f ( x 0 ) x 1 − x 0 x 2 − x 1 = \frac{\frac{f(x_2)-f(x_0)}{x_2-x_0}-\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_1} =x2−x1x2−x0f(x2)−f(x0)−x1−x0f(x1)−f(x0)
可以看到 a 2 a_2 a2实际上是差商的差商。
由此扩展到更大的n，要求的n次多项式应该有如下形式： P n ( x ) = a 0 + a 1 ( x − x 0 ) + P_n(x) = a_0 + a_1(x-x_0) + Pn(x)=a0+a1(x−x0)+ a 2 ( x − x 0 ) ( x 2 − x 1 ) + ⋅ ⋅ ⋅ ⋅ + a n ( x − x 0 ) ⋅ ⋅ ⋅ ( x − x n − 1 ) a_2(x-x_0)(x_2-x_1)+ ···· + a_n(x-x_0)···(x-x_{n-1}) a2(x−x0)(x2−x1)+⋅⋅⋅⋅+an(x−x0)⋅⋅⋅(x−xn−1)
问题的关键在于如何快速高效地求解这些系数。
为此引入均差的概念：
设 f [ x 0 , x k ] = f ( x k ) − f ( x 0 ) x k − x 0 f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0} f[x0,xk]=xk−x0f(xk)−f(x0)为函数f(x)关于点 x 0 x_0 x0和 x k x_k xk的一阶均差。 f [ x 0 , x 1 , x k ] = f [ x 0 , x k ] − f [ x 0 , x 1 ] x k − x 1 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−x1f[x0,xk]−f[x0,x1]为函数f(x)的二阶均差。一般地，称 f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k ] = f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k − 2 , x k ] − f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k − 1 ] x k − x k − 1 f[x_0,x_1,···,x_k] = \frac{f[x_0,x_1,···,x_{k-2},x_k]-f[x_0,x_1,···,x_{k-1}]}{x_k-x_{k-1}} f[x0,x1,⋅⋅⋅,xk]=xk−xk−1f[x0,x1,⋅⋅⋅,xk−2,xk]−f[x0,x1,⋅⋅⋅,xk−1]
均差具有如下的性质：
(1)K阶均差可以表为函数值 f ( x 0 ) , f ( x 1 ) ⋅ ⋅ ⋅ , f ( x k ) f(x_0),f(x_1)···,f(x_k) f(x0),f(x1)⋅⋅⋅,f(xk)的线性组合，即 f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k ] = ∑ j = 0 k f ( x j ) ( x j − x 0 ) ⋅ ⋅ ⋅ ( x j − x j − 1 ) ( x j − x j + 1 ) ⋅ ⋅ ⋅ ( x j − x k ) f[x_0,x_1,···,x_k]=\sum_{j=0}^{k}\frac{f(x_j)}{(x_j-x_0)···(x_j-x_{j-1})(x_j-x_{j+1})···(x_j-x_k)} f[x0,x1,⋅⋅⋅,xk]=j=0∑k(xj−x0)⋅⋅⋅(xj−xj−1)(xj−xj+1)⋅⋅⋅(xj−xk)f(xj)

证：一阶均差显然成立。
假设对于k-1阶均差均成立，则对于k阶均差： f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k ] = f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k − 2 , x k ] − f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k − 1 ] x k − x k − 1 f[x_0,x_1,···,x_k] = \frac{f[x_0,x_1,···,x_{k-2},x_k]-f[x_0,x_1,···,x_{k-1}]}{x_k-x_{k-1}} f[x0,x1,⋅⋅⋅,xk]=xk−xk−1f[x0,x1,⋅⋅⋅,xk−2,xk]−f[x0,x1,⋅⋅⋅,xk−1] = ∑ j = 0 , j ≠ k − 1 k f ( x j ) ( x j − x 0 ) ⋅ ⋅ ⋅ ( x j − x j − 1 ) ( x j − x j + 1 ) ⋅ ⋅ ⋅ ( x k − x k ) =\sum_{j=0,j≠k-1}^{k}\frac{f(x_j)}{(x_j-x_0)···(x_j-x_{j-1})(x_j-x_{j+1})···(x_k-x_{k})} =j=0,j̸=k−1∑k(xj−x0)⋅⋅⋅(xj−xj−1)(xj−xj+1)⋅⋅⋅(xk−xk)f(xj) − ∑ j = 0 , j ≠ k k f ( x j ) ( x j − x 0 ) ⋅ ⋅ ⋅ ( x j − x j − 1 ) ( x j − x j + 1 ) ⋅ ⋅ ⋅ ( x k − x k − 1 ) -\sum_{j=0,j≠k}^{k}\frac{f(x_j)}{(x_j-x_0)···(x_j-x_{j-1})(x_j-x_{j+1})···(x_k-x_{k-1})} −j=0,j̸=k∑k(xj−x0)⋅⋅⋅(xj−xj−1)(xj−xj+1)⋅⋅⋅(xk−xk−1)f(xj) = ∑ j = 0 k f ( x j ) ( x j − x k − 1 ) ( x j − x 0 ) ⋅ ⋅ ⋅ ( x j − x j − 1 ) ( x j − x j + 1 ) ⋅ ⋅ ⋅ ( x j − x k ) ( x k − x k − 1 ) =\sum_{j=0}^{k}\frac{f(x_j)(x_j-x_{k-1})}{(x_j-x_0)···(x_j-x_{j-1})(x_j-x_{j+1})···(x_j-x_k)(x_k-x_{k-1})} =j=0∑k(xj−x0)⋅⋅⋅(xj−xj−1)(xj−xj+1)⋅⋅⋅(xj−xk)(xk−xk−1)f(xj)(xj−xk−1) − ∑ j = 0 k f ( x j ) ( x j − x k ) ( x j − x 0 ) ⋅ ⋅ ⋅ ( x j − x j − 1 ) ( x j − x j + 1 ) ⋅ ⋅ ⋅ ( x j − x k ) ( x k − x k − 1 ) -\sum_{j=0}^{k}\frac{f(x_j)(x_j-x_{k})}{(x_j-x_0)···(x_j-x_{j-1})(x_j-x_{j+1})···(x_j-x_{k})(x_k-x_{k-1})} −j=0∑k(xj−x0)⋅⋅⋅(xj−xj−1)(xj−xj+1)⋅⋅⋅(xj−xk)(xk−xk−1)f(xj)(xj−xk)
= ∑ j = 0 k f ( x j ) ( x k − x k − 1 ) ( x j − x 0 ) ⋅ ⋅ ⋅ ( x j − x j − 1 ) ( x j − x j + 1 ) ⋅ ⋅ ⋅ ( x j − x k ) =\sum_{j=0}^{k}\frac{f(x_j)(x_k-x_{k-1})}{(x_j-x_0)···(x_j-x_{j-1})(x_j-x_{j+1})···(x_j-x_{k})} =j=0∑k(xj−x0)⋅⋅⋅(xj−xj−1)(xj−xj+1)⋅⋅⋅(xj−xk)f(xj)(xk−xk−1)
这个性质表明，均差与[]内部节点排列顺序无关
(2)由性质(1)可得： f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k ] = f [ x 1 , x 2 , ⋅ ⋅ ⋅ , x k ] − f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x k − 1 ] x k − x 0 f[x_0,x_1,···,x_k]=\frac{f[x_1,x_2,···,x_k]-f[x_0,x_1,···,x_{k-1}]}{x_k-x_0} f[x0,x1,⋅⋅⋅,xk]=xk−x0f[x1,x2,⋅⋅⋅,xk]−f[x0,x1,⋅⋅⋅,xk−1]
(3) f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x n ] = f ( n ) ( c ) n ! , c ∈ [ a , b ] f[x_0,x_1,···,x_n]=\frac{f^{(n)}(c)}{n!},c∈[a,b] f[x0,x1,⋅⋅⋅,xn]=n!f(n)(c),c∈[a,b]（拉格朗日中值定理可证）

根据上述均差的定义，有 f ( x ) = f ( x 0 ) + f [ x , x 0 ] ( x − x 0 ) f(x)=f(x_0)+f[x,x_0](x-x_0) f(x)=f(x0)+f[x,x0](x−x0) f [ x , x 0 ] = f [ x 0 , x 1 ] + f [ x , x 0 , x 1 ] ( x − x 1 ) f[x,x_0]=f[x_0,x_1]+f[x,x_0,x_1](x-x_1) f[x,x0]=f[x0,x1]+f[x,x0,x1](x−x1) ⋅ · ⋅ ⋅ · ⋅ ⋅ · ⋅ f [ x , x 0 , ⋅ ⋅ ⋅ , x n − 1 ] = f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x n ] + f [ x , x 0 , ⋅ ⋅ ⋅ , x n ] ( x − x n ) f[x,x_0,···,x_{n-1}]=f[x_0,x_1,···,x_n]+f[x,x_0,···,x_n](x-x_n) f[x,x0,⋅⋅⋅,xn−1]=f[x0,x1,⋅⋅⋅,xn]+f[x,x0,⋅⋅⋅,xn](x−xn)
依次将后一项代入前一项即可得： f ( x ) = f ( x 0 ) + f [ x 0 , x 1 ] ( x − x 0 ) + f [ x 0 , x 1 , x 2 ] ( x − x 0 ) ( x − x 1 ) + ⋅ ⋅ ⋅ f(x) = f(x_0) + f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+··· f(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋅⋅⋅ + f [ x 0 , x 1 , ⋅ ⋅ ⋅ , x n ] ( x − x 0 ) ⋅ ⋅ ⋅ ( x − x n − 1 ) + f [ x , x 0 , x 1 , ⋅ ⋅ ⋅ , x n ] w n + 1 ( x ) +f[x_0,x_1,···,x_n](x-x_0)···(x-x_{n-1})+f[x,x_0,x_1,···,x_n]w_{n+1}(x) +f[x0,x1,⋅⋅⋅,xn](x−x0)⋅⋅⋅(x−xn−1)+f[x,x0,x1,⋅⋅⋅,xn]wn+1(x) = P n ( x ) + R n ( x ) =P_n(x)+R_n(x) =Pn(x)+Rn(x)
其中 P n ( x ) = f ( x ) − f [ x , x 0 , x 1 , ⋅ ⋅ ⋅ , x n ] w n + 1 ( x ) P_n(x) = f(x) - f[x,x_0,x_1,···,x_n]w_{n+1}(x) Pn(x)=f(x)−f[x,x0,x1,⋅⋅⋅,xn]wn+1(x)为最终的插值多项式，而 R n ( x ) = f [ x , x 0 , x 1 , ⋅ ⋅ ⋅ , x n ] w n + 1 ( x ) R_n(x) = f[x,x_0,x_1,···,x_n]w_{n+1}(x) Rn(x)=f[x,x0,x1,⋅⋅⋅,xn]wn+1(x)则为误差。
以上便是牛顿插值多项式的计算过程，它的程序实现较为简便，只需逐阶次计算均差即可，复杂度为 O ( n 2 ) O(n^2) O(n2)