c语言中spline函数 c++ spline

急求二维矩阵三次样条插值的程序,要C或C++语言。望有注释！谢谢！ 邮箱lsmann9@hotmail.com

/*********************下面程序是C语言程序（标准C）******************/

超过10多年行业经验，技术领先，服务至上的经营模式，全靠网络和口碑获得客户，为自己降低成本，也就是为客户降低成本。到目前业务范围包括了：网站建设、网站制作，成都网站推广，成都网站优化，整体网络托管，小程序制作，微信开发，手机APP定制开发，同时也可以让客户的网站和网络营销和我们一样获得订单和生意！

/* 计算给定M0,Mn值的三次样条插值多项式 */

/*给定离散点(1.1,0.4),(1.2,0.8),(1.4,1.65),(1.5,1.8),M0=Mn=0,*/

/*用M关系式构造三次样条插值多项式S(x),计算S(1.25)。 */

/*************************************************************/

#include stdio.h

#define Max_N 20

main()

{int i,k,n;

double h[Max_N+1],b[Max_N+1],c[Max_N+1],d[Max_N+1],M[Max_N+1];

double u[Max_N+1],v[Max_N+1],yy[Max_N+1],x[Max_N+1],y[Max_N+1];

double xx,p,q,S;

printf("\nPlease input n value:"); /*输入插值点数n*/

do

{ scanf("%d",n);

if(nMax_N)

printf("\nplease re-input n value:");

}

while(nMax_N||n=0);

printf("Input x[i],i=0,...%d:\n",n-1);

for(i=0;in;i++) scanf("%lf",x[i]);

printf("Input y[i],i=0,...%d:\n",n-1);

for(i=0;in;i++) scanf("%lf",y[i]);

printf("\nInput the M0,Mn value:");

scanf("%lf%lf",M[0],M[n]);

printf("\nInput the x value:"); /*输入计算三次样条插值函数的x值*/

scanf("%lf",xx);

if((xxx[n-1])||(xxx[0]))

{printf("Please input a number between %f and %f.\n",x[0],x[n-1]);

return;

}

/*计算M关系式中各参数的值*/

h[0]=x[1]-x[0];

for(i=1;in;i++)

{h[i]=x[i+1]-x[i];

b[i]=h[i]/(h[i]+h[i-1]);

c[i]=1-b[i];

d[i]=6*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1])/(h[i]+h[i-1]);

}

/*用追赶法计算Mi,i=1,...,n-1*/

d[1]-=c[1]*M[0];

d[n-1]-=b[n-1]*M[n];

b[n-1]=0; c[1]=0; v[0]=0;

for(i=1;in;i++)

{u[i]=2-c[i]*v[i-1];

v[i]=b[i]/u[i];

yy[i]=(d[i]-c[i]*y[i-1])/u[i];

}

for(i=1;in;i++)

{M[n-i]=yy[n-i]-v[n-i]*M[n-i+1];

}

/*计算三次样条插值函数在x处的值*/

k=0;

while(xx=x[k]) k++;

k=k-1;

p=x[k+1]-xx;

q=xx-x[k];

S=(p*p*p*M[k]+q*q*q*M[k+1])/(6*h[k])+(p*y[k]+q*y[k+1])/h[k]-h[k]*(p*M[k]+q*M[k+1])/6;

printf("S(%f)=%f\n",xx,S); /*输出*/

getch();

}

/*----------------------------------- End of file ------------------------------------*/

/*程序输入输出：

Please input n value:4

Input x[i],i=0,...3:

1.1 1.2 1.4 1.5

Input y[i],i=0,...3:

0.4 0.8 1.65 1.8

Input the M0,Mn value: 0 0

Input the x value:1.25

S(1.250000)=1.033171

*/

谁有MATLAB的插值函数interp1（）的C程序，或者非扭结边界条件下的三次样条插值，

spline函数可以实现三次样条插值

x = 0:10;

y = sin(x);

xx = 0:.25:10;

yy = spline(x,y,xx);

plot(x,y,'o',xx,yy)

另外fnplt csapi这两个函数也是三次样条插值函数，具体你可以help一下！

已知函数值C语言中用三次样条插值求某个点值程序

void SPL(int n, double *x, double *y, int ni, double *xi, double *yi)； 是你所要。

已知 n 个点 x,y; x 必须已按顺序排好。要插值 ni 点，横坐标 xi[], 输出 yi[]。

程序里用double 型，保证计算精度。

SPL调用现成的程序。

现成的程序很多。端点处理方法不同，结果会有不同。想同matlab比较，你需 尝试 调用 spline（）函数 时，令 end1 为 1, 设 slope1 的值，令 end2 为 1 设 slope2 的值。

#include

#include

int spline (int n, int end1, int end2,

double slope1, double slope2,

double x[], double y[],

double b[], double c[], double d[],

int *iflag)

