牛顿插值法C语言实现

牛顿插值法C语言实现

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

一、插值法

设函数$y=f(x)$在区间$\left[a,b\right]$上有定义，且已知在点$a\le x_0<x_1<...<x_n\le b$上的值$y_0,y_1,...,y_n$,如果构造出来一个简单函数$P(x)$，使得$P(x_i)=y_i,i=0,1,2,...n$那么称$P(x)$为$f(x)$的插值函数，点$x_0,x_1...x_n$为插值节点，$\left[a,b\right]$为差值区间，求差值函数$P(x)$的方法称为插值法。

如果令$P(x)$为多项式，即
$$P(x)=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^n$$
其中$a_i$为实数，就称$P(x)$为插值多项式，相应的插值法称为多项式插值。

二、均差

2.1定义

称$f\left[x_0,x_k\right]=\frac{f\left(x_k\right)-f\left(x_0\right)}{x_k-x_0}$为函数$f(x)$关于点$x_0,x_k$的 一阶均差 ，$f\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(x)$的二阶均差，一般的，称
$$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}}(\ast )$$
为$f(x)$的$k$阶均差

2.2 性质

(1)$k$阶均差克表示为函数值$f(x_0),f(x_1),...,f(x_k)$的线性组合，即
$$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\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]$$
(2)由性质1和$(\ast )$式可得

$$f\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}$$
(3)若$f(x)$在$[a,b]$上存在$n$阶导数，且结点$x_0,x_1,...,x_n\in [a,b]$,则$n$阶均差与导数的关系为
$$f\left[x_0,x_1,...,x_n\right]=\frac{f^{(n)}(\xi )}{n!},\xi \in [a,b].$$

2.3均差表

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

$x_k$	$f(x_k)$	一阶均差	二阶均差	三阶均差	四阶均差
$x_0$	$f(x_0)$
$x_1$	$f(x_1)$	$f\left[x_0,x_1\right]$
$x_2$	$f(x_2)$	$f\left[x_1,x_2\right]$	$f\left[x_0,x_1,x_2\right]$
$x_3$	$f(x_3)$	$f\left[x_2,x_3\right]$	$f\left[x_1,x_2,x_3\right]$	$f\left[x_0,x_1,x_2,x_3\right]$
$x_4$	$f(x_4)$	$f\left[x_3,x_4\right]$	$f\left[x_2,x_3,x_4\right]$	$f\left[x_1,x_2,x_3,x_4\right]$	$f\left[x_0,x_1,x_2,x_3,x_4\right]$
…	…	…	…	…	…

三、牛顿插值法

3.1 差值多项式的逐次生成

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

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

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

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

$f(x)=f(x_0)+f[x,x_0](x-x_0),$

$f[x,x_0]=f[x_0,x_1]+f[x_0,x_1,x_2](x-x_1),$

$...$

$f[x,x_0,...,x_{n-1}]=f[x_0,x_1,...,x_n]+f[x,x_0,...,x_n](x-x_n).$

通过迭代得到
$$\begin{gathered} 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) \end{gathered}$$
其中
$$\begin{array}{rl}{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{array}$$

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

4.1 整体思路

我们已知了$n+1$个点$(x_0,f(x_0)),(x_1,f(x_1))...(x_n,f(x_n))$通过牛顿差值多项式的形式$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})$

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

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)=\frac{n(n-1)}{2}$个空间浪费，所以我们使用一维数组来存储。我们以上表为例

$x_k$	$f(x_k)$	一阶均差	二阶均差	三阶均差	四阶均差
$x_0$	$f(x_0)$
$x_1$	$f(x_1)$	$f\left[x_0,x_1\right]$
$x_2$	$f(x_2)$	$f\left[x_1,x_2\right]$	$f\left[x_0,x_1,x_2\right]$
$x_3$	$f(x_3)$	$f\left[x_2,x_3\right]$	$f\left[x_1,x_2,x_3\right]$	$f\left[x_0,x_1,x_2,x_3\right]$
$x_4$	$f(x_4)$	$f\left[x_3,x_4\right]$	$f\left[x_2,x_3,x_4\right]$	$f\left[x_1,x_2,x_3,x_4\right]$	$f\left[x_0,x_1,x_2,x_3,x_4\right]$
…	…	…	…	…	…

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

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

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

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

我们根据差商定义有
$$f[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[x_k,x_{k+1},x_{k+2}]$是第$k$个元素，则 $f[x_{k+1},x_k]$是第$k-n$个元素，$f[x_{k+2},x_{k+1}]$是第$k-n+1$个元素；

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

归纳得

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

所以我们可以通过一个二重for循环来得出所有差商表的数，循环变量$i$代表差商阶数减1，循环变量$j$代表本阶差商的第几个数，$i+1$阶差商有$n-(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 输出均差表

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

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

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

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

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 根据牛顿插值多项式预测下一个节点的值

先在程序里得到预测的横坐标$a$，然后我们在创建一个有$n$个元素的数组来存$P_n(x)$除了差商和$f(x_0)$的部分
$$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})$$

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

先算$(a-x_0)$，$(a-x_0)(a-x_1)$，…，$(a-x_0)(a-x_1)...(a-x_{n-2})$

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

然后再算$P_n(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 x$为例子

得到的$\ln 0.54\approx -0.616143$