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

最新推荐文章于 2024-09-12 11:26:15 发布

frboom

最新推荐文章于 2024-09-12 11:26:15 发布

阅读量1.9k

点赞数 1

分类专栏：数值计算

本文链接：https://blog.csdn.net/frboom/article/details/108537399

版权

数值计算专栏收录该内容

2 篇文章

订阅专栏

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

在拉格朗日插值法中，

一般插值公式为： $L_{n}\left ( x \right )=\sum^{_{n}}_{k=0}\left [ y_k\prod_{j=0;j\neq k}^{n} \left ( \frac{x-x_{j}}{x_k-x_j} \right )\right ]$

其中 $l_{n}\left ( x \right )=\prod_{j=0;j\neq k}^{n} \left ( \frac{x-x_{j}}{x_k-x_j} \right )$ 的每次计算都依赖于全部的插值结点，在增加或减少结点时，必须全部重新计算。（不是特别理解）为了克服这个缺点从而引入了牛顿插值法。

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

$\begin{matrix} N_n\left ( x \right )= \ f\left [ x_0 \right ] +f[x_0,x_1](x-x_0) \\ +f\left [x_0,x_1,x_2 \right ](x-x_0)(x-x_1)+... \\ \ +f\left [x_0,x_1,...,x_n \right ](x-x_0)(x-x_1)\cdot ...\cdot (x-x_{n-1}) \end{matrix}$

其中定义：

$f(x_0)$ 为零阶差商

$f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_j-x_i}$ 为一阶差商

$f[x_0,x_1,x_2]=\frac{f[x_j,x_k]-f[x_i,x_j]}{x_j-x_i}$ 为二阶差商

以此类推。。。

另有公式：

$f[x_0,x_1,...,x_m]=\sum_{j=0}^{m} \left ( \frac{f(x_j)}{(x_j-x_0)\cdot ...\cdot (x_j-x_{j-1})(x_j-x_{j+1})\cdot ...\cdot (x_j-x_m)} \right )$

在编程实现的过程中就需要考虑，怎么计算 $f[x_0,x_1,...,x_m]$ ，

可分为定义法和公式法。

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

这样设计必然需要设定一个数组用于存储中间 $f[x_i,x_{i+1},...,x_j]$ 。

公式法：利用公式直接计算得到 $f[x_0,x_1,...,x_m]$

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

下面是代码阶段：

在代码中用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;
}

测试样例：

参考：