How to write the calculation code for gauss seidel and SOR method ?

[LindseyCLS]

Recently, I am given an assignment on writing a c++ code for both gauss seidel and sor method. I want to use the same calculation function for both so I need to include OMEGA in the calculation part where gauss seidel(omega=1)and SOR (omega=0). But i tried for a week but still the program won't run probably.

#include <stdio.h>
#include <iostream>
#include <cmath>
#include<iomanip>
using namespace std;

void gaussSOR(){
int i,j,count=0,flag=0;
cout.setf(ios::fixed);
int A[10][10]={{-2,1,0,0,0,0,0,0,0,0},{1,-2,1,0,0,0,0,0,0,0},{0,1,-2,1,0,0,0,0,0,0},{0,0,1,-2,1,0,0,0,0,0},{0,0,0,1,-2,1,0,0,0,0},{0,0,0,0,1,-2,1,0,0,0},{0,0,0,0,0,1,-2,1,0,0},{0,0,0,0,0,0,1,-2,1,0},{0,0,0,0,0,0,0,1,-2,1},{0,0,0,0,0,0,0,0,1,-2}};
int B[10][1]={{10},{20},{30},{40},{50},{60},{70},{80},{90},{100}};
double x[i],OMEGA,y[i],eps=0.0001;


do{
for(i=0;i<10;i++){
y[i]=x[i];
x[i]=B[10][1];
for(j=0;j<10;j++){
if(j!=i){
x[i]=x[i]-A[i][j]*x[j];
}

}
x[i]=x[i]/A[i][i];
if (abs(x[i]-y[i])<=eps) //Compare the new value with the last value
flag++;
cout<<x[i]<<setw(18);
}
cout<<"\n";
count++;
}while(flag<10);

cout<<"\n The solution is as follows:\n";
for (i=0;i<10;i++)
cout<<"x"<<i<<" = "<<x[i]<<endl;

}



int main(){
int i,j,k;
int A[10][10]={{-2,1,0,0,0,0,0,0,0,0},{1,-2,1,0,0,0,0,0,0,0},{0,1,-2,1,0,0,0,0,0,0},{0,0,1,-2,1,0,0,0,0,0},{0,0,0,1,-2,1,0,0,0,0},{0,0,0,0,1,-2,1,0,0,0},{0,0,0,0,0,1,-2,1,0,0},{0,0,0,0,0,0,1,-2,1,0},{0,0,0,0,0,0,0,1,-2,1},{0,0,0,0,0,0,0,0,1,-2}};
int B[10][1]={{10},{20},{30},{40},{50},{60},{70},{80},{90},{100}};
double x[i],OMEGA;
cout<<"This is a program to calculate the 10x10 tridiagonal linear system (Ax=B).";
cout<<"\nA is a 10x10 tridiagonal matrix given as : \n";
for(i=0;i<10; i++){
for ( j=0;j<10;j++){
cout<< A[i][j]<<" ";
}
cout<<endl;
}
cout<<"\n\n B matrix is given as :\n";
for( i=0;i<10;i++){
for ( j=0;j<1;j++){
cout<<B[i][j]<<" ";
}
cout<<endl;
}
cout<<"\n\n x is the 10x1 matrics of variables that need to solve.\n";
cout<<"\n\n Please enter the value for OMEGA.";
cout<<"\n Let OMEGA equal to 1 if you want to use GAUSS-SEIDEL method.";
cout<<"\nFor SOR method, please enter the weightage value for the OMEGA. ";
cout<<"\n**noted that the weightage value should larger than 1 but smaller than 2.";
cout<<"\n OMEGA= ";
cin>>OMEGA;

if (OMEGA==1){
cout<<"You had chosen GAUSS SEIDEL method.";
gaussSOR();
}
else {
cout<<"You had chosen SOR method.";
gaussSOR();

}
return 0;
}

[lastchance]

@LindseyCLS,

In future, please put your code in code tags (first item in the format menu). Then it will be possible to see the indentation and structure of your code.

Get basic Gauss-Seidel working first, then worry about over-relaxing it.

On average, half of your x[j] won't have been initialised before you run the line
x[i]=x[i]-A[i][j]*x[j];
I suggest that before you start iterating they are explicitly set to 0.

It should be
x[i]=B[i][1];
NOT
x[i]=B[10][1];
B[10][1] is actually out of bounds anyway. In fact it would make far more sense to have B[] as a one-dimensional array.

You don't reset flag to 0 at the start of each iteration.

You should pass arrays A and B as arguments, not try to declare them separately in main() and your gaussSOR() routine.

Declaration of arrays x and y is wrong:
double x[i],y[i];
as i hasn't been set here.

