三次样条差值C语言程序，工程常用算法

/* 三次样条插值计算算法 */ #include "math.h" #include "stdio.h" #include "stdlib.h" /* N：已知节点数N+1 R：欲求插值点数R+1 x，y为给定函数f(x)的节点值{x(i)} (x(i)<x(i+1)) ,以及相应的函数值{f(i)} 0<=i<=N 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<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); 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); } Top