数值分析的两种插值
数值分析的两种插值：拉格朗日插值和牛顿插值
拉格朗日四次插值
python代码：
1 def lagrange(x,y,n,xk,Nk): 2 f=[0.0,0.0,0.0,0.0,0.0] 3 for i in range(0,n): 4 f[i]=y[i] 5 k=0 6 for j in range(0,5): 7 if(i!=j): 8 f[i]=f[i]*(xk-x[j])/(x[i]-x[j]) 9 k=k+1 10 Nk=Nk+f[i] 11 print("当x={}时，f(x)={}".format(xk,Nk)) 12 return 0 13 14 def main(): 15 insertNumber = int(input("请输入插值点的数:")) 16 17 predictNumber = int(input("请输入预估点的数:")) 18 19 array=input("以空格为间隔输入插值点的x值：") 20 insertNumberX=[float(n) for n in array.split()] 21 22 array=input("以空格为间隔输入插值点的y值：") 23 insertNumberY=[float(n) for n in array.split()] 24 25 array=input("以空格为间隔输入预估点的x值：") 26 predictNumberX=[float(n) for n in array.split()] 27 28 py=0.0 29 for i in range(0,predictNumber): 30 lagrange(insertNumberX,insertNumberY,insertNumber,predictNumberX[i],py) 31 32 if __name__ == '__main__': 33 main()
实验结果：
备注：
基函数代码：
f[i]=f[i]*(xk-x[j])/(x[i]-x[j])
C语言代码：
#include <stdio.h> #define MAX 20 //输入点的结构 typedef struct stPoint { double x; double y; } Point; int main() { int n; int i,j; Point points[MAX]; double x,tmp,lagrange=0;//这个x是你将要计算的f(x)插值点，tmp是拉格朗日基函数，larange是根据拉格朗日函数得出f(x)的值 printf("请输入被插值点的个数:(它是从0开始的，所以输入3代表4个点）"); scanf("%d",&n); if(n>MAX) { printf("您输入的个数过多．"); return 1; } if(n<=0) { printf("您输入的点数太少．"); return 1; } //输入插值点的x值和y值 printf("请输入插值点的x值和y值：\n"); for(i=0;i<=n;i++) scanf("%lf%lf",&points[i].x,&points[i].y); //输入计算拉格朗日插值多项式的x值 for(int a=0;a<=3;a++){ printf("\n请输入计算拉格朗日插值多项式的x值:"); scanf("%lf",&x); //利用拉格朗日插值公式计算函数x值的对应f(x) for(i=0;i<=n;i++) { for(j=0,tmp=1;j<=n;j++) { if(j==i) //去掉xi与xj相等的情况 continue; //范德蒙行列式下标就是j!=k,相等分母为０就没意义了 tmp=tmp*(x-points[j].x)/(points[i].x-points[j].x);//这个就是套公式，没什么难度 //tmp是拉格朗日基函数 } lagrange=lagrange+tmp*points[i].y; //最后计算基函数*y,全部加起来，就是该x项的拉格朗日函数了 } //拉格朗日函数计算完毕，代入所求函数x的值，求解就ok了 printf("\n拉格朗日函数f(%lf)=%lf\n",x,lagrange); } return 0; }
实验结果：
核心代码：
tmp=tmp*(x-points[j].x)/(points[i].x-points[j].x);
-----------------------------------------------------------------------------------------------------------------------------------------------------
以上是两种代码算一题，答案有些出入，代码是没什么问题的，结果不一样我这代码小白就搞不清楚了！
牛顿四次插值：
#include <stdio.h> #include <stdlib.h> void data(double* x, double* y, int n); //x-横坐标，y-纵坐标，f-插值系数，n插值节点个数 void newton(double* x, double* y, double* f, int n); void printnew(double* x, double* y, double* f, int n); void newvalue(double* x, double* y, double* f, int n); int main(void) { int n; double *x, *y, *f; printf("输入要插值节点的个数:"); scanf("%d", &n); x = (double*)malloc(sizeof(double) * n); y = (double*)malloc(sizeof(double) * n); f = (double*)malloc(sizeof(double) * (n - 1) * n / 2); data(x, y, n); newton(x, y, f, n - 1); printnew(x, y, f, n); do { newvalue(x, y, f, n); } while (1); return 0; } void data(double* x, double* y, int n) { //读取初始数据 int i = 0; while (i < n) { printf("x[%d]:", i); scanf("%lf", &x[i]); printf("y[%d]:", i); scanf("%lf", &y[i]); i++; } } void newton(double* x, double* y, double* f, int n) { //建立牛顿插值多项式的系数 int i = 0, j, k = 0; for (i = 0; i < n; i++) for (j = 0; j < n - i; j++) { if (i == 0) f[k++] = (y[j + 1] - y[j]) / (x[j + 1] - x[j]); else { f[k] = (f[k + i - n ] - f[k + i - n - 1]) / (x[j + i + 1] - x[j]); k++; } } } void printnew(double* x, double* y, double* f, int n) { //输出差商表 int i, j, k = 0; printf("差商表:\n"); printf("x\t "); for (i = 0; i < n; i++) printf("f(x%d)\t\t", i); printf("\n"); for (i = 0; i < n; i++) printf("----------------"); printf("\n"); for (i = 0; i < n; i++) { printf("%-10g %-10g", x[i], y[i]); k = i - 1; for (j = 0; j < i; j++) { printf(" %-10g", f[k]); k += n - 2 - j; } if (j == i) printf("\n"); } } void newvalue(double* x, double* y, double* f, int n) { //根据牛顿插值多项式预测下一个节点的值 double a, *b; int i, k = 0; b = (double*)malloc(sizeof(double) * n); printf("输入要要插入节点的X的值:"); scanf("%lf", &a); b[0] = 1.0; for (i = 0; i < n - 1; i++) b[i + 1] = b[i] * (a - x[i]); for (i = 0; i < n; i++) { if (i == 0) a = y[0]; else { a += b[i] * f[k]; k += n - i; } } printf("插值节点对应的Y的值:%g\n", a); }
