最小二乘法程序.docx
- 文档编号:26715260
- 上传时间:2023-06-22
- 格式:DOCX
- 页数:17
- 大小:20.45KB
最小二乘法程序.docx
《最小二乘法程序.docx》由会员分享,可在线阅读,更多相关《最小二乘法程序.docx(17页珍藏版)》请在冰豆网上搜索。
最小二乘法程序
最小二乘法程序
include
#ineludevconio.h〉
include
#inelude
#defineN5//N个点
#defineT3//T次拟合
#defineW1〃权函数
#definePRECISION0.00001
floatpow_refloata,intn)
{
inti;
if(n==0)
retum
(1);
floatres=a;
for(i=1;i { res*=a; } retum(res); } voidmutiple(floata[][N],floatb|][T+1]ffloatcQ[T+1]) { floatres=0; inti,j,k; for(i=0;i for(j=0;j { res=0; for(k=0;k { res+=a[i][k]*b[k][j]; c[i][j]=res; } } } voidmatrix_trans(floata[][T+1],floatb[][N]) { inti,j; for(i=0;i } } } voidinit(floatx_y[][2],intn) { inti; printf(“请输入%(1个己知点: \n”,N); for(i=0;i { printf("(x%dy%d): “,i,i); scanf(H%f%f”,&x_y[i][0],&x_y[i][1]); } } voidget_A(floatmatrix_A[][T+1],floatx_y[][2],intn) { inti,j; for(i=0;i { forO=0;j { matrix_A[i][j]=W*pow_n(x_y[i][0],j); } } } voidprint_array(floatarray[][T+1],intn) { inti,j; for(i=0;i { forO=0;j { printf("%-g",array[i]O]); } printfCAn"); } } voidconvert(floatargu[][T+2],intn)inti,j,k,p,t;floatrate,temp; for(i=1;i { if(argu[i-1][i-1]==0) { for(p=i;p { if(argu[p][i-1]! =0) break; } if(p==n) { printfC方程组无解! \n“); exit(O); } for(t=0;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 { arguU][k]-=argu[i-1][k]*rate; if(fabs(argu[j][k])v二PRECISION)argu[j][k]=O; } } } } voidcompute(floatargu[][T+2],intn,floatrootQ) { inti,j; floattemp; 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]; } } voidget_y(floattrans_A[][N],floatx_y[][2],floaty[],intn) { inti,j; floattemp; for(i=0;i { temp=0; for(j=0;j { temp+=trans_A[i][jrx_y[j]⑴; } y[i]=temp; } } voidcons_formula(floatcoef_A[][T+1],floaty[],floatcoef_form[][T+2]) { inti,j; for(i=0;i { for(j=0;j { if(j==T+1) coef_form[i][j]=y[i]; else coef_form[i][j]=coef_A[i][j]; } } } voidprint_root(floata[],intn) { inti,j; printf("%d个点的%(1次拟合的多项式系数为: \nu,N,T); for(i=0;i { printf("a[%d]=%g,,',i+1,a[i]); } printf("\n"); printf('拟合曲线方程为: \ny(x)=%g",a[O]); for(i=1;i { printf("+%g“,a[i]); for(j=0;j { printf("*X"); } } pnntf("\n"); } voidprocess() { float x_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef」ormu[T+1 ][T+2],y[T+1],a[T+1]; init(x_y,N); get_A(matrix_A,x_y,N); printfC1矩阵A为: \nH); print_array(matrix_A,N); matrix_trans(matrix_A,t「ans_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); } voidmain() { process(); } ]]> 23: 57 vranknum>user1 vContent> [CDATA[ 你可以改一下 不从终端输入,直接在程序中给出参数 请输入5个己知点: (xOy0): -2-0.1 (x1y1): -10.1 (x2y2): 00.4 (x3y3): 10.9 (x4y4): 21.6 矩阵A为: 1 ・2 4 -8 1 1 -1 1 0 0 0 1 1 1 1 1 2 4 8 法矩阵为: 5 0 10 0 0 10 0 34 10 0 34 0 0 34 0 13 5个点的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 ]]> 26: 11 vranknum>user1v/「anknum> vContent> [CDATA[include voidmain() { intnum,i; floatx,y,l,m,n,p,a,b; i=1; l=0.0; m=0.0; n=0.0; p=0.0; printf("请输入你想计算的x,y的个数: “); scanf("%d",&num); if(num>=1) { while(i<=num); {printf("i#输入x的值”); scant(“%lf”,&x); printf("请输入y的值”); scant(“%lf”,&y); l+=x; m+=y; n+=x*y; p+=x*x; i++; } a=(numF-l*m)/(num*p-ri); b=(p*m-n*l)/(num*p-l*l); printf(”最小二乘法所算得的斜率和截距分别为%彳和%f\n”ab); } elseprintfCmun"输入有误! ); ftinclude #include llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll ////////////////////// 〃矩阵结构体 structMatrix { intm,n;//m为行数,n为列数 double林pm;〃指向矩阵二维数组的指针 }; 〃初始化矩阵mt,并置矩阵的行为m,列为n voidInitMatrix(structMatrix*mt,intm,intn) { irrti; 〃指定矩阵的行和列 mt->n=n; 〃为矩阵分配内存 mt->pm=newdouble*[m]; for(i=0;i { mt->pm[i]=newdouble[n]; } } 〃用0初始化矩阵mt,并置矩阵的行为m,列为n voidInitMatrixO(structMatrix*mt,intm,intn) InitMatrix(mt,m,n);for(i=0;i } 〃销毁矩阵mt voidDestroyMatrix(structMatrix*mt) { inti; 〃释放矩阵内存 for(i=0;i { delete[]mt->pm[i]; } delete[]mt->pm; 〃矩阵相乘mt3=4ntl*mt2 //成功返回1,失败返回0 intMatrixMul(Matrix*mtl,Matrix*mt2,Matrix*mt3) { inti,j,k; doubles; 〃检査行列号是否匹配 if(mtl->n! =mt2->m||mtl->m! Ft3->m||mt2->n! =mt3->ii)return0; for(i=0;i for(j=0;j { s=0.0; for(k=0;k } return1; } 〃矩阵Wimt2=T(mtl) //成功返回1,失败返回0 intMatrixTran(Matrix*mtl,Matrix*mt2) { inti,j; 〃检査行列号是否匹配 if(mtl->m! =mt2->n||mtl->n! ^nt2->m)return0; for(i=0;i for(j=0;j return1; } 〃控制台显示矩阵mt voidConsoleShowMatrix(Matrix*mt) { inti,j; for(i=0;i { printf("\n"); for(j=0;j } printf("\n"); }llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll////////////////////// //Guass列主元素消元法求解方程Ax=b,a=(A,b) intGuassl(double**a,double*x,intn) inti,j,k,numl,*h,t; double*1,max; l=newdouble[n]; h=newint[n]; for(i=0;i for(i=l;i { max=fabs(a[h[i-1]][iT]); numl=i-l; 〃列元的最大值for(j=i;j { if(fabs(a[h[j]][i-l])>max) {numl=h[j]; max=fabs(a[h[j]][i-1]); } } if(max<0.00000000001)return0; 〃交换行号 if(numl>i-l) { t=h[i]; h[i]=h[numl]; h[numl]=t; } for(j=i;j for(k=i;k a[h[j]][k]=a[h[j]][k]-l[j]*a[h[i-l]][k]; } for(i=n-l;i>=0;i—) { x[i]=a[h[i]][n]; for(j=i+l;j } 〃清除临时数组内存 delete[]1; delete[]h; return1; Illlllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll /////////////////////// //最小二乘法求解矩阵Ax=b intMinMul(MatrixA,Matrixb,double*x) { inti,j; if(b.n! =l)return0; if(A«m! =b・m)return0; MatrixTranA;//定义A的转置 InitMatrixO(ATranA,A.n,A.m); MatrixTran(&A,&TranA); MatrixTranA_A;〃定义A的转置与A的乘积矩阵 InitMatrixO(&TranA_A,A.n,A.n); MatrixMul(ATranA,&A,&TranA_A);//A的转置与A的乘积 MatrixTranA_b;〃定义A的转置与b的乘积矩阵 InitMatrixO(&TranA_b,A.n,1); MatrixMul(&TranA,&b,&TranA_b);//A的转置与b的乘积 DestroyMatrix(&TranA);〃释放A的转置的内存 MatrixTranA_A_b;〃定义增广矩阵 InitMatrixO(&TranA_A_b,TranA_A.m,TranA_A.m+1); 〃增广矩阵赋值 for(i=0;i { for(j=0;j TranA_A_b.pm[i][TranA_A_b.m]=TranA_b.pm[i][0]; } DestroyMatrix(&TranA_A); DestroyMatrix(&TranA_b); //Guass列主消元法求麻 Guassl(TranApm,x,TranA_A_b.m); DestroyMatrix(&TranA_A_b); return1; } intMinMul(double**A,double*b,intm,intn,double*x) { intr,i; MatrixAl,bl; Al.pm=newdouble*[m]; Al.m=m; Al.n=n; InitMatrix(&bl,m,1); for(i=0;i bl.pm[i][O]=b[i];}r=MinMul(Al,bl,x);delete[]A1.pm;DestroyMatrix(&bl);returnr; }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 最小二乘法 程序