华剑缘
一切都在这个过程中获得,将那些目标埋藏于心中
posts - 19,comments - 20,trackbacks - 0
/////////////////////////////////////// /
//                                     //
//                                    //
//数值分析  直接三角分解法求方程组    //
//                                    //
//  在VC2005下编译通过,由于VC6.0的   //
//   不标准会有错误                   //
//                                    //
//   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
<<"正在析构中~~~~"<<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
<<"请输入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
<<"请输入b:"<<endl;
 
for(int j=0;j<size;j++){
  cout
<<""<<j+1<<"个:"<<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
<<"计算U得:"<<endl;
 U.Disp();
 cout
<<"计算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;
}


试验总结:将高斯消去法改写为紧凑形式,可以直接从矩阵 A 的元素的导计算 L U 元素的递推公式,而不需任何中间步骤,一旦实现了矩阵 A U L 分解那么就等价于求解两个三角形方程组。

  注意: 编成语言中的数组以 0’ 为首元素,数组的一位偏移最容易出错;

             注意变量的作用域;

posted on 2006-05-30 21:22 华剑缘 阅读(1290) 评论(0)  编辑 收藏 引用

只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理