Newton 插值

最新推荐文章于 2024-04-16 15:18:30 发布

zephyr_john

最新推荐文章于 2024-04-16 15:18:30 发布

阅读量363

点赞数 2

分类专栏：算法

本文链接：https://blog.csdn.net/zephyr_john/article/details/72356941

版权

算法专栏收录该内容

3 篇文章 0 订阅

订阅专栏

define MAXTIMES 10 //最高构造10次插值函数

typedef struct {
double x;
double y;
}Point; //插值节点类型

int Times; //实际要构造的次数
double DCList[(MAXTIMES + 1)*(MAXTIMES + 2) / 2]; //差商列表

Point knotList[MAXTIMES + 1];

void initKnotList(); //插值节点数据输入
void initDCList(); //构造差商表
void displayDCList();

double Newton(double u); //Newton 插值函数，计算点x处的函数值

int main()
{
double u; //计算该点的插值函数值

initKnotList();
initDCList();

printf("%f \n", Newton(u));


while (getchar() != '\n');
getchar();
return 0;

}

void initKnotList()
{
int i;

printf("input Times:");
scanf_s("%d", &Times, sizeof(int));

for (i = 0; i <= Times; i++)
{
    scanf_s("%lf", &knotList[i].x, sizeof(double));
    scanf_s("%lf", &knotList[i].y, sizeof(double));
}

}

void initDCList()
{
int i, j;
int prePosi, currPosi;

currPosi = 0;
for (i = 0; i <= Times; i++)                 //第i阶差商
{
    for (j = 0; j <= Times - i; ++j)         //第i阶的第j个差商
    {
        if (i == 0)                          //0阶差商
        {
            DCList[currPosi + j] = knotList[j].y;
        }
        else
        {
            DCList[currPosi + j] = 
            (DCList[prePosi + j] - DCList[prePosi + j + 1]) / 
            (knotList[j].x - knotList[j + i].x);
        }
    }
    prePosi = currPosi;
    currPosi = (i + 1)*Times + 1;
}

}

void displayDCList()
{
int i, j;
int currPosi;

currPosi = 0;
for (i = 0; i <= Times; i++)
{
    for (j = 0; j <= Times - i; ++j)
    {
        printf("%f ", DCList[currPosi+j]);
    }
    currPosi = ((i+1)*Times) + 1;
    printf("\n");
}

}

double Newton(double u)
{
int i, j;
int DCNumber=0;
double sum = 0;
double w;

for (i = 0; i <= Times; i++)
{
    w = 1;
    for (j = 0; j < i; ++j)
    {
        w *= (u - knotList[j].x);
    }

    sum += DCList[DCNumber] * w;
    DCNumber = ((i + 1)*Times + 1);
}
return sum;

}