• Newton插值的C++实现


    Newton(牛顿)插值法具有递推性,这决定其性能要好于Lagrange(拉格朗日)插值法。其重点在于差商(Divided Difference)表的求解。

    步骤1. 求解差商表,这里采用非递归法(看着挺复杂挺乱,这里就要自己动笔推一推了,闲了补上其思路),这样,其返回的数组(指针)就是差商表了,

    /*
     * 根据插值节点及其函数值获得差商表
     * 根据公式非递归地求解差商表
     * x: 插值节点数组
     * y: 插值节点处的函数值数组
     * lenX: 插值节点的个数
     * return: double类型数组
    */
    double * getDividedDifferenceTable(double x[], double y[], int lenX) {
        double *result = new double[lenX*(lenX - 1) / 2];
        for (int i = 0; i < lenX - 1; i++) {
            result[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
        }
        int step = lenX - 1;  // 增加步长
        int yindex = lenX - 1;  // 分子的基准值,线性减速递增
        int xgap = 2;  // 分母的间距,每次+1
        int xstep = 2;
        while (step >= 1) {
            for (int i = 0; i < step - 1; i++) {
                result[yindex + i] = (result[yindex - step + i + 1] - result[yindex - step + i]) / (x[xstep + i] - x[xstep + i - xgap]);
            }
            yindex += (step - 1);
            xstep++;
            step--;
            xgap++;
        }
    
        return result;
    }

    步骤 2. 取得差商表每一列的第一个作为基函数的系数

    /**
     * 从差商表中获取一定阶数的某个差商
     * dd: 差商表
     * len: 差商表的长度
     * rank: 所取差商的阶数
     * index: rank阶差商的第index个差商
     * return: 返回double类型所求差商
     */
    double getDividedDifference(double dd[], int len, int rank, int index) {
        // 根据差商表的长度求解差商的最高阶数
        // 由于差商表是三角形的,因此有规律可循,解方程即可
        int rankNum = (int)(0.5 + sqrt(2 * len + 0.25) - 1);
        printf("%d 阶差商 
    ", rank);
        int pos = 0;
        int r = 1;
        // 根据n+(n+1)+...+(n-rank)求得rank阶差商在差商表数组中的索引
        while (rank > 1)
        {
            //     printf("geting dd's index, waiting...
    ");
            pos += rankNum;
            rank--;
            rankNum--;
        }
        pos += index;  // 然后加上偏移量,意为rank阶差商的第index个差商
    
        return dd[pos];
    }

    步骤3. Newton 插值主流程

    /*
     * Newton Interpolating 牛顿插值主要代码
     * x: 插值节点数组
     * y: 插值节点处的函数值数组
     * lenX: 插值节点个数(数组x的长度)
     * newX: 待求节点的数组
     * lennx: 待求节点的个数(数组lenX的长度)
     */
    double *newtonInterpolation(double x[], double y[], int lenX, double newX[], int lennx) {
        // 计算差商表
        double *dividedDifferences = getDividedDifferenceTable(x, y, lenX);
        // 
        double *result = new double[lennx];
    
        // 求差商表长度
        int ddlength = 0;
        for (int i = 1; i < lenX; i++)
        {
            ddlength += i;
        }
        // 打印差商表信息
        printf("=======================差商表的信息=======================
    ");
        printf("差商个数有:%d, 待插值节点个数:%d
     =======================差商表=======================", ddlength, lennx);
        int ifnextrow = 0;  // 控制打印换行
        for (int i = 0; i < ddlength; i++) {
            if (ifnextrow == (i)*(i + 1) / 2)
                printf("
    ");
            printf("%lf, ", dividedDifferences[i]);
            ifnextrow++;
        }
        printf("
    ");
        // 计算Newton Interpolating 插值
        double ddi = 0;
        double coef = 1;
        for (int i = 0; i < lennx; i += 2) {
            coef = (double)1;
            result[i] = y[i];
            printf("=============打印计算过程=============
    ");
            for (int j = 1; j < lenX -1; j++) {
                coef *= (newX[i] - x[j - 1]);
                ddi = getDividedDifference(dividedDifferences, ddlength, j, 0);
                printf("取得差商: %lf
    ", ddi);
                result[i] = result[i] + coef * ddi;
                printf("计算result[%d] + %lf * %lf", i, coef, ddi);
                printf("获得结果result %d: %lf
    ", i, result[i]);
            }
            printf("result[%d] = %lf
    ", i, result[i]);
    
            // =======================选做题:求误差==========================
            printf("求解截断误差 所需差商:%lf
    ", getDividedDifference(dividedDifferences, ddlength, lenX-1, 0));
            // printf("求解截断误差 所需多项式:%lf
    ", coef);
            printf("求解截断误差 所需多项式的最后一项:%lf - %lf
    ", newX[i] , x[lenX - 2]);
            result[i + 1] = getDividedDifference(dividedDifferences, ddlength, lenX - 1, 0)*coef*(newX[i] - x[lenX - 2]);
            // ===============================================================
        }
        if (dividedDifferences) {
            delete[] dividedDifferences;
        }
        return result;
    }

    步骤4. 主函数示例

    int main()
    {
        std::cout << "Hello World!
    ";
        printf("Newton插值!
    ");
        double x[] = {0.40, 0.55, 0.65, 0.80, 0.90, 1.05};
        //double x[] = { 1.5, 1.6, 1.7 };
        int lenX = sizeof(x) / sizeof(x[0]);
        double y[] = {0.41075, 0.57815, 0.69675, 0.88811, 1.02652, 1.25382};
        //double y[] = { 0.99749, 0.99957, 0.99166 };
        double newx[] = {0.596};
        //double newx[] = { 1.609 };
        // 计算牛顿插值结果和截断误差
        double *result = newtonInterpolation(x, y, lenX, newx, 1);
        printf("=====================最终结果=====================
    ");
        printf("求得sin(1.609)的近似值为:%lf
    ", *result);
        printf("截断误差:%.10lf
    ", *(result + 1));
        if (result) {
            delete[] result;
        }
        return 0;
    }
  • 相关阅读:
    java基础5 (一维)数组和二维数组
    Java 内部类
    Java 抽象类和接口
    Java 多态
    Java 继承
    Java 包(package)
    Java String类和StringBuffer类
    Java 数组
    Java 封装与类
    Java 概述和编程基础
  • 原文地址:https://www.cnblogs.com/zhaoke271828/p/13891817.html
Copyright © 2020-2023  润新知