牛顿插值法C语言实现

请注意：在数学推导过程和解释中的插值法都是用的n+1n+1n+1个结点来推导，但是程序中是以nnn个结点来编写，数学推导中的nnn是多项式的次数，程序中的nnn是结点的个数

一、插值法

设函数y=f(x)y=f(x)y=f(x)在区间[a,b]\left[a,b\right][a,b]上有定义，且已知在点a≤x0<x1<...<xn≤ba\le x_0<x_1<...<x_n\le ba≤x0<x1<...<xn≤b上的值y0,y1,...,yny_0,y_1,...,y_ny0,y1,...,yn,如果构造出来一个简单函数P(x)P(x)P(x)，使得P(xi)=yi,i=0,1,2,...nP(x_i)=y_i,i=0,1,2,...nP(xi)=yi,i=0,1,2,...n那么称P(x)P(x)P(x)为f(x)f(x)f(x)的插值函数，点x0,x1...xnx_0,x_1...x_nx0,x1...xn为插值节点，[a,b]\left[a,b\right][a,b]为差值区间，求差值函数P(x)P(x)P(x)的方法称为插值法。

如果令P(x)P(x)P(x)为多项式，即
P(x)=a0+a1x+a2x2+...+anxnP(x)=a_0+a_1x+a_2x^2+...+a_nx^n P(x)=a0+a1x+a2x2+...+anxn
其中aia_iai为实数，就称P(x)P(x)P(x)为插值多项式，相应的插值法称为多项式插值。

二、均差

2.1定义

称f[x0,xk]=f(xk)−f(x0)xk−x0f \left [x_0,x_k \right] = \frac{f \left( x_k \right) - f \left( x_0 \right)}{x_k-x_0}f[x0,xk]=xk−x0f(xk)−f(x0)为函数f(x)f(x)f(x)关于点x0,xkx_0,x_kx0,xk的 一阶均差 ，f[x0,x1,xk]=f[x0,xk]−f[x1,xk]xk−x1f\left[x_0,x_1,x_k\right]=\frac {f\left[x_0,x_k\right] - f\left[x_1,x_k\right]} {x_k-x_1}f[x0,x1,xk]=xk−x1f[x0,xk]−f[x1,xk]为f(x)f(x)f(x)的二阶均差，一般的，称
f[x0,x1,...,xk]=f[x0,x1,...,xk]−f[x0,x1,...,xk−1]xk−xk−1(∗)f \left[ x_0,x_1,...,x_k \right]=\frac {f \left[ x_0,x_1,...,x_k \right] - f \left[ x_0,x_1,...,x_{k-1} \right]} {x_k-x_{k-1}} \ \ \ \ \ \ \ \ \ \ (*) f[x0,x1,...,xk]=xk−xk−1f[x0,x1,...,xk]−f[x0,x1,...,xk−1] (∗)
为f(x)f(x)f(x)的kkk阶均差

2.2 性质

(1)kkk阶均差克表示为函数值f(x0),f(x1),...,f(xk)f(x_0),f(x_1),...,f(x_k)f(x0),f(x1),...,f(xk)的线性组合，即
f[x0,x1,...,xk]=∑j=0kf(xj)(xj−x0)...(xj−xj−1)(xj−xj+1)...(xj−xk)f \left[ x_0,x_1,...,x_k \right]=\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)
这个性质表明均差与结点的排列次序无关，称为均差的对称性，即
f[x0,x1,...,xk]=f[x1,x0,...,xk]=...=f[x1,x2,...,xk,x0]f \left[ x_0,x_1,...,x_k \right]=f \left[ x_1,x_0,...,x_k \right]=...=f \left[ x_1,x_2,...,x_k,x_0 \right] f[x0,x1,...,xk]=f[x1,x0,...,xk]=...=f[x1,x2,...,xk,x0]
(2)由性质1和(∗)(*)(∗)式可得

