c语言实现,用最小二乘法求解方程数多于未知变量的线性方程组的最适解(即矛盾方程组)

  一、代码:

/***************************************************************************
 *            ArgMin.h
 *  Author 
 *  Fri Aug  5 17:18:47 2005
 *  Copyright  2005 
 *  Email 
 ****************************************************************************/


#include "stdio.h"
#include "math.h"

/* Interface */
//solve the linear equations using "argmin" method
int ArgMin(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);
//Solve a series of linear equations
int SolveLinearEqts(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);
//Elimination of unknowns
int ReduceUnknowns(double *mtrx_tmp,int COLUMN,int rowNo,int colNo,int rowNum,int colNum);


/* Begin Coding */
//ArgMin using "solvelineareqts". 

int ArgMin(double *mtrx_tmp,int COLUMN,int rowNum,int colNum,double *solution)
{
  int k,l,j;  
  double *eqt;
  int *q;
  q=&colNum;
  eqt=(double *)malloc(colNum*colNum*sizeof(double));
  //printf("%d",*q);
 
  for(k=0;k<*q-1;k++)
    for(l=0;l<*q-1;l++)
    {
      eqt[k*COLUMN+l]=0;
      for(j=0;j<rowNum;j++)
        eqt[k*COLUMN+l]+=*(mtrx_tmp+COLUMN*j+k)*(*(mtrx_tmp+COLUMN*j+l));
    }
  for(k=0;k<*q-1;k++)
  {
    eqt[k*COLUMN+*q-1]=0;
    for(j=0;j<rowNum;j++)
      eqt[k*COLUMN+*q-1]+=*(mtrx_tmp+COLUMN*j+*q-1)*(*(mtrx_tmp+COLUMN*j+k));
  }
 
    /* show the contents of eqts */
    /*int m,n;
    printf("eqt:/n");
    for(m=0;m<*q-1;m++)
    {
      for(n=0;n<*q;n++)  
        printf("%.2lf/t",eqt[m*COLUMN+n]);
      printf("/n");
    }*/

  if(!SolveLinearEqts(eqt,COLUMN,(*q-1),*q,solution)) {
    return 0;
  }
 
  return 1;
}

//Solve the solution of a series of linear equations
int SolveLinearEqts(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution)
{
  int i,j;
  double tmpSum;
  if(rowNum!=(colNum-1)) {
    printf("Can't solve the equations because equation number ");
    printf("is not the same as unknow parameters!/n");
    return 0;
  }
  //reduce unknown parameters
  for(i=0;i<colNum-2;i++)
  {  
    if(!ReduceUnknowns(inMtrx,COLUMN,i,i,rowNum,colNum)) {
      printf("/nNeed more Equations to solve elements in matrix A/n");
      printf("(Tip: You can try setting a different /"theta/" value or ");
      printf("checking the data introduced to function/"CalculateCoeff/"!)/n");
      return 0;
    }
  }
  //Calculate the equation at the bottom to acquire value of the first
  //variable, then Substitute the solved variables to the equations in
  //order to solve more variables:
  for(i=rowNum-1;i>=0;i--)
  {
    tmpSum=0;
    for(j=i+1;j<rowNum;j++)
      tmpSum+=*(solution+j)*(*(inMtrx+i*COLUMN+j));
      
    if(fabs(inMtrx[i*COLUMN+i])<0.000001) {
      printf("/nNeed more Equations to solve elements in matrix A/n");
      printf("(Tip: You can try setting a different /"theta/" value, /n or ");
      printf("checking the data introduced to function/"CalculateCoeff/"!)/n");
      return 0;      
    }
    *(solution+i)=(*(inMtrx+i*COLUMN+colNum-1)-tmpSum)/(*(inMtrx+i*COLUMN+i));
  }
  return 1;
}

/*Column number of Array is COLUMN_A, but actually the number of figures
in each row is "colNum"! */
int ReduceUnknowns(double *mtrx_tmp,int COLUMN,int rowNo,int colNo,int rowNum,int colNum)
{
  int i,j,tmpInt;
  double tmpDbl1,tmpDbl2;
  double rowMax,colMax;
  double *p,*mp;

  /*Adjust the rowNo whose value is zero in matrix_tmp[rowNo][colNo] down to
  a suitable site*/
  p=mtrx_tmp+rowNo*COLUMN+colNo;
  mp=mtrx_tmp+rowNo*COLUMN+colNo;
  //select main variable:
  /*rowMax=fabs(*p);
  for(j=colNo+1;j<colNum;j++)
  {
    p++;
    if(rowMax<fabs(*p)){
      rowMax=fabs(*p);
    }
  }    
  if(rowMax==0) return 0;
  p+=COLUMN-(colNum-colNo-1);
  colMax=fabs(*(mtrx_tmp+rowNo*COLUMN+colNo))/rowMax;
  tmpInt=rowNo;  
  for(i=rowNo+1;i<rowNum;i++)
  {  
    rowMax=fabs(*p);
    for(j=colNo+1;j<colNum;j++)
    {
      p++;
      if(rowMax<fabs(*p)){
        rowMax=fabs(*p);
      }
    }    
    p+=COLUMN-(colNum-colNo-1);
    if(rowMax==0) return 0;
    tmpDbl1=fabs(*(mtrx_tmp+i*COLUMN+colNo))/rowMax;
    if(colMax<tmpDbl1) {
      colMax=tmpDbl1;
      tmpInt=i;
    }
  }
  //change the whole row of "tmpInt" to "rowNo"
  for(j=colNo;j<colNum;j++)
  {
    tmpDbl1=*(mtrx_tmp+rowNo*COLUMN+j);
    *(mtrx_tmp+rowNo*COLUMN+j)=*(mtrx_tmp+tmpInt*COLUMN+j);
    *(mtrx_tmp+tmpInt*COLUMN+j)=tmpDbl1;
  }*/
    
  /*Transform the matrix_tmp using linear calculation methods*/
  tmpDbl1=*p;
  p+=COLUMN;
  for(i=rowNo+1;i<rowNum;i++)
  {
    if(fabs(*p)<0.000001) continue;
    tmpDbl2=*p;
    for(j=colNo;j<colNum;j++)
    {
      *p=*p-*mp*tmpDbl2/tmpDbl1;
      p++;
      mp++;
    }
    p+=COLUMN-(colNum-colNo);
    mp-=colNum-colNo;
  }
  return 1;
}

 

