完备三次样条插值.cpp

数值计算的一些经典算法

#include<iostream.h>
#include<math.h>
#define F(x) (1/(1+25*(x)*(x)))
#define F1(x) -50*x/((1+25*x*x)*(1+25*x*x))
#define N 32
#define M 200
#define TOL 1.0e-8 
void dca1(double h[N],double d[N+1],double a[N+1],double c[N]);
void pursue(double d[N+1],double c[N],double a[N+1],double p[N+1],double q[N]);
void back(double x[N+1],double b[N+1],double p[N+1],double q[N],double a[N+1]);
double g(double x[N+1],double y,double b[N+1],double h[N],double m[N+1]);
void error(double a[]);
void main()
{
   double x[33],f[M+1],e[N+1],b[N+1],h[N],d[N+1],a[N+1],c[N+1],p[N+1],q[N+1],m[N+1];
	for(int i=0;i<=N;i++)
		x[i]=-1+(i)*2./N;
	for(i=0;i<N;i++) h[i]=2./N;
    for(i=0;i<N;i++)
        e[i]=(F(x[i+1])-F(x[i]))/h[i];
    b[0]=6*(e[0]-F1(x[0]));
	b[N]=6*(F1(x[N])-e[N-1]);
	for(i=1;i<N;i++)
	b[i]=6*(e[i]-e[i-1]);
	dca1(h,d,a,c);
	pursue(d,c,a,p,q);
    back(m,b,p,q,a);
	for(int j=0;j<=M;j++)
	{
		double t=g(x,(-1+j*2./M),e,h,m);
        f[j]=fabs(F(-1+j*2./M)-t);
	}
    error(f);
}
void dca1(double h[N],double d[N+1],double a[N+1],double c[N])//构造主对角线元素和次对角线元素 
{ d[0]=2*h[0],d[N]=2*h[N-1],a[N]=h[N-1],c[0]=h[0]; 
   for(int i=1;i<N;i++) 
   { d[i]=2*(h[i-1]+h[i]); 
     c[i]=h[i]; 
     a[i]=h[i-1]; 
   } 
} 
void pursue(double d[N+1],double c[N],double a[N+1],double p[N+1],double q[N])//Crout分解 
{ if(fabs(d[0])<TOL) {cout<<"Method failed"; return;} 
  p[0]=d[0]; 
  q[0]=c[0]/d[0]; 
  for(int k=1;k<N;k++) 
   { p[k]=d[k]-a[k]*q[k-1]; 
    if(fabs(p[k])<TOL) {cout<<"Method failed"; return;} 
    q[k]=c[k]/p[k]; 
   } 
   p[N]=d[N]-a[N]*q[N-1]; 
   if(fabs(p[N])<TOL) 
  {cout<<"Method failed"; return;} 
} 
void back(double x[N+1],double b[N+1],double p[N+1],double q[N],double a[N+1])//回代过程 
{ double y[N+1]; 
  y[0]=b[0]/p[0];
  for(int k=1;k<=N;k++) 
    { y[k]=(b[k]-a[k]*y[k-1])/p[k]; }
   x[N]=y[N]; 
   for(k=N-1;k>=0;k--) x[k]=y[k]-q[k]*x[k+1]; 
} 
double g(double x[N+1],double y,double b[N+1],double h[N],double m[N+1])
{   int c;
    for(int i=0;i<=N;i+=1)
		if(y>=x[i]&&y<=x[i+1]) {c=i; break;}
	double t=F(x[c])+(b[c]-(h[c]/6)*(m[c+1]+2*m[c]))*(y-x[c])+\
		     (m[c]/2)*(y-x[c])*(y-x[c])+((m[c+1]-m[c])/(6*h[c]))*(y-x[c])*(y-x[c])*(y-x[c]);
	cout<<y<<","<<t<<endl;
	return t;   
}
void error(double a[])
{   double max=a[0];
    int t;
	for(int i=1;i<=M;i++)
		if(a[i]>=max) {max=a[i],t=i;}
	cout<<"N= "<<N<<'\t';
	cout<<"x= "<<-1+(t)*2./M<<'\t';
	cout<<"errorMax= "<<max<<endl;
}