Help with Gaussian Elimination

[mcook228]

I have an assignment to write code to calculate a cubic spline. The values for my coefficients are correct, but the Gaussian elimination part of my program has me really confused. I was given code in C by my instructor, but am struggling in making it work with my code.

If anyone has any suggestions I would GREATLY appreciate it!

#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <fstream.h>

void main() {
    int n = 4;
    double x[4];
    double y[4];
    double matrix[4][4];
    double a[4][4];
    double b[4];
    int i, j, m;
    double ratio, temp;
    int count;

    for (m=0; m<n; m++) {
        for(i=0; i<n; i++) {
            a[m][i]=0.0;
        }
    }

    x[0]=1;
    y[0]=1;
    x[1]=2;
    y[1]=.5;
    x[2]=3;
    y[2]=1/3;
    x[3]=4;
    y[3]=.25;

// this portion works fine and the results check out....

    j=1;

    for(m = 1; m<n-1; m++) {
        a[m][j-1] = (x[j]-x[j-1])/6;
        a[m][j] = (x[j+1]-x[j-1])/3;
        a[m][j+1] = (x[j+1]-x[j])/6;
        b[m] = ((y[j+1]-y[j])/(x[j+1]-x[j]))-((y[j]-y[j-1])/(x[j]-x[j-1]));

        cout << a[m][j-1] << " M[" << m << "][" << j-1 << "] + ";
        cout << a[m][j] << " M[" << m << "][" << j << "] + ";
        cout << a[m][j+1] << " M[" << m << "][" << j+1 << "] ";
        cout << " = " << b[m] << endl << endl;

        j++;
    }


// Here's where the problems begin...  I just get garbage data...
    
    /* Gaussian Elimination (provided by instructor) */

    for(i=1;i<(n-1);i++){
        for(j=(i+1);j<n;j++){
            ratio = a[j][i] / a[i][i];
            cout << "ratio " << a[j][i] << "/" << a[i][i] << " = " << ratio << endl;
            for(count=i;count<n;count++) {
                a[j][count] -= (ratio * matrix[i][count]);
                cout << "a[" << j << "][" << count << "] = " << a[j][count] << endl;
            }
            b[j] -= (ratio * b[i]);
            cout << "b[" << j << "] = " << b[j] << endl;
        }
    }
  
    for (i=0;i<=n-1;i++){
        for(j=0;j<n;j++){
           cout << matrix[i][j] << endl;
        }
    }

    /* Back substitution */

    x[n-1] = b[n-1] / matrix[n-1][n-1];
    for(i=(n-2);i>=0;i--){
        temp = b[i];
        for(j=(i+1);j<n;j++){
            temp -= (matrix[i][j] * x[j]);
        }
        x[i] = temp / matrix[i][i];
    }
    cout << "Answer:" << endl;
    for(i=0;i<n;i++){
        cout << "x[" << i << "] = " << x[i] << endl;

    }
    return;
}

[Salem]

> #include <iostream.h>
These are old-style headers.
New C++ compilers should support
#include <iostream>
using namespace std; // so we can do cout rather than std::cout

> #include <math.h>
The way to include a 'C' header in new C++ is to prefix the name of the header with a 'c', so
#include <cmath>

> void main()
main returns an int.

> a[j][count] -= (ratio * matrix[count]);
As far as I can tell, matrix is uninitialised at this point, so you just introduce garbage into your calculations.

[iamthwee]

Yup Salem is right.

Your way to do Gaussian elimination is piss poor or you copied it wrong.

Here's one I did earlier. I've used a trivial example to make it clear what's going on.

Firstly, let us assume x,y,z \in \mathbb{R}

and da question...

5x  + 2y  + 1z = 17
1x  + 4y  + 2z = 16
2x  + 1y  + 4z = 11

Solution 
x=2
y=3
z=1

Then the code for guassian elimination would be as follows:

/* Program to demonstrate gaussian elimination
   on a set of linear simultaneous equations
 */

#include <iostream>
#include <cmath>

using namespace std;

const int n=3;
const double eps = 1.e-15;

/*Preliminary pivoting strategy
  Pivoting function
 */
 double pivot(double a[][n],double b[],int i)
 {
     int j=i;
     double t=0;
     
     for(int k=i; k<n; k+=1)
     {
         double aki = fabs(a[k][i]);
         if(aki>t)
         {
             t=aki;
             j=k;
         }
     }
     if(j>i)
     {
         double dummy;
         for(int L=0; L<n; L+=1)
         {
             dummy = a[i][L];
             a[i][L]= a[j][L];
             a[j][L]= dummy;
         }
         double temp = b[j];
         b[i]=b[j];
         b[j]=temp;
     }
     return a[i][i];
 }        
 
/* Forward elimination */
void triang(double a[][n],double b[])
{
    for(int i=0; i<n-1; i+=1)
    {
        double diag = pivot(a,b,i);
        if(fabs(diag)<eps)
        {
            cout<<"zero det"<<endl;
            return;
        }
        for(int j=i+1; j<n; j+=1)
        {
            double mult = a[j][i]/diag;
            for(int k = i+1; k<n; k+=1)
            {
                a[j][k]-=mult*a[i][k];
            }
            b[j]-=mult*b[i];
        }
    }
}

/*
DOT PRODUCT OF TWO VECTORS
*/
double dotProd(double u[], double v[], int k1, int k2)
{
  double sum = 0;
  for(int i = k1; i <= k2; i += 1)
  {
     sum += u[i] * v[i];
  }
  return sum;
}

/*
BACK SUBSTITUTION STEP
*/
void backSubst(double a[][n], double b[], double x[])
{
  for(int i =  n-1; i >= 0; i -= 1)
  {
    x[i] = (b[i] - dotProd(a[i], x, i + 1,  n-1))/ a[i][i];
  }
}
/*
REFINED GAUSSIAN ELIMINATION PROCEDURE
*/
void gauss(double a[][n ], double b[], double x[])
{
   triang(a, b);
   backSubst(a, b, x);
}                               

// EXAMPLE MAIN PROGRAM
int main()
{
  double a[n][n];
  double b[n];
  double x[n];

  //hard code matrix for clarity
  a[0][0] = 5.0;
  a[0][1] = 2.0;
  a[0][2] = 1.0;
  
  a[1][0] = 1.0;
  a[1][1] = 4.0;
  a[1][2] = 2.0;
  
  a[2][0] = 2.0;
  a[2][1] = 1.0;
  a[2][2] = 4.0;

  b[0] = 17.0;
  b[1] = 16.0;
  b[2] = 11.0;

  gauss(a, b, x);
  for( int i = 0; i < n; i += 1)
  {
     cout << x[i] << endl;
  }
  std::cin.get();
   return 0;
}

Get it?

Recall also that due to rounding corrections, errors propagate easily throughout the system...so you have to watch out for that.

[mcook228]

Thank you so much! I appreciate the sample code, it definitely helped clarify some things. I'm going to make some changes and try it again!

Thanks, again!

[iamthwee]

Sorry, actually something is wrong in the code I have given. That will teach me to double check the code with more than just one test example!

Rather strangely, it works ok for the one example I have given but fails with others!

Rather than spend ages debugging that I have provided you with another code snippet.

This SHOULD do the job properly...

#include <iostream>
#include <cmath>
using namespace std;

void gauss(int N, // number of unknowns
		   float A [20] [21], // coefficients and constants
		   float result[20],
		   bool& err)
