牛顿插值多项式算法

NewtonDifference.h

#pragma once

class NewtonDifference
{
public:
	NewtonDifference(void);
	virtual ~NewtonDifference(void);

	// 三阶插值函数
	float f3(float x);

	// 牛顿插值多项式算法
	void NewInt(double* x, double* y
		, int n, double& xi, double* yint, double* ea);
};

NewtonDifference.cpp

#include "NewtonDifference.h"
#include <string.h>
#include <stdio.h>

NewtonDifference::NewtonDifference(void)
{
    //f3(2);

    // 开始牛顿插值
    double *x = new double[4];
    x[0] = 1.0f;
    x[1] = 4.0f;
    x[2] = 6.0f;
    x[3] = 5.0f;

    double *y = new double[4];
    y[0] = 0.0f;
    y[1] = 1.386294f;
    y[2] = 1.791759f;
    y[3] = 1.609438f;

    int n = 3;
    double * yint = new double[n];
    double * ea = new double[n];

    // 插出一个值
    double xi = 2.0f;
    //double y = 0.0f;
    //double ea = 0.0f;
    NewInt(x, y, n, xi, yint, ea);
}

NewtonDifference::~NewtonDifference(void)
{
    
}

#pragma warning(disable: 4244)
float NewtonDifference::f3(float x)
{
    float f3 = 0.0f;
    f3 = 0.0 + 0.4620981f * (x-1)
        + (-0.05187311f)* (x-1) * (x-4)
        + 0.007865529f * (x - 1) * (x - 4)
        * (x - 6);

    return f3;
}

// x, y ,采样的固定点
// n, 多项式插值的阶数
// xi, 要估计的自变量值
// yint, 不同阶数下插值出的y值
// ea, 不同阶数下的估计误差值
// 牛顿插值多项式算法
void NewtonDifference::NewInt(double* x, double* y
            , int n, double& xi, double* yint, double* ea)
{

    double *fdd = new double[n*n];
    ::memset(fdd, 0, sizeof(float)*n*n);

    // 0阶差商
    for(int i = 0; i <= n; i++)
    {
        fdd[i*n] = y[i];

        printf("x[%d]=%f, y[%d]=%f\n", i, x[i], i, y[i]);
    }

    // 1阶及以后的差商
    for (int j = 1; j <= n; j++)
    {
        printf("\n%d order difference devide:\n", j);
        for (int i = 0; i <= n -j; i++)
        {
            // 当前阶数下，各采样点的差商
            fdd[i*n+j] = (fdd[(i+1)*n + j-1] - fdd[i*n + j-1])
                / (x[i+j] - x[i]);

            printf("fdd[%d][%d]=%f\n\n", i,j, fdd[i*n+j]);
        }
    }

    double xterm = 1.0f;
    double yterm = 0.0f;
    yint[0] = fdd[0*n + 0];

    // 处理得到的各阶差商
    for(int order = 1; order <= n; order++)
    {
        // 多项式中,x项的处理
        xterm = xterm * (xi - x[order-1]);

        // fdd[][]是一个二维数组
        yterm = yint[order - 1]
                + fdd[0*n + order] * xterm;

        // 当前阶数下的估值误差
        // Rn = fn+1(x) - fn(x);
        ea[order] = yterm - yint[order - 1];

        // 可观察不同阶数下的估值精度
        yint[order] = yterm;

        printf("yint[%d]=%f, ea[%d]=%f\n", order, yint[order]
        ,order, ea[order]);
    }
}