often.txt

几个常用的算法

☉ 设为首页 | 加入收藏 | 广告服务. 
    
  
 
网站首页 C语言 C++语言 Windows编程 接口&驱动 汇编语言 编程资料 软 件 其它文章 等级考试辅导专栏 
 
 
 
　 当前位置：网站首页>>C++语言>>C++基础>>学习园地 双击自动滚屏   
 
 
 
矩阵直接三角分解法的c++语言实现  
 

 
发表日期：2003年10月28日  出处：原创  作者：萧何  已经有2045位读者读过此文 
 


   矩阵的计算是非常有用的，也有很多的算法，其中直接三角分解法是比较高效而且误差相对较小的一种算法，我用面向对象的设计思想，将直接三角分解法用C++实现，我的程序还有很多的不足，希望各位给予指教。（直接三角分解法的具体算法 清参照《计算机数值算法》）

#include<iostream.h>

class GS2
{
public:
  GS2 (int i,int j,double e)               //构造矩阵
  {
   x=i;
   y=j;
         eps=e;
   a=new double[x*(x+y)];
  }
  void build();
  void return_x();                   //回代求 x[i]
  bool judge (int);              //判断是否|P|<=EPS
  bool Doolittle ();              //doolittle直接三角分解法的主体程序
  void changcemain (int);       //选主元
  void show();                  //打印x
  void print();
  ~GS2()
  {
   delete []a;
  }
protected:
 double *a,eps;
    int x,y;
};

void GS2::build()
{
 for(int i=0;i<x;i++)
  for(int j=0;j<x+y;j++)
  {
   cout<<"x["<<i+1<<"]["<<j+1<<"]=";

   cin>>a[j+(x+y)*i];
  }
}
//
void GS2::print()
{
 for(int i=0;i<x;i++)
  for(int j=0;j<x+y;j++)
  {
   cout<<"x["<<i+1<<"]["<<j+1<<"]="

   <<a[j+(x+y)*i]<<endl;
  }
}
//

bool GS2::Doolittle()
{
 double temp;
    

 for(int r=0;r<x-1;r++)
 {  
 
  if(r==0)
  {
   changcemain(r);
        if (!judge(r))
               return false;
        else
     for(int i=r+1;i<x;i++)
       a[r+(x+y)*i]/=a[r+(x+y)*r];
  }  
        else
  {
      for(int i=r;i<x;i++)
   {
    temp=0;
     for(int k=0;k<r;k++)
      temp+=a[k+(x+y)*i]*a[r+(x+y)*k];
             
     a[r+(x+y)*i]-=temp;
   }
  
   changcemain(r);
   
   if (!judge(r))
           return false;

   for( i=r+1;i<x;i++)
   {
    a[r+(x+y)*i]/=a[r+(x+y)*r];
   }
   for(i=r+1;i<x+y;i++)
   {
    temp=0;
    for(int k=0;k<r;k++)
    {
     temp+=a[k+(x+y)*r]*a[i+(x+y)*k];
    }
         a[i+(x+y)*r]-=temp;
   }
  }  //end else

 }

 for(int i=x-1;i<x+y;i++)
 {
  temp=0;
  for(int k=0;k<x-1;k++)
  {
   temp+=a[k+(x+y)*(x-1)]*a[i+(x+y)*k];
  }
  a[i+(x+y)*(x-1)]-=temp;
 }
     
    if (!judge(x-1))
           return false;
    
 return true;
}

void GS2::return_x()
{
 double temp;
 for(int j=0;j<y;j++)
 {

    for(int i=x-1;i>=0;i--)
 {
  temp=0;
  for(int k=x-1;k>i;k--)
  {
          temp+=a[k+(x+y)*i]*a[(x+j)+(x+y)*k];
  }
    a[(x+j)+(x+y)*i]-=temp;
 
 a[(x+j)+(x+y)*i]/=a[i+(x+y)*i];
 }
 }
}

void GS2::show()
{
  for(int j=0;j<y;j++)
  {
   cout<<endl;
   for(int i=0;i<x;i++)
   {
    cout<<"x["<<j+1<<"]["<<i+1<<"]="<<a[(x+j)+(x+y)*i]<<'\t';
   }
  }
}

void GS2::changcemain(int r)
{
 double temp=a[r+(x+y)*r];
    int t=r;
 for(int i=r+1;i<x;i++)
 {
  if( a[r+(x+y)*i] > temp )
  {  
   temp=a[r+(x+y)*i];
       t=i;
  }
 }
 if (t!=r)
  for(i=0;i<x+y;i++)
  {
   temp=a[i+(x+y)*r];
   a[i+(x+y)*r]=a[i+(x+y)*t];
   a[i+(x+y)*t]=temp;
  }
}

bool GS2::judge(int r)
{ 
 if (a[r+(x+y)*r]>0)
   return ((a[r+(x+y)*r]>eps )?true:false);
 else
      return (((-a[r+(x+y)*r])>eps )?true:false);
}
   