// Solve system of N linear equations with N unknowns
// using Gaussian elimination with scaled partial pivoting
// First N rows and N+1 columns of A contain the system
// with right-hand sides of equations in column N+1
// err returns true if process fails; false if it is successful
// original contents of A are destroyed
// solution appears in column N
{
	int indx[20];
	float scale[20];
	float maxRatio;
	int maxIndx;
	int tmpIndx;
	float ratio;
	float sum;
    
	for (int i = 0; i < N; i++) indx[i] = i;	// index array initialization

	// determine scale factors

	for (int row = 0; row < N; row++)
	{
		scale[row] = abs(A[row][0]);
		for (int col = 1; col < N; col++)
		{
			if (abs(A[row][col]) > scale[row]) scale[row] = abs(A[row][col]);
		}
	}

// forward elimination

	for (int k = 0; k < N; k++)
	{
		// determine index of pivot row
		maxRatio = abs(A[indx[k]] [k])/scale[indx[k]];
		maxIndx = k;
		for (int i = k+1; i < N; i++)
		{
			if (abs(A[indx[i]] [k])/scale[indx[i]] > maxRatio)
			{
				maxRatio = abs(A[indx[i]] [k])/scale[indx[i]];
				maxIndx = i;
			}
		}
		if (maxRatio == 0) // no pivot available
		{
			err = true;
			return;
		}
		tmpIndx =indx[k]; indx[k]=indx[maxIndx]; indx[maxIndx] = tmpIndx;

		// use pivot row to eliminate kth variable in "lower" rows
		for (int i = k+1; i < N; i++)
		{
			ratio = -A[indx[i]] [k]/A[indx[k]] [k];
			for (int col = k; col <= N; col++)
			{
				A[indx[i]] [col] += ratio*A[indx[k]] [col];
			}
		}
	}

	// back substitution

	for (int k = N-1; k >= 0; k--)
	{
		sum = 0;
		for (int col = k+1; col < N; col++)
		{
			sum += A[indx[k]] [col] * A[indx[col]] [N];
		}
		A[indx[k]] [N] = (A[indx[k]] [N] - sum)/A[indx[k]] [k];
	}

/*
	cout << endl;
	for (int r=0; r<N; r++)
	{
		cout << indx[r];
		for (int c=0; c<=N; c++) cout<<"  " << A[r][c];
		cout << endl;
	}
*/
	for (int k = 0; k < N; k++) result[k] = A[indx[k]] [N];
}



int main()
{
	float A[20][21];
	float X[20];
	bool err;
	A[0][0] = 5; 
    A[0][1] = 2; 
    A[0][2] = 1; 
    A[0][3] = 17;
     
	A[1][0] = 1; 
    A[1][1] = 4; 
    A[1][2] = 2; 
    A[1][3] = 16; 
    
	A[2][0] = 2; 
    A[2][1] = 1; 
    A[2][2] = 4; 
    A[2][3] = 11; 
	
	gauss(3, A, X, err);
	for (int i=0; i<3; i++) cout << X[i] << "  ";
	std::cin.get();
	return 0;
}

[mcook228]

Again, thank you SO much!!! I really appreciate the sample code and your comments! ;)

[tformed]

I modified the program so i can accept user input for the x y z and b values, but how do modify it to accept the number of unknowns?

[iamthwee]

>but how do modify it to accept the number of unknowns?

Increase the size of your array, the code should be able to handle arbitary square matrix values.

[tformed]

>but how do modify it to accept the number of unknowns?

Increase the size of your array, the code should be able to handle arbitary square matrix values.

I am sure I can do it within the code, but during runtime I am not sure how.

This is what professor asks us to do,

Program to solve linear equation with n as demission.
AX=B

a11X1 + a12X2 + .......A1nXn = b1
a21X1 + a22X2 + .......A2nXn = b2
-----------------------------------------
an1X1 + an2X2 + .......AnnXn = bn

[iamthwee]

Well it's pretty simple, instead of hard coding it prompt the user for a value.

Ever heard of cin?

[tformed]

Yes I have heard of that lol, and plus i appreciate your help greatly.
What I have so far:

