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:

- The accuracy checking is not correct
- In calculation of the solution, a term
`b[i]`

is missing - 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; }