插值（三）——Newton插值（C++）

差商

对于一系列互不相等的$x_0,x_1,\cdots ,x_n$，有差商定义如下：
$$f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}(i\ne j,i,j=0,1,\cdots ,n)$$
称此为$f(x)$在$x_i,x_j$处的一阶差商，又称
$$f(x_i,x_j,x_k)=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}$$
为$f(x)$在$x_i,x_j,x_k$处的二阶差商。
因此可以得到n阶差商的定义：
$$f[x_0,x_1,\cdots ,x_n]=\frac{f[x_0,x_1,\cdots ,x_{n-1}]-f[x_1,x_2,\cdots ,x_n]}{x_0-x_n}$$
差商的计算方法如下：

牛顿插值法

与拉格朗日插值法不同，如果需要增加插值点，拉格朗日插值法需要重新计算所有插值点的值，重新计算基函数，而牛顿插值法只需要重附上一项即可，不需要重新计算。
牛顿插值法需要利用到差商，其插值公式如下：
$$\begin{array}{rl}{N}_n(x)=& f(x_0)+f[x_0,x_1]\ast (x-x_0)+f[x_0,x_1,x_2]\ast (x-x_0)(x-x_1)+\\ & f[x_0,x_1,x_2,x_3]\ast (x-x_0)(x-x_1)(x-x_2)+\cdots +\\ & f[x_0,x_1,\cdots ,x_n]\ast (x-x_0)(x-x_1)\cdots (x-x_{n-1})\end{array}$$

代码

//牛顿插值法
#include<iostream>
#include<cmath>
//求差商
void poor_quotient(double *xi,double *yi,double *Ni,int n)
{
    double **a_temp=new double*[n];
    for (int i=0; i < n;i++)
    {
        a_temp[i]=new double[n];
        a_temp[i][0]=yi[i];
    }
    for (int j = 1; j < n;j++)
    {
        for (int i=j; i < n;i++)
        {
            a_temp[i][j]=(a_temp[i][j-1]-a_temp[i-1][j-1])/(xi[i]-xi[i-j]);
        }
    }
        for (int i = 0; i < n; i++)
        {
            Ni[i] = a_temp[i][i];
        }
    for (int i = 0; i < n;i++)
    {
        delete[] a_temp[i];
    }
    delete[] a_temp;
}
double Newton(double *Ni,double *xi,double x,int n)
{
    double result =Ni[0];
    double temp = 1.0f;
    for (int i = 1;i<n;i++)
    {
        temp *= x - xi[i-1];
        result += Ni[i] * temp;
    }
    return result;
}
double fx(double x)
{
    //这里填入被插值函数
    //如果插值区间包括分母为0的情况需要手动调整
    return sqrt(x + 3);
}
int main()
{
    int n;//插值点数
    std::cout<<"请输入插值点数：";
    std::cin>>n;
    double error[3]={0.0f,0.0f,0.0f};//误差评价
    double *x = new double [n];
    double *y = new double [n];
    double *Ni = new double [n];
    double a=3, b=10;//定义插值区间
    //生成插值数据
    for (int i = 0;i<n;i++)
    {
        x[i] = a + (b-a)*i/(n-1);
        y[i] = fx(x[i]);
    }

    //计算差商
    poor_quotient(x, y, Ni, n);
    //std::cout << std::endl;
    // 验证误差，抽取区间内的100个点进行误差检测
    for (int i = 0; i < 100; i++)
    {
        double x_temp = a + (b - a) * i / (100 - 1);
        double y_temp = fx(x_temp);
        double y_temp2 = Newton(Ni, x, x_temp, n);
        if (fabs(y_temp - y_temp2) > error[0])
        {
            error[0] = fabs(y_temp - y_temp2);
            // std::cout << y_temp2 << std::endl;
        }
        error[1] += fabs((y_temp - y_temp2) / y_temp);
        error[2] += ((y_temp - y_temp2)) * ((y_temp - y_temp2));
    }
    //err[0]得到的是在区间内最大的误差
    //err[1]得到的是平均最大相对误差
    //err[2]是均方根误差
    error[1] = error[1]/100;
    error[2] = sqrt(error[2]/100);
    for(int i = 0;i < 3;i++)
    {
        std::cout<<"误差"<<i+1<<"："<<error[i]<<std::endl;
    }
    delete[] x;
    delete[] y;
    delete[] Ni;
    return 0;
}

结果分析

与多项式插值、拉格朗日插值相比，对于相同的测试函数，其具有更高的精度，并且其龙格现象也更加不明显。