新三源码x3ym

1. Akima 插值和样条插值的C语言源代码，要有注释。

Akima插值

附带的图片为运行结果

#include"stdio.h"

#include"math.h"

#include"interpolation.h"

voidinterpolation_akima(AKINTEPap){

intnum,k,kk,m,l;

doublepio,*mtr,*x,*y,u[5],p,q;

num=ap->n;k=ap->k;

pio=ap->t;mtr=ap->s;

x=ap->x;y=ap->y;

if(num<1){

return;

}

elseif(num==1){

mtr[0]=mtr[4]=y[0];

return;

}

elseif(num==2){

mtr[0]=y[0];

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

if(k<0)

mtr[4]=(y[0]*(pio-x[1])-y[1]*(pio-x[0]))/(x[0]-x[1]);

return;

}

if(k<0){

if(pio<=x[1])kk=0;

elseif(pio>=x[num-1])kk=num-2;

else{

kk=1;m=num;

while(((kk-m)!=1)&&((kk-m)!=-1)){

l=(kk+m)/2;

if(pio<x[l-1])m=l;

elsekk=l;

}

kk--;

}

}

elsekk=k;

if(kk>=num-1)kk=num-2;

u[2]=(y[kk+1]-y[kk])/(x[kk+1]-x[kk]);

if(num==3){

if(kk==0){

u[3]=(y[2]-y[1])/(x[2]-x[1]);

u[4]=2.0*u[3]-u[2];

u[1]=2.0*u[2]-u[3];

u[0]=2.0*u[1]-u[2];

}

else{

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

u[0]=2.0*u[1]-u[2];

u[3]=2.0*u[2]-u[1];

u[4]=2.0*u[3]-u[2];

}

}

else{

if(kk<=1){

u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);

if(kk==1){

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

u[0]=2.0*u[1]-u[2];

if(num==4)u[4]=2.0*u[3]-u[2];

elseu[4]=(y[4]-y[3])/(x[4]-x[3]);

}

else{

u[1]=2.0*u[2]-u[3];

u[0]=2.0*u[1]-u[2];

u[4]=(y[3]-y[2])/(x[3]-x[2]);

}

}

elseif(kk>=(num-3)){

u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);

if(kk==(num-3)){

u[3]=(y[num-1]-y[num-2])/(x[num-1]-x[num-2]);

u[4]=2.0*u[3]-u[2];

if(num==4)u[0]=2.0*u[1]-u[2];

elseu[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);

}

else{

u[3]=2.0*u[2]-u[1];

u[4]=2.0*u[3]-u[2];

u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);

}

}

else{

u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);

u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);

u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);

u[4]=(y[kk+3]-y[kk+2])/(x[kk+3]-x[kk+2]);

}

}

mtr[0]=fabs(u[3]-u[2]);

mtr[1]=fabs(u[0]-u[1]);

if((fabs(mtr[0])<0.0000001)&&(fabs(mtr[1])<0.0000001))

p=(u[1]+u[2])/2.0;

elsep=(mtr[0]*u[1]+mtr[1]*u[2])/(mtr[0]+mtr[1]);

mtr[0]=fabs(u[3]-u[4]);

mtr[1]=fabs(u[2]-u[1]);

if((fabs(mtr[0])<0.0000001)&&(fabs(mtr[1])<0.0000001))

q=(u[2]+u[3])/2.0;

elseq=(mtr[0]*u[2]+mtr[1]*u[3])/(mtr[0]+mtr[1]);

mtr[0]=y[kk];

mtr[1]=p;

mtr[3]=x[kk+1]-x[kk];

mtr[2]=(3.0*u[2]-2.0*p-q)/mtr[3];

mtr[3]=(q+p-2.0*u[2])/(mtr[3]*mtr[3]);

if(k<0){

p=pio-x[kk];

mtr[4]=mtr[0]+mtr[1]*p+mtr[2]*p*p+mtr[3]*p*p*p;

}

return;

}

main()

{

doublex[11]={3.0,5.0,8.0,13.0,17.0,25.0,27.0,29.0,31.0,35.0,39.0};

doubley[11]={7.0,10.0,11.0,17.0,23.0,18.0,13.0,6.0,3.0,1.0,0.0};

AKINTEaa={11,x,y,-1,14.0,{0}};

AKINTEab={11,x,y,-1,28.0,{0}};

printf(" ");

interpolation_akima(&aa);

printf("x=%6.3f,f(x)=%e ",aa.t,aa.s[4]);

printf("mtr0=%e,mtr1=%e,mtr2=%e,mtr3=%e ",aa.s[0],aa.s[1],aa.s[2],aa.s[3]);

printf(" ");

interpolation_akima(&ab);

printf("x=%6.3f,f(x)=%e ",ab.t,ab.s[4]);

printf("mtr0=%e,mtr1=%e,mtr2=%e,mtr3=%e ",ab.s[0],ab.s[1],ab.s[2],ab.s[3]);

printf(" ");

}

