Newton 插值

来源:互联网 发布:淘宝怎样延迟久点收货 编辑:程序博客网 时间:2024/05/16 18:18

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;

}

原创粉丝点击