二、测试

/******************************************************************************/
     (1)int ArgMin(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);
       函数功能:用最小二乘法求解方程数多于未知变量的线性方程组的最适解。
       测试代码:
          #define Q 5
          #define ROWNUM 6
          int main()
          {
            double solution[Q-1],eqt[Q][Q]={0};
            double a[ROWNUM][Q]={{3,-2,4,6,-11},{4,3,2,9,-2},{2,6,8,3,4},{2,4,5,3,3},
                                  {0,0,8,6,-20},{3,4,5,6,0}};
            ArgMin(&a[0][0],Q,ROWNUM,Q,&solution[0]);
            for(i=0;i<Q-1;i++) printf("%.14lf/n",solution[i]);
            return 0;
          }
       测试结果:
          eqts:
          42.000000000    38.000000000    61.000000000    84.000000000    -27.000000000
          38.000000000    81.000000000    86.000000000    69.000000000    52.000000000
          61.000000000    86.000000000    198.000000000   159.000000000   -161.000000000
          84.000000000    69.000000000    159.000000000   207.000000000   -183.000000000
          3.00000000000000
          2.00000000000000
          -1.00000000000000
          -2.00000000000000

        说明:a, 待求解的方程组数据:
                {3,-2,4,6,-11},
                {4,3,2,9,-2},
                {2,6,8,3,4},
                {2,4,5,3,3},
                {0,0,8,6,-20},
                {3,4,5,6,0}
                以上数据中共有四个变量,六个方程;
              b, 测试显示的数据中,“eqts”为通过最小二乘法得出的待求解的线性方程组增广矩阵;下面四行为计算求得的解,与实际完全符合。

     (2)int SolveLinearEqts(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);
       函数功能:用消元法解N×N的线性方程组(即未知数和方程数相同的方程组)。
       测试代码:
          double b[4][5]={{3,-2,4,6,-11},{4,3,2,9,-2},{2,6,8,3,4},{2,4,5,3,3}};
          double c[5];
          SolveLinearEqts(&b[0][0],5,4,5,&c[0]);
          for(i=0;i<4;i++) printf("%.14lf/n",c[i]);    
       测试结果:
          3.00000000000000
          2.00000000000000
          -1.00000000000000
          -2.00000000000000

       说明:求得的解与实际完全符合。


/******************************************************************************/
  (c)int ArgMin(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);
  功能说明:用最小二乘法求解方程数多于未知变量的线性方程组的最适解。
  参数说明:
    (1)inMtrx: 线性方程组以增广矩阵的形式用此数组传入方程;
    (2)COLUMN: 数组inMtrx的列数为COLUMN;
    (3)rowNum: 数组inMtrx中数据的行数;
    (4)colNum: 数组inMtrx中数据的列数;
    (5)solution: 求得方程的解从此参数传出;
    (6)返回值:0,表示方程组无唯一解;1,表示函数正常结束。

  (d)int SolveLinearEqts(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);
  功能说明:解N×N的线性方程组(即未知数和方程数相同的方程组)。
  参数说明:
    (1)inMtrx: 线性方程组以增广矩阵的形式用此数组传入方程;
    (2)COLUMN: 数组inMtrx的列数为COLUMN;
    (3)rowNum: 数组inMtrx中数据的行数;
    (4)colNum: 数组inMtrx中数据的列数;
    (5)solution: 求得方程的解从此参数传出;
    (6)返回值:0,表示方程组无唯一解;1,表示函数正常结束。

  (e)int ReduceUnknowns(double *mtrx_tmp,int COLUMN,int rowNo,int colNo,int rowNum,int colNum);
  功能说明:将储存在mtrx_tmp中、表示待求解线性方程组的增广矩阵变换成阶梯形矩阵。每调用函数一次,通过变换矩阵会将其中元素(rowNo,colNo)对应的那一列第rowNo行以下(不包括第rowNo行)的元素全变为0。
  参数说明:
    (1)mtrx_tmp: 储存待求解线性方程组增广矩阵的数组以此参数传入其首地址;
    (2)COLUMN: 储存数据的数组mtrx_tmp列数为COLUMN;
    (3)rowNo: 元素(rowNo,colNo)的行号;
    (4)colNo: 元素(rowNo,colNo)的列号;
    (5)rowNum:  mtrx_tmp指向的数组中数据的行数;
    (6)colNum:  mtrx_tmp指向的数组中数据的列数;
    (7)返回值:0,表示元素(rowNo,colNo)对应的那一列第rowNo行及以下(包括第rowNo行)的元素已经都是0,此时线性方程组没有唯一解;1,表示函数正常结束。

  • 2
    点赞
  • 17
    收藏
    觉得还不错? 一键收藏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值