int main()
{
 int n,m;
 double e;
 cout<<"please input n , m and EPS :";
 cin>>n>>m>>e;

 GS2 gs2(n,m,e);
 gs2.build();
   
 if(!gs2.Doolittle())
 {
  cout<<"此方程无解!"<<endl;
  return 1;
 }
 gs2.print();
 gs2.return_x();
 gs2.show();
 return 0;
}

     
几个常用的算法，包括消元，迭代，三角分解等，都是解普通线性方程组的   
  void   gauss(double   **a,double   *b,double   *x,int   n)   
  //高斯列主元消去法求解线性方程组   
  {   
  int   i,j,k,ik;   
  double   mik,temp;   
  for(k=0;k<n;k++)   
  {   
  mik=-1;   
  for(i=k;i<n;i++)   
  if(ABS(a[i][k])>mik)   
  {   
  mik=ABS(a[i][k]);   
  ik=i;   
  }   
  for(j=k;j<n;j++)   
  SWAP(a[ik][j],a[k][j]);   
  SWAP(b[k],b[ik]);   
  b[k]/=a[k][k];   
  for(i=n-1;i>=k;i--)   
  a[k][i]/=a[k][k];   
  for(i=k+1;i<n;i++)   
  {   
  b[i]-=a[i][k]*b[k];   
  for(j=n-1;j>=k;j--)   
  a[i][j]-=a[i][k]*a[k][j];   
  }   
  }   
  for(i=n-1;i>=0;i--)   
  {   
  x[i]=b[i];   
  for(j=i+1;j<n;j++)   
  x[i]-=a[i][j]*x[j];   
  }   
  }   
    
  void   jacobi(double   **a,double   *b,double   *x,int   n,int   Times)   
  //雅可比迭代求线性方程组   
  {   
  int   i,j,t;   
  double   *xx=new   double[n];   
  for(t=0;t<Times;t++)   
  {   
  for(i=0;i<n;i++)   
  {   
  xx[i]=b[i];   
  for(j=0;j<n;j++)   
  if(j!=i)   
  xx[i]-=a[i][j]*x[j];   
  xx[i]/=a[i][i];   
  }   
  for(i=0;i<n;i++)   
  x[i]=xx[i];   
  }   
  delete   xx;   
  }   
    
  void   gs(double   **a,double   *b,double   *x,int   n,int   Times)   
  //高斯-赛德尔迭代法求解线性方程组   
  {   
  int   i,j,t;   
  for(t=0;t<Times;t++)   
  {   
  for(i=0;i<n;i++)   
  {   
  x[i]=b[i];   
  for(j=0;j<n;j++)   
  if(j!=i)   
  x[i]-=a[i][j]*x[j];   
  x[i]/=a[i][i];   
  }   
  }   
  }   
    
  void   doolittle_deco(double   **a,int   n)   
  //Doolittle三角分解，分解后的矩阵保存在a中   
  {   
  int   i,j,k,m;   
  for(k=0;k<n;k++)   
  {   
  for(j=k;j<n;j++)   
  for(m=0;m<k;m++)   
  a[k][j]-=a[k][m]*a[m][j];   
  for(i=k+1;i<n;i++)   
  {   
  for(m=0;m<k;m++)   
  a[i][k]-=a[i][m]*a[m][k];   
  a[i][k]/=a[k][k];   
  }   
  }   
  }   
    
  void   doolittle_back(double   **a,double   *b,double   *x,int   n)   
  //Doolittle三角分解法求解线性方程组的回代过程   
  {   
  int   i,j;   
  double   *y;   
  y=new   double[n];   
  for(i=0;i<n;i++)   
  {   
  y[i]=b[i];   
  for(j=0;j<i;j++)   
  y[i]-=a[i][j]*y[j];   
  }   
  for(i=n-1;i>=0;i--)   
  {   
  x[i]=y[i];   
  for(j=i+1;j<n;j++)   
  x[i]-=a[i][j]*x[j];   
  x[i]/=a[i][i];   
  }   
  delete   y;   
  }   
    
  void   doolittle(double   **a,double   *b,double   *x,int   n)   
  //Doolittle三角分解法求解线性方程组   
  {   
  doolittle_deco(a,n);   
  doolittle_back(a,b,x,n);   
  }
Top

11 楼saint001（saint001）回复于 2003-12-17 21:34:37 得分 14前面加上   
  #include   "stdio.h"   
  #include   "math.h"   
  #define   ABS(x)   (x)>0?(x):-(x)   
  #define   SWAP(a,b)   {temp=(a);(a)=(b);(b)=temp;}   
    
  对它们的调用   
  int   main(int   argc,   char*   argv[])   
  {   
  int   i,j,n;   
  double   **a,*b,*x;   
  FILE   *fp=fopen("d:\\data.txt","r");   
  fscanf(fp,"%d",&n);   
  a=new   double*[n];   
  b=new   double[n];   
  x=new   double[n];   
  for(i=0;i<n;i++)   
  {   
  a[i]=new   double[n];   
  for(j=0;j<n;j++)   
  fscanf(fp,"%lf",a[i]+j);   
  fscanf(fp,"%lf",b+i);   
  }   
  fclose(fp);   
    
  gauss(a,b,x,n);//Gauss主元消去法   
  for(i=0;i<n;i++)   
  printf("%f\t",x[i]);   
  printf("\n");   
    
  for(i=0;i<n;i++)   
  x[i]=1;   
  jacobi(a,b,x,3,10);//Jacobi迭代法   
  for(i=0;i<n;i++)   
  printf("%f\t",x[i]);   
  printf("\n");   
    
  for(i=0;i<n;i++)   
  x[i]=1;   
  gs(a,b,x,n,10);//GS迭代法   
  for(i=0;i<n;i++)   
  printf("%f\t",x[i]);   
  printf("\n");   
    
  doolittle(a,b,x,n);//Doolittle三角分解   
  for(i=0;i<n;i++)   
  printf("%f\t",x[i]);   
  printf("\n");   
    
  for(i=0;i<n;i++)   
  delete   a[i];   
  delete   a,b,x;   
  return   0;   
  }