#include <iostream>
#include <cmath>
using namespace std;
void gauss(int N, // number of unknowns
float A [20] [21], // coefficients and constants
float result[20],
bool& err)
// Solve system of N linear equations with N unknowns
// using Gaussian elimination with scaled partial pivoting
// First N rows and N+1 columns of A contain the system
// with right-hand sides of equations in column N+1
// err returns true if process fails; false if it is successful
// original contents of A are destroyed
// solution appears in column N
{
int indx[20];
float scale[20];
float maxRatio;
int maxIndx;
int tmpIndx;
float ratio;
float sum;

for (int i = 0; i < N; i++) indx[i] = i; // index array initialization
// determine scale factors
for (int row = 0; row < N; row++)
{
scale[row] = abs(A[row][0]);
for (int col = 1; col < N; col++)
{
if (abs(A[row][col]) > scale[row]) scale[row] = abs(A[row][col]);
}
}
// forward elimination
for (int k = 0; k < N; k++)
{
// determine index of pivot row
maxRatio = abs(A[indx[k]] [k])/scale[indx[k]];
maxIndx = k;
for (int i = k+1; i < N; i++)
{
if (abs(A[indx[i]] [k])/scale[indx[i]] > maxRatio)
{
maxRatio = abs(A[indx[i]] [k])/scale[indx[i]];
maxIndx = i;
}
}
if (maxRatio == 0) // no pivot available
{
err = true;
return;
}
tmpIndx =indx[k]; indx[k]=indx[maxIndx]; indx[maxIndx] = tmpIndx;
// use pivot row to eliminate kth variable in "lower" rows
for (int i = k+1; i < N; i++)
{
ratio = -A[indx[i]] [k]/A[indx[k]] [k];
for (int col = k; col <= N; col++)
{
A[indx[i]] [col] += ratio*A[indx[k]] [col];
}
}
}
// back substitution
for (int k = N-1; k >= 0; k--)
{
sum = 0;
for (int col = k+1; col < N; col++)
{
sum += A[indx[k]] [col] * A[indx[col]] [N];
}
A[indx[k]] [N] = (A[indx[k]] [N] - sum)/A[indx[k]] [k];
}
 
for (int k = 0; k < N; k++) result[k] = A[indx[k]] [N];
}
int main()
{
float A[20][21];
float X[20];
bool err;
int N;
float x1,x2,x3,y1,y2,y3,z1,z2,z3,b1,b2,b3;
cout << "Enter X1: ";
cin>> x1;
A[0][0] = x1; 
cout << "Enter Y1: ";
cin>> y1;
A[0][1] = y1;
cout << "Enter Z1: ";
cin>> z1;
A[0][2] = z1;
cout << "Enter B1: ";
cin>> b1;
A[0][3] = b1;
//----------
cout << "Enter X2: ";
cin>> x2;
A[1][0] = x2;
cout << "Enter Y2: ";
cin>> y2;
A[1][1] = y2;
cout << "Enter Z2: ";
cin>> z2;
A[1][2] = z2;
cout << "Enter B2: ";
cin>> b2;
A[1][3] = b2; 

//-----------
cout << "Enter X3: ";
cin>> x3;
A[2][0] = x3;
cout << "Enter Y3: ";
cin>> y3;
A[2][1] = y3;
cout << "Enter Z3: ";
cin>> z3;
A[2][2] = z3;
cout << "Enter B3: ";
cin>> b3;
A[2][3] = b3; 

gauss(3, A, X, err);
for (int i=0; i<3; i++) cout << X[i] << " ";
std::cin.get();
return 0;
}

But where do I place the N as demission?

[iamthwee]

