-1

this is a piece of code for a simple iteration method for solving systems of linear algebraic equations:

double* iter(double** a, double* y, int n, int& iter)
{
    double* res = new double[n];
    int i, j;
    
    for (i = 0; i < n; i++)
    {
        res[i] = y[i] / a[i][i];
    }

    double eps = 0.0001;
    double* Xn = new double[n];

    do {
        iter++;

        for (i = 0; i < n; i++) {
            Xn[i] = y[i] / a[i][i];
            for (j = 0; j < n; j++) {
                if (i == j)
                    continue;
                else {
                    Xn[i] -= a[i][j] / a[i][i] * res[j];
                }
            }
        }

        bool flag = true;
        for (i = 0; i < n - 1; i++) {
            if (fabs(Xn[i] - res[i]) > eps) {
                flag = false;
                break;
            }
        }

        for (i = 0; i < n; i++) {
            res[i] = Xn[i];
        }

        if (flag)
            break;
    } while (1);

    return res;
}

and formula for it: enter image description here

but I would like to implement the seidel method.and slightly changed the code according to the formula below

for (i = 0; i < n; i++) {
            Xn[i] = y[i] / a[i][i];
            for (j = 0; j < i-1; j++) {
                    Xn[i] -= a[i][j] / a[i][i] * Xn[j];
            }
            for (j = i+1; j < n; j++){
                Xn[i] -= a[i][j] / a[i][i] * res[j];
            }
        }

enter image description here

but I'm not getting exactly what I expected:

enter image description here

I would be grateful if you could tell me where I made a mistake. thank you in advance for your answers.

enter image description here

Vs_De_S
  • 155
  • 1
  • 6
  • 2
    *I would be grateful if you could tell me where I made a mistake* -- Do a hand calculation, then [use a debugger](https://stackoverflow.com/questions/25385173/what-is-a-debugger-and-how-can-it-help-me-diagnose-problems) and see where your program diverges from the hand calculation. – PaulMcKenzie Oct 31 '21 at 08:48
  • 1
    BTW, you would be far better off using `std::vector` instead of `new double[]` and `double *`. The code you have now leaks memory. – PaulMcKenzie Oct 31 '21 at 08:51
  • in the implement of second formula: for(j = i+1; j < n; ++j) should change to for(j = i; j < n; ++j). since in your formula, j is in range [i+1, n], but in your code j is in range [i+1, n). :) – che.wang Oct 31 '21 at 09:04
  • @che.wang if I do as you said, it will turn out what is in the screenshot (I attached to the question) – Vs_De_S Oct 31 '21 at 09:10
  • @Vs_De_S an example input. so we can check it ourself. – che.wang Oct 31 '21 at 09:34
  • @che.wang I kind of managed to solve the problem, but I didn't figure out how myself. I haven't changed anything at all...but there was some other problem, haha.. – Vs_De_S Oct 31 '21 at 10:03

1 Answers1

0

Your mistake lies in the new implementation.

The first sum of the Seidel method sums up to the element before the diagonal, while your for loop goes up to two elements before the diagonal.

Instead of for(j = 0; j < i-1; j++) you should have for(j = 0; j < i; j++)

Note that Gauss Seidel method is applicable if the elements on the diagonal are non-zero.

kiszkot
  • 19
  • 2