C语言实现三次样条插值的子函数

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

创新互联公司自2013年起，是专业互联网技术服务公司，拥有项目网站设计、成都网站设计网站策划，项目实施与项目整合能力。我们以让每一个梦想脱颖而出为使命，1280元井陉做网站,已为上家服务,为井陉各地企业和个人服务,联系电话:18980820575

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

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

SPL调用现成的程序。

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

#include stdio.h

#include math.h

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;

}

分段线性插值问题

#includestdio.h

#includemath.h

double Lagrange1(double *x, double *y, double xx) //拉格郎日插值

{

int i,j;

double *a,yy=0.000;

a=new double[6];

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

{

a[i]=y[i];

for(j=0;j 6;j++)

if(j!=i)

a[i]*=(xx-x[j])/(x[i]-x[j]);

yy+=a[i];

}

delete a;

return yy;

}

double Lagrange2(double *x, double *y, double input) //分段线性插值

{

double output;

int i;

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

{

if (x[i] = input x[i+1] = input)

{

output=y[i] +(y[i+1]-y[i])*(input-x[i])/(x[i+1]-x[i]);

break;

}

}

return output;

}

double Lagrange3(double *x,double *y,double u) //分段二次插值

{

int i,k=0;

double v;

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

{

if(ux[1])

{

k=0;

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]));

}

if((x[i]uu=x[i+1])(fabs(u-x[i])=fabs(u-x[i+1])))

{

k=i-1;

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]));

}

if ((x[i]uu=x[i+1])fabs(u-x[i])fabs(u-x[i+1]))

{

k=i;

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]));

}

if(ux[4])

{

k=3;

v=y[k]*(u-x[k+1])*(u-x[k+2])/((x[k]-x[k+1])*(x[k]-x[k+2]))+y[k+1]*(u-x[k])*(u-x[k+2])/((x[k+1]-x[k])*(x[k+1]-x[k+2]))+y[k+2]*(u-x[k])*(u-x[k+1])/((x[k+2]-x[k])*(x[k+2]-x[k+1]));

}

}

return v;

}

void main()

{

double x[6] = {0.0, 0.1, 0.195, 0.3, 0.401, 0.5},y[6] = {0.39894,0.39695,0.39142,0.38138,0.36812,0.35206};

double u;

scanf("%lf",u);

printf("%f\n",Lagrange1(x,y,u)); //拉格郎日插值

printf("%f\n",Lagrange2(x,y,u)); //分段线性插值

printf("%f\n",Lagrange3(x,y,u)); //分段二次插值

}

怎么用c语言编程一个分段函数？

#include

int main()

{

int x,y;

scanf("%d",x);

if(0xx10) y=3*x+2;

else

{if(x=0) y=0;

else

{if (x0) y=x*x;

else printf("go die\n");

}

}

printf("%d",y);

return 0;

}该程序的分段函数如下:

f(x)=3x+2  (0x10)

f(x)=1         (x=0)

f(x) = x*x    (x0)

#include stdio.h

#include math.h

void main()

{

float x;

double y;

printf("Please input the value of x:");

scanf("%f",x);

if(x=-10x=4)

{

y=fabs(x-2);

printf("y=%.2f\n",y);

}

else if(x=5x=7)

{

y=x+10;

printf("y=%.2f\n",y);

}

else if(x=8x=12)

{

y=pow(x,4);

printf("y=%.2f\n",y);

}

else

printf("No answer\n");

}

C 语言编程 分段抛物线插值

我这里有2个程序,第一个用了2个函数

第二个用了1个函数,感觉误差小些

说下用法吧

先输入数据 两个两个的输,中间用空格隔开,比如

please input data1: 11 11.08然后回车

依次输入完五组数据

完了你想查某个温度的溶解度,他要求你输入温度你输入11.5

它就输出对应的溶解度,然后他提示你是否继续查溶解度,是就输入y,想结束程序

就输入n.

#includestdio.h

#includestdlib.h

#define NUMBER 5

typedef struct

{

double x;

double y;

}Point;

double * parabola(Point*, Point*, Point*);

double calculate(double* , double );

int main()

