目前位置: VCer资源中心 >>> VCer代码 >>> C++/MFC基础

[本帖已阅读2768次 分值80 回复3次] 张贴资源 发回信箱 控制面板

解高斯方程组算法(自己写的和清华教材上的比较)

提供者:barco 张贴时间:2004-05-17 23:36:27.0 出处:vcer.net 作者:不祥

解高斯方程组算法(自己写的和清华教材上的比较)(2004-05-17 23:36:27.0)


巴壳


 
级别: VCer小兵
头衔: VCer会员

经验: 462
作品: 15
分会: 华中分会
注册: 2004-05-13 22:18:24.0
登录: 2004-07-22 21:24:52.0

重点是方程组类Linequ的方法 Solve(),解法一是清华教材《C++语言程序设计》郑莉等编的算法,好像有错(反正测试计算出来试错的main()为测试码),且不怎么明白!于是自己写了个,自己命名为逐行消去法!请大家看看!呵呵写得不好请帮忙改进!

//   Linqu.h

//Begin of file Linequ.h

#ifndef LINQU_H

#define LINQU_H

#include <iostream.h>

#include <math.h>

class Matrix                   //The definition of basic class

{

public:                        //Interface for outside

 Matrix(int dim=2);         //Constructor

 ~Matrix();                 //Destructor

 void setMatrix(double* rmatr);   //Initialize Matrix

 void printM();             //Display Matrix

protected:

 int index;                 //Dimation of Matrix

 double* MatrixA;

};

class Linequ:public Matrix     //Definition of public externs class

{

public:

 Linequ(int dim=2);

 ~Linequ();

 void setLinequ(double* a,double* b);   //Initialize the equation

 void printL();                         //Display the equation

 int Solve();                             //Solution of gaosi equation

 void showX();                          //Display the result

private:

 double* sum;

 double* solu;                         //result of equation

};

#endif

//End of file Linequ.h

//  Linequ.cpp

//Begin of file Linequ.cpp

#include "Linqu.h"

//Realization of basic class

void Matrix::setMatrix(double* rmatr)         //set Matrix

{

 for(int i=0;i<index*index;i++)

  *(MatrixA+i)=rmatr[i];

}

Matrix::Matrix(int dim)                   //Constructor

{

 index=dim;                              //Give protected member value

 MatrixA=new double[index*index];        //Distribution of dynamic storage

}

Matrix::~Matrix()

{

 delete[] MatrixA;                      //Release storage

}

void Matrix::printM()                     //Display members

{

 cout<<"The Matrix is:"<<endl;

 for(int i=0;i<index;i++)

 {

  for(int j;j<index;j++)

   cout<<*(MatrixA+i*index+j)<<"  ";

  cout<<endl;

 }

}

//Realization of derive class

Linequ::Linequ(int dim):Matrix(dim)    //Constructor of class Linequ

{

 sum=new double[dim];                  //Dynamic destribution

 solu=new double[dim];

}

Linequ::~Linequ()

{

 delete[] sum;

 delete[] solu; 

}

void Linequ::setLinequ(double* a,double* b)    //Set equation

{

 setMatrix(a);                            //Call basic function

 for(int i=0;i<index;i++)

  sum[i]=b[i];

}

void Linequ::printL()               //Display equation

{

 cout<<"The Line equation is:"<<endl;

 for(int i=0;i<index;i++)

 {

  for(int j=0;j<index;j++)

   cout<<*(MatrixA+i*index+j)<<"  ";

  cout<<"      "<<sum[i]<<endl;

 }

}

void Linequ::showX()              //display result

{

 cout<<"The Result is:"<<endl;

 for(int i=0;i<index;i++)

  cout<<"X["<<i<<"]="<<solu[i]<<endl;

}

//解法 一 :全选主元高斯消去法

/*int Linequ::Solve()           //Solution of Line equation

{

 int *js,l,k,i,j,is,p,q;

 double d,t;

 js=new int[index];

 l=1;

 for(k=0;k<=index-2;k++)   //The process of expunction

 {

  d=0.0;

  js[k] = k;

  is = k;

  for(i=k;i<=index-1;i++)

   for(j=k;j<=index-1;j++)

   {

    t=fabs(MatrixA[i*index+j]);

    if(t>d){d=t;js[k]=j;is=i;}

   }

   if(d+1.0==1.0) l=0;

   else

   {

    if(js[k]!=k)

     for(i=0;i<=index-1;i++)

     {

      p=i*index+k; q=i*index+js[k];

      t=MatrixA[p];

      MatrixA[p]=MatrixA[q];

      MatrixA[q]=t;

     }

    if(is!=k)

    {

     for(j=k;j<=index-1;j++)

     {

      p=k*index+j; q=is*index+j;

      t=MatrixA[p];

      MatrixA[p]=MatrixA[q];

      MatrixA[q]=t;

     }

     t=sum[k]; sum[k]=sum[is]; sum[is]=t;

    }

   }

   if(l==0)

   {

    delete[] js; cout<<"fail"<<endl;

    return 0;

   }

   d=MatrixA[k*index+k];

   for(j=k+1;j<=index-1;j++)      //****************

   {

    p=k*index+j;

    if(d!=0)MatrixA[p]=MatrixA[p]/d;

   }

   sum[k]=sum[k]/d;

   for(i=k+1;i<=index-1;i++)        //****************

   {

    for(j=k+1;j<=index-1;j++)

    {

     p=i*index+j;

     MatrixA[p]=MatrixA[p]-MatrixA[i*index+k]*MatrixA[k*index+j];

    }

    sum[i]=sum[i]-MatrixA[i*index+k]*sum[k];

   }

 }

 d=MatrixA[(index-1)*index+index-1];

 if(fabs(d)+1.0==1.0)

 {

  delete[] js;

  cout<<"fail"<<endl;

  return 0;

 }

 solu[index-1]=sum[index-1]/d;

 for(i=index-2;i>=0;i--)

 {

  t=0.0;

  for(j=i+1;j<=index-1;j++)

   t=t+MatrixA[i*index+j]*solu[j];

  solu[i]=solu[i]-t;

 }

 js[index-1]=index-1;

 for(k=index-1;k>=0;k--)

  if(js[k]!=k)

  {

   t=solu[k]; solu[k]=solu[js[k]]; solu[js[k]]=t;

  }

 delete[] js;

 return 1;

}*/

