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

4 篇文章 0 订阅

一、代码：

/***************************************************************************
*            ArgMin.h
*  Author
*  Fri Aug  5 17:18:47 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
评论
03-09 2万+
05-29 3万+
12-08 1万+
08-13
06-13
10-01

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

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

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