Repeat: fix basic Gauss-Seidel first; worry about over-relaxing it later - that would be a trivial change of
x[i]=x[i]/A[i][i];
to
x[i]=(1-OMEGA) * y[i] + OMEGA * x[i]/A[i][i];
OMEGA should also be passed as an argument of function gaussSOR().

[LindseyCLS]

thanks for the reply! For the B array, since it is B[10][1] , if i declare it as one dimensional array, will it not appear properly ?? If not , then I should write B as B[n] or B[10] or B[1] ?
Sorry for a lot of questions , because I'm really weak in programming.

[lastchance]

You can declare it as a one-dimensional array
double B[10];

The corresponding value for x[i] at the start of each iteration then becomes
x[i] = B[i];

Your equation system is tri-diagonal, BTW. There are better solvers than G-S for tri-diagonal systems.

[LindseyCLS]

thanks for the advice ! Too bad this is what my lecturer wants . He wants gs and sor method to solve this. I will try my best to re-do my program! Later , I really hop that you can help me check and point out my mistakes .

[LindseyCLS]

For the array passing as arguements,
i tried
void gaussSOR(int A[][],int B[][])
but it seems like not so correct

[lastchance]

Include the dimensions in the statement. Also, neither A or B are int, and we have already suggested that B should be a one-dimensional array.

[LindseyCLS]

#include <stdio.h>
#include <iostream>
#include <cmath>
#include<iomanip>
using namespace std;

//function prototype
void gaussSOR(int A[10][10], double B[10]);

void gaussSOR(int A[10][10], double B[10]){
int i,j,count=0,flag=0;
cout.setf(ios::fixed);
double x[i],OMEGA,y[i],eps=0.0001;
do{

for(i=0;i<10;i++){
y[i]=x[i];
x[i]=B[i];
for(j=0;j<10;j++){
if(j!=i){
x[i]=x[i]-A[i][j]*x[j];
}

}
x[i]=x[i]/A[i][i];
if (abs(x[i]-y[i])<=eps) //Compare the new value with the last value
flag++;
cout<<x[i]<<setw(18);
}
cout<<"\n";
count++;
}while(flag<10);

cout<<"\n The solution is as follows:\n";
for (i=0;i<10;i++)
cout<<"x"<<i<<" = "<<x[i]<<endl;

}



int main(){
int i,j,k;
int A[10][10]={{-2,1,0,0,0,0,0,0,0,0},{1,-2,1,0,0,0,0,0,0,0},{0,1,-2,1,0,0,0,0,0,0},{0,0,1,-2,1,0,0,0,0,0},{0,0,0,1,-2,1,0,0,0,0},{0,0,0,0,1,-2,1,0,0,0},{0,0,0,0,0,1,-2,1,0,0},{0,0,0,0,0,0,1,-2,1,0},{0,0,0,0,0,0,0,1,-2,1},{0,0,0,0,0,0,0,0,1,-2}};
double B[10]={10,20,30,40,50,60,70,80,90,100};
double x[i],OMEGA;
cout<<"This is a program to calculate the 10x10 tridiagonal linear system (Ax=B).";
cout<<"\nA is a 10x10 tridiagonal matrix given as : \n";
for(i=0;i<10; i++){
for ( j=0;j<10;j++){
cout<< A[i][j]<<" ";
}
cout<<endl;
}
cout<<"\n\n B matrix is given as :\n";
for( i=0;i<10;i++){
cout<<B[i]<<" ";
}
cout<<endl;
cout<<"\n\n x is the 10x1 matrics of variables that need to solve.\n";
//declare initial guesses as zero
for (i=0;i<10;i++){
x[i]=0;
}
cout<<"\n\n Please enter the value for OMEGA.";
cout<<"\n Let OMEGA equal to 1 if you want to use GAUSS-SEIDEL method.";
cout<<"\nFor SOR method, please enter the weightage value for the OMEGA. ";
cout<<"\n**noted that the weightage value should larger than 1 but smaller than 2.";
cout<<"\n OMEGA= ";
cin>>OMEGA;

if (OMEGA==1){
cout<<"You had chosen GAUSS SEIDEL method.";
gaussSOR(A,B);
}
else {
cout<<"You had chosen SOR method.";
gaussSOR(A,B);

}
return 0;
}

[LindseyCLS]

my latest updates. I think the passing arrays as argument part still got mistakes

[lastchance]

@LindseyCLS,

Please use code tags (first element in the format menu). It is impossible to see the structure of your code as it is at the moment - look at the posted version!

You have ignored most of what I said in previous posts.


