数值计算（二）之插值法与线性回归（拉格朗日插值法，牛顿插值法，赫米特插值法，最小二乘法）

插值法

插值法主要解决的问题就是，用一个多项式函数来逼近原函数，或者用多项式函数来拟合离散数据，在计算机图形的处理中，插值法应用广泛

插值法基本的有两种，拉格朗日插值法和牛顿插值法

还有一种要求更为严格的差值方法赫米特（Hiemite）插值法

插值法用多项式来逼近原函数，可以证明给定N个插值节点，只能构造唯一的，最高次不高于n的差值多项式

一般我们可以使用待定系数法，列出若干个方程，来求解系数，得到差值多项式。

拉格朗日插值法则构造一个差值基函数，保证这一部分在X取X0的时候等于1，在其他的时候等于0，不让其他的点干扰

基本形式为

插值基函数为

牛顿插值法则是引入差商的概念

关于差商的计算，我们可以使用递推的方式

最后牛顿插值法的基本形式为

拉格朗日插值法不具有承袭性，当增加了一个新的插值节点，我们需要全部重新计算

而牛顿插值法，改进了这一点，以上情况出现时，只需要在后面增加多一项即可。

拉格朗日插值法

double Lagrange(double *p,int n,double x)
{//P是插值节点数组，n是插值节点格式，x是自变量的取值
    int i,j,k;
    double numerator,denominator,ans;
    ans=0;//numerator分子，denominator分母
    for(k=0;k<n;k++)
    {
        numerator=1;
        denominator=1;
        for(i=0;i<n;i++)
        {
            if(i!=k)
            {
                numerator*=(x-(*(p+i*2)));
                denominator*=( (*(p+k*2)) - (*(p+i*2)) );
            }//求插值基函数值
            else
                continue;
        }//End for-i
        ans=ans+( (numerator/denominator) * (*(p+k*2+1)) );//插值基函数*函数值Yi
    }//End for-k
    return ans;
}

牛顿插值法

double Newton(double *p,int n,double x)
{//P是插值节点数组，n是插值节点格式，x是自变量的取值
    int i,j,k;
    double numerator,denominator,ans=0;
    double *np=(double *)malloc(n*sizeof(double));//用以保存差商
    //使用滚动数组保存插值，节约存储空间
    double ny=(x-(*p));//x-x0
    for(i=0;i<n;i++)//差商数组初始化
        *(np+i)= *(p+i*2+1);
    ans=(*(np+n-1));//一阶差商即为函数值
    for(k=1;k<n;k++)
    {
        for(i=n-1;i>=k;i--)
        {//从后往前求差商，保存在np数组内
            numerator=( *(np+i-1) ) - ( *(np+i));
            denominator=( *(p+(i-k)*2)) - ( *(p+i*2) );
            *(np+i)=(numerator/denominator);
        }//End for-i
        ans+=(*(np+n-1))*ny;
        ny*=(x-(*(p+k*2) ));//累乘(x-xi)的结果
    }//End for-k
    return ans;
}

Hiemite插值法

赫米特插值法在前两种插值法的基础上，不仅要求在插值节点处函数值相同，而且导数值也相同，比较常见的就是二点三次插值，即给定两个插值节点和函数值和一个点的导数值

基本形式为

赫米特插值法

double Hiemite(double *p,int n,double x)
{//p表示插值节点，n表示插值节点个数，x是自变量的取值
 //p是一个三维数组，第一维为X，第二维度为y，第三维度为导数值y′
    int i,j,k;
    double ans=0;
    double numerator,denominator;//前者为分子，后者为分母
    double *lx=(double *)malloc(n*sizeof(double));//记录插值基函数的平方
    double *ax=(double *)malloc(n*sizeof(double));//记录系数α
    double *bx=(double *)malloc(n*sizeof(double));//记录系数β
    for(k=0;k<n;k++)
    {
        numerator=1;    denominator=1;
        for(i=0;i<n;i++)
        {
            if(i!=k){
                numerator*=(x- (*(p+i*3)));
                denominator*=(*(p+k*3) - *(p+i*3));
            }//End if
        }//End for-i
        *(lx+k)=(numerator)/(denominator);
        *(lx+k)=pow( (*(lx+k)) ,2);//求插值基函数的平方

        *(ax+k)=1-2*( (x-(*(p+k*3))) / (*(p+k*3) - *(p+((k+1)%2)*3) ) );
        *(bx+k)=(x-( *(p+k*3) ));
        ans+=(*(p+k*3+1)) * (*(ax+k)) * (*(lx+k));//Yi*αi*Li^2
        ans+=(*(p+k*3+2)) * (*(bx+k)) * (*(lx+k));//Y′i*βi*Li^2
    }//End for-k
    return ans;
}

最小二乘法

二乘法是为了在一系列离散点中找到一条曲线，最好地拟合这些点

为了达到最好的拟合效果，我们使其曲线的残差平方和最小

最小二乘法

double Fitting(double *p,int n)
{
    int i,j,k;
    double b,a;
    double ans[4]={0};
    double *m=(double *)malloc(2*3*sizeof(double));

    for(i=0;i<n;i++){//用X,Y,x*x,x*y的累加和
        ans[0]+=(*(p+i*2));//保存Xi求和结果
        ans[1]+=(*(p+i*2+1));//保存Yi求和结果
        ans[2]+=( (*(p+i*2)) * (*(p+i*2)) );//保存Xi*Xi求和结果
        ans[3]+=( (*(p+i*2)) * (*(p+i*2+1)) );//保存Xi*Yi求和结果
    }//End for-i
    *(m+0)=n;   *(m+1)=ans[0];   *(m+2)=ans[1];//第一个方程组的系数
    *(m+3)=ans[0];  *(m+4)=ans[2];  *(m+5)=ans[3];//第二个方程组的系数
    //以下过程模拟解方程的过程
    for(i=0;i<3;i++){//i循环两个方程组的三个系数
        *(m+i)/=n;//第一个方程组左右两边除以N
        *(m+i+3)-=( ans[0] * ( *(m+i) ) );//模拟第一个方程组乘以Xi的累加和后消去a的过程
    }

    b= ( *(m+5) )/( *(m+4) );//a消去以后直接求b
    a= ( (*(m+2))-( b*(*(m+1))) ) /(*m);//带入第一个方程组中求出a
    printf("b=%lf a=%lf\n",b,a);
    printf("线性拟合的结果为Y=%lfX",b);
    a<0?printf("%lf\n",a):printf("+%lf\n",a);
}