C++ my Gauss–Seidel method doesn’t work wery well

I want to solve system of equations (3×3) in C++ useing Gauss–Seidel method. I read some of pages about this: Wiki Page2 and Page3. By pattern oneself on those pages I wrote this code:

#include <iostream>
#include <cmath>
#define epsilon 0.001

int main()
{
    //declaration and initialization of variables
    long double a[3][3], b[3], x[3], xo[3];

    a[0][0] = 3;   a[0][1] = -6.5; a[0][2] = 3.33; b[0] = 7;
    a[1][0] = 4.7; a[1][1] = 8;    a[1][2] = -17;  b[1] = 5;
    a[2][0] = 6;   a[2][1] = 2;    a[2][2] = 8;    b[2] = -5.2;
    //By calculator: x=0,5063353908 y=-1,215147148 z=-0,725964756

    int n = 3;

    //initialization of x array
    for (int i = 0; i < n; i++)
    {
        x[i] = b[i] / a[i][i];
    }

    //start Gauss-Seidel method
    do
    {
        for (int i = 0; i < n; i++)
        {
            //xo is useing to chceck accuracy
            xo[i] = x[i];
            for (int j = 0; j < n; j++)
            {
                if (i != j) x[i] -= (a[i][j] * x[j]) / a[i][i];
            }
        }
      //chceck the accuracy
    } while (fabs(x[0] - xo[0]) < epsilon &&
             fabs(x[1] - xo[1]) < epsilon &&
             fabs(x[2] - xo[2]) < epsilon);

    std::cout << "x = " << x[0] << std::endl;
    std::cout << "y = " << x[1] << std::endl;
    std::cout << "z = " << x[2] << std::endl;

    return 0;
}

This code doesn’t work very well – the result of the program isn’t the same like solving by calculator. I really don’t have an idea what is going on.

Answer

There are several issues in your code:

  1. The accuracy checking is not correct
  2. In calculation of the solution, a term b[i] is missing
  3. The matrix is not diagonal dominant, far from it: the algorithm cannot converge

In addition, to avoid getting nan when the algorithm doesn’t converge, I have added a condition to limit the number of iterations.

I have also added a final check, i.e. calculation of ax, to see if a good result is effectively obtained.

#include <iostream>
#include <cmath>
#define epsilon 0.001

int main()
{
    //declaration and initialization of variables
    long double a[3][3], b[3], x[3], xo[3];

    a[0][0] = 3;   a[0][1] = -6.5; a[0][2] = 3.33; b[0] = 7;
    a[1][0] = 4.7; a[1][1] = 8;    a[1][2] = -17;  b[1] = 5;
    a[2][0] = 6;   a[2][1] = 2;    a[2][2] = 8;    b[2] = -5.2;
    //By calculator: x=0,5063353908 y=-1,215147148 z=-0,725964756
    
    a[0][0] = 12;       // to make matrix diadonal dominant
    a[1][2] = -1;

    int n = 3;

    //initialization of x array
    for (int i = 0; i < n; i++) {
        x[i] = b[i] / a[i][i];
    }

    //start Gauss-Seidel method
    int iter = 0;
    int n_iter_max = 20;
    do {
        for (int i = 0; i < n; i++) {
            //xo is used to check accuracy
            xo[i] = x[i];
            x[i] = b[i];
            for (int j = 0; j < n; j++) {
                if (i != j) x[i] -= (a[i][j] * x[j]);
            }
            x[i] /= a[i][i];
        }
        iter++;
      //check the accuracy
    } while ((std::abs(x[0] - xo[0]) > epsilon ||
             std::abs(x[1] - xo[1]) > epsilon ||
             std::abs(x[2] - xo[2]) > epsilon) && (iter <= n_iter_max));

    std::cout << "x = " << x[0] << std::endl;
    std::cout << "y = " << x[1] << std::endl;
    std::cout << "z = " << x[2] << std::endl;
    
//  Final check

    double b_cal[3];
    for (int i = 0; i < n; ++i) {
        b_cal[i] = 0.0;
        for (int j = 0; j < n; j++) {
            b_cal[i] += a[i][j] * x[j];
        }
    }
    std::cout << "A*x = ";
    for (int i = 0; i < n; ++i) {
        std::cout << b_cal[i] << " (" << b_cal[i] - b[i] << ")  ";
    }
    std::cout << "n";

    return 0;
}