数值分析的两种插值

数值分析的两种插值：拉格朗日插值和牛顿插值

拉格朗日四次插值

 

python代码：

def lagrange(x,y,n,xk,Nk):
    f=[0.0,0.0,0.0,0.0,0.0]
    for i in range(0,n):
        f[i]=y[i]
        k=0
        for j in range(0,5):
            if(i!=j):
                f[i]=f[i]*(xk-x[j])/(x[i]-x[j])
                k=k+1
        Nk=Nk+f[i]
    print("当x={}时，f(x)={}".format(xk,Nk))
    return 0

def main():
    insertNumber = int(input("请输入插值点的数:"))
    
    predictNumber = int(input("请输入预估点的数:"))
    
    array=input("以空格为间隔输入插值点的x值：")
    insertNumberX=[float(n) for n in array.split()]

    array=input("以空格为间隔输入插值点的y值：")
    insertNumberY=[float(n) for n in array.split()]

    array=input("以空格为间隔输入预估点的x值：")
    predictNumberX=[float(n) for n in array.split()]

    py=0.0
    for i in range(0,predictNumber):
        lagrange(insertNumberX,insertNumberY,insertNumber,predictNumberX[i],py)

if __name__ == '__main__':
    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值：
");
   for(i=0;i<=n;i++)
      scanf("%lf%lf",&points[i].x,&points[i].y);
    //输入计算拉格朗日插值多项式的x值
   for(int a=0;a<=3;a++){
       printf("
请输入计算拉格朗日插值多项式的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("
拉格朗日函数f(%lf)=%lf
",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("差商表:
");
    printf("x	  ");
    for (i = 0; i < n; i++)
        printf("f(x%d)		", i);
    printf("
");
    for (i = 0; i < n; i++)
        printf("----------------");
    printf("
");
    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("
");
    }
}
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
", a);
}