gauss-seigel-1.cpp

线性方程解法算法

/*************************************************************
*这个程序采用高斯—赛德尔迭代法解线性方程组。
*程序设计要求：设计测试主程序
*测试线性方程：
*	10*x1+x2+2*x3=44
*	2*x1+10*x2+x3=51
*	x1+2*x2+10*x3=61
*    解：x1=3,x2=4,x3=5
**************************************************************/
#include <iostream.h>
#include <math.h>
#include <conio.h>
#include <iomanip.h>

/********************************************
*                                           *
*            Class Matrix                   *
*                                           *
********************************************/
class matrix//高斯—赛德尔矩阵算法类
{
friend void operator<<(ostream &,matrix &);
friend void operator>>(istream &,matrix &);

protected:

int row,column;

//'mat' is my actual matrix
double **mat;

//size of 'variable' always =(column-1)
int varnum;

//contains the solution set
double *variable;

//counts number of iterations
int itercount;

//performs a single iteration on all the equations
//the argument is the relaxation coefficient
void iteration(double);

//checks the allowable error
bool epsilon(double*,double*,int,double);

public:
matrix(int,int);

//small function for asking the user for the number of rows and columns
static void initialize(int&,int&);

//prepares function for reduction
void rearrange();

//solves by gauss-seigel method,
//perfomes multiple iterations until a solution is found
//the argument is the relaxation-coefficient, and this is "past" to 'iteration'
void solve(double);

//prints the solution
void show_answer();

};
//END OF CLASS MATRIX



/***************************
*      构造函数             *
***************************/
matrix::matrix(int r=1, int c=1):row(r),column(c)
{
//creates my matrix dynamically, but there are still no values in the positions
mat=new double*[row];

	for(int i=0;i<row;i++)
	{
	mat[i]=new double[column];
	}
}
//end of constructor

/**********************************
*    矩阵行列设置                   
*    ROWS:线性方程数
*    COLUMNS:系数及常数的个数 	
**********************************/
void matrix::initialize(int &i,int &j)
{
cout<<"How many ROWS?>";
cin>>i;
cout<<"How many COLUMNS?>";
cin>>j;
//clrscr();
}

/**********************************
*       高斯—赛德尔矩阵变换函数     *
**********************************/
//rearranges the system so that it can be solved by gauss-seigel method
//must always rearange before solving
void matrix::rearrange()
{

varnum=column-1;
//'variable' will contain the solution set
//they will get initialized with the first guess in the 'solve' function
variable=new double[varnum];


  	/*divides all terms by the desired variables
	(the variable which is being solved for) coefficient*/
	for(int i=0;i<row;i++)
	{
   double coefficient=mat[i][i];
		for(int j=0;j<column;j++)
   	{
   	mat[i][j]/=coefficient;
      }
	}

   //all variables except the diagonal are being brought to the other side
	for( i=0;i<row;i++)
   {
   	for(int j=0;j<column-1;j++)
      {
      mat[i][j]*=-1;
      }
   mat[i][i]=0;
   }
}

//END OF 'REARRANGE' FUNCTION


/***********************************
*       迭代函数		          *
***********************************/
//performes a single iteration on the system, using passed assumed values
void matrix::iteration(double lambda)
{

//'last' is for the relaxation equation
double last;

	for(int i=0;i<varnum;i++)
   {
   last=variable[i];
   variable[i]=0;
   	for(int j=0;j<varnum;j++)
      {
      variable[i]+=mat[i][j]*variable[j];
      }
   variable[i]+=mat[i][column-1];

   //new value after relaxation
   variable[i]=last+lambda*(variable[i]-last);//lambda是松弛因子，可以控制迭代速度
   }
}

/*********************************
*        求解函数 	         *
*********************************/
void matrix::solve(double lambda)
{
//DECLARATIONS AND INITIALIZATIONS
   //initializes first guess
	for (int i=0;i<varnum;i++)
	{
	variable[i]=0;
	}


itercount=0;

//this is the allowable error this value can be changed to suit the problem
double critereon=0.0001;

double *newest=new double[varnum];
double *last=new double[varnum];

for( i=0;i<varnum;i++)
{
newest[i]=variable[i];
}
//END OF DECLARATIONS AND INITIALIZATIONS

//MAIN BODY OF THE FUNCTION STARTS HERE-
//while('condition not met'){perform another iteration}
  	do
	{
      for(int i=0;i<varnum;i++)
   	{
      last[i]=newest[i];
      }

   //this is the most important part,
   //everything else in this loop is only to check that the critereon is met
	iteration(lambda);

      for( i=0;i<varnum;i++)
      {
		newest[i]=variable[i];
      }

	itercount++;
   }
	while(epsilon(newest,last,varnum,critereon));

}

/**************************************
*        Epsilon Critereon            *
**************************************/
bool matrix::epsilon(double *newest,double *last,int size,double critereon)
{
	for(int i=0;i<size;i++)
   {
      //if(it has not met the critereon)
   	if((fabs(newest[i]-last[i])/newest[i])>critereon)
      {
      //then (return 1)=(condition not met) and the loop is repeated
      return 1;
      }
   }

//critereon has been met
return 0;

}

/**************************************
*        输出求解结果函数		      *
**************************************/
//'solve' function must be executed before 'show_answer' function
void matrix::show_answer()
{

	for(int i=0;i<varnum;i++)
	{
	cout<<"X"<<(i+1)<<"="<<variable[i]<<endl;
	}

cout<<endl<<endl<<"Number of Iterations="<<itercount<<endl;
}

/**************************************
*        矩阵输入、输出流重载函数        *
**************************************/
void operator<<(ostream& out,matrix& m)
{
for (int i=0;i<m.row;i++)
	{
		for (int j=0;j<m.column;j++)
		{
		out<<m.mat[i][j]<<setw(10);
		}
   out<<endl<<endl;
	}
}

void operator>>(istream& in,matrix& m)
{
for(int i=0;i<m.row;i++)
	{
   	for(int j=0;j<m.column;j++)
      {
      cout<<"value of indice ["<<(i+1)<<"]["<<(j+1)<<"]";
      in>>m.mat[i][j];
      }
	}

//clrscr();
}
//end of stream operators


/***********************
*                      *
*        MAIN          *
*                      *
***********************/
void main()
{
cout
<<"This Program solves linear systems by the Gauss-Seigel Method"
<<endl
<<"Try to input the matrix so that there is Diagonal Dominance"
<<endl<<endl;


}