高斯赛德尔迭代算法.cpp

高斯-赛德尔迭代法--源代码及关键源代码。包括了高斯-赛德尔矩阵算法类

/*******************************************
*            Class Matrix                   *
*******************************************/
#include <iostream.h>
#include <math.h>
#include <conio.h>
#include <iomanip.h>

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 funtion 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
/**************************
* constructor--构造函数       *
**************************/
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:系数及常数的个数 	
*********************************/
/*********************************
*    static initialization        *
*********************************/
void matrix::initialize(int &i,int &j)
{
cout<<"How many ROWS?>";
cin>>i;
cout<<"How many COLUMNS?>";
cin>>j;
//clrscr();
}
/*********************************
*       高斯-赛德尔矩阵变换函数  *
*       rearrange function           *
*********************************/
//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
/**********************************
*       iteration funtion--迭代函数     *
**********************************/
//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);
   }
}

/********************************
*        solve function--求解函数    *
********************************/
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 criterion=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 criterion is met
	iteration(lambda);
      for( i=0;i<varnum;i++)
      {
		newest[i]=variable[i];
      }
	itercount++;
   }
	while(epsilon(newest,last,varnum,criterion));
}
/*******************************************
*        Epsilon Criterion-Epsilon 精度要求   *
*******************************************/
bool matrix::epsilon(double *newest,double *last,int size,double criterion)
{
	for(int i=0;i<size;i++)
   {
      //if(it has not met the criterion)
   	if((fabs(newest[i]-last[i])/newest[i])>criterion)
      {
      //then (return 1)=(condition not met) and the loop is repeated
      return 1;
      }
   }
//criterion has been met
return 0;
}
/************************************************
*   show answer function-- 输出求解结果函数         *
*************************************************/
//'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;
}
/***********************************************************
*        stream operators--矩阵输入、输出流重载函数            *
***********************************************************/
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;
int i,j;
matrix::initialize(i,j);
matrix one(i,j);
cin>>one;
char X;
cout<<"Do you want to see the matrix?:(y/n)>";
cin>>X;
if(X=='y'){
cout<<one;}
one.rearrange();
double relax_coef;
for(char x='y';x=='y' && x!='n';)
{
cout<<endl<<"Relaxation Coefficient(1=no relaxation):>";
cin>>relax_coef;
one.solve(relax_coef);
one.show_answer();
cout<<endl<<"try another Relaxation Coefficient? :(y/n)>";
cin>>x;
}
cin.ignore(128,'\n');
cin.ignore(128,'\n');
}