最小二乘法程序.doc

最小二乘法程序
#include
#include
#include
#include
#define N 5//N个点
#define T 3 //T次拟合
#define W 1//权函数
#define PRECISION 0.00001
float pow_n(float a,int n)
{
int i;
if(n==0)
return(1);
float res=a;
for(i=1;iid mutiple(float a[][N],float b[][T+1],float c[][T+1])
{
float res=0;
int i,j,k;
for(i=0;i)
{
b[j][i]=a[i][j];
}
}
}
void init(float x_y[][2],int n)
{
int i;
printf("请输入%d个已知点：\n",N);
for(i=0;i1;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;ii++)
{
for(j=i;jk++)
{
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)
{
for(j=0;j
for(i=1;irix_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();
}
]]>

2007-4-19 19:23:57

一级(初级)
user1
100
40389872
5478010
1526752
jiangxc2004
0

2 4 8
法矩阵为：
5 0 10 0
0 10 0 34
10 0 34 0
0 34 0 130
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
]]>

2007-4-19 19:26:11

一级(初级)
user1
100
40390406
5478010
1526752
jiangxc2004
0

void main ()
{
int num,i;
float x,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("请输入x的值");
scanf ("%lf",&x);
printf("请输入y的值");
scanf ("%lf",&y);
l+=x;
m+=y;
n+=x*y;
p+=x*x;
i++;
}
a=(num*n-l*m)/(num*p-l*l);
b=(p*m-n*l)/(num*p-l*l);
printf("最小二乘法所算得的斜率和截距分别为%f和%f\n",a,b);
}
else printf("mun"输入有误！);
}

#include
#include
//////////////////////////////////////////////////////////////////////////////////////////
//矩阵结构体
struct Matrix
{
 int m,n;//m为行数,n为列数
 double **pm;//指向矩阵二维数组的指针
};
//初始化矩阵mt，并置矩阵的行为m，列为n
void InitMatrix(struct Matrix *mt,int m,int n)
{
 int i;
 //指定矩阵的行和列
 mt->m=m;
 mt->n=n;
 //为矩阵分配内存
 mt->pm=new double *[m];
 for(i=0;ipm[i]=new double[n];
 }
}
//用0初始化矩阵mt，并置矩阵的行为m，列为n
void InitMatrix0(struct Matrix *mt,int m,int n)
{
 int i,j;
    InitMatrix(mt,m,n);
 for(i=0;i)
  for(j=0;jpm[i][j]=0.0;
}
//销毁矩阵mt
void DestroyMatrix(struct Matrix *mt)
{
 int i;
 //释放矩阵内存
 for(i=0;im;i++)
 {
  delete []mt->pm[i];
 }
 delete []mt->pm;
}
//矩阵相乘mt3=mt1*mt2
//成功返回1，失败返回0
int MatrixMul(Matrix *mt1,Matrix *mt2,Matrix *mt3)
{
 int i,j,k;
 double s;
 //检查行列号是否匹配
 if(mt1->n!=mt2->m||mt1->m!=mt3->m||mt2->n!=mt3->n) return 0;
 for(i=0;im;i++)
  for(j=0;jn;j++)
  {
            s=0.0;
   for(k=0;kn;k++) s=s+mt1->pm[i][k]*mt2->pm[k][j];
   mt3->pm[i][j]=s;
  }
  return 1;
}
//矩阵转置mt2=T(mt1)
//成功返回1，失败返回0
int MatrixTran(Matrix *mt1,Matrix *mt2)
{
 int i,j;
 //检查行列号是否匹配
 if(mt1->m!=mt2->n||mt1->n!=mt2->m) return 0;
 for(i=0;im;i++)
  for(j=0;jn;j++) mt2->pm[j][i]=mt1->pm[i][j];
  return 1;
}
//控制台显示矩阵mt
void ConsoleShowMatrix(Matrix *mt)
{
 int i,j;
 for(i=0;im;i++)
 {
  printf"\n");
  for(j=0;jn;j++) printf("%2.4f ",mt->pm[i][j]);
 }
 printf("\n");
}
//////////////////////////////////////////////////////////////////////////////////////////
//Guass列主元素消元法求解方程Ax=b,a=(A,b)
int Guassl(double **a,double *x,int n)
{
 int i,j,k,numl,*h,t;
    double *l,max;
 l=new double[n];
 h=new int[n];
 for(i=0;imax)
   {
    numl=h[j];
    max=fabs(a[h[j]][i-1]);
   }
  }
  if(max<0.00000000001) return 0;
  //交换行号
  if(numl>i-1)
  {
   t=h[i];
   h[i]=h[numl];
   h[numl]=t;
  }
  for(j=i;j=0;i--)
 {
        x[i]=a[h[i]][n];
        for(j=i+1;j<n;j++) x[i]=x[i]-a[h[i]][j]*x[j];
        x[i]=x[i]/a[h[i]][i];
 }

 //清除临时数组内存
 delete []l;
 delete []h;
 return 1;
}
///////////////////////////////////////////////////////////////////////////////////////////
//最小二乘法求解矩阵Ax=b
int MinMul(Matrix A,Matrix b,double *x)
{
 int i,j;
 if(b.n!=1) return 0;
 if(A.m!=b.m) return 0;
 Matrix TranA;//定义A的转置
 InitMatrix0(&TranA,A.n,A.m);
    MatrixTran(&A,&TranA);
 Matrix TranA_A;//定义A的转置与A的乘积矩阵
 InitMatrix0(&TranA_A,A.n,A.n);
    MatrixMul(&TranA,&A,&TranA_A);//A的转置与A的乘积
 Matrix TranA_b;//定义A的转置与b的乘积矩阵
 InitMatrix0(&TranA_b,A.n,1);
 MatrixMul(&TranA,&b,&TranA_b);//A的转置与b的乘积
 DestroyMatrix(&TranA);//释放A的转置的内存
 Matrix TranA_A_b;//定义增广矩阵
 InitMatrix0(&TranA_A_b,TranA_A.m,TranA_A.m+1);
 //增广矩阵赋值
 for(i=0;iA_b.m]=TranA_b.pm[i][0];
 }
 DestroyMatrix(&TranA_A);
 DestroyMatrix(&TranA_b);
 //Guass列主消元法求解
 Guassl(TranA_A_b.pm,x,TranA_A_b.m);
 DestroyMatrix(&TranA_A_b);
 return 1;
}
int MinMul(double **A,double *b,int m,int n,double *x)
{
 int r,i;
 Matrix Al,bl;
 Al.pm=new double *[m];
 Al.m=m;
 Al.n=n;
 InitMatrix(&bl,m,1);
 for(i=0;i