正文

用c++写成的最小二乘法的源代码2009-03-30 20:55:00

【评论】 【打印】 【字体: 】 本文链接:http://blog.pfan.cn/weiyutian/41872.html

分享到:

  #include <stdio.h>#include <conio.h>#include <math.h>#include <process.h>#define N 5//N个点#define T 3 //T次拟合#define W 1//权函数#define PRECISION 0.00001float pow_n(float a,int n){int i;if(n==0)return(1);float res=a;for(i=1;i<n;i++){res*=a;}return(res);}void mutiple(float a[][N],float b[][T+1],float c[][T+1]){float res=0;int i,j,k;for(i=0;i<T+1;i++)for(j=0;j<T+1;j++){res=0;for(k=0;k<N;k++){res+=a[i][k]*b[k][j];c[i][j]=res;}}}void matrix_trans(float a[][T+1],float b[][N]){int i,j;for(i=0;i<N;i++){for(j=0;j<T+1;j++){b[j][i]=a[i][j];}}}void init(float x_y[][2],int n){int i;printf("请输入%d个已知点:\n",N);for(i=0;i<n;i++){printf("(x%d y%d):",i,i);scanf("%f %f",&x_y[i][0],&x_y[i][1]);}}void get_A(float matrix_A[][T+1],float x_y[][2],int n){int i,j;for(i=0;i<N;i++){for(j=0;j<T+1;j++){matrix_A[i][j]=W*pow_n(x_y[i][0],j);}}}void print_array(float array[][T+1],int n){int i,j;for(i=0;i<n;i++){for(j=0;j<T+1;j++){printf("%-g",array[i][j]);}printf("\n");}}void convert(float argu[][T+2],int n){int i,j,k,p,t;float rate,temp;for(i=1;i<n;i++){for(j=i;j<n;j++){if(argu[i-1][i-1]==0){for(p=i;p<n;p++){if(argu[p][i-1]!=0)break;}if(p==n){printf("方程组无解!\n");exit(0);}for(t=0;t<n+1;t++){temp=argu[i-1][t];argu[i-1][t]=argu[p][t];argu[p][t]=temp;}}rate=argu[j][i-1]/argu[i-1][i-1];for(k=i-1;k<n+1;k++){argu[j][k]-=argu[i-1][k]*rate;if(fabs(argu[j][k])<=PRECISION)argu[j][k]=0;}}}}void compute(float argu[][T+2],int n,float root[]){int i,j;float temp;for(i=n-1;i>=0;i--){temp=argu[i][n];for(j=n-1;j>i;j--){temp-=argu[i][j]*root[j];}root[i]=temp/argu[i][i];}}void get_y(float trans_A[][N],float x_y[][2],float y[],int n){int i,j;float temp;for(i=0;i<n;i++){temp=0;for(j=0;j<N;j++){temp+=trans_A[i][j]*x_y[j][1];}y[i]=temp;}}void cons_formula(float coef_A[][T+1],float y[],float coef_form[][T+2]){int i,j;for(i=0;i<T+1;i++){for(j=0;j<T+2;j++){if(j==T+1)coef_form[i][j]=y[i];elsecoef_form[i][j]=coef_A[i][j];}}}void print_root(float a[],int n){int i,j;printf("%d个点的%d次拟合的多项式系数为:\n",N,T);for(i=0;i<n;i++){printf("a[%d]=%g,",i+1,a[i]);}printf("\n");printf("拟合曲线方程为:\ny(x)=%g",a[0]);for(i=1;i<n;i++){printf(" + %g",a[i]);for(j=0;j<i;j++){printf("*X");}}printf("\n");}void process(){float x_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef_formu[T+1][T+2],y[T+1],a[T+1];init(x_y,N);get_A(matrix_A,x_y,N);printf("矩阵A为:\n");print_array(matrix_A,N);matrix_trans(matrix_A,trans_A);mutiple(trans_A,matrix_A,coef_A);printf("法矩阵为:\n");print_array(coef_A,T+1);get_y(trans_A,x_y,y,T+1);cons_formula(coef_A,y,coef_formu);convert(coef_formu,T+1);compute(coef_formu,T+1,a);print_root(a,T+1);}void main(){process();}]]></Content><PostDateTime>2007-4-19 19:23:57</PostDateTime></Reply><Reply><PostUserNickName></PostUserNickName><rank>一级(初级)</rank><ranknum>user1</ranknum><credit>100</credit><ReplyID>40389872</ReplyID><TopicID>5478010</TopicID><PostUserId>1526752</PostUserId><PostUserName>jiangxc2004</PostUserName><Point>0</Point><Content><![CDATA[你可以改一下不从终端输入,直接在程序中给出参数 请输入5个已知点:(x0 y0):-2 -0.1(x1 y1):-1 0.1(x2 y2):0 0.4(x3 y3):1 0.9(x4 y4):2 1.6矩阵A为:1         -2        4         -81         -1        1         -11         0         0         01         1         1         11         2         4         8法矩阵为:5         0         10        00         10        0         3410        0         34        00         34        0         1305个点的3次拟合的多项式系数为:a[1]=0.408571,    a[2]=0.391667,    a[3]=0.0857143, a[4]=0.00833333,拟合曲线方程为:y(x)=0.408571 + 0.391667*X + 0.0857143*X*X + 0.00833333*X*X*X ]]></Content><PostDateTime>2007-4-19 19:26:11</PostDateTime></Reply><Reply><PostUserNickName></PostUserNickName><rank>一级(初级)</rank><ranknum>user1</ranknum><credit>100</credit><ReplyID>40390406</ReplyID><TopicID>5478010</TopicID><PostUserId>1526752</PostUserId><PostUserName>jiangxc2004</PostUserName><Point>0</Point><Content><![CDATA[这样就可以直接调用process()函数了! 二次拟合的话就把宏 T 成2;拟合点的数目 N 也可以修改!也可以去到注释的部分进行返回值的调用!

阅读(759) | 评论(0)


版权声明:编程爱好者网站为此博客服务提供商,如本文牵涉到版权问题,编程爱好者网站不承担相关责任,如有版权问题请直接与本文作者联系解决。谢谢!

评论

暂无评论
您需要登录后才能评论,请 登录 或者 注册