三次样条插值的实现

1、程序比较简单的：

#include<iostream>

#include<iomanip>

usingnamespacestd;

constintMAX=50;

floatx[MAX],y[MAX],h[MAX];

floatc[MAX],a[MAX],fxym[MAX];

floatf(intx1,intx2,intx3){

floata=(y[x3]-y[x2])/(x[x3]-x[x2]);

floatb=(y[x2]-y[x1])/(x[x2]-x[x1]);

return(a-b)/(x[x3]-x[x1]);

}//求差分

voidcal_m(intn){//用追赶法求解出弯矩向量M……

floatB[MAX];

B[0]=c[0]/2;

for(inti=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];

}

voidprintout(intn);

intmain(){

intn,i;charch;

do{

cout<<"Pleaseputinthenumberofthedots:";

cin>>n;

for(i=0;i<=n;i++){

cout<<"PleaseputinX"<<i<<':';

cin>>x[i];//cout<<endl;

cout<<"PleaseputinY"<<i<<':';

cin>>y[i];//cout<<endl;

}

for(i=0;i<n;i++)//求步长

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

cout<<"Please输入边界条件 1:已知两端的一阶导数 2:两端的二阶导数已知 默认:自然边界条件 ";

intt;

floatf0,f1;

cin>>t;

switch(t){

case1:cout<<"PleaseputinY0'Y"<<n<<"' ";

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;

case2:cout<<"PleaseputinY0"Y"<<n<<"" ";

cin>>f0>>f1;

c[0]=a[n]=0;

fxym[0]=2*f0;fxym[n]=2*f1;

break;

default:cout<<"不可用 ";//待定

};//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<<" 输出三次样条插值函数： ";

printout(n);

cout<<"Doyoutohaveanthertry?y/n:";

cin>>ch;

}while(ch=='y'||ch=='Y');

return0;

}

voidprintout(intn){

cout<<setprecision(6);

for(inti=0;i<n;i++){

cout<<i+1<<":["<<x[i]<<","<<x[i+1]<<"] "<<" ";

/*

cout<<fxym[i]/(6*h[i])<<"*("<<x[i+1]<<"-x)^3+"<<<<"*(x-"<<x[i]<<")^3+"

<<(y[i]-fxym[i]*h[i]*h[i]/6)/h[i]<<"*("<<x[i+1]<<"-x)+"

<<(y[i+1]-fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x-"<<x[i]<<") ";

cout<<endl;*/

floatt=fxym[i]/(6*h[i]);

if(t>0)cout<<t<<"*("<<x[i+1]<<"-x)^3";

elsecout<<-t<<"*(x-"<<x[i+1]<<")^3";

t=fxym[i+1]/(6*h[i]);

if(t>0)cout<<"+"<<t<<"*(x-"<<x[i]<<")^3";

elsecout<<"-"<<-t<<"*(x-"<<x[i]<<")^3";

cout<<" ";

t=(y[i]-fxym[i]*h[i]*h[i]/6)/h[i];

if(t>0)cout<<"+"<<t<<"*("<<x[i+1]<<"-x)";

elsecout<<"-"<<-t<<"*("<<x[i+1]<<"-x)";

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

if(t>0)cout<<"+"<<t<<"*(x-"<<x[i]<<")";

elsecout<<"-"<<-t<<"*(x-"<<x[i]<<")";

cout<<endl<<endl;

}

cout<<endl;

}

2、程序比较复杂的：

（程序前面的01.，02.，03.等等为语句编号，实际应用时请一一删除）01./*=======================================================================*/

02.#include<stdio.h>

03.////////////////////////////////////////////////////////////////////////////////

04.#defineMAXNUM50//定义样条数据区间个数最多为50个

05.typedefstructSPLINE//定义样条结构体，用于存储一条样条的所有信息

06.{//初始化数据输入

07.floatx[MAXNUM+1];//存储样条上的点的x坐标，最多51个点

08.floaty[MAXNUM+1];//存储样条上的点的y坐标，最多51个点

09.unsignedintpoint_num;//存储样条上的实际的点的个数

10.floatbegin_k1;//开始点的一阶导数信息

11.floatend_k1;//终止点的一阶导数信息

12.//floatbegin_k2;//开始点的二阶导数信息

13.//floatend_k2;//终止点的二阶导数信息

14.//计算所得的样条函数S(x)

15.floatk1[MAXNUM+1];//所有点的一阶导数信息

16.floatk2[MAXNUM+1];//所有点的二阶导数信息

17.//51个点之间有50个段，func[]存储每段的函数系数

18.floata3[MAXNUM],a1[MAXNUM];

19.floatb3[MAXNUM],b1[MAXNUM];

20.//分段函数的形式为Si(x)=a3[i]*{x(i+1)-x}^3+a1[i]*{x(i+1)-x}+

21.//b3[i]*{x-x(i)}^3+b1[i]*{x-x(i)}

22.//xi为x[i]的值，xi_1为x[i+1]的值

23.}SPLINE,*pSPLINE;

