/*第6题 高斯-赛德尔迭代法--源代码及关键源代码注解如下：*/
/*******************************************
* 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');
}
