三次样条.cpp

最新版三次样条插值源程序

#include <stdio.h> 
#include <math.h> 
#include <memory.h> 

void swap(double & a, double & b) 
{ 
double temp; 
temp=a; 
a=b; 
b=temp; 
} 

bool GuassSolve(double *A,double *B,double *X,int n) 
{ 
int i,j,k,row; 
for(k=0;k<n-1;k++) 
{ 
double max=0.0; 
for(i=k;i<n;i++) 
{ 
if(fabs(A[i*n+k])>max) 
{ 
max=fabs(A[i*n+k]); 
row=i; 
} 
} 
if(max==0.0) return false; 
if(row!=k) 
{ 
for(j=0;j<n;j++) 
swap(A[row*n+j],A[k*n+j]); 
swap(B[row],B[k]); 
} 
for(i=k+1;i<n;i++) 
{ 
double m=A[i*n+k]/A[k*n+k]; 
for(j=k+1;j<n;j++) 
{ 
A[i*n+j]=A[i*n+j]-m*A[k*n+j]; 
} 
B[i]=B[i]-m*B[k]; 
} 
} 
X[n-1]=B[n-1]/A[n*n-1]; 
for(k=n-2;k>=0;k--) 
{ 
double sum=0.0; 
for(j=k+1;j<n;j++) 
sum+=A[k*n+j]*X[j]; 
X[k]=(B[k]-sum)/A[k*n+k]; 
} 
return true; 
} 

void Apend(double miu[],double rmd[],double d[],double x[],double y[],double h 
[], 
  double pf1,double pf2, 
  double ppf1,double ppf2, 
  int n,int flag) 
{ 
switch(flag) 
{ 
case 0: 
miu[n-1]=1; 
rmd[0]=1; 
d[0]=6*((y[1]-y[0])/(x[1]-x[0])-pf1)/h[0]; 
d[n-1]=6*(pf2-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/h[n-2]; 
break; 
case 1: 
rmd[0]=miu[n-1]=0; 
d[0]=2*ppf1; 
d[n-1]=2*ppf2; 
break; 
case 2: 
rmd[0]=rmd[n-1]=h[0]/(h[n-1]+h[0]); 
miu[0]=miu[n-1]=1-rmd[n-1]; 
d[0]=d[n-1]=6*((y[1]-y[0])/(x[1]-x[0])-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/(h[0 
]+h[n-2]); 
break; 
} 
} 

void T(double x[],double y[],int n) 
{ 
int flag=0; 
double   pf1=3.0; 
double   pf2=-4.0; 
double   ppf1=0; 
double   ppf2=0; 
double * h=new double [n]; 
double * miu=new double [n]; 
double * rmd=new double [n]; 
double * d=new double[n]; 

memset(d,0,n*sizeof(double)); 
memset(rmd,0,n*sizeof(double)); 
memset(miu,0,n*sizeof(double)); 

for(int i=1;i<n;i++) 
{ 
h[i-1]=x[i]-x[i-1]; 
} 

for(i=1;i<n-1;i++) 
{ 
rmd[i]=h[i]/(h[i-1]+h[i]); 
miu[i]=h[i-1]/(h[i-1]+h[i]); 
d[i]=6*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/(h[i-1]+h[i 
]); 
} 

Apend(miu,rmd,d,x,y,h,pf1,pf2,ppf1,ppf2,n,0); 

double * A=new double [n*n],* M=new double [n]; 

memset(A,0,n*n*sizeof(double)); 

for(i=1;i<n-1;i++) 
{ 
A[i*n+i]=2; 

A[i*n+i-1]=miu[i]; 

A[i*n+i+1]=rmd[i]; 
} 

if(n==2) 
{ 
A[0]=2; 
A[n-2]=miu[0]; 
A[1]=rmd[0]; 
A[(n-1)*n+n-1]=2; 
A[(n-1)*n+n-2]=miu[n-1]; 
A[(n-1)*n+1]=rmd[n-1]; 
} 
else 
{ 
A[0]=2; 
A[1]=rmd[0]; 
A[(n-1)*n+n-1]=2; 
A[(n-1)*n+n-2]=miu[n-1]; 
} 

GuassSolve(A,d,M,n); 

for(i=0;i<n;i++) 
{ 
printf("M[%d]=%f\n",i+1,M[i]); 
} 

} 

void Initiate(double * x,double * y) 
{ 
x[0]=27.7; 
y[0]=4.1; 

x[1]=28; 
y[1]=4.3; 

x[2]=29; 
y[2]=4.1; 

x[3]=30; 
y[3]=3.0; 
} 

void main() 
{ 
int n=4; 
double * x=new double[n]; 
double * y=new double[n]; 

Initiate(x,y); 

T(x,y,n); 
}