三次样条插值计算算法.docVIP

/* 三次样条插值计算算法 */#include math.h #include stdio.h #include stdlib.h /*N：已知节点数N+1R：欲求插值点数R+1x，y为给定函数f(x)的节点值{x(i)} (x(i) x(i+1)) ,以及相应的函数值{f(i)} 0 =i =NP0=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 N;k++)h[k]=x[k+1]-x[k];for(k=1;k N;k++)a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k N;k++)c[k]=1-a[k];for(k=1;k N;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 N;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[0] || u[i] x[N]) {/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m)