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

1、  一、代码：/* *            ArgMin.h *  Author   *  Fri Aug  5 17:18:47 2005 *  Copyright  2005   *  Email   */#include "stdio.h"#include "math.h"

2、;/* Interface */solve the linear equations using "argmin" methodint ArgMin(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);/Solve a series of linear equationsint SolveLinearEqts(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);/Elimination of unknownsint

3、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=&

4、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+)          eqtk*COLUMN+l=0;      for(j=0;j<rowNum;j+) &

5、#160;      eqtk*COLUMN+l+=*(mtrx_tmp+COLUMN*j+k)*(*(mtrx_tmp+COLUMN*j+l);      for(k=0;k<*q-1;k+)      eqtk*COLUMN+*q-1=0;    for(j=0;j<rowNum;j+)      eqtk*COLUMN+*q-1+=*(mtrx_tmp+COLUMN*j

6、+*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+)    &#

7、160;     printf("%.2lf/t",eqtm*COLUMN+n);      printf("/n");    */  if(!SolveLinearEqts(eqt,COLUMN,(*q-1),*q,solution)     return 0;      return 1;/Solve the solution of a series of lin

8、ear equationsint 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

9、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 eleme

10、nts in matrix A/n");      printf("(Tip: You can try setting a different /"theta/" value or ");      printf("checking the data introduced to function/"CalculateCoeff/"!)/n");      r

11、eturn 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; 

12、0;  for(j=i+1;j<rowNum;j+)      tmpSum+=*(solution+j)*(*(inMtrx+i*COLUMN+j);          if(fabs(inMtrxi*COLUMN+i)<0.000001)       printf("/nNeed more Equations to solve elements in matrix A/n"

13、;);      printf("(Tip: You can try setting a different /"theta/" value, /n or ");      printf("checking the data introduced to function/"CalculateCoeff/"!)/n");      return 0; 

14、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

15、,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_tmprowNocolNo down to   a suitable site*/  p=mtrx_tmp+rowNo*COLUMN+colNo; 

16、 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=

17、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+)       &

18、#160;  p+;      if(rowMax<fabs(*p)        rowMax=fabs(*p);                  p+=COLUMN-(colNum-colNo-1);    if(rowMax=0) return 0;   

19、 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<colN

20、um;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*/  t

21、mpDbl1=*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;    &

22、#160; p+;      mp+;        p+=COLUMN-(colNum-colNo);    mp-=colNum-colNo;    return 1; 二、测试/*/     (1)int ArgMin(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);  

23、0;    函数功能：用最小二乘法求解方程数多于未知变量的线性方程组的最适解。       测试代码：          #define Q 5          #define ROWNUM 6          int main()&#

24、160;                     double solutionQ-1,eqtQQ=0;            double aROWNUMQ=3,-2,4,6,-11,4,3,2,9,-2,2,6,8,3,4,2,4,5,3,3,    &#

25、160;                             0,0,8,6,-20,3,4,5,6,0;            ArgMin(&a00,Q,ROWNUM,Q,&soluti

26、on0);            for(i=0;i<Q-1;i+) printf("%.14lf/n",solutioni);            return 0;                

27、测试结果：          eqts:          42.000000000    38.000000000    61.000000000    84.000000000    -27.000000000       &

28、#160;  38.000000000    81.000000000    86.000000000    69.000000000    52.000000000          61.000000000    86.000000000    198.000000000   159.00000

29、0000   -161.000000000          84.000000000    69.000000000    159.000000000   207.000000000   -183.000000000          3.00000000000000  

30、60;       2.00000000000000          -1.00000000000000          -2.00000000000000        说明：a, 待求解的方程组数据：       

31、;         3,-2,4,6,-11,                4,3,2,9,-2,                2,6,8,3,4,     

32、0;          2,4,5,3,3,                0,0,8,6,-20,                3,4,5,6,0     

33、           以上数据中共有四个变量，六个方程；              b, 测试显示的数据中，“eqts”为通过最小二乘法得出的待求解的线性方程组增广矩阵；下面四行为计算求得的解，与实际完全符合。     (2)int SolveLinearEqts(double *inMtrx,int COLUMN,int ro

34、wNum,int colNum,double *solution);       函数功能：用消元法解N×N的线性方程组（即未知数和方程数相同的方程组）。       测试代码：          double b45=3,-2,4,6,-11,4,3,2,9,-2,2,6,8,3,4,2,4,5,3,3;      

35、    double c5;          SolveLinearEqts(&b00,5,4,5,&c0);          for(i=0;i<4;i+) printf("%.14lf/n",ci);           测试结果：

36、0;         3.00000000000000          2.00000000000000          -1.00000000000000          -2.00000000000000   

37、    说明：求得的解与实际完全符合。/*/  (c)int ArgMin(double *inMtrx,int COLUMN,int rowNum,int colNum,double *solution);  功能说明：用最小二乘法求解方程数多于未知变量的线性方程组的最适解。  参数说明：    (1)inMtrx: 线性方程组以增广矩阵的形式用此数组传入方程；    (2)COLUMN: 数组inMtrx的列数为COLUMN；    (3)rowNum: 数组inMtrx中数据的行数