|
牛頓插值的關(guān)鍵在于差商表的計算，差商表第一行是y值，為了配合計算，在該矩陣上方配上節(jié)點x0、x1、x2……xnf[x0,x1]=[f(x1)?f(x0)]/(x1?x0)
f[x0,x1,x2]=[f[x1,x2]?f[x0,x1]]/(x2?x0)……
所以只要計算矩陣內(nèi)上三角值即可。
#include <stdio.h>
#include <stdlib.h>
int main()
{
float table(int n,float a1[10],float a2[10],float a3[10][10]);
float newton(int n,float a4[10][10],float a5[10]);
float arrX[10],arrY[10],arrL[10][10];
int num,i;
printf("請輸入插值節(jié)點的個數(shù)（個數(shù)應(yīng)小于10）：");
scanf("%d",&num);
printf("請輸入各個插值節(jié)點的值：\n");
for(i=0; i<num; i++)
{
printf("請輸入X%d值：",i+1);
scanf("%f",&arrX[i]);
printf("請輸入Y%d值：",i+1);
scanf("%f",&arrY[i]);
}
table(num,arrX,arrY,arrL);
newton(num,arrL,arrX);
return 0;
}
float table(int n,float a1[10],float a2[10],float a3[10][10])
{
int i,j;
for(i=0; i<n; i++)
{
a3[0][i]=a2[i];//第一行初始化為y值
}
for(i=0; i<n; i++)
for(j=n-1; j>i; j--)//從一行最后往前循環(huán)，到i=j為止，即上三角全部計算賦值
{
a3[i+1][j]=(a3[i][j]-a3[i][j-1])/(a1[j]-a1[j-1*(i+1)]);//差商表計算，最后一項arrX[j-1*(i+1)]腳標(biāo)是計算步長
}
for(i=1; i<n; i++)//下三角未計算賦值，系統(tǒng)會隨機分配值給下三角的每一位，故賦值0會使差商表輸出更整齊
for(j=0; j<i; j++)
{
a3[i][j]=0.0;
}
printf("差商表為：\n");
printf("----------------------------------------------------------\n");
for(i=0; i<n; i++)
{
for(j=0; j<n; j++)
{
if(j%n==0)//num個數(shù)一行輸出
printf("\n");
printf("%f\t",a3[i][j]);
}
}
printf("\n");
printf("----------------------------------------------------------\n");
return 0;
}
float newton(int n,float a4[10][10],float a5[10])
{
int i;
float x,y,t1=1.0;
while(1)
{
printf("請輸入要插入節(jié)點的X值：");
scanf("%f",&x);
y=a4[0][0];
for(i=1; i<n; i++)
{
t1=t1*(x-a5[i-1]);
y=y+a4[i][i]*t1;//差商表對角線值依次乘以(x-x0)(x-x1)……
}
printf("插值結(jié)果為：%f",y);
printf("\n");
}
return 0;
}
運行結(jié)果如下：
|
