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; // 分母的间距，每次+1int 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 阶差商 \n", rank);int pos = 0;int r = 1;// 根据n+(n+1)+...+(n-rank)求得rank阶差商在差商表数组中的索引while (rank > 1){// printf("geting dd's index, waiting...\n");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("=======================差商表的信息=======================\n");printf("差商个数有：%d, 待插值节点个数：%d\n =======================差商表=======================", ddlength, lennx); int ifnextrow = 0; // 控制打印换⾏for (int i = 0; i < ddlength; i++) {if (ifnextrow == (i)*(i + 1) / 2)printf("\n");printf("%lf, ", dividedDifferences[i]);ifnextrow++;}printf("\n");// 计算Newton Interpolating 插值double ddi = 0;double coef = 1;for (int i = 0; i < lennx; i += 2) {coef = (double)1;result[i] = y[i];printf("=============打印计算过程=============\n");for (int j = 1; j < lenX -1; j++) {coef *= (newX[i] - x[j - 1]);ddi = getDividedDifference(dividedDifferences, ddlength, j, 0);printf("取得差商： %lf\n", ddi);result[i] = result[i] + coef * ddi;printf("计算result[%d] + %lf * %lf", i, coef, ddi);printf("获得结果result %d: %lf\n", i, result[i]);}printf("result[%d] = %lf\n", i, result[i]);// =======================选做题：求误差==========================printf("求解截断误差所需差商:%lf\n", getDividedDifference(dividedDifferences, ddlength, lenX-1, 0));// printf("求解截断误差所需多项式：%lf\n", coef);printf("求解截断误差所需多项式的最后⼀项：%lf - %lf\n", 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!\n";printf("Newton插值!\n");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("=====================最终结果=====================\n");printf("求得sin(1.609)的近似值为：%lf\n", *result);printf("截断误差：%.10lf\n", *(result + 1));if (result) {delete[] result;}return0;}。