f[x0,x1,...,xk]=f[x1,x2,...,xk]−f[x0,x1,...,xk−1]xk−x0f \left[ x_0,x_1,...,x_k \right]=\frac {f\left[ x_1,x_2,...,x_k \right] - f\left[ x_0,x_1,...,x_{k-1} \right]} {x_k-x_0} f[x0,x1,...,xk]=xk−x0f[x1,x2,...,xk]−f[x0,x1,...,xk−1]
(3)若f(x)f(x)f(x)在[a,b][a,b][a,b]上存在nnn阶导数，且结点x0,x1,...,xn∈[a,b]x_0,x_1,...,x_n \in [a,b]x0,x1,...,xn∈[a,b],则nnn阶均差与导数的关系为
f[x0,x1,...,xn]=f(n)(ξ)n!,ξ∈[a,b].f\left[ x_0,x_1,...,x_n \right]=\frac {f^{(n)}(\xi)}{n!},\ \ \ \xi\in[a,b]. f[x0,x1,...,xn]=n!f(n)(ξ), ξ∈[a,b].

2.3均差表

均差的计算可以使用均差表

xkx_kxk	f(xk)f(x_k)f(xk)	一阶均差	二阶均差	三阶均差	四阶均差
x0x_0x0	f(x0)f(x_0)f(x0)
x1x_1x1	f(x1)f(x_1)f(x1)	f[x0,x1]f\left[x_0,x_1 \right]f[x0,x1]
x2x_2x2	f(x2)f(x_2)f(x2)	f[x1,x2]f\left[x_1,x_2 \right]f[x1,x2]	f[x0,x1,x2]f\left[ x_0,x_1,x_2 \right]f[x0,x1,x2]
x3x_3x3	f(x3)f(x_3)f(x3)	f[x2,x3]f\left[x_2,x_3 \right]f[x2,x3]	f[x1,x2,x3]f\left[ x_1,x_2,x_3 \right]f[x1,x2,x3]	f[x0,x1,x2,x3]f\left[ x_0,x_1,x_2,x_3 \right]f[x0,x1,x2,x3]
x4x_4x4	f(x4)f(x_4)f(x4)	f[x3,x4]f\left[x_3,x_4 \right]f[x3,x4]	f[x2,x3,x4]f\left[ x_2,x_3,x_4 \right]f[x2,x3,x4]	f[x1,x2,x3,x4]f\left[ x_1,x_2,x_3,x_4 \right]f[x1,x2,x3,x4]	f[x0,x1,x2,x3,x4]f\left[x_0,x_1,x_2,x_3,x_4\right]f[x0,x1,x2,x3,x4]
…	…	…	…	…	…

三、牛顿插值法

3.1 差值多项式的逐次生成

当n=1n=1n=1时，有两个点，这时的插值多项式P(x)P(x)P(x)满足P(x0)=f(x0),P(x1)=f(x1)P(x_0)=f(x_0),P(x_1)=f(x_1)P(x0)=f(x0),P(x1)=f(x1)，且P(x)=a0+a1xP(x)=a_0+a_1xP(x)=a0+a1x，带入得
P1(x)=f(x0)+f(x1)−f(x0)x1−x0(x−x0)P_1(x)=f(x_0)+\frac{f(x_1)-f(x_0)} {x_1-x_0}\left(x-x_0\right) P1(x)=f(x0)+x1−x0f(x1)−f(x0)(x−x0)
这可以看成P0(x)=f(x0)P_0(x)=f(x_0)P0(x)=f(x0)的修正，即
P1(x)=P0(x)+a1(x−x0)P_1(x)=P_0(x)+a_1(x-x_0) P1(x)=P0(x)+a1(x−x0)
其中a1=f(x1)−f(x0)x1−x0a_1=\frac{f(x_1)-f(x_0)}{x_1-x_0}a1=x1−x0f(x1)−f(x0)是函数f(x)f(x)f(x)的一阶差商。

当n=2n=2n=2时，有三个已知点，插值多项式P(x)P(x)P(x)满足P(x0)=f(x0),P(x1)=f(x1),P(x2)=f(x2)P(x_0)=f(x_0),P(x_1)=f(x_1),P(x_2)=f(x_2)P(x0)=f(x0),P(x1)=f(x1),P(x2)=f(x2)