int main()
{
   float A[20][21];
   float X[20];
   bool err;
   int N;
   float x1, x2, x3, y1, y2, y3, z1, z2, z3, b1, b2, b3;
   cout << "Enter X1: ";
   cin >> x1;
   A[0][0] = x1;
   cout << "Enter Y1: ";
   cin >> y1;
   A[0][1] = y1;
   cout << "Enter Z1: ";
   cin >> z1;
   A[0][2] = z1;
   cout << "Enter B1: ";
   cin >> b1;
   A[0][3] = b1;

   cout << "Enter X2: ";
   cin >> x2;
   A[1][0] = x2;
   cout << "Enter Y2: ";
   cin >> y2;
   A[1][1] = y2;
   cout << "Enter Z2: ";
   cin >> z2;
   A[1][2] = z2;
   cout << "Enter B2: ";
   cin >> b2;
   A[1][3] = b2;

   cout << "Enter X3: ";
   cin >> x3;
   A[2][0] = x3;
   cout << "Enter Y3: ";
   cin >> y3;
   A[2][1] = y3;
   cout << "Enter Z3: ";
   cin >> z3;
   A[2][2] = z3;
   cout << "Enter B3: ";
   cin >> b3;
   A[2][3] = b3;

   gauss ( 3, A, X, err );
   for ( int i = 0; i < 3; i++ ) cout << X[i] << " ";
   std::cin.get();
   return 0;
}

Um ever heard of a nested for loop?

>But where do I place the N as demission?

Um maybe before you start entering stuff into the array?

[tformed]

I have heard of the nested for loop.

Unfortunately I do no know how to apply it with it showing the output of X1, X2, X3 etc.

