// 函数功能: 求解拉格朗日微分矩阵在tk处的值在第i个配置点上的分量 // 输入: point:n * 1维矩阵,是插值点的时间维度,n为插值节点数目(升序) // k0:为当前点的标号 // j:以1开始的标号 // 输出: f当前点的导数值 // 方法: 拉格朗日微分矩阵的基本公式 // 编写: 李兆亭,2019 / 11 / 12 // 李兆亭, 2019 / 12 / 23修改程序说明 // 李兆亭, 2019 / 12 / 23修改tk为k0,适用于配置点的编程计算,原式更适用于一般情况 #include "GPM.h" int DLagrange_point(mat* point1, int *k0, int *j, double *f) { // 初始化 int row_len = (size(*point1))(0); double b = 1.0, a_sum = 0; double tj, ti, a, tk; // 分母计算, 分母与标号L无关 for (int L = 0; L < row_len; L++) { if (L != *j) { // (*point1).print("(*point1)"); tj = (*point1)(L, 0); ti = (*point1)(*j, 0); b = b*(ti - tj); } } // 分子计算,分子与标号有关 for (int i0 = 0; i0 < row_len; i0++) { if (i0 == *j) { a = 0; } else { a = 1; for (int L = 0; L < row_len; L++) { if (L == *j || L == i0) { } else { tj = (*point1)(L, 0); tk = (*point1)(*k0, 0); a = a*(tk - tj); } } } a_sum = a_sum + a; } // 相除得结果 *f = a_sum / b; return 1; }