//解法 二  逐行消去法 自创

int Linequ::Solve()

{

 int *js,k,i,j,l,flag,p,q;

 double d,t;

 js=new int[index];

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

  js[i]=0;

 for(i=0;i<index;i++)                //纪录系数矩阵各行行首连零个数

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

  {

   if(MatrixA[i*index+j]==0) js[i]++;

   else break;

  }

 do{

   flag=0;                //设置循环和判断标志

   

   for(i=0;i<index-1;i++)

   {

    if(js[i] > js[i+1])   //如果前一行的零个数比后一行多,便进行交换

    {

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

     {  //整行交换

      p=i*index+j;  q=(i+1)*index+j;

      t=MatrixA[p]; MatrixA[p]=MatrixA[q]; MatrixA[q]=t;

     }

     t=sum[i]; sum[i]=sum[i+1]; sum[i+1]=t;//对应和向量交换

     l = js[i];  js[i] = js[i+1]; js[i+1] = l;//零计数 交换

    }

   }

   

   for(i=0;i<index;i++)  //各行元素分别除以第一非零元素

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

     if(MatrixA[i*index+j]!=0)

     {

      d=MatrixA[i*index+j]; //找到第一个非零元素

      for(k=j;k<index;k++)

       MatrixA[i*index+k]=MatrixA[i*index+k]/d;

      sum[i]=sum[i]/d;

      break;  //跳出该行

     }

   

   for(i=0;i<index-1;i++)

    if(js[i]==js[i+1])  //若存在非零元素个数相等的行(肯定是相零行)flag=1

    {

     flag=1;

     break;

    }

    

   if(flag == 1) //有非零元素个数相等的行便进行相减

   {

    for(i=index-1;i>0;i--)

     if(js[i] == js[i-1])

     {

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

       MatrixA[i*index+j] -= MatrixA[(i-1)*index+j];

      sum[i] -= sum[i-1];

      js[i] = 0;

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

       if(MatrixA[i*index+j] == 0) js[i]++; //重新纪录行零个数

       else break;

     }

   }

   else if(js[index-1] >= index-1) //若不存在零数相等的行,

    flag = 2;                   //那么看最后一行是否只有一个非零或全零,

  }                                   //是这样就结束while循环

  while(flag!=2);

  

  if(js[index-1]==index)  //若最后一行全零,则无解或有无穷组解

  {

   cout<<"There infinit solutions or no solution!"<<endl;

   delete[] js;

   return 0;

  }

  else{                  //若最后一行只有一个非零,则有唯一解

   solu[index-1]=sum[index-1]; //因为最后一行一列元素为1,解等于对应的和

   for(i=index-2;i>=0;i--)

   {

    solu[i]=sum[i];       //重倒数第二行开始逐行回代

    for(j=i+1;j<index;j++)

     solu[i] -= MatrixA[i*index+j]*solu[j];

   }

   delete[] js;

   return 1;

  }

}

//End of file Linequ.cpp

/*************************************************************************

/  test.cpp

**************************************************************************/

#include "Linqu.h"

void main()

{

 double a[]=

 {0.2368,0.2471,0.2568,1.2671,

 0.1968,0.2071,1.2168,0.2271,

 0.1581,1.1675,0.1768,0.1871,

 1.1161,0.1254,0.1397,0.1490};

 double b[]={1.8471,1.7471,1.6471,1.5471};

 Linequ equl(4);

 equl.setLinequ(a,b);

 equl.printL();

 if(equl.Solve())

  equl.showX();

 else

  cout<<"Fail"<<endl;

 int t;

 cin>>t;

}

/////////////////////////////////////////////////////////////////////////

注:转载文章需注明来源:VCer.net 文章地址:http://vcer.net/2225.html

  如果你觉得VCer.net不错,而且你愿意为VCer.net捐赠一元钱,那么点击后面的捐赠按钮吧:) vcer.net捐赠

[回复该贴] [加入个人书签]
[投票结果]

A: 评分 10 0% (0 票)
B: 评分 5 0% (0 票)
C: 评分 0 0% (0 票)
D: 评分 -5 0% (0 票)
E: 评分 -10 0% (0 票)

 


re:解高斯方程组算法(自己写的和清华教材上的比较)
这个Matrix类的析构函数~Matrix()应该是为 virtual ~Matrix()吧?

kingwind 于 2004-05-20 19:23:24.0 编辑 [回复该贴]

re:解高斯方程组算法(自己写的和清华教材上的比较)
呵呵,不好意思,我直接从我的VC里COPY过来的,没注意!

barco 于 2004-05-18 11:52:08.0 编辑 [回复该贴]

re:解高斯方程组算法(自己写的和清华教材上的比较)
可是没有排版好代码呢~

bluejoe 于 2004-05-18 08:08:25.0 编辑 [回复该贴]