当前位置:文档之家› (C语言)最小二乘法的曲线拟合

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

/*最小二乘法的曲线拟合*/
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define max100
void main()
{
int i,j,k,m,N,mi;
float mx,temp;
float X[max][max],Y[max],x[max],y[max],a[max];
FILE*fp1;
if((fp1=fopen("in1.txt","r"))==NULL)/*输入拟合曲线的次数m以及已知的数据组数N*/
{
printf("Can't open this file!\n");
exit(0);
}
for(i=0;i<2;i++)
fscanf(fp1,"%d%d",&m,&N);
fclose(fp1);
FILE*fp2;
if((fp2=fopen("in2.txt","r"))==NULL)/*输入已知的N组数据坐标*/
{
printf("Can't open this file!\n");
exit(0);
}
for(i=0;i<N;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;k<N;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;k<N;k++)
temp=temp+y[k]*pow(x[k],i);
Y[i]=temp;
}
for(j=0;j<m;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(j<mi)
{
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);
}。

相关主题