可表示为
P2(x)=P1(x)+a2(x−x0)(x−x1)P_2(x)=P_1(x)+a_2(x-x_0)(x-x_1) P2(x)=P1(x)+a2(x−x0)(x−x1)
显然P2(x)P_2(x)P2(x)满足上述条件，将P(x2)=f(x2)P(x_2)=f(x_2)P(x2)=f(x2)带入得到
a2=P2(x2)−P1(x2)(x2−x0)(x2−x1)=f(x2)−f(x0)x2−x0−f(x1)−f(x0)x1−x0x2−x1a_2=\frac {P_2(x_2)-P_1(x_2)}{(x_2-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} a2=(x2−x0)(x2−x1)P2(x2)−P1(x2)=x2−x1x2−x0f(x2)−f(x0)−x1−x0f(x1)−f(x0)
所以a2=f[x2,x0]−f[x1,x0]x2−x1a_2=\frac{f[x_2,x_0]-f[x_1,x_0]}{x_2-x_1}a2=x2−x1f[x2,x0]−f[x1,x0]是函数f(x)f(x)f(x)的二阶差商

实际上，根据均差定义，将xxx看成[a,b][a,b][a,b]上的一个点，可得

f(x)=f(x0)+f[x,x0](x−x0),f(x)=f(x_0)+f[x,x_0](x-x_0),f(x)=f(x0)+f[x,x0](x−x0),

f[x,x0]=f[x0,x1]+f[x0,x1,x2](x−x1),f[x,x_0]=f[x_0,x_1]+f[x_0,x_1,x_2](x-x_1),f[x,x0]=f[x0,x1]+f[x0,x1,x2](x−x1),

.........

f[x,x0,...,xn−1]=f[x0,x1,...,xn]+f[x,x0,...,xn](x−xn).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(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)+f[x,x0,...,xn]ωn+1(x)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_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})+f[x,x_0,...,x_n]\omega_{n+1}(x) f(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)+f[x,x0,...,xn]ωn+1(x)
其中
Pn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)R(x)=f(x)−Pn(x)=f[x,x0,...,xn]ωn+1(x)\begin{aligned} P_n(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_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1}) \\ R_{\ }(x) &= f(x)-P_n(x)=f[x,x_0,...,x_n]\omega_{n+1}(x) \end{aligned} Pn(x)R (x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)=f(x)−Pn(x)=f[x,x0,...,xn]ωn+1(x)

四、牛顿插值法的C语言实现

4.1 整体思路

我们已知了n+1n+1n+1个点(x0,f(x0)),(x1,f(x1))...(xn,f(xn))(x_0,f(x_0)),(x_1,f(x_1))...(x_n,f(x_n))(x0,f(x0)),(x1,f(x1))...(xn,f(xn))通过牛顿差值多项式的形式Pn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)P_n(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_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})Pn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)

我们未知的数据只有差商值，所以我们输入数据之后将各阶均差算出来就可以得到牛顿差值多项式。

4.2 初始数据

4.2.1变量说明

double* x;   //自变量值的集合
double* y;  //函数值的集合
double* f;  //插值系数，
int n;      //插值节点个数

4.2.2 输入各个变量

首先我们要确定插值节点的个数，所以先在main函数中输入变量n的值

printf("输入要插值节点的个数:");scanf("%d", &n);

然后我们要为x指针和y指针开辟内存空间，x和y指针开n个内存空间

x = (double*)malloc(sizeof(double) * n );
y = (double*)malloc(sizeof(double) * n );

再输入各个节点的值，在这里我们创建一个data函数，里面采用for循环

void data(double* x, double* y, int n){for(int i=0;i<n+1;i++) {cout << "x[" << i << "]:";cin >> x[i];cout << "y[" << i << "]:";cin >> y[i];}
}

现在我们就得到了三个变量（x,y,n）并且都已赋值

4.3 计算各阶均差和数据的存储

这里先介绍差商的存储，因为差商可以表示成一个表，所以用矩阵或二维数组存储是最方便的，但是考虑到如果有n个结点，生成的矩阵就有1+2+...+(n−1)=n(n−1)21+2+...+(n-1)=\frac{n(n-1)}21+2+...+(n−1)=2n(n−1)个空间浪费，所以我们使用一维数组来存储。我们以上表为例

