2

I made a program to solve linear equations using Gauss Jordan elimination method. The program is working correctly but in some cases instead of giving the answer as 0, it returns a very small value.

#include<iostream>
#include<iomanip>
#include<cassert>
#define N 10
using namespace std;

//printing out the array
void print(float x[][N], int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<=n;j++)
            cout << setprecision(5) << setw(15) <<x[i][j];
        cout << endl;
    }
    cout << endl;
}

//to normalise the leading entries to 1
void normalize(float x[][N], int n, int i, int j){
    float fac = x[i][j];
    for(int k=0;k<=n;k++)
        x[i][k] /= fac;
}

//check if the leading entry is a zero
bool chk_zero(float x[][N], int n, int i){
    int j, k, c{0};
    if(x[i][i]==0){
    for(j=i;j<n-1;j++){
        c=1;
        while(x[j+c][i]==0 && i+c<n){
            c++;
            if(i+c==n-1)
                assert(i+c==n-1 && "Equation has no solution");
                return false;
            }

        for(k=0;k<=n;k++){
            swap(x[i][k], x[i+c][k]);
            }
        return true;

        }
    }
    return true;

}

//Gauss Jordan elimination method
void GaussJordan(float x[][N], int n){
    int i, j, k, c;
    float rat;

    for(i=0;i<n;i++){
        //not taking the zero in pivot column case
        chk_zero(x, n, i);
        normalize(x, n, i, i);

        for(j=0;j<n;j++){
            if (i != j){
               float fac{x[j][i]};
               for(k=0;k<=n;k++)
                   x[j][k] = x[j][k]-fac*x[i][k];
            }
        }
    }
}



int main(){
    float arr[][N] = { {0, 5, 1, 2},
                       {2, 11, 5, 3},
                       {1, 0, 0, 0} };

    int n = sizeof(arr)/sizeof(arr[0]);

    print(arr, n);

    GaussJordan(arr, n);

    cout << n << endl;

    print(arr, n);
    return 0;
}

the output I get is:

 1              0              0     2.3842e-08
 0              1              0            0.5
-0             -0              1           -0.5

The output I should get is:

 1              0              0              0
 0              1              0            0.5
-0             -0              1           -0.5

The value 2.3842e-08 should be a zero. Is it because of the precision of the floating point in C++ ?

If so, what should I do in order to round such low values to 0 without loss of data ?

Also, why is there a "-0" instead of a 0.

Circuit_Breaker0.7
  • 125
  • 1
  • 5
  • Have you tried to step through the code statement by statement in a debugger while monitoring all variables and their values? – Some programmer dude Sep 15 '21 at 16:07
  • I also suggest you reindent your code consistently (decent editor have functionality to do it for you) as that will make it easier to see some possible bugs. Other possible bugs include statements which will never be reached. – Some programmer dude Sep 15 '21 at 16:12
  • The program works fine in almost all cases, but sometimes instead of returning 0, it prints a very small value close to 0. The only error I encountered was a value tending to zero but not zero. – Circuit_Breaker0.7 Sep 15 '21 at 16:17
  • 1
    *what should I do in order to round such low values to 0 without loss of data* - You can't do anything. Floating point math is inherently inexact. What you can do is to improve precision by using some tricks (like partial pivoting). – Evg Sep 15 '21 at 16:43

2 Answers2

1

The value 2.3842e-08 should be a zero. Is it because of the precision of the floating point in C++ ?

Yes, it is.

If so, what should I do in order to round such low values to 0 without loss of data ?

You can use an epsilon value to decide whether a float value is considered to be 0:

const float epsilon = 0.000001; // Depends on the use case of your application
if (std::abs(some_float) < epsilon) {
    // Treat some_float as 0
}

However, this approach has some cons. Read here for more.

Also, why is there a "-0" instead of a 0.

This is because of how float numbers are represented. Most C++ compilers use the IEEE 754 standard to implement floating point arithmetics.

Few things to add (unrelated to the question):

  • You should avoid using using namespace std as it pollutes the global namespace. Especially for large codebases.
  • Use std::array<> instead of raw arrays as they are more safe.
Zakk
  • 1,935
  • 1
  • 6
  • 17
  • I understood what you said, but how do other "libraries/programs" for linear algebra overcome this problem ? The epsilon solution is not foolproof, or there are other algorithms to evaluate the matrix ? – Circuit_Breaker0.7 Sep 16 '21 at 07:51
  • @Circuit_Breaker0.7 Take a look at https://www.ttmath.org/ , and https://en.wikipedia.org/wiki/List_of_C%2B%2B_multiple_precision_arithmetic_libraries if you want more choices. – Zakk Sep 16 '21 at 09:29
1

The value 2.3842e-08 should be a zero. Is it because of the precision of the floating point in C++ ?

No. It is a feature of float point number representation in binary number system according IEEE-754 standard. It is implemented in hardware and is used by almost all the languages. All the floating point calculations are performed with errrors so you should set right rounding mode for FPU. See <fenv.h> for details.

Bwana The Master
  • 121
  • 2