「数学」- 手撸 Newton插值法 C++

具体思路现在，没空补，看代码之前，建议推一下newton插值公式。

 

#include <iostream>
#include <cstdio>
#include <queue>
#include <stdlib.h>
#include <ctime>
/// Newton 插值法的实现方法
#define NAN 0x3f3f3f3f
class Newdon{
    // （x-x0）(x-x1).....(x-xn) *a
    private:
        double  fac[1000]={0}; // o ~ n 的u多项式系数, 系数从0到n
        double  fac_ploy[1000]={0}; // 新增一项的（x-x0）(x-x1).....(x-xn)的各个多项式的系数，用于构造新的多项式系数
        int fac_len = 0;  //有效系数的长度 从0开始

        double get_fac_ploy(double);// （x-x0）(x-x1).....(x-xn)的值
    public:
        Newdon(double *,double *,int);
        void add_node(double,double);
        double get_value(double x);
        void test_print();
};

void Newdon::test_print()
{
    for(int i=0;i<fac_len;i++)
        printf(" %lf ",fac[i]);
    printf("\n");
    for(int i=0;i<fac_len;i++)
        printf(" %lf ",fac_ploy[i]);
    printf("\n");
}
Newdon::Newdon(double *x_value,double *y_value, int len)
{
    if(len<1) printf("len is too small to creat Newton\n");
    // 系数构造方式 比较巧妙
    //动态构造 新建新newton项式的系数， 当然newton项式之间的系数有递推关系
    for(int i=0;i<len;i++)
    {
        //test_print();
        add_node(x_value[i],y_value[i]);
    }
}

// 已有一个节点的基础上 添加一个新节点
void Newdon::add_node(double x,double y)
{
    //第一个节点
    if(fac_len == 0)
    {
        fac[0] = y;
        fac_ploy[0] = -x; fac_ploy[1] = 1;//下一个多项式的（x-x0）
        fac_len ++;
        return;
    }
    // 之后的节点
    double tmp1 = get_value(x);
    fac_len++;// 有效的fac_ploy 比当前fac多一个
    double tmp2 = get_fac_ploy(x);
    double ploy_a = (y - tmp1) / tmp2; // 求出系数a，之后将系数相乘相加即可
    for(int i=0;i<fac_len;i++)
    {
        fac[i] += fac_ploy[i] * ploy_a;
    }
    //随之，增加有效的ploy系数
    for(int i=fac_len-1;i>=0;i--)
    {
        fac_ploy[i+1] += fac_ploy[i];
        fac_ploy[i] = fac_ploy[i] * (-x);
    }
}
double Newdon::get_fac_ploy(double x)
{
    if(fac_len <= 0){
        printf("fac_len is <= 0\n");
        return 0;
    }
    double res = fac_ploy[fac_len-1];
    for(int i=fac_len-2;i>=0;i--)
        res  =  res * x + fac_ploy[i];
    return res;
}
double Newdon::get_value(double x)
{
    if(fac_len <= 0){
        printf("fac_len is <= 0\n");
        return 0;
    }
    double res = fac[fac_len-1];
    for(int i=fac_len-2;i>=0;i--)
        res  =  res * x + fac[i];
    return res;
}

int main(){
    srand(time(0));
    double y[10] = {10,20,49,50,40,0,-40,-30,-10,20};
    double x[10] = {1,2,3,4,5,6,7,8,9,10};
    int len = 10;
   /* for(int i=0;i<len;i++)
    {
        x[i] = rand()%10000;
        y[i] = rand()%20000;
    }
    */
    Newdon first = Newdon(x,y,len);
    first.test_print();
    int a = 11, b =15;
    for(int i=a;i<b;i++)
    printf("c(%d) = %.2lf\n",i,first.get_value(i));
}

可用MATLAB的Newton的插值法源代码进行验证！

代码转自

%newton.m
%求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
%输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量，Y是纵坐标向量，
%x是以向量形式输入的m个插值点，M在[a,b]上满足｜f~(n+1)(x)｜≤M
%注：f~(n+1)(x)表示f(x)的n+1阶导数
%输出的量：向量y是向量x处的插值，误差限R，n次牛顿插值多项式L及其系数向量C，
%差商的矩阵A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
    z = x(t);
    A = zeros(n,n);
    A(:,1) = Y';
    s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
        for j = 2 : n
            for i = j : n
                A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
            end
            q1 = abs(q1*(z-X(j-1)));
            c1 = c1 * j;
        end
        C = A(n, n); q1 = abs(q1*(z-X(n)));
        for k = (n-1):-1:1
            C = conv(C, poly(X(k)));
            d = length(C);
            C(d) = C(d) + A(k,k);%在最后一维，也就是常数项加上新的差商
        end
        y(t) = polyval(C,z);
        R(t) = M * q1 / c1;
end
L = poly2sym(C);

运行

M  = 10000000000000

X = [1 2 3 4 5 6 7 8 9 10 ]
Y = [ 10 20 49 50  40  0 -40 -30  -10 20 ]
x = [11,12,13,14]
 newton(X,Y,x,M)