xkx_kxk	f(xk)f(x_k)f(xk)	一阶均差	二阶均差	三阶均差	四阶均差
x0x_0x0	f(x0)f(x_0)f(x0)
x1x_1x1	f(x1)f(x_1)f(x1)	f[x0,x1]f\left[x_0,x_1 \right]f[x0,x1]
x2x_2x2	f(x2)f(x_2)f(x2)	f[x1,x2]f\left[x_1,x_2 \right]f[x1,x2]	f[x0,x1,x2]f\left[ x_0,x_1,x_2 \right]f[x0,x1,x2]
x3x_3x3	f(x3)f(x_3)f(x3)	f[x2,x3]f\left[x_2,x_3 \right]f[x2,x3]	f[x1,x2,x3]f\left[ x_1,x_2,x_3 \right]f[x1,x2,x3]	f[x0,x1,x2,x3]f\left[ x_0,x_1,x_2,x_3 \right]f[x0,x1,x2,x3]
x4x_4x4	f(x4)f(x_4)f(x4)	f[x3,x4]f\left[x_3,x_4 \right]f[x3,x4]	f[x2,x3,x4]f\left[ x_2,x_3,x_4 \right]f[x2,x3,x4]	f[x1,x2,x3,x4]f\left[ x_1,x_2,x_3,x_4 \right]f[x1,x2,x3,x4]	f[x0,x1,x2,x3,x4]f\left[x_0,x_1,x_2,x_3,x_4\right]f[x0,x1,x2,x3,x4]
…	…	…	…	…	…

从一阶均差开始，f[x0,x1]f\left[x_0,x_1 \right]f[x0,x1]作为f[0]f[\ 0\ ]f[ 0 ]，f[x1,x2]f\left[x_1,x_2 \right]f[x1,x2]作为f[1]f[\ 1\ ]f[ 1 ] …一阶均差有n-1个数，二阶均差有n-2个数 … n-1阶均差有1个数。所以一共需要1+2+...+(n−1)=n(n−1)21+2+...+(n-1)=\frac {n(n-1)}21+2+...+(n−1)=2n(n−1)个空间，所以我们申请n(n−1)2\frac{n(n-1)}22n(n−1)个空间

f = (double*)malloc(sizeof(double) * (n - 1) * n / 2);

我们现在有了n(n−1)2\frac {n(n-1)}22n(n−1)个内存空间，前n-1个为一阶差商，所以我们可以通过一个for循环来算出来所有一阶差商

for(j=0,k=0;i<n;j++){f[k++]=( y[j + 1] - y[j] ) / ( x[j + 1] - x[j])
}

我们根据差商定义有
f[xk,xk+1,xk+2]=f[xk+2,xk+1]−f[xk+1,xk]xk+2−xkf[x_k,x_{k+1},x_{k+2}]=\frac {f[x_{k+2},x_{k+1}]-f[x_{k+1},x_k]} {x_{k+2}-x_k} f[xk,xk+1,xk+2]=xk+2−xkf[xk+2,xk+1]−f[xk+1,xk]
假设f[xk,xk+1,xk+2]f[x_k,x_{k+1},x_{k+2}]f[xk,xk+1,xk+2]是第kkk个元素，则 f[xk+1,xk]f[x_{k+1},x_{k}]f[xk+1,xk]是第k−nk-nk−n个元素，f[xk+2,xk+1]f[x_{k+2},x_{k+1}]f[xk+2,xk+1]是第k−n+1k-n+1k−n+1个元素；

假设f[xk,xk+1,xk+2,xk+3]f[x_k,x_{k+1},x_{k+2},x_{k+3}]f[xk,xk+1,xk+2,xk+3]是第kkk个元素，则f[xk,xk+1,xk+2]f[x_k,x_{k+1},x_{k+2}]f[xk,xk+1,xk+2]是第k−n+1k-n+1k−n+1个元素f[xk+1,xk+2,xk+3]f[x_{k+1},x_{k+2},x_{k+3}]f[xk+1,xk+2,xk+3]是第k−n+2k-n+2k−n+2个元素

归纳得

假设i+1i+1i+1阶差商第一个元素f[xk,xk+1,...xk+i]f[x_k,x_{k+1},...x_{k+i}]f[xk,xk+1,...xk+i]是第kkk个元素，则f[xk,xk+1,...xk+i−1]f[x_k,x_{k+1},...x_{k+i-1}]f[xk,xk+1,...xk+i−1]是第k−n+(i+1)−2=k−n+i−1k-n+(i+1)-2=k-n+i-1k−n+(i+1)−2=k−n+i−1个元素，f[xk+1,xk+2,...xk+i]f[x_{k+1},x_{k+2},...x_{k+i}]f[xk+1,xk+2,...xk+i]是第k−n+ik-n+ik−n+i个元素

所以我们可以通过一个二重for循环来得出所有差商表的数，循环变量iii代表差商阶数减1，循环变量jjj代表本阶差商的第几个数，i+1i+1i+1阶差商有n−(i+1)=n−i−1n-(i+1)=n-i-1n−(i+1)=n−i−1个数，所以循环可以写出。

