Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > C++ Programming with NR

Reply
 
Thread Tools Display Modes
  #1  
Old 03-12-2006, 09:00 AM
Cloud Cloud is offline
Registered User
 
Join Date: Mar 2006
Posts: 1
Gauss-Seidel method in C++ program

I must add to may program GS algorithm for soloving linear equation (i already hav gauss elimination and LU algorithm).
I make function (but it not works, please corect me :

Code:
where:
n -dimensin of main matrix
tab - main matrix
itmax - max number of iteration
prec- precisin
start - initial for xold and x vector


void GS_solve(type **tab,type *b, type *x, const int &n,double prec, int itmax, type start){
	double suma,m;
	double error=1.;
	type *xold=new type[n];
	for(int i=0;i<n;++i){ 
                xold[i]=start;
                x[i]=start;
        }
	for(int it=1;it<itmax;++it){
		if(error>prec){
			for(int i=0;i<n;++i){
				suma=0;
				for(int j=0;j<i;++j){
					suma+=tab[i][j]*x[j];
				}
				for(int j=i+1;j<n;++j){
					suma+=tab[i][j]*xold[j];
				}
				x[i]=(b[i]-suma)/tab[i][i];
			}
		}else break;
		error=0.;
		for(int k=0;k<n;++k){
			m=fabs(x[k]-xold[k]);
			if(error<m) error=m;
			xold[k]=x[k];
		}
	}
	delete [] xo;
}
PS.Sorry for my bad englisch
Reply With Quote
  #2  
Old 06-01-2008, 04:17 PM
Flying_Runner Flying_Runner is offline
Registered User
 
Join Date: Jun 2008
Posts: 1
Here's what I did for the Gauss Seidel. You can get away with only using one array for your solution Vector (X) and therefore you only have one loop, which makes your code look cleaner.

For the system Ax = b


double * x = new double [rows];
double sum = 0;

for (int l = 0; l < rows; l++)
x[l] = 0;

while(iterations--)
{
for(int i = 1; i <= rows; i++)
{
sum = b[i-1];
for(int j = 1; j <= rows; j++)
{
if( i != j )
{
sum = sum - (el[i][j] * x[j-1]);
}
}
x[i-1] = sum/el[i][i];
}
}


Background on my code: el is my two-dimensional array (Matrix A) and it is 1 based, to make my matrix class nice (or at least I think it is). However, the X and b vectors are zero based, which is why I have the j-1 and i-1 for those. Also, I know a lot of people will have the while loop condition statement check between previous values. I just gave it a set number of iterations to go, but it wouldn't be hard to change that. This also doesn't have relaxations either. That wouldn't be hard to add either. I was mainly trying to get the bare bones for the algorithm while keeping everything really clean. Hope this helps!

Last edited by Flying_Runner; 06-01-2008 at 04:18 PM. Reason: Add some explanations
Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 03:27 AM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2013, Jelsoft Enterprises Ltd.