牛顿差商多项式的理解与C++实现

这篇接着上一篇，对牛顿基本插值多项式的理解与代码实现

在拉格朗日插值法中，

一般插值公式为：

其中的每次计算都依赖于全部的插值结点，在增加或减少结点时，必须全部重新计算。（不是特别理解）为了克服这个缺点从而引入了牛顿插值法。

牛顿基本插值多项式的一般公式为：

其中定义：

为零阶差商

为一阶差商

为二阶差商

以此类推。。。

另有公式：

在编程实现的过程中就需要考虑，怎么计算，

可分为定义法和公式法。

• 定义法：基本思路便是先算出一阶差商，然后算出二阶差商，。。。

这样设计必然需要设定一个数组用于存储中间。

• 公式法：利用公式直接计算得到

这样设计的后果很明显，在公式中会有大量乘法运算，速度会比较慢。

---

下面是代码阶段：

在代码中用quotient_array用于存储所有差商，

quotient_list用于存储从x0开始的差商。

#include<iostream>
#define maxsize 100
using namespace std;
double quotient_array[maxsize][maxsize];
void D_creat_quotient_list(int n,double x[maxsize],double y[maxsize],double quotient_list[maxsize])
{//构造差商表
    //定义法 占空间但是省时间
    for(int k=0;k<n;k++)
    {
        for(int i=0;i<n-k;i++)
        {
            if(k==0)
            {
                quotient_array[i][0]=y[i];
            }
            else
            {
                quotient_array[i][k]=(quotient_array[i+1][k-1]-quotient_array[i][k-1])/(x[i+k]-x[i]);
            }
        }
        quotient_list[k]=quotient_array[0][k];
    }
}
void F_creat_quotient_list(int n,double x[maxsize],double y[maxsize],double quotient_list[maxsize])
{//构造差商表
    //公式法 省空间 但是费时间
    quotient_list[0]=y[0];
    double temp;
    for(int k=1;k<n;k++)
    {
        quotient_list[k]=0;
        for(int j=0;j<=k;j++)
        {
            temp = 1;
            for(int i=0;i<=k;i++)
            {
                if(i!=j)
                {
                    temp*=(x[j]-x[i]);
                }
            }
            temp = y[j]/temp;
            quotient_list[k]+=temp ;
        }
    }
}

double F_interpolitation(int n,double x[maxsize],double y[maxsize],double vx,double result,double quotient_list[maxsize])
{//差商牛顿插值法
    //两种方法二选一
    D_creat_quotient_list(n,x,y,quotient_list);//定义法
    //F_creat_quotient_list(n,x,y,quotient_list);//公式法
    double temp;
    result = quotient_list[0];
    for(int k=1;k<=n;k++)
    {
        temp=1;
        for(int i=0;i<k;i++)
        {
            temp*=(vx-x[i]);
        }
        result+=quotient_list[k]*temp;
    }
    return result;
}
int main()
{
    int n;
    double x[maxsize],y[maxsize],result=0,vx,temp;
    double quotient_list[maxsize];
    cin>>n>>vx;
    for(int i=0;i<n;i++)
    {
        cin >> x[i] >> y[i];
    }
    result=F_interpolitation(n,x,y,vx,result,quotient_list);
    cout<<result<<endl;
    return 0;
}

---

测试样例：

参考：

x	100	121	144	169
	10	11	12	13

分别用定义法和公式法求的值

定义法：

4
115
100 10
121 11
144 12
169 13

---

公式法：

4
115
100 10
121 11
144 12
169 13

对比也可以看出速度上是有差距的，一般都用定义法求值。

其中两个函数得到的quotient_list的值相同。

quotient_list

10

0.047619

-0.00009411

0.0000003138