牛顿插值 三次样条插值 错误 为何?
#include <math.h>#include <stdio.h>
#define n 6
double YY(double x)
{
double y;
y=1.0/(1.0+25.0*x*x);
return y;
}
double NN(double x[],double y[],double xx)
{
int i,k;
double yy;
for(k=1;k<n;k++)
{
for(i=n-1;i>=k;i--)
y[i]=(y[i]-y[i-1])/(x[i]-x[i-k]);
}
yy=y[n-1];
for(i=n-2;i>=0;i--)
yy=yy*(xx-x[i])+y[i];
return yy;
}
double SS(double x[],double y[],double xx)
{
double M[n],h[n-1],c [n-1],a[n-1],b[n],d[n],u[n],mm[n],l[n-1],m[n],h1,x1,x2,yy;
int i,k,K=1;
for(i=0;i<n;i++)
M[i]=y[i];
for(k=1;k<3;k++)
{for(i=n-1;i>=k;i--)
M[i]=(M[i]-M[i-1])/(x[i]-x[i-k]);
}
h[0]=x[1]-x[0];
for(i=1;i<n-1;i++)
{h[i]=x[i+1]-x[i];
c[i]=h[i]/(h[i-1]+h[i]);
a[i-1]=1-c[i];
b[i]=2.0;
M[i]=6*M[i+1];
}
M[0]=0.0;M[n-1]=0.0;c[0]=0.0;b[0]=2.0;a[n-2]=0.0;b[n-1]=2.0;
for(i=0;i<n;i++)
d[i]=M[i];
u[0]=b[0];
mm[0]=d[0];
for(k=1;k<n;k++)
{
l[k-1]=a[k-1]/u[k-1];
u[k]=b[k]-l[k-1]*c[k-1];
mm[k]=d[k]-l[k-1]*mm[k-1];
}
m[n-1]=mm[n-1]/u[n-1];
for(k=n-2;k>=0;k--)
m[k]=(mm[k]-c[k]*m[k+1])/u[k];
for(i=0;i<n-1;i++)
{if(xx<=x[i])
{K=i;break;}
else K=i+1;
}
h1=x[K]-x[K-1];
x1=x[K]-xx;
x2=xx-x[K-1];
yy=((m[K-1]*x1*x1*x1)/6+(m[K]*x2)/6+y[K-1]*x1-(m[K-1]*x1*h1*h1)/6+y[K]*x2-(m[K]*x2*h1*h1)/6)/h1;
return yy;
}
main()
{int i,k=99;
double x[n],y[n],X[99],Y[99],N[99],S[99],EN,ES;
for(i=0;i<n;i++)
{x[i]=-1.0+(2.0*i)/(n-1);
y[i]=YY(x[i]);
}
for(i=0;i<n;i++)
printf("x[%d]=%.4f ",i,x[i]);
printf("\n");
for(i=0;i<n;i++)
printf("y[%d]=%.4f ",i,y[i]);
printf("\n");
for(i=0;i<k;i++)
{X[i]=-1.0+0.02*(i+1);
Y[i]=YY(X[i]);
N[i]=NN(x,y,X[i]);
S[i]=SS(x,y,X[i]);
}
EN=fabs(Y[0]-N[0]);
ES=fabs(Y[0]-S[0]);
for(i=1;i<k;i++)
{if(fabs(Y[i]-N[i])>EN)
EN=fabs(Y[i]-N[i]);
if(fabs(Y[i]-S[i])>ES)
ES=fabs(Y[i]-S[i]);
}
printf("X Y N S \n");
for(i=0;i<=10;i++)
printf("%f %f %f %f\n",X[i],Y[i],N[i],S[i]);
printf("EN=%f ES=%f\n",EN,ES);
getch();}
为何运行结果:越来越大。没搞明白。。。。。。
是newton插值和三次样条插值
请教高手 先谢了