/**/ /////////////////////////////////////// /
//?????????????????????????????????????//
//????????????????????????????????????//
//數(shù)值分析??直接三角分解法求方程組????//
//????????????????????????????????????//
//??在VC2005下編譯通過(guò)，由于VC6.0的???//
//???不標(biāo)準(zhǔn)會(huì)有錯(cuò)誤???????????????????//
//????????????????????????????????????//
//???2006.5.30???v0.1?????????????????//
/**/////////////////////////////////////////

#include?"iostream"
using?namespace?std;

class?Matrix
{
private:
?double**?A;??????//矩陣A
?double?*b;???????//向量b
public:
?int?size;
?Matrix(int?);
?~Matrix();
friend?double*?Dooli(Matrix&?);
?void?Input();
?void?Disp();

};

Matrix::Matrix(int?x)
{
?size=x;

?//為向量b分配空間并初始化為0
?b=new?double?[x];
?for(int?j=0;j<x;j++)
??b[j]=0;

?//為向量A分配空間并初始化為0
?A=new?double*?[x];
?for(int?i=0;i<x;i++)
??A[i]=new?double?[x];
?for(int?m=0;m<x;m++)
??for(int?n=0;n<x;n++)
???A[m][n]=0;
}

Matrix::~Matrix()
{
????cout<<"正在析構(gòu)中~~~~"<<endl;
????delete?b;
????for(int?i=0;i<size;i++)
????????delete?A[i];
????delete?A;
}

void?Matrix::Disp()
{
?for(int?i=0;i<size;i++)
?{
??for(int?j=0;j<size;j++)
???cout<<A[i][j]<<"??";
??cout<<endl;
?}
}

void?Matrix::Input()
{
?cout<<"請(qǐng)輸入A:"<<endl;
?for(int?i=0;i<size;i++)
??for(int?j=0;j<size;j++){
???cout<<"第"<<i+1<<"行"<<"第"<<j+1<<"列:"<<endl;
??cin>>A[i][j];
??}
???cout<<"請(qǐng)輸入b:"<<endl;
?for(int?j=0;j<size;j++){
??cout<<"第"<<j+1<<"個(gè):"<<endl;
???cin>>b[j];
?}
?
}
?
double*?Dooli(Matrix&?A)
{
?double?*Xn=new?double?[A.size];
?Matrix?L(A.size),U(A.size);

?//分別求得U，L的第一行與第一列
???for(int?i=0;i<A.size;i++)
??????U.A[0][i]=A.A[0][i];
???for(int?j=1;j<A.size;j++)
??????L.A[j][0]=A.A[j][0]/U.A[0][0];

//分別求得U，L的第r行，第r列
?????double?temp1=0,temp2=0;
?for(int?r=1;r<A.size;r++){
?????//U
?????for(int?i=r;i<A.size;i++){
?????????for(int?k=0;k<r-1;k++)
????????????temp1=temp1+L.A[r][k]*U.A[k][i];?
????????????U.A[r][i]=A.A[r][i]-temp1;
?????}
?????//L
?????for(int?i=r+1;i<A.size;i++){
??????????for(int?k=0;k<r-1;k++)
????????????temp2=temp2+L.A[i][k]*U.A[k][r];
?????????????L.A[i][r]=(A.A[i][r]-temp2)/U.A[r][r];
?????}
?}
?cout<<"計(jì)算U得:"<<endl;
?U.Disp();
?cout<<"計(jì)算L的:"<<endl;
?L.Disp();
?
?double?*Y=new?double?[A.size];

?Y[0]=A.b[0];
?for(int?i=1;i<A.size;i++?){
?????double?temp3=0;
?????for(int?k=0;k<i-1;k++)
?????????temp3=temp3+L.A[i][k]*Y[k];
?????Y[i]=A.b[i]-temp3;
?}

?Xn[A.size-1]=Y[A.size-1]/U.A[A.size-1][A.size-1];
?for(int?i=A.size-1;i>=0;i--){
?????double?temp4=0;
?????for(int?k=i+1;k<A.size;k++)
?????????temp4=temp4+U.A[i][k]*Xn[k];
?????Xn[i]=(Y[i]-temp4)/U.A[i][i];
?}
?
return?Xn;
}
?
int?main()
{
?Matrix?B(4);
?B.Input();
?double?*X;
?X=Dooli(B);
?cout<<"~~~~解得:"<<endl;
?for(int?i=0;i<B.size;i++)
?????cout<<"X["<<i<<"]:"<<X[i]<<"?";
?cout<<endl<<"呵呵呵呵呵";
?return?0;
}

試驗(yàn)總結(jié)：將高斯消去法改寫為緊湊形式，可以直接從矩陣 A 的元素的導(dǎo)計(jì)算 L ， U 元素的遞推公式，而不需任何中間步驟，一旦實(shí)現(xiàn)了矩陣 A 的 U ， L 分解那么就等價(jià)于求解兩個(gè)三角形方程組。

? 注意： 編成語(yǔ)言中的數(shù)組以 ’0’ 為首元素，數(shù)組的一位偏移最容易出錯(cuò)；

???????????? 注意變量的作用域；