(C语言)最小二乘法的曲线拟合

/*最小二乘法的曲线拟合*/#includestdio.h#includemath.h#includestdlib.h#definemax100voidmain(){inti,j,k,m,N,mi;floatmx,temp;floatX[max][max],Y[max],x[max],y[max],a[max];FILE*fp1;if((fp1=fopen(in1.txt,r))==NULL)/*输入拟合曲线的次数m以及已知的数据组数N*/{printf(Can'topenthisfile!\n);exit(0);}for(i=0;i2;i++)fscanf(fp1,%d%d,&m,&N);fclose(fp1);FILE*fp2;if((fp2=fopen(in2.txt,r))==NULL)/*输入已知的N组数据坐标*/{printf(Can'topenthisfile!\n);exit(0);}for(i=0;iN;i++)fscanf(fp2,%f%f,&x[i],&y[i]);fclose(fp2);for(i=0;i=m;i++)/*由给定的点得系数矩阵*/for(j=i;j=m;j++){temp=0;for(k=0;kN;k++)temp=temp+pow(x[k],(i+j));X[i][j]=temp;X[j][i]=X[i][j];}for(i=0;i=m;i++)/*由给定的点得右端矩阵列向量*/{temp=0;for(k=0;kN;k++)temp=temp+y[k]*pow(x[k],i);Y[i]=temp;}for(j=0;jm;j++)/*列主元素消去法解该线性方程组，得拟合曲线多项式的系数*/{for(i=j+1,mi=j,mx=fabs(X[j][j]);i=m;i++)if(fabs(X[i][j])mx){mi=i;mx=fabs(X[i][j]);}if(jmi){temp=Y[j];Y[j]=Y[mi];Y[mi]=temp;for(k=j;k=m;k++){temp=X[j][k];X[j][k]=X[mi][k];X[mi][k]=temp;}}for(i=j+1;i=m;i++){temp=-X[i][j]/X[j][j];Y[i]+=Y[j]*temp;for(k=j;k=m;k++)X[i][k]+=X[j][k]*temp;}}a[m]=Y[m]/X[m][m];for(i=m-1;i=0;i--){a[i]=Y[i];for(j=i+1;j=m;j++)a[i]-=X[i][j]*a[j];a[i]/=X[i][i];}FILE*fp3;/*输出拟合所得的曲线方程*/fp3=fopen(out.txt,w);fprintf(fp3,P(x)=(%f),a[0]);for(i=1;i=m;i++)fprintf(fp3,+(%f)*x^%d,a[i],i);fclose(fp3);}