void newton(double* x, double* y, double* f, int n)
{ int i = 0, j, k = 0;for (i = 0; i < n; i++)for (j = 0; j < n - i; j++) {if (i == 0)f[k++] = (y[j + 1] - y[j]) / (x[j + 1] - x[j]);else {f[k] = (f[k + i - n] - f[k + i - n - 1]) / (x[j + i + 1] - x[j]);k++;}}
}

4.4 输出均差表

fff的前n−1n-1n−1个元素是一阶差商，接下来的n−2n-2n−2个元素是二阶差商，…

所以我们只需要将xxx指针中的元素循环一次，每次循环输出一行数据。

第一行的数据有x[0],y[0]x[0],y[0]x[0],y[0]，第二行数据有x[1],y[1],f[0]x[1],y[1],f[0]x[1],y[1],f[0]。

所以第i+1i+1i+1行数据有i+2i+2i+2个数，有iii个差商的数据，所以我们可以用一个nnn次的循环，先输出x[i]x[i]x[i]和y[i]y[i]y[i]，再用一个iii次循环来输出差商数据。

void printnew(double* x, double* y, double* f, int n)
{ /*表头*/int i, j, k = 0;printf("差商表:\n");printf("x\t  ");for (i = 0; i < n; i++)printf("f(x%d)\t\t", i);printf("\n");for (i = 0; i < n; i++)printf("----------------");printf("\n");/*表头部分结束*/for (i = 0; i < n; i++) {printf("%-10g  %-10g", x[i], y[i]);k = i - 1;for (j = 0; j < i; j++) {printf("     %-10g", f[k]);k += n - 2 - j;}if (j == i)printf("\n");}
}

4.5 根据牛顿插值多项式预测下一个节点的值

先在程序里得到预测的横坐标aaa，然后我们在创建一个有nnn个元素的数组来存Pn(x)P_n(x)Pn(x)除了差商和f(x0)f(x_0)f(x0)的部分
Pn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)P_n(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_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1}) Pn(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+...+f[x0,x1,...,xn](x−x0)(x−x1)...(x−xn−1)

b = (double*)malloc(sizeof(double) * n);

先算(a−x0)(a-x_0)(a−x0)，(a−x0)(a−x1)(a-x_0)(a-x_1)(a−x0)(a−x1)，…，(a−x0)(a−x1)...(a−xn−2)(a-x_0)(a-x_1)...(a-x_{n-2})(a−x0)(a−x1)...(a−xn−2)

b[0]=1;
for (i = 0; i < n - 1; i++)b[i + 1] = b[i] * (a - x[i]);

然后再算Pn(a)P_n(a)Pn(a)

for (i = 0; i < n; i++) {if (i == 0)a = y[0];else {a += b[i] * f[k];k += n - i;}
}

最后输出

printf("插值节点对应的Y的值:%g\n", a);

所以我们得到了完整的函数体

void getValue(double* x, double* y, double* f, int n)
{ double a, * b;int i, k = 0;b = (double*)malloc(sizeof(double) * n);printf("输入要要插入节点的X的值:");cin >> a;b[0] = 1.0;for (i = 0; i < n - 1; i++)b[i + 1] = b[i] * (a - x[i]);for (i = 0; i < n; i++) {if (i == 0)a = y[0];else {a += b[i] * f[k];k += n - i;}}printf("插值节点对应的Y的值:%g\n", a);
}

4.6编写main函数，完成程序

我们现在有四个函数

void data(double* x, double* y, int n);
void newton(double* x, double* y, double* f, int n);
void printnew(double* x, double* y, double* f, int n);
void newvalue(double* x, double* y, double* f, int n);

他们的功能分别为：输入数据，计算系数，输出差商表，预测结点的值。

所以我们可以通过这几个函数来编写main函数

int main()
{int n;double* x, * y, * f;printf("输入要插值节点的个数:");scanf("%d", &n);x = (double*)malloc(sizeof(double) * n);y = (double*)malloc(sizeof(double) * n);f = (double*)malloc(sizeof(double) * (n - 1) * n / 2);data(x, y, n);newton(x, y, f, n - 1);printnew(x, y, f, n);do {newvalue(x, y, f, n);} while (1);return 0;
}

4.7测试程序

我们运行程序，用y=ln⁡xy=\ln xy=lnx为例子

得到的ln⁡0.54≈−0.616143\ln 0.54\approx-0.616143ln0.54≈−0.616143