double x[i],OMEGA,y[i];
This is (a) illegal in standard C++ (because you can't declare variable size arrays of this type) and (b) i has no value at this point anyway. You don't need a whole array y[].


matrix A[][] should be double, not int.


x[i]=x[i]-A[i][j]*x[j];
Most of those x[j] will not have been initialised before they are used on the first iteration. You should have (e.g.) zeroed them before starting the iteration.


You still aren't re-setting flag to 0 at the start of each iteration. This isn't a terribly good way of checking convergence anyway.


You should call gaussSOR in the same way, irrespective of the value of OMEGA. You don't need an if...else block. Also, you should pass OMEGA as an argument of gaussSOR as well as the arrays.


I suggest that you declare array size as a global variable at the start:
const int SIZE = 10;
Then you can have the following function prototype for your solver:
void gaussSOR( double A[SIZE][SIZE], double B[SIZE], double X[SIZE], double OMEGA );
You will need to declare arrays of these sizes in main() and initialise A[][] and B[]. You should also do all your output in main().


This is what I got for your problem. (I reduced the tolerance to 1e-6 rather than 1e-4).

Converged to tolerance 1e-06 with OMEGA = 1.57 in 41 iterations.
Solution:
              X             AX              B
       -200.000         10.000         10.000
       -390.000         20.000         20.000
       -560.000         30.000         30.000
       -700.000         40.000         40.000
       -800.000         50.000         50.000
       -850.000         60.000         60.000
       -840.000         70.000         70.000
       -760.000         80.000         80.000
       -600.000         90.000         90.000
       -350.000        100.000        100.000

[LindseyCLS]

i thought i set the flag to zero already or I set it in a wrong way ?

int i,j,count=0,flag=0;

[lastchance]

Please read what I wrote.

You set flag to 0 at the start. However, you didn't RESET it to 0 for subsequent iterations of your do loop (i.e. just after the do{ statement). So flag just keeps growing and you will incorrectly believe the system is converged after only a few loops.

Counting the number of elements whose change is small enough on one iteration is less common than simply finding the maximum change on that iteration and comparing it with eps. Again, the variable used to find the maximum needs to be reset at the start of the do loop.

[LindseyCLS]

I did read every word you write and I already tried my best to fix the problem. But I am just a beginner so I still got a lot to learn. Sorry if my dumbness in c++ offended you and I really appreciate everything you taught me.
Back to the flag, how can I reset it to zero in the do loop ?
I tried to find in my books and google and I still can't find a way to do it .

[lastchance]

Think of the order in which things are arranged. Sometimes it may be better to type in the loop structure without filling it in.

   int flag;     // it is irrelevant whether you initialise it here
    // ...
   do {
      flag = 0;      // Must be reset to 0 at start of loop
      // ... update all x[i] and, where necessary, flag
   } while ( flag < 10 );  

Note that 10 is a "magic number" - you may want to change the size of your array. It would be far better declaring const int SIZE = 10; at the start of your program and testing as while ( flag < SIZE ); This would also allow you to pass arrays as arguments much more easily.

[LindseyCLS]

#include <stdio.h>
#include <iostream>
#include <cmath>
#include<iomanip>
using namespace std;

//Global constants
const int SIZE=10;

//function prototype
void gaussSOR( double A[SIZE][SIZE], double B[SIZE], double x[SIZE], double OMEGA  );

int main(){
	int i,j;
	double A[SIZE][SIZE]={{-2,1,0,0,0,0,0,0,0,0},{1,-2,1,0,0,0,0,0,0,0},{0,1,-2,1,0,0,0,0,0,0},{0,0,1,-2,1,0,0,0,0,0},{0,0,0,1,-2,1,0,0,0,0},{0,0,0,0,1,-2,1,0,0,0},{0,0,0,0,0,1,-2,1,0,0},{0,0,0,0,0,0,1,-2,1,0},{0,0,0,0,0,0,0,1,-2,1},{0,0,0,0,0,0,0,0,1,-2}};
	double B[SIZE]={10,20,30,40,50,60,70,80,90,100};
	double x[SIZE],OMEGA;
	cout<<"This is a program to calculate the 10x10 tridiagonal linear system (Ax=B).";
	cout<<"\nA is a 10x10 tridiagonal matrix given as : \n";
		for(i=0;i<SIZE; i++){
			for ( j=0;j<SIZE;j++){
			cout<< A[i][j]<<" ";
		}
		cout<<endl;
		}
	cout<<"\n\n B matrix is given as :\n";
	for( i=0;i<SIZE;i++){
		cout<<B[i]<<" ";
				}
		cout<<endl;
	cout<<"\n\n x is the 10x1 matrics of variables that need to solve.\n"; 
	//declare initial guesses as zero
	for (i=0;i<SIZE;i++){
			x[i]=0;
	}
cout<<"\n\n Please enter the value for OMEGA.";
cout<<"\n Let OMEGA equal to 1 if you want to use GAUSS-SEIDEL method.";
cout<<"\nFor SOR method, please enter the weightage value for the OMEGA. ";
cout<<"\n**noted that the weightage value should larger than 1 but smaller than 2.";
cout<<"\n OMEGA= ";
cin>>OMEGA;

//calling the iterations function
gaussSOR(A,B,x,OMEGA );

return 0; 
}

void gaussSOR( double A[SIZE][SIZE], double B[SIZE], double x[SIZE], double OMEGA ){
int i,j,count=0,flag;
cout.setf(ios::fixed);
double y,eps=0.000001;
	do{
	flag=0;
	//declare the initial guesses as zero
	for ( i=0;i<SIZE;i++){
		x[i]=0;
	}
		for(i=0;i<SIZE;i++){
		y=x[i];
		x[i]=B[i];
			for(j=0;j<SIZE;j++){
					if(j!=i){
						x[i]=x[i]-A[i][j]*x[j];
						}

			}
			x[i]=x[i]/A[i][i];
			if (abs(x[i]-y)<=eps) //Compare the new value with the last value
			flag++;
			cout<<x[i]<<setw(18);
			}
	cout<<"\n";
	count++; 
}while(flag<SIZE); 

cout<<"\n The solution is as follows:\n";
for (i=0;i<10;i++)
cout<<"x"<<i<<" = "<<x[i]<<endl; 

}

[LindseyCLS]

my latest updates.
1. I'm not sure whether I correct in calling the function
2.for the x[ i ] which is the initial guesses I declared it as zero, but from your last post you say to zero the x [ j ] before iterations. So , i'm quite confuse here whether i should zeroed x[i] or x[ j ] or both ?

[lastchance]

WOW! Almost there!

(1) Your function call is fine.

(2) Take your lines 56-58

for ( i=0;i<SIZE;i++){
		x[i]=0;
	}

and put them BEFORE the do loop. (Actually, you have set them in main() anyway, so you could simply remove them.) If you leave them in the loop then you would keep re-setting x[] to 0 and you would loop for ever.

For my sanity, please comment out your output lines cout<<x[i]<<setw(18); and cout<<"\n"; inside the do-loop. (They are fine outside the loop.) It takes sufficient iterations that you don't want to see the whole solution vector on every pass of the do-loop. Just output count (the number of iterations) and your final solution before you return from the routine. Also, unless you want to see a lot of ".999999" in your solution you may want to add
cout << setprecision(4);
before outputting your final solution.

This then works!!!


To include your over-relaxation parameter OMEGA, simply change your line
x[i]=x[i]/A[i][i];
to
x[i] = ( 1.0 - OMEGA ) * y + OMEGA * x[i] / A[i][i];

Your indentation is mildly inconsistent, which makes it quite hard to read your code in places.

[LindseyCLS]

Thanks !! I corrected my code and finally I got the solutions ! But just one more question, is is possible for me to calculate how many iterations in total for both of the method ? Just to see which method is faster to calculate the solution?
Also , I will try my best to fix the identation problem !

[lastchance]

Hello @LindseyCLS.

First of all, well done! Especially for surviving my sharpness of reply!

You have a variable count which is correctly tracking the number of loops. I would print out its value before (or even with) your output line
cout<<"\n The solution is as follows:\n";


The line change to incorporate varying OMEGA is in my previous post. When I tried your problem I found that the optimal value (in terms of number of iteration loops) for this particular system of equations was about 1.57

However, this optimal value is very case-specific. It would be different for every matrix and my own observation is that it tends to reduce toward 1.0 for large matrices, so doing away with any advantages gained. For nonlinear problems (my own field is CFD) we actually tend to stabilise solutions by employing under-relaxation (OMEGA < 1), a stable numerical solution being difficult to achieve otherwise. For your linear system you don't need this.


Some other thoughts. If your matrix is going to be tri-diagonal then you might care to investigate a tri-diagonal matrix algorithm (TDMA) solver, which would be direct and much faster than Gauss-Seidel. For systems with other diagonals (often arising from 2-d problems) another useful technique is line-Gauss-Seidel ("alternating-direction-implicit", ADI).

[LindseyCLS]

Thanks for the extra knowledge about the TDMA , I will try to google and learn about it !
So back to the coding , if I want to print out the number of iterations , I need to print out the number of count first right before
cout<<"\n The solution is as follows:\n";

Does I take your word correctly ??