C语言最小二乘法.doc

函数逼近与曲线拟合，用最小二乘法进行曲线拟合的C或C+编写的完整程序！已知x 0 5 10 15 20 25 30 35 40 45 50 55 y 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 近似解析表达式为y=at+bt2+ct3求a,b,c曲线拟合： #include #include #include #include Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3); void main() int i ,n ,m ; double *x,*y,*a,dt1,dt2,dt3,b; n = 12;/ 12个样点 m = 4; /3次多项式拟合 b = 0; /x的初值为0 /*分别为x,y,a分配存贮空间*/ x = (double *)calloc(n,sizeof(double); if(x = NULL) printf(内存分配失败n); exit (0); y = (double *)calloc(n,sizeof(double); if(y = NULL) printf(内存分配失败n); exit (0); a = (double *)calloc(n,sizeof(double); if(a = NULL) printf(内存分配失败n); exit (0); for(i=1;i=n;i+) xi-1=b+(i-1)*5; /*每隔5取一个点，这样连续取12个点*/ y0=0; y1=1.27; y2=2.16; y3=2.86; y4=3.44; y5=3.87; y6=4.15; y7=4.37; y8=4.51; y9=4.58; y10=4.02; y11=4.64; /*xi-1点对应的y值是拟合已知值*/ Smooth(x,y,a,n,m,&dt1,&dt2,&dt3); /*调用拟合函数*/ for(i=1;i=m;i+) printf(a%d = %.10fn,(i-1),ai-1); printf(拟合多项式与数据点偏差的平方和为：n); printf(%.10en,dt1); printf(拟合多项式与数据点偏差的绝对值之和为：n); printf(%.10en,dt2); printf(拟合多项式与数据点偏差的绝对值最大值为：n); printf(%.10en,dt3); free(x); /*释放存储空间*/ free(y); /*释放存储空间*/ free(a); /*释放存储空间*/ Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3)/(x,y,a,n,m,dt1,dt2,dt3 ) /double *x; /*实型一维数组，输入参数，存放节点的xi值*/ /double *y; /*实型一维数组，输入参数，存放节点的yi值*/ /double *a; /*双精度实型一维数组，长度为m。返回m一1次拟合多项式的m个系数*/ /int n; /*整型变量，输入参数，给定数据点的个数*/ /int m; /*整型变量，输入参数，拟合多项式的项数*/ /double *dt1; /*实型变量，输出参数，拟合多项式与数据点偏差的平方和*/ /double *dt2; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值之和*/ /double *dt3; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值最大值*/ int i ,j ,k ; double *s,*t,*b,z,d1,p,c,d2,g,q,dt; /*分别为s ,t ,b分配存贮空间*/ s = (double *)calloc(n,sizeof(double); if(s = NULL) printf(内存分配失败n); exit (0); t = (double *)calloc(n,sizeof(double); if(t = NULL) printf(内存分配失败n); exit (0); b = (double *)calloc(n,sizeof(double); if(b = NULL) printf(内存分配失败n); exit (0); z = 0; for(i=1;i=n;i+) z=z+xi-1/n; /*z为各个x的平均值*/ b0=1; d1=n; p=0; c=0; for(i=1;i1) t1=1; t0=-p; d2=0; c=0; g=0; for(i=1;i=n;i+) q=xi-1-z-p; d2=d2+q*q; c=yi-1*q+c; g=(xi-1-z)*q*q+g; c=c/d2; p=g/d2; q=d2/d1; d1=d2; a1=c*t1; a0=c*t0+a0; for(j=3;j=4) for(k=j-2;k=2;k-) sk-1=-p*tk-1+tk-2-q*bk-1; s0=-p*t0-q*b0; d2=0; c=0; g=0; for(i=1;i=1;k-) q=q*(xi-1-z)+sk-1; d2=d2+q*q; c=yi-1*q+c; g=(xi-1-z)*q*q+g; c=c/d2; p=g/d2; q=d2/d1; d1=d2; aj-1=c*sj-1; tj-1=sj-1; for(k=j-1;k=1;k-) ak-1=c*sk-1+ak-1; bk-1=tk-1; tk-1=sk-1; *dt1=0; *dt2=0; *dt3=0; for(i=1;i=1;k-) q=q*(xi-1-z)+ak-1; d