計算機程序?qū)嵺`報告

計算機程序設計

實踐報告

學院：資源加工與生物工程學院

班級：生物工程類1001班

指導教師：楊英杰

學號：_______0306100103_______

姓名：_________王蕾________

完成日期：2011年6月24日

1某選礦車間某班取樣化驗數(shù)據(jù)如表1,利用礦量平衡和金屬平衡可以列

出含四個方程的線性方程組，請利用行線性變換編程求出各精礦及尾礦的

理論產(chǎn)率，并填入表1相應位置。

表1各礦流的化驗品位

礦產(chǎn)化驗品位％

率

流%PbZnS

原礦100.003.0811.9214.58

鉛精礦X157.164.0821.19

鋅精礦X20.7956.5621.71

硫精礦

X32.062.9040.24

尾礦X40.250.734.93

解題思路：首先，本題主要的就是設計一個求解線性方程組的程序，而求解

線性方程組主要是通過矩陣的線性行變換來求解。由系數(shù)矩陣與常數(shù)項向量組成

增廣矩陣，通過線性行變換，當系數(shù)矩陣變?yōu)閱挝痪仃嚂r，那么變化后對應的常

數(shù)量就是對應的方程組地解。所以，解線性方程組的解法布驟歸結(jié)起來就是矩陣

的線性行變換。程序運行成功后，將其添加到界面程序中。

利用礦量和三種組分的進出平衡可列出四個線性方程組：

X1+x=100

x2+x3+4

57.16xi+0.79x2+2.06X3+0.25x4=3.08

4.08xi+56.56X2+2.90x3+0.73x4=11.92

21.19xt+31.17X2+40.24x3+4.93x4=14.58

獨立方程個數(shù)多于未知數(shù)個數(shù)

X=X]D=1111Y=100

x257.160.792.060.253.08

2.90

X34.0856.560.7311.92

xt21.1931.1740.244.9314.58

則方程組可以用矩陣方程表示為：

D*X=Y

程序如下：

#include<iostream.h>

classMatrix

(

intRowMax,ColMax;

public:

float**Mat;

Matrix(intrm,intcm)〃構(gòu)造函數(shù)

(

RowMax=rm;

ColMax=cm;

Mat=newfloat*[rm];

for(inti=0;i<rm;i++)

Mat[i]=newfloat[cm];

)

~Matrix(){}〃析構(gòu)函數(shù)

voidInput()〃輸入

(

cout?〃請輸入矩陣：\n〃；

for(inti=0;i<RowMax;i++)

for(intj=0;j<ColMax;j++)

cin?Mat[i][j];

}

intGetRowMaxO{returnRowMax;}〃取矩陣的行數(shù)

intGetColMaxO{returnColMax;}//取矩陣的列數(shù)

voidMatrix::GetRow(constintn)〃取矩陣某一行的全部元素

{

cout〈〈〃Matrix(:,”<<n<X〃)=\n"<<'\t’;

for(inti=0;i<ColMax;i++)

cout<<Mat[n-l][i]?,\t9;

cout?,\n';

}

voidMatrix::GetCol(constintn)〃取矩陣某一列的全部元素

(

cout?”Matrix(:,”<<n<<")=\n";

for(inti=0;KRowMax;i++)

cout?,\tf?Mat[i][n-l]<<，\n;

)

voidPrint()〃顯示矩陣

cout?z，Matrix=\n，z;

for(inti=0;i<RowMax;i++)

for(intj=0;j<ColMax;j++)

cout?Mat[i][j]<<，\tJ;

cout<<，\n，;

)

)

Matrixoperator*(Matrix&m)//矩陣乘法；可直接寫Matrix

tl(2,3),t2(3,2),t3(2,2);t3=tl*t2

(

Matrixt(RowMax,m.ColMax);

for(inti=0;i<RowMax;i++)

for(intj=0;j<m.ColMax;j++)

(

floatsum=0;

for(intk=0;k<ColMax;k++)

sum+=Mat[i]Mat[k][j];

t.Mat[i][j]=sum;

)

returnt;

)

MatrixTrans();〃矩陣轉(zhuǎn)置

MatrixInv();〃矩陣求逆

}；

MatrixMatrix::Trans()

(

constintrm=RowMax;

constintcm=ColMax;

Matrixtmat(cm,rm);

for(inti=0;i<cm;i++)

for(intj=0;j<rm;j++)

tmat.Mat[i][j]=Mat[j][i];

cout??zMatrix，二、n〃；

for(i=0;i<cm;i++)

(

cout?，\V;

for(intj=0;j<rm;j++)

cout?tmat.Mat[i][j]\t*;

cout?,\n';

}

returntmat;

MatrixMatrix::Inv()

(

constintrm=RowMax;

constintcm=ColMax;

Matriximat(rm,cm);

Matrixmat(rm,cm);

for(inti=0;i<rm;i++)

for(intj=0;j<rm;j++)

mat.Mat[i][j]=Mat[i][j];

for(i=0;i<rm;i++)

for(intj=0;j<rm;j++)

imat.Mat[i][j]=0;

for(i=0;i<rm;i++)

imat.Mat[i][i]=l;

for(i=0;i<rm;i++)

(

〃若對角線上元素為0,則進行換行運算，使得對

角線上元素非零

(

intnpi;〃最大值的行下標

for(intj=0;j〈rm;j++)〃找到最大值的下標

if(mat.Mat[j][i]>mat.Mat[m][i])

(

m=j;

)

for(j=0;j<rm;j++)

(

floatn,x;

n=mat.Mat[i][j];

x=imat.Mat[i][j];

mat.Mat[i][j]=mat.Mat[m][j];

mat.Mat[m][j]=n;

imat.Mat[i][j]=imat.Mat[m][j];

imat.Mat[m][j]=x;

}

)

if(mat.[i]!=1)

(

floatt=mat.Mat[i][i];

for(intj=0;j<rm;j++)

mat.Mat[i][j]=mat.Mat[i][j]/t;

imat.Mat[i][j]=imat.Mat[i][j]/t;

)

)

for(intj=0;j<rm;j++)

(

if(i!=j)

(

floatt=mat.Mat[j][i]/mat.Mat[i][i];

for(intk=0;k<rm;k++)

(

mat.Mat[j][k]-=mat.Mat[i][k]*t;

imat.Mat[j][k]-=imat.Mat[i][k]*t;

)

)

)

}

cout?，，Matrix=\nzz;

for(i=0;i<rm;i++)

(

for(intj=0;j<rm;j++)

cout?mat.Mat[i][j]?>\tf;

cout?J\n;

)

cout?，zInverse=\n，z;

for(i=0;i<rm;i++)

(

for(intj=0;j<rm;j++)

cout?imat.Mat[i][j]<<'\t';

cout?,\n';

)

returnimat;

)

voidmain()//main函數(shù)中可自行更改自己想寫的程序段，我這邊寫出部分只是為了舉例

(

Matrixmatl(4,4),mat2(4,1),mat3(4,4);

matl.Input();

mat2.Input();

mat3=matl.Inv()*mat2;

mat3.Print();

)

運行結(jié)果截圖:

2因變量y是自變量x的非線性函數(shù)，已得到實測數(shù)據(jù)見表2,假定y可

以用x的3次多項式近似表示，即y對x的多項式回歸模型為

23

y=bo+b?x+b2x+b3x

請求出模型蠶數(shù)b。、h、bz和b3,計算模型值ym與相對誤差r%并填入表2

中。

表2實測數(shù)據(jù)

xyymr%

-4-210

-2-20

01.1

250

5590

解題思路：將表2中實測數(shù)據(jù)帶入三次多項式，得出包含4個未知數(shù)的線性方程組:

b0+b|(-4)+b2(-4)~+b3(-4)=-210即1-4匕+16b2-64b3=-210

bo+b](_2)+ba(12)~+b(-2)*=—20-_=-

3bo2b|+4b28b320

b+b,O+bO2+bO3=l.1

o23b0=1.1

b+b2+b22+b23=50

0I23bo+2bi+4b2+8b3=5O

bo+M+b^+bj'=590鳳+5片25b2+125b：；=590

只有4個未知數(shù)，而有5個獨立方程,顯然這是一個矛盾方程組，

令b=b0D=1-4-64Y=-210

bi4I-8-20

ba0I01.1

ba24850

1525125590

b稱為模型參數(shù)向量，D稱為結(jié)構(gòu)矩陣,Y稱為觀測值向量，則線性方程組可以用矩陣

方程表示為：

Db=Y

這個矛盾方程組可以用最小二乘法求出最小二乘解。

程序設計：

#include<iostream.h>

classMatrix

intRowMax,ColMax;

public:

float**Mat;

Matrix(intrm,intcm)

RowMax=rm;

ColMax=cm;

Mat二newfloat*[rm];

for(inti=0;i<rm;i++)

Mat[i]=newfloat[cm];

)

"Matrix(){}

voidInput()

(

cout<〈〃請輸入實測數(shù)據(jù)(矩陣形式):\n〃;

for(inti=0;i<RowMax;i++)

for(intj=0;j<ColMax;j++)

cin?Mat[i][j];

intGetRowMaxO{returnRowMax;}

intGetColMaxO(returnColMax;}

voidMatrix::GetRow(constintn)

(

cout?z*\nMatrix=\nz，;

for(inti=0;i<ColMax;i++)

cout?Mat[n-l][i]?,\tf;

cout?,\n';

)

voidMatrix::GetCol(constintn)

(

constintrm=RowMax;

constintcm=0;

Matrixmatl(rm,cm);

for(inti=0;i<rm;i++)

matl.Mat[i][0]=Mat[i][n-1];

cout?/znMatrix=\n/，;

for(int1=0;KRowMax;1++)

(

cout?matl.Mat[1][0]<<?\tJ;

cout?，\n;

)

)

voidPrint()

(

cout?〃Matrix=\n〃；

for(inti=0;i〈RowMax;i++)

(

for(intj=0;j<ColMax;j++)

cout?Mat[i][j]<<，\t';

cout<<，\n;

)

)

Matrixoperator+(Matrix&m)

(

Matrixt(RowMax,ColMax);

for(inti=0;i<RowMax;i++)

for(intj=0;j<ColMax;j++)

t.Mat[i][j]=Mat[i][j]+m.Mat[i][j];

returnt;

)

Matrixoperator-(Matrix&m)

(

Matrixt(RowMax,ColMax);

for(inti=0;i<RowMax;i++)

for(intj=0;j<ColMax;j++)

t.Mat[i][j]=Mat[i]Mat[i][j];

returnt;

}

Matrixoperator*(Matrix&m)

(

Matrixt(RowMax,m.ColMax);

for(inti=0;i<RowMax;i++)

for(intj=0;j<m.ColMax;j++)

(

floatsum=0;

for(intk=0;k<ColMax;k++)

sum+=Mat[i][k]*m.Mat[k][j];

t.Mat[i][j]=sum;

)

returnt;

)

MatrixTrans();

MatrixInv();

}；

MatrixMatrix::Trans()

(

constintrm=RowMax;

constintcm=ColMax;

Matrixtmat(cm-l,rm);

for(inti=0;i<cm-l;i++)

for(intj=0;j<rm;j++)

tmat.Mat[i][j]=Mat[j][i];

cout?，，Matrix,=\n〃；

for(i=0;i<cm-l;i++)

(

cout?,\tJ;

for(intj=0;j<rm;j++)

cout?tmat.Mat[i][j]<<，\t，;

cout?,\n';

}

returntmat;

MatrixMatrix::Inv()

(

constintrm=RowMax;

constintcm=ColMax;

Matriximat(rm,cm-l);

Matrixmat(rm,cm);

Matrixmatl(rm,cm);

for(inti=0;i<rm;i++)

for(intj=0;j<rm;j++)

mat.Mat[i][j]=Mat[i][j];

for(int1=0;Krm;l++)

matl.Mat[1][0]=Mat[1][cm-1];

for(i=0;i<rm;i++)

for(intj=0;j<rm;j++)

imat.Mat[i][j]=0;

for(i=0;i<rm;i++)

imat.Mat[i][i]=l;

inth,p,q,j,k;

floattemp,temp2;

for(i=0;i<rm;i++)

(

for(k=i;k<rm;k++)

(

if(mat.Mat[k][i]!=0)

(

h=k;

k=rm;

)

)

if(k==rm)

(

cout<<〃該矩陣沒有相應的逆矩陣。〃Gendl;

i=rm;

)

else

(

if(h!=i)

(

for(j=0;j<rm;j++)

mat.Mat[i][j]=mat.Mat[i][j]+mat.Mat[h][j];

matl.Mat[i][j]=matl.Mat[i][j]+matl.Mat[h][j];

imat.Mat[i][j]=imat.Mat[i][j]+imat.Mat[h][j];

)

)

temp=mat.Mat[i][i];

for(j=0;j<rm;j++)

(

mat.Mat[i][j]=mat.Mat[i][j]/temp;

matl.Mat[i][j]=matl.Mat[i][j]/temp;

imat.Mat[i][j]=imat.Mat[i][j]/temp;

}

for(p=0;p<rm;p++)

(

if(i!=p)

(

temp2=mat.Mat[p][i];

for(q=0;q<rm;q++)

(

mat.Mat[p][q]=mat.Mat[p][q]-temp2*mat.Mat[i][q];

matl.Mat[p][q]=matl.Mat[p][q]-temp2*matl.Mat[i][q];

imat.Mat[p][q]=imat.Mat[p][q]-temp2*imat.Mat[i][q];

)

)

)

}

)

if(k!=rm)

(

cout?/zanswer=\n，z;

for(l=0;l<rm;l++)

(

cout?matl.Mat[1][0]?>\t';

cout?,\n;

)

)

returnimat;

)

voidmain()

{floatbO,bl,b2,b3;

inti,j;

floatym[5],r[5];

Matrixmatl(6,2),mat3(4,6),mat4(4,4),mat5(4,5),mat6(4,1),mat7(4,1),mat8(5,5);

matl.Input();

for(i=0;i<5;i++)

(

mat8.Mat[i][0]=l;

mat8.Mat[i][l]=matl.Mat[i][0];

mat8.Mat[i][2]=matl.Mat[i][0]*matl.Mat[i][0];

mat8.Mat[i][3]=matl.Mat[i][0]*matl.Mat[i][0]*matl.Mat[i][0];

mat8.Mat[i][4]=matl.Mat[i][1];

)

cout<<“變換后："<<endl;

for(i=0;i<5;i++)

(

for(j=0;j<4;j++)

cout<<mat8.Mat[i][j]?>W;

cout?>\n';

)

mat3=mat8.Trans();

mat5=mat3*mat8;

for(i=0;i<4;i++)

(

mat6.Mat[i][0]=mat5.Mat[i][4];

)

mat4=mat5.Inv();

mat7=mat4*mat6;

b0=mat7.Mat[0][0];

bl=mat7.Mat[1][0];

b2=mat7.Mat[2][0];

b3=mat7.Mat[3][0];

cout<〈〃回歸模型是：y=〃?b(K<〃+〃<<bl<<〃*x+〃<Xb2<<〃*x2+〃<<b3<<〃*x3〃〈<endl;

ym[0]=bO+bl*mat8.Mat[0][l]+b2*mat8.Mat[0][2]+b3*mat8.Mat[0][3];

ym[1]=bO+b1*mat8.Mat[1][1]+b2*mat8.Mat[1][2]+b3*mat8.[3];

ym[2]=bO+b1*mat8.Mat[2][l]+b2*mat8.Mat[2][2]+b3*mat8.Mat[2][3];

ym[3]=b0+bl*mat8.Mat[3][1]+b2*mat8.Mat[3][2]+b3*mat8.Mat[3][3];

ym[4]=b0+bl*mat8.Mat[4][1]+b2*mat8.Mat[4][2]+b3*mat8.Mat[4][3];

cout<<“模型值"<<endl;

for(i=0;i<4;i++)

(

cout?，/ym/，<<(i+l)?，，=z/<<ym[i]<<endl;

)

cout<<"相對誤差"<<endl;

for(i=0;i<4;i++)

rLiJ=(ymLiJ_matl.MatLiJLU)/matl.MatLiJLU；

(i+1)?/，=w?r[i]<<endl;

)

cout<<“回歸模型合理"<<endl;

)

運行結(jié)果截圖：

*E:\swl00103\T2\Debug\T2.exe*日回口

▲

1-416-64

1-24-8

1000

1248

1525125

Matrixj=

11111

-4-2025

1640425

-64-808125

answer=

1.817

1.4884

3.18194

4.01003

回歸模型是：y=l-81699+1.48845*x+3.18193*x2+4.01002*x3

模型推

yml=-209.867

yn2=-20.5124

yn3=1.81699

ym4=49.6018

相對誤差

rl=-0.000631496

r2=0.0256183

r3=0.651807

r4=-0.00796394

回歸模型合理

Pressanykeytocontinue

3請用行線性變換求下面三元一次線性方程組的解和系數(shù)矩陣的逆矩陣。

-x1+x2+x3=13

x1-x2+x3=9

x1+x2-x3=1

解題思路：本題可利用矩陣的行線性變換求解該線性方程組與逆矩陣。已知當用線性

變換把系數(shù)矩陣變換成單位矩陣時，放常數(shù)項向量的列就是方程組的解，放單位矩陣

的位置就是系數(shù)矩陣的逆矩陣。

本題中令X=X1A=-111B=13

Xz1-119

X311-11

X稱為未知數(shù)向量，A稱為系數(shù)矩陣，B稱為常數(shù)項向量。

則線性方程組可以用矩陣方程表示為:

A*X=B

用A的逆矩陣A'左乘上面矩陣方程兩邊，左邊的就變換成了單位矩陣和解的乘積，

故有

X=A'B

所以，X是方程組的解，A)是該系數(shù)矩陣的逆矩陣。

此題根據(jù)類運算做出，最先定義了一個矩陣類，然后根據(jù)其中定義的Input,Print,

相乘,矩陣求逆等函數(shù)書寫主函數(shù)，先輸入系數(shù)矩陣matl,再輸入等式右邊的矩陣mat3,

將matl求逆賦予mat2再左乘mat3,即得到所要求的值xl,x2,x3.同時輸出matl的逆

矩陣mat2.

程序如下：

#include<iostream.h>

classMatrix

(

intRowMax,ColMax;

public:

float**Mat;

Matrix(intrm,intcm)〃構(gòu)造函數(shù)

(

RowMax=rm;

ColMax=cm;

Mat二newfloat*[rm];

for(inti=0;i<rm;i++)

Mat⑴=newfloat[cm];

)

~Matrix(){}〃析構(gòu)函數(shù)

voidInput()〃輸入

(

coutG〃請輸入矩陣：\n〃；

for(inti=0;i<RowMax;i++)

for(intj=0;j<ColMax;j++)

cin?Mat[i][j];

)

intGetRowMaxO{returnRowMax;}〃取矩陣的行數(shù)

intGetColMaxO{returnColMax;}//取矩陣的列數(shù)

voidMatrix::GetRow(constintn)〃取矩陣某一行的全部元素

(

cout?，zMatrix(:,/，<<n<<，/)=\n，，?,\t';

for(inti=0;i<ColMax;i++)

cout?Mat[n-l][i]?，\t);

cout?,\n;

}

voidMatrix::GetCol(constintn)〃取矩陣某一列的全部元素

cout<X"Matrix(:,”<<n<<")=\n";

for(inti=0;i<RowMax;i++)

cout<<J\t*?Mat[i][n-l]<<，\n;

)

voidPrint()〃顯示矩陣

(

cout?〃Matrix=\n〃；

for(inti=0;i<RowMax;i++)

(

for(intj=0;j<ColMax;j++)

cout?Mat[i][j]<<?\t*;

cout<<，\n';

}

)

Matrixoperator*(Matrix&m)〃矩陣乘法；可直接寫Matrix11(2,3),t2(3,2),t3(2,2);t3=tl*t2

(

Matrixt(RowMax,m.ColMax);

for(inti=0;i<RowMax;i++)

for(intj=0;j<m.ColMax;j++)

(

floatsum=0;

for(intk=0;k<ColMax;k++)

sum+=Mat[i][k]*m.Mat[k][j];

t.Mat[i][j]=sum;

}

returnt;

)

MatrixTrans();〃矩陣轉(zhuǎn)置

MatrixInv();〃矩陣求逆

}；

MatrixMatrix::Trans()

(

constintrm=RowMax;

constintcm=ColMax;

Matrixtmat(cm,rm);

for(inti=0;i<cm;i++)

for(intj=0;j<rm;j++)

tmat.Mat[i][j]=Mat[j][i];

cout?/zMatrix,=\n〃;

for(i=0;i<cm;i++)

cout?，\t?;

for(intj=0;j<rm;j++)

cout<<tmat.Mat[i][j]?，\tJ;

cout?，\n;

}

returntmat;

)

MatrixMatrix::Inv()

(

constintrm=RowMax;

constintcm=ColMax;

Matriximat(rm,cm);

Matrixmat(rm,cm);

for(inti=0;i<rm;i++)

for(intj=0;j<rm;j++)

mat.Mat[i][j];

for(i=O;i<rm;i++)

for(intj=0;j<rm;j++)

imat.Mat[i][j]=0;

for(i=0;i<rm;i++)

imat.Mat[i][i]=l;

for(i=0;i<rm;i++)

(

1￡(0=1^1.\^江月[H)〃若對角線上元素為0,則進行換行運算，使得對角線上

元素非零

(

intm二i;〃最大值的行下標

for(intj=0;j++)〃找到最大值的下標

if(mat.Mat[j][i]>mat.Mat[m][i])

(

m=j;

)

for(j=0;j<rm;j++)

(

floatn,x;

n二mat.Mat[i][j];

x=imat.Mat[i][j];

mat.Mat[i][j]=mat.Mat[m][j];

mat.Mat[m][j]=n;

imat.Mat[i][j]=imat.Mat[m][j];

imat.Mat[m][j]=x;

if(mat.Mat[i][i]!=1)

floatt=mat.Mat[i][i];

for(intj=0;j<rm;j++)

(

mat.Mat[i][j]=mat.Mat[i]

imat.Mat[i][Mat[i]

)

)

for(intj=0;j<rm;j++)

(

if(i!=j)

(

floatt=mat.Mat[j][i]/mat.Mat[i][i];

for(intk=0;k<rm;k++)

(

mat.Mat[j][k]-=mat.Mat[i][k]*t;

imat.Mat[j][k]-=imat.Mat[i]

)

)

)

)

cout?/zMatrix=\n，z;

for(i=0;i<rm;i++)

(

for(intj=0;j<rm;j++)

cout?mat.Mat[i][j]\tf;

cout?，\n;

)

cout<<〃Inverse=\n〃；

for(i=0;i<rm;i++)

(

for(intj=0;j<rm;j++)

cout?imat.Mat[i][j]?,\tf;

cout?,\n;

)

returnimat;

)

voidmain()〃main函數(shù)中可自行更改自己想寫的程序段，我這邊寫出部分只是為了舉例

{

Matrixmatl(3,3),mat2(3,3),mat3(3,1),mat4(3,3);

matl.Input();

mat3.Input();

mat2=matl.Inv();

mat4=mat2:4:mat3;

mat2.Print();

mat4.Print();

運行結(jié)果界面截圖

4在表1中增加兩種組分的化驗品位如表3,可以列出六個線性方程，顯

然是矛盾的方程組，請用最小二乘法求精/尾礦產(chǎn)率的最小二乘解，并填

入表3相應位置。（礦物專業(yè)必做，生物專業(yè)可選）

表3各礦流的化驗品位

礦產(chǎn)化驗品位%

率

流%PbZnSFeCd

原礦100.003.0811.9214.589.320.0162

鉛精礦Xi57.164.0821.1910.260.0075

鋅精礦

X20.7956.5631.174.750.0880

硫精礦

X32.062.9040.2441.230.0045

尾礦X40.250.734.935.190.0004

解題思路：利用礦量和五種組分的進出平衡可列出六個線性方程組:

X1+X+X+

23x4=100

57.16xt+0.79X2+2.06X3+0.25x4=3.08

4.08x,+56.56x2+2.90X3+0.73x4=11.92

21.19xt+31.17X2+40.24x3+4.93x.,=14.58

10.26xi+4.75x2+41.23x3+5.19X4=9.32

0.0075xi+0.0880x2+0.0045x3+0.0004x,,=0.0162

由于獨立方程個數(shù)多于未知數(shù)個數(shù)，屬于矛盾方程組,

令

X=XiD=1111Y=100

X257.160.792.060.253.08

X34.0856.562.900.7311.92

X421.1931.1740.244.9314.58

10.264.7541.235.199.32

0.00750.08800.00450.00040.0162

X是未知數(shù)向量，D是系數(shù)矩陣，Y是常數(shù)項向量（它們的元素則用帶下標的對應小寫

字母

表示），則方程組可以用矩陣方程表示為：

D*X=Y

當獨立方程個數(shù)多于未知數(shù)時，就構(gòu)成矛盾方程組，矛盾方程組是無解的。

用最小二乘法可以合理地解決這類矛盾方程組的求解問題。最小二乘法的原理是

尋

找一組合適的解，使得方程組等號兩邊差值的平方之和達到最小值，這組解稱為最小

乘解，最小二乘解是唯一的。

對于一般的含有p未知數(shù)n個線性方程的矛盾方程組，可以構(gòu)造平方和函數(shù)Q如

nP

Q(xi,x2,...,Xp)=￡(y-dZdij+Xj))~2

i=lj=l

平方和Q是未知數(shù)x“X2,...,Xp的函數(shù)，若Q有極小值，則它對Xi,X2,...,X),的p

個偏導數(shù)應該等于零，由此得出的p個線性方程就能確定一組唯一的解，即矛盾方程

組的最小二乘解。

用D'表示D的轉(zhuǎn)置矩陣，令A=D'*D和B=D'*Y,最小二乘解X滿足下面的矩陣方

程：

A*X=B

這是有唯一解的線性方程組，可以用前面介紹過的線性方程組求解方法求解。A稱為

信

息矩陣，是個pxp的方陣，行列數(shù)等于未知數(shù)的個數(shù)。A的逆矩陣稱為相關矩陣，相

關矩

陣在回歸分析中檢驗各因子的顯著性時是很有用的。

矩陣A=D'*D和B=D'*Y可按式計算

n

==

aij=Xdki*dkj(il>???,P>jl,??-,p)

k=l

n

bi=￡dki*yk(i=l,...,p)

k=l

程序如下：

#include<iostream.h>

classMatrix

(

intRowMax,ColMax;

public:

float**Mat;

Matrix(intrm,intcm)〃構(gòu)造函數(shù)

(

RowMax=rm;

ColMax=cm;

Mat=newfloat*[rm];

for(inti=0;i<rm;i++)

Mat[i]=newfloat[cm];

)

'Matrix(){}〃析構(gòu)函數(shù)

voidInput()〃輸入

(

cout<<〃請輸入矩陣：\n〃；

for(inti=0;KRowMax;i++)

for(intj=0;j<ColMax;j++)

cin?Mat[i][j];

)

intGetRowMaxO{returnRowMax;}〃取矩陣的行數(shù)

intGetColMaxO(returnColMax;}//取矩陣的列數(shù)

voidMatrix::GetRow(constintn)//取矩陣某一行的全部元素

cout?"Matrix(:,”<<n<<")=\n"<<'\t’;

for(inti=0;i<ColMax;i++)

cout<<Mat[n-l][i]?，\tJ;

cout?,\n';

)

voidMatrix::GetCol(constintn)〃取矩陣某一列的全部元素

(

cout<<"Matrix(:,"<<n<<")=\n";

for(inti=0;i<RowMax;i++)

cout?，\t*?Mat[i][n-l]<<，\n;

)

voidPrint()〃顯示矩陣

(

cout<X〃Matrix=\n”;

for(inti=0;KRowMax;i++)

(

for(intj=0;j<ColMax;j++)

cout?Mat[i][j]?>\t*;

cout<<,\n;

}

}

Matrixoperator*(Matrix&m)〃矩陣乘法；可直寫Matrix

tl(2,3),t2(3,2),t3⑵2);t3=tl*t2

(

Matrixt(RowMax,m.ColMax);

for(inti=0;i<RowMax;i++)

for(intj=0;j<m.ColMax;j++)

(

floatsum=0;

for(intkz:0;k<ColMax;k++)

sum+=Mat[i][k]*m.Mat[k][j];

t.Mat[i][j]=sum;

}

returnt;

)

MatrixTrans();〃矩陣轉(zhuǎn)置

MatrixInv();〃矩陣求逆

)；

MatrixMatrix::Trans()

{

constintrm=RowMax;

constintcm=ColMax;

Matrixtmat(cm,rm);

for(inti=0;i<cm;i++)

for(intj=0;j<rm;j++)

tmat.Mat[i][j]=Mat[j][i];

cout?，，Matrix，=\n〃;

for(i=0;i<cm;i++)

(

cout?,\t*;

for(intj=0;j<rm;j++)

cout<<tmat.Mat[i][j]?'\t';

cout?，\n，;

)

returntmat;

)

MatrixMatrix::Inv()

(

constintrm=RowMax;

constintcm=ColMax;

Matriximat(rm,cm);

Matrixmat(rm,cm);

for(inti=0;i<rm;i++)

for(intj=0;j<rm;j++)

mat.Mat[i][j]=Mat[i][j];

for(i=0;i<rm;i++)

for(intj=0;j<rm;j++)

imat.Mat[i][j]=0;

for(i=0;i<rm;i++)

imat.Mat[i][i]=l;

for(i=0;i<rm;i++)

(

for(intj=0;j<rm;j++)

(

if(i!=j)

(

floatt=mat.Mat[j][i]/mat.Mat[i][i];

for(intk=0;k<rm;k++)

(

mat.Mat[j][k]-=mat.Mat[i][k]*t;

imat.Mat[j][k]-=imat.Mat[i][k]*t;

}

}

}

)

for(i=0;i<rm;i++)

if(mat.Mat[i][i]!=1)

(

floatt=mat.Mat[i][i];

for(intj=0;j<rm;j++)

(

mat.Mat[i][j]=mat.Mat[i][j]/t;

imat.[j]=imat.Mat[i][j]/t;

)

)

cout?，zMatrix=\nz/;

for(i=0;i<rm;i++)

(

for(intj=0;j<rm;j++)

cout?mat.Mat[i][j]<<，\t);

cout?,\n;

)

cout?〃Inverse=\n〃;

for(i=0;i<rm;i++)

(

for(intj=0;j<rm;j++)

cout<<imat.Mat[i][j]?>\t);

cout?，\n;

)

returnimat;

)

voidmain()〃main函數(shù)中可.自行更改自己想寫的程序段，我這邊寫出部分只是為了舉例

(

Matrixmatl(6,4),mat2(4,6),mat3(4,4),mat4(6,1),mat5(4,1),mat6(4,1);

matl.Input();

mat4.Input();

mat2=matl.Trans();

mat3=mat2*matl;

mat5=mat2*mat4;

mat6=mat3.Inv()*mat5;

mat6.Print();

}

運行結(jié)果界面截圖:

不"D:\s"100103\T4\Debug\I4.exe"

k0.264.7541.235.19

e.00750.00800.00450.0004

情輸入矩陣：

1003081192：L4589321.62

Matrix*=

157.164.0821.1910.260.0075

10.7956.5631.714.750.008

12.062.940.2441.230.0045

10.250.734.935.190.0004

f!atrix=

bl1.30491e-009-1.04254e-0095.54498e-012

b.91?19e-0081-6.68624e-0089.95294e-012

b.96474e-0073.18114e-00?1