I might understand the concept of it providing the amount of A[#][#] but after that I am lost once again. :(

I hope I am disturbing you in any way.

[iamthwee]

Hint:

for ( int i = 0; i < 4; i++ )
{
   for ( int j = 0; j < 4; j++ )
   {
      cin >> A[i][j];
   }
}

[tformed]

rofl, now it is asking me for the variables an infinite amount of time.

it just repeats the same statement :(

[iamthwee]

Show us the code...

[tformed]

Don't laugh, as I am a newb to this and trying to learn. And the help from the professor does not help much since English is not his primary language.

int main()
{
float A[20][21];
float X[20];
bool err;
int n;
int i,j;
float x,y,z;


for ( int i = 0; i < 4; i++ )
{
for ( int j = 0; j < 4; j++ )
{
cin >> A[i][j];
cout << "Enter X's: ";
cin >> x;
A[i][j] = x;
}
for ( int j = 0; j < 4; j++ )
{
cin >> A[i][j];
cout << "Enter Y's: ";
cin >> y;
A[i][j] = y;
}
for ( int j = 0; j < 4; j++ )
{
cin >> A[i][j];
cout << "Enter Z's: ";
cin >> z;
A[i][j] = z;
}
}
gauss ( 3, A, X, err );
for ( int i = 0; i < 3; i++ ) cout << X[i] << " ";
std::cin.get();
return 0;
}

[iamthwee]

Hint

for ( int i = 0; i < 4-1; i++ )
{
for ( int j = 0; j < 4; j++ )
{
cin >> A[j];
}
}

This is the only thing you need, just this, once only. (change the value '4' for variables)

[tformed]

Something like this:

int main()
{
float A[20][21];
float X[20];
bool err;
int i,j;
float x,y,z,b;
 
 
for ( int i = 0; i < 4; i++ )
{
for ( int j = 0; j < 4; j++ )
{
cout<< "Enter demission of n: ";
cin >> A[i][j];
A[i][j] = x;
cout << "Enter X: ";
cin >> y;
A[i][j] = y;
cout << "Enter Y: ";
cin >> z;
A[i][j] = z;
cout << "Enter Z: ";
cin >> b;
A[i][j] = b;
}
}
gauss ( 3, A, X, err );
for ( int i = 0; i < 3; i++ ) cout << X[i] << " ";
std::cin.get();
return 0;
}

now it repeats xyz over again...

[iamthwee]

Um I have a feeling I could be here for a long time. How about:

string crap = "xyzb";
  
for ( int i = 0; i < 4-1; i++ )
{
for ( int j = 0; j < 4; j++ )
{

 cout <<"enter "<<crap[j]<<i+1<<":";
 cin>>A[i][j];
 cout <<"\n";
}
cout<<"\n";
}

[tformed]

but all it does it state the x1 and increments it.

so it give me this

x1:
y2:
z3:

x1:
y2:
z3:

I do not understand how n comes into play.

[iamthwee]

but all it does it state the x1 and increments it.

so it give me this

x1:
y2:
z3:

x1:
y2:
z3:

Nope you must have changed something by accident in the code I posted.

I do not understand how n comes into play

for ( int i = 0; i < N-1; i++ )
{
for ( int j = 0; j < N; j++ )
{

cout <<"enter "<<crap[j]<<i+1<<":";
cin>>A[j];
cout <<"\n";
}
cout<<"\n";
}

[tformed]

Ok this is what I get now:

int main()
{
float A[20][21];
float X[20];
bool err;
int i,j;
int x,y,z,b;
int n;
cout<< "Enter n: ";
cin>> n;
string some = "XYZB";
for ( int i = 0; i < n-1; i++ )
{
for ( int j = 0; j < n; j++ )
{
cout <<"enter "<<some[j]<<i+1<<":";
cin>>A[j];
cout <<"\n";
}
cout<<"\n";
}

gauss ( 3, A, X, err );
for ( int i = 0; i < 3; i++ ) cout << X << " ";
std::cin.get();
return 0;
}

output is:

when I press n=2
x1 = 5
y1 =2
z1 =1
x2 = 1
y2 =4
z2 =2

1.56e013 0 -7.819

i am still confounded on what n does. I see it in the program, but what is the motive behind having n as demission. When i place a different n it gives me the correct answers but sometimes, the program crashes when i place 13

[iamthwee]

>When i place a different n it gives me the correct answers but sometimes, the program crashes when i place 13

Of course! Don't you also think N has a part to play in the below expression too?

gauss ( 3, A, X, err );
for ( int i = 0; i < 3; i++ ) cout << X[i] << " ";
std::cin.get();
return 0;

Or were you expecting me to explain that as well?

[tformed]

explaining would be great mr.iamthwee :)

I placed n in the code below, still shows same.
don't get frustrated with me please. :'(

Like this:

gauss ( 3, A, X, err );
for ( int i = 0; i < n; i++ ) cout << X[i] << " ";
std::cin.get();
return 0;

[iamthwee]

when I press n=2
x1 = 5
y1 =2
z1 =1
x2 = 1
y2 =4
z2 =2

For the above n should be 3 shouldn't it?

[iamthwee]

Like this:

gauss ( 3, A, X, err );
for ( int i = 0; i < n; i++ ) cout << X << " ";
std::cin.get();
return 0;

Hint: Something is still not right.

[tformed]

of I understand how this is working little by little, n is the number of variables within the equation including XYZB

but when the number is greater as in lets say 6, it crashes.

[iamthwee]

>but when the number is greater as in lets say 6, it crashes
Then allocate more space for the initial array!

float A [100] [100], // coefficients and constants
float result[100],

Ok look, I don't mean to be mean, but that code I've shown you should have put you at least a week in front of the rest of your class. What you are trying to do now is "making it look nice" crapola. I'm sure you have enough time to figure this out.

[On a side note]
Whether you choose to use my code example as your actual assignment is your decision, if you can't explain what it is doing to your teacher then expect a fail (this is highly likely cos you're asking me how to do the simplest of things right now!). If you're clever you'll try to understand how it's working the magic.

In any case, you need to step back and take a long hard think, cos I'm just going round in circles here. Good night! :)

[tformed]

you are right, I have bugged you too much. It's that the professor does not teach but expects great things out of us. This is our final project, yet we only made it as far as using switch statements.
I will try to make this thing work, as I must.

Again thanks for your help, and sorry for the inconvenience.