{

double x = 0, y = 0;

double* f = NULL;

int i = 0, n = 0;//n为要插的中间那个点的位置

char c = 'y';

Point p[NUMBER];

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

{

printf("Please input the point%d:", i+1);

scanf("%lf%lf", p[i].x, p[i].y);

while(10 != getchar())

{

continue;

}

}

while('n' != c)

{

printf("Please input the temperature:");

scanf("%lf", x);

while(10 != getchar())

{

continue;

}

//计算插值点的位子

if(2 * x = (p[1].x + p[2].x)) n = 1;

else if(2 * x = (p[NUMBER-2].x + p[NUMBER-1].x)) n = NUMBER - 2;

else n = (int)NUMBER / 2;

printf("%d\n", n);

f = parabola(p[n-1],p[n],p[n+1]);//计算抛物线方程

y = calculate(f,x);//计算溶解度

printf("The solubility at this temperature is %lf\n", y);

printf("Continue?(y/n) ");

scanf("%c", c);

while(10 != getchar())

{

continue;

}

//释放内存

free(f);

}

return 0;

}

//下面为计算抛物线方程的函数

//设抛物线函数为y = a * x^2 + b * x + c

//以下f[0]为a,f[1]为b,f[2]为c

double * parabola(Point* p1, Point* p2, Point* p3)

{

double temp1 = 0, temp2 = 0;

double * f = NULL;

f = (double*)calloc(3,sizeof(double));

if(NULL == f)

{

printf("Calloc failed!\n");

return NULL;

}

temp1 = (p2-y - p1-y)/(p2-x - p1-x);

temp2 = (p3-y - p2-y)/(p3-x - p2-x);

f[0] = (temp2 - temp1)/(p3-x - p1-x);

f[1] = temp1 - f[0] * (p1-x + p2-x);

f[2] = p1-y - p1-x * (p1-x * f[0] + f[1]);

return f;

}

//根据温度计算溶解度

double calculate(double* f, double x)

{

double y;

y = x * (f[0] * x + f[1]) + f[2];

return y;

}

//这里是第二个程序

#includestdio.h

#includestdlib.h

#define NUMBER 5

typedef struct

{

double x;

double y;

}Point;

double parabola(Point*, Point*, Point*, double);

int main()

{

double x = 0, y = 0;

int i = 0, n = 0;//n为要插的中间那个点的位置

char c = 'y';

Point p[NUMBER];

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

{

printf("Please input the data%d:", i+1);

scanf("%lf%lf", p[i].x, p[i].y);

while(10 != getchar())

{

continue;

}

}

while('n' != c)

{

printf("Please input the temperature:");

scanf("%lf", x);

while(10 != getchar())

{

continue;

}

//计算插值点的位子

if(2 * x (p[1].x + p[2].x)) n = 1;

else if(2 * x (p[NUMBER-2].x + p[NUMBER-1].x)) n = NUMBER - 2;

else n = (int)NUMBER / 2;

//printf("%d\n", n);

y = parabola(p[n-1],p[n],p[n+1],x);

printf("The solubility at this temperature is %lf\n", y);

printf("Continue?(y/n) ");

scanf("%c", c);

while(10 != getchar())

{

continue;

}

}

return 0;

}

//直接套用公式计算溶解度

double parabola(Point* p1, Point* p2, Point* p3, double x)

{

double temp1 = 0, temp2 = 0, temp3 = 0, y = 0;

double a = 0, b = 0, c = 0;

temp1 = (p2-y - p1-y) * (x - p1-x);

temp2 = (p3-y - p1-y) * (p2-x - p1-x) - (p2-y - p1-y) * (p3-x - p1-x);

temp3 = (p3-x - p1-x) * (p3-x - p1-x) * (p2-x - p1-x);

y = p1-y + temp1 / (p2-x - p1-x) + temp2 * (x - p1-x) * (x - p2-x) / temp3;

return y;

}

该程序误差绝对小 但结果和你不是完全符合

你的数据好象有点问题 14度的溶解度不太正常

当前名称：c语言分段插值的子函数 c语言分段函数的计算
标题URL：http://scyanting.com/article/doeidpd.html