三次样条插值函数c语言程序

#include
#include
using namespace std;

const int MAX = 50;
float x[MAX], y[MAX], h[MAX];
float c[MAX], a[MAX], fxym[MAX];

float f(int x1, int x2, int x3){
float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);
float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);
return (a - b)/(x[x3] - x[x1]);
} //求差分
void cal_m(int n){ //用追赶法求解出弯矩向量M……
float B[MAX];
B[0] = c[0] / 2;
for(int i = 1; i < n; i++)
B[i] = c[i] / (2 - a[i]*B[i-1]);

fxym[0] = fxym[0] / 2;
for(i = 1; i <= n; i++)
fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);
for(i = n-1; i >= 0; i--)
fxym[i] = fxym[i] - B[i]*fxym[i+1];
}
void printout(int n);

int main(){
int n,i; char ch;
do{
cout<<"Please put in the number of the dots:";
cin>>n;
for(i = 0; i <= n; i++){
cout<<"Please put in X"<cin>>x[i]; //cout<cout<<"Please put in Y"<cin>>y[i]; //cout<}
for(i = 0; i < n; i++) //求 步长
h[i] = x[i+1] - x[i];
cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";
int t;
float f0, f1;
cin>>t;
switch(t){
case 1:cout<<"Please put in Y0\' Y"<cin>>f0>>f1;
c[0] = 1; a[n] = 1;
fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];
fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];
break;
case 2:cout<<"Please put in Y0\" Y"<cin>>f0>>f1;
c[0] = a[n] = 0;
fxym[0] = 2*f0; fxym[n] = 2*f1;
break;
default:cout<<"不可用\n";//待定
};//switch
for(i = 1; i < n; i++)
fxym[i] = 6 * f(i-1, i, i+1);
for(i = 1; i < n; i++){
a[i] = h[i-1] / (h[i] + h[i-1]);
c[i] = 1 - a[i];
}
a[n] = h[n-1] / (h[n-1] + h[n]);
cal_m(n);
cout<<"\n输出三次样条插值函数：\n";
printout(n);

cout<<"Do you to have anther try ? y/n :";
cin>>ch;
}while(ch == 'y' || ch == 'Y');
return 0;
}
void printout(int n){
cout<for(int i = 0; i < n; i++){
cout</*
cout<<<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<cout<float t = fxym[i]/(6*h[i]);
if(t > 0)cout<else cout<<-t<<"*(x - "<t = fxym[i+1]/(6*h[i]);
if

(t > 0)cout<<" + "<else cout<<" - "<<-t<<"*(x - "<cout<<"\n\t";
t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];
if(t > 0)cout<<"+ "<else cout<<"- "<<-t<<"*("<t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];
if(t > 0)cout<<" + "<else cout<<" - "<<-t<<"*(x - "<cout<}
cout<}


























#include "math.h "
#include "stdio.h "
#include "stdlib.h "


/*
N：已知节点数N+1
R：欲求插值点数R+1
x，y为给定函数f(x)的节点值{x(i)} (x(i) P0=f(x0)的二阶导数；Pn=f(xn)的二阶导数
u:存插值点{u(i)} 0 <=i <=R
求得的结果s(ui)放入s[R+1] 0 <=i <=R
返回0表示成功，1表示失败
*/
int SPL(int N,int R,double x[],double y[],double P0,double Pn,double u[],double s[])
{
/*声明局部变量*/
double *h; /*存放步长:{hi} 0 <=i <=N-1 */
double *a; /*存放系数矩阵{ai} 1 <=i <=N ； 分量0没有利用 */
double *c; /*先存放系数矩阵{ci} 后存放{Bi} 0 <=i <=N-1 */
double *g; /*先存放方程组右端项{gi} 后存放求解中间结果{yi} 0 <=i <=N */
double *af; /*存放系数矩阵{a(f)i} 1 <=i <=N ； */
double *ba; /*存放中间结果 0 <=i <=N-1*/
double *m; /*存放方程组的解{m(i)} 0 <=i <=N ； */

int i,k;
double p1,p2,p3,p4;

/*分配空间*/
if(!(h=(double*)malloc(N*sizeof(double)))) exit(1);
if(!(a=(double*)malloc((N+1)*sizeof(double)))) exit(1);
if(!(c=(double*)malloc(N*sizeof(double)))) exit(1);
if(!(g=(double*)malloc((N+1)*sizeof(double)))) exit(1);
if(!(af=(double*)malloc((N+1)*sizeof(double)))) exit(1);
if(!(ba=(double*)malloc((N)*sizeof(double)))) exit(1);
if(!(m=(double*)malloc((N+1)*sizeof(double)))) exit(1);

/*第一步：计算方程组的系数*/
for(k=0;k h[k]=x[k+1]-x[k];
for(k=1;k a[k]=h[k]/(h[k]+h[k-1]);
for(k=1;k c[k]=1-a[k];
for(k=1;k g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);
c[0]=a[N]=1;
g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;
g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;

/*第二步：用追赶法解方程组求{m(i)} */
ba[0]=c[0]/2;
g[0]=g[0]/2;
for(i=1;i {
af[i]=2-a[i]*ba[i-1];
g[i]=(g[i]-a[i]*g[i-1])/af[i];
ba[i]=c[i]/af[i];
}
af[N]=2-a[N]*ba[N-1];
g[N]=(g[N]-a[N]*g[N-1])/af[N];

m[N]=g[N]; /*P110 公式：6.32*/
for(i=N-1;i> =0;i--)
m[i]=g[i]-ba[i]*m[i+1];

/*第三步：求值*/
for(i=0;i <=R;i++)
{
/*判断u(i)属于哪一个子区间，即确定k */
if(u[i] x[N])
{

/*释放空间*/
free(h);
free(a);
free(c);
free(g);
free(af);
free(ba);
free(m);
return 1;
}
k=0;
while(u[i]> x[k+1])
k++;


p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3);
p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);
p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);
p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);
s[i]=p1+p2+p3+p4;
}

/*释放空间*/
free(h);
free(a);
free(c);
free(g);
free(af);
free(ba);
free(m);

return 0;
}


void main()
{
int N,R;
double *x,*y,*u,*s;
double P0,Pn;
int i;

/*验证算法：*/
N=7;
R=6;

/*分配空间*/
if(!(x=(double*)malloc((N+1)*sizeof(double))))
{
printf( "malloc error!\n ");
exit(1);
}
if(!(y=(double*)malloc((N+1)*sizeof(double))))
{
printf( "malloc error!\n ");
exit(1);
}
if(!(u=(double*)malloc((R+1)*sizeof(double))))
{
printf( "malloc error!\n ");
exit(1);
}
if(!(s=(double*)malloc((R+1)*sizeof(double))))
{
printf( "malloc error!\n ");
exit(1);
}

x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;
y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0.9463;
u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;

P0=-0.4794;
Pn=-0.9463;

if(!SPL( N, R, x, y, P0, Pn, u, s))
{
/*打印结果*/
printf( "\nx= ");
for(i=0;i <=N;i++)
printf( "%8.1f ",x[i]);

printf( "\ny= ");
for(i=0;i <=N;i++)
printf( "%8.4f ",y[i]);

printf( "\n\nu= ");
for(i=0;i <=R;i++)
printf( "%9.2f ",u[i]);
printf( "\ns= ");
for(i=0;i <=R;i++)
printf( "%9.5f ",s[i]);
printf( "\nsin= ");
for(i=0;i <=R;i++)
printf( "%9.5f ",sin(u[i]));
}

/*释放空间*/
free(x);
free(y);
free(u);
free(s);