{

int nm1, ib, i, ascend;

double t;

nm1 = n - 1;

*iflag = 0;

if (n 2)

{ /* no possible interpolation */

*iflag = 1;

goto LeaveSpline;

}

ascend = 1;

for (i = 1; i n; ++i) if (x[i] = x[i-1]) ascend = 0;

if (!ascend)

{

*iflag = 2;

goto LeaveSpline;

}

if (n = 3)

{

d[0] = x[1] - x[0];

c[1] = (y[1] - y[0]) / d[0];

for (i = 1; i nm1; ++i)

{

d[i] = x[i+1] - x[i];

b[i] = 2.0 * (d[i-1] + d[i]);

c[i+1] = (y[i+1] - y[i]) / d[i];

c[i] = c[i+1] - c[i];

}

/* ---- Default End conditions */

b[0] = -d[0];

b[nm1] = -d[n-2];

c[0] = 0.0;

c[nm1] = 0.0;

if (n != 3)

{

c[0] = c[2] / (x[3] - x[1]) - c[1] / (x[2] - x[0]);

c[nm1] = c[n-2] / (x[nm1] - x[n-3]) - c[n-3] / (x[n-2] - x[n-4]);

c[0] = c[0] * d[0] * d[0] / (x[3] - x[0]);

c[nm1] = -c[nm1] * d[n-2] * d[n-2] / (x[nm1] - x[n-4]);

}

/* Alternative end conditions -- known slopes */

if (end1 == 1)

{

b[0] = 2.0 * (x[1] - x[0]);

c[0] = (y[1] - y[0]) / (x[1] - x[0]) - slope1;

}

if (end2 == 1)

{

b[nm1] = 2.0 * (x[nm1] - x[n-2]);

c[nm1] = slope2 - (y[nm1] - y[n-2]) / (x[nm1] - x[n-2]);

}

/* Forward elimination */

for (i = 1; i n; ++i)

{

t = d[i-1] / b[i-1];

b[i] = b[i] - t * d[i-1];

c[i] = c[i] - t * c[i-1];

}

/* Back substitution */

c[nm1] = c[nm1] / b[nm1];

for (ib = 0; ib nm1; ++ib)

{

i = n - ib - 2;

c[i] = (c[i] - d[i] * c[i+1]) / b[i];

}

b[nm1] = (y[nm1] - y[n-2]) / d[n-2] + d[n-2] * (c[n-2] + 2.0 * c[nm1]);

for (i = 0; i nm1; ++i)

{

b[i] = (y[i+1] - y[i]) / d[i] - d[i] * (c[i+1] + 2.0 * c[i]);

d[i] = (c[i+1] - c[i]) / d[i];

c[i] = 3.0 * c[i];

}

c[nm1] = 3.0 * c[nm1];

d[nm1] = d[n-2];

}

else

{

b[0] = (y[1] - y[0]) / (x[1] - x[0]);

c[0] = 0.0;

d[0] = 0.0;

b[1] = b[0];

c[1] = 0.0;

d[1] = 0.0;

}

LeaveSpline:

return 0;

}

double seval (int n, double u,

double x[], double y[],

double b[], double c[], double d[],

int *last)

{

int i, j, k;

double w;

i = *last;

if (i = n-1) i = 0;

if (i 0) i = 0;

if ((x[i] u) || (x[i+1] u))

{

i = 0;

j = n;

do

{

k = (i + j) / 2;

if (u x[k]) j = k;

if (u = x[k]) i = k;

}

while (j i+1);

}

*last = i;

w = u - x[i];

w = y[i] + w * (b[i] + w * (c[i] + w * d[i]));

return (w);

}

void SPL(int n, double *x, double *y, int ni, double *xi, double *yi)

{

double *b, *c, *d;

int iflag,last,i;

b = (double *) malloc(sizeof(double) * n);

c = (double *)malloc(sizeof(double) * n);

d = (double *)malloc(sizeof(double) * n);

if (!d) { printf("no enough memory for b,c,d\n");}

else {

spline (n,0,0,0,0,x,y,b,c,d,iflag);

if (iflag==0) printf("I got coef b,c,d now\n"); else printf("x not in order or other error\n");

for (i=0;ini;i++) yi[i] = seval(ni,xi[i],x,y,b,c,d,last);

free(b);free(c);free(d);

};

}

main(){

double x[6]={0.,1.,2.,3.,4.,5};

double y[6]={0.,0.5,2.0,1.6,0.5,0.0};

double u[8]={0.5,1,1.5,2,2.5,3,3.5,4};

double s[8];

int i;

SPL(6, x,y, 8, u, s);

for (i=0;i8;i++) printf("%lf %lf \n",u[i],s[i]);

return 0;

}