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; }