24.typedefintRESULT;//返回函数执行的结果状态，下面为具体的返回选项

25.#ifndefTRUE

26.#defineTRUE1

27.#endif

28.#ifndefFALSE

29.#defineFALSE-1

30.#endif

31.#ifndefNULL

32.#defineNULL0

33.#endif

34.#ifndefERR

35.#defineERR-2

36.#endif

37.//////////////////////////////////////////////////////////////////////////////////

38./*===============================================================================

39.***函数名称：Spline3()

40.***功能说明：完成三次样条差值，其中使用追赶法求解M矩阵

41.***入口参数：(pSPLINE)pLine样条结构体指针pLine中的x[],y[],num,begin_k1,end_k1

42.***出口参数：(pSPLINE)pLine样条结构体指针pLine中的函数参数

43.***返回参数：返回程序执行结果的状态TRUEorFALSE

44.================================================================================*/

45.RESULTSpline3(pSPLINEpLine)

46.{

47.floatH[MAXNUM]={0};//小区间的步长

48.floatFi[MAXNUM]={0};//中间量

49.floatU[MAXNUM+1]={0};//中间量

50.floatA[MAXNUM+1]={0};//中间量

51.floatD[MAXNUM+1]={0};//中间量

52.floatM[MAXNUM+1]={0};//M矩阵

53.floatB[MAXNUM+1]={0};//追赶法中间量

54.floatY[MAXNUM+1]={0};//追赶法中间变量

55.inti=0;

56.////////////////////////////////////////计算中间参数

57.if((pLine->point_num<3)||(pLine->point_num>MAXNUM+1))

58.{

59.returnERR;//输入数据点个数太少或太多

60.}

61.for(i=0;i<=pLine->point_num-2;i++)

62.{//求H[i]

63.H[i]=pLine->x[i+1]-pLine->x[i];

64.Fi[i]=(pLine->y[i+1]-pLine->y[i])/H[i];//求F[x(i),x(i+1)]

65.}

66.for(i=1;i<=pLine->point_num-2;i++)

67.{//求U[i]和A[i]和D[i]

68.U[i]=H[i-1]/(H[i-1]+H[i]);

69.A[i]=H[i]/(H[i-1]+H[i]);

70.D[i]=6*(Fi[i]-Fi[i-1])/(H[i-1]+H[i]);

71.}

72.//若边界条件为1号条件，则

73.U[i]=1;

74.A[0]=1;

75.D[0]=6*(Fi[0]-pLine->begin_k1)/H[0];

76.D[i]=6*(pLine->end_k1-Fi[i-1])/H[i-1];

77.//若边界条件为2号条件，则

78.//U[i]=0;

79.//A[0]=0;

80.//D[0]=2*begin_k2;

81.//D[i]=2*end_k2;

82./////////////////////////////////////////追赶法求解M矩阵

83.B[0]=A[0]/2;

84.for(i=1;i<=pLine->point_num-2;i++)

85.{

86.B[i]=A[i]/(2-U[i]*B[i-1]);

87.}

88.Y[0]=D[0]/2;

89.for(i=1;i<=pLine->point_num-1;i++)

90.{

91.Y[i]=(D[i]-U[i]*Y[i-1])/(2-U[i]*B[i-1]);

92.}

93.M[pLine->point_num-1]=Y[pLine->point_num-1];

94.for(i=pLine->point_num-1;i>0;i--)

95.{

96.M[i-1]=Y[i-1]-B[i-1]*M[i];

97.}

98.//////////////////////////////////////////计算方程组最终结果

99.for(i=0;i<=pLine->point_num-2;i++)

100.{

101.pLine->a3[i]=M[i]/(6*H[i]);

102.pLine->a1[i]=(pLine->y[i]-M[i]*H[i]*H[i]/6)/H[i];

103.pLine->b3[i]=M[i+1]/(6*H[i]);

104.pLine->b1[i]=(pLine->y[i+1]-M[i+1]*H[i]*H[i]/6)/H[i];

105.}

106.returnTRUE;

107.}

108.//////////////////////////////////////////////////////////////////////////////////

109.SPLINEline1;

110.pSPLINEpLine1=&line1;

111.//////////////////////////////////////////////////////////////////////////////////

112.main()

113.{

114.line1.x[0]=27.7;

115.line1.x[1]=28;

116.line1.x[2]=29;

117.line1.x[3]=30;

118.line1.y[0]=4.1;

119.line1.y[1]=4.3;

120.line1.y[2]=4.1;

121.line1.y[3]=3.0;

122.line1.point_num=4;

123.line1.begin_k1=3.0;

124.line1.end_k1=-4.0;

125.Spline3(pLine1);

126.return0;

127.}

128.//////////////////////////////////////////////////////////////////////////////////