拉格朗日插值数学原理
给定平面上 n + 1 n + 1 n+1个不同的数据点 ( x k , f ( x k ) ) (x_{k},f(x_{k})) (xk,f(xk)), k = 0 , 1 , ⋯ , n k = 0,1,\cdots,n k=0,1,⋯,n, x i ≠ x j x_{i} \neq x_{j} xi=xj, i ≠ j i \neq j i=j;则满足条件
P n ( x k ) = f ( x k ) , k = 0 , 1 , ⋯ , n P_{n}\left( x_{k} \right) = f\left( x_{k} \right),\ \ k = 0,1,\cdots,n Pn(xk)=f(xk), k=0,1,⋯,n
的 n n n次拉格朗日插值多项式
P n ( x ) = ∑ k = 0 n f ( x k ) l k ( x ) P_{n}(x) = \sum_{k = 0}^{n}{f(x_{k})}l_{k}(x) Pn(x)=k=0∑nf(xk)lk(x)
是存在唯一的。若 x k ∈ [ a , b ] , k = 0 , 1 , ⋯ , n x_{k} \in \lbrack a,b\rbrack,k = 0,1,\cdots,n xk∈[a,b],k=0,1,⋯,n,且函数 f ( x ) f(x) f(x)充分光滑,则当 x ∈ [ a , b ] x \in \lbrack a,b\rbrack x∈[a,b]时,有误差估计式
f ( x ) − P n ( x ) = f ( n + 1 ) ( ξ ) ( n + 1 ) ! ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n ) , ξ ∈ [ a , b ] f\left( x \right) - P_{n}\left( x \right) = \frac{f^{\left( n + 1 \right)}\left( \xi \right)}{\left( n + 1 \right)!}\left( x - x_{0} \right)\left( x - x_{1} \right)\cdots\left( x - x_{n} \right),\ \ \xi \in \lbrack a,b\rbrack f(x)−Pn(x)=(n+1)!f(n+1)(ξ)(x−x0)(x−x1)⋯(x−xn), ξ∈[a,b]
程序设计
核心代码
double x, y = 0.0;scanf("%lf", &x);double a[N + 1], b[N + 1];int n = 0;while (scanf("%lf%lf", &a[n], &b[n]) >= 2) n++;n--;for (int k = 0; k <= n; k++) {double l = 1.0;for (int j = 0; j <= n; j++) {if (j != k) l *= (x - a[j]) / (a[k] - a[j]);}y += l * b[k];}printf("x = %.3lf\ny = %.3lf", x, y);
作为真正的程序员!一定要会设计框架
!!!!!!!
so 我设计了这么一个Lagrange测试框架 只要修改参数即可!!!!
框架示例
可以修改:
- n的数量
- x的数量
- 最大的n值
- 若干n的值
- 若干x的值
- 左右区间值
- 插值点取值方法
X()
- y计算方法
Y()
#include <cmath>
#include <cstdio>
#define N1 3 // n amount
#define N2 4 // x amount
#define N3 20 // n maxint Ns[N1] = {5, 10, 20};
double x[N2] = {-0.95, -0.05, 0.05, 0.95};
double l = -1.0;
double r = 1.0;double X(int k, int n) {double h = (r - l) / n;return l + k * h;
}double Y(double x) { return 1 / (1 + x * x); }int main() {for (int i = 0; i < N2; i++) printf("\tx=%.2lf", x[i]);printf("\n");for (int i = 0; i < N1; i++) {double a[N3 + 1], b[N3 + 1];int n = Ns[i];for (int k = 0; k <= n; k++) {a[k] = X(k, n); // xb[k] = Y(a[k]); // y}printf("n=%d", n);for (int p = 0; p < N2; p++) {double y = 0.0;for (int k = 0; k <= n; k++) {double l = 1.0;for (int j = 0; j <= n; j++) {if (j != k) l *= (x[p] - a[j]) / (a[k] - a[j]);}y += l * b[k];}printf("\t%.6lf", y);}printf("\n");}printf("Actual");for (int p = 0; p < N2; p++) printf("\t%.6lf", Y(x[p]));return 0;
}
输出
x=0.75 x=1.75 x=2.75 x=3.75 x=4.75
n=5 0.528974 0.373325 0.153733 -0.025954 -0.015738
n=10 0.678990 0.190580 0.215592 -0.231462 1.923631
n=20 0.636755 0.238446 0.080660 -0.447052 -39.952449
Actual 0.640000 0.246154 0.116788 0.066390 0.042440
还没完!!!!!!得出的数据直接复制到Word里,全选+点击“转换文本为表格”
即可立刻变成表格!!!!!
- 逗号分隔,导出.csv文件,用excel打开复制……那也行???