插值（一）——多项式插值（C++）

插值

插值的作用是可以将原本比较难计算的函数转换为误差在一定范围内的多项式，比如在单片机中直接计算 $\sqrt{x}、\log_2x$ 之类的函数是比较麻烦的，但是使用插值的方法就可以将其转换为误差可控的只有乘法和加减法的多项式，从而简化计算。

多项式插值

多项式插值是利用多项式来拟合一系列离散的数据点，从而达到简化计算的目的。本文主要介绍最“暴力”的插值方法。

设定所需构造的插值多项式：
$D_n(x)=a_0+a_1x+a_2x^2+\cdots +a_nx^n$
由插值条件可得：
$D_n(x_i)=y_i,i=0,1,\cdots n$
由此可以得到对应的方程组：
$\begin{cases} 1a_0+a_1x_0+a_2x_0^2+\cdots +a_nx_0^n=y_0\\ 1a_0+a_1x_1+a_2x_1^2+\cdots +a_nx_1^n=y_1\\ \cdots\\ 1a_0+a_1x_n+a_2x_n^2+\cdots +a_nx_n^n=y_n \end{cases}$
用于插值规定了输入的 $x_i$ 要互不相同，因此得到的行列式D具有唯一性，且为范德蒙德行列式。
$D=\begin{vmatrix} 1&x_0&x_0^2&\cdots &x_0^n\\ 1&x_1&x_1^2&\cdots &x_1^n\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&x_n&x_n^2&\cdots &x_n^n \end{vmatrix}=\prod_{0\le j\le i \le n}(x_i-x_j)$
不难看出这个行列式唯一且不为0。
只需要解出方程组即可，比如克拉默法则：
$D_i=\begin{vmatrix} 1&x_0&\cdots &x_0^{i-1}& y_0 &x_0^{i+1}&\cdots &x_0^n\\ 1&x_1&\cdots &x_1^{i-1}& y_1 &x_1^{i+1}&\cdots &x_1^n\\ \vdots&\vdots&\cdots&\vdots&\vdots&\vdots&\cdots&\vdots\\ 1&x_n&\cdots &x_n^{i-1}& y_n &x_n^{i+1}&\cdots &x_n^n\\ \end{vmatrix}$
然后 $a_i=\frac{D_i}{D}$ 即可解出所有的系数

代码

代码里面有部分求行列式值的算法是参考此博客的，使用了里面的使用代数余子式分解的方法求解对应的值，因为本文只是讨论本方法的插值效果，故没有考虑运行速度，仅是为了简化非主要部分的代码才采用代数余子式求行列式的值的方法，如果需要提升运行速度，可以选择该博客里的另一个方法进行替代。

//多项式插值
#include<iostream>
#include<cmath>
//使用代数余子式进行求解
double determinant_value(double **D,int n)
{
    //递归终点
    if(n==1)
    {
        return  D[0][0];
    }
    else if(n==2)
    {
        return D[1][1]*D[0][0]-D[0][1]*D[1][0];
    }
    else{
        double D_value=0;
        for(int k=0;k<n;k++)
        {
            double **D_temp=new double *[n-1];
            int i2=1;//由于是根据第0行第j列展开，所以原本的行列式直接从第一行开始拷贝
            for(int i1=0;i1<n-1;i1++)
            {
                D_temp[i1]=new double[n-1];
                int j2=0;
                for(int j1=0;j1<n-1;j1++)
                {
                    if(j2==k)
                    {
                        j2++;
                    }//去除第j列
                    D_temp[i1][j1]=D[i2][j2++];
                }
                i2++;
            }
            D_value+=(((k+2)&0x0001)?(-1):1)*D[0][k]*determinant_value(D_temp,n-1);//计算代数余子式与对应项的乘积
            for(int i=0;i<n-1;i++)
            {
                delete[] D_temp[i];
            }
            delete[] D_temp;
        }
        return D_value;
    }
}
//多项式插值系数计算，输入插值坐标x,y，存储系数Di的数组和插值点数n
void polynomial_interpolation(double *x,double *y,double *D_i,int n)
{
    double **D=new double*[n];
    double D_value;
    double **D_temp=new double*[n];
    //其实这里求D可以直接利用范德蒙德行列式的性质求解
    for(int i = 0;i < n;i++)
    {
        D[i]=new double[n];
        D_temp[i]=new double[n];
        for(int j = 0;j < n;j++)
        {
            if(j==0)
            {
                D[i][j]=1;
            }
            else{
                D[i][j]=D[i][j-1]*x[i];
            }
        }
    }
    D_value=determinant_value(D,n);
    //下面是为了求D_i，每次只需要修改两列数据
    for (int i = 0; i < n;i++)
    {
        if(i==0)
        {
            for (int i1 = 0; i1 < n;i1++)
            {
                D_temp[i1][0]=y[i1];
                for(int j1 = 1;j1 < n;j1++)
                {
                    D_temp[i1][j1]=D[i1][j1];
                }
            }
        }
        else{
            for(int i2 = 0;i2 < n;i2++)
            {
                D_temp[i2][i-1]=D[i2][i-1];
                D_temp[i2][i]=y[i2];
            }
        }
        //求多项式系数
        D_i[i] = determinant_value(D_temp, n) / D_value;
    }

    for (int i = 0; i < n; i++)
    {
        delete[] D[i];
        delete[] D_temp[i];
    }
    delete[] D;
    delete[] D_temp;
}
double fx(double x)
{
    //这里填入被插值函数
    //如果插值区间包括分母为0的情况需要手动调整
    return sqrt(x + 3);
}
//利用得到的多项式系数计算函数值
double Dx(double x,int n,double *D_i)
{
    double x_temp = 1.0;
    double y_temp = 0;
    for(int j = 0;j < n;j++)
    {
        y_temp += D_i[j]*x_temp;
        x_temp*= x;
    }
    return y_temp;
}
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 *D_i = 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]);
    }
    // std::cout<<"请输入插值点坐标x："<<std::endl;
    // for(int i = 0;i < n;i++)
    // {
    //     std::cin>>x[i];
    // }
    // std::cout<<"请输入插值点坐标y："<<std::endl;
    // for(int i = 0;i < n;i++)
    // {
    //     std::cin>>y[i];
    // }
    polynomial_interpolation(x, y, D_i, n);
    std::cout<<"插值多项式为：y(x)="<<D_i[0];
    for(int i = 1;i < n;i++)
    {
        if(D_i[i]>0)
        {
            std::cout << "+"<<D_i[i] <<"x^"<<i;
        }
        else{
            std::cout << D_i[i] <<"x^"<<i;
        }
        
    }
    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 = Dx(x_temp, n, D_i);
        if(fabs(y_temp-y_temp2)>error[0])
        {
            error[0] = fabs(y_temp-y_temp2);
        }
        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;
    }
    //std::cout<<fx(3.4)<<"  "<<Dx(3.4, n, D_i)<<std::endl;
    delete[] x;
    delete[] y;
    delete[] D_i;
    return 0;
}

结果与分析

运行上面的代码如下：
在这里插入图片描述

在这里插入图片描述

这里的误差1为区间内最大的误差，误差2是平均相对误差，误差3是均方根误差，这里的误差评判标准仅作参考。
不难发现当插值点数为3和5时误差都非常小，可以根据精度需求进行选取，但是当插值点数为10时，误差将会非常大，所以选取合适的插值点数是十分必要的，否则会因为龙格现象，导致插值函数在插值点附近的值波动很大，从而导致误差很大。
这种插值方法是非常不建议使用的，因为随着阶数增长其计算量会非常大，在我的电脑上3和5都是瞬间计算完成，但是到了10就需要几秒了，到了12就算很久都算不出来