Gauss-Jacobi iteration method

Question

I'm trying to write a programm that solves system of equations Ax=B using Gauss-Jacobi iteration method.

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

int main(void) {
    double **a, *b, *x, *f, eps = 1.e-2, c;  
    int n = 3, m = 3, i, j, bool = 1, d = 3;
    /* printf("n=") ; scanf("%d", &n);
       printf("m=") ; scanf("%d", &n) */
   
   
    a =malloc(n * sizeof *a);
    for (i = 0; i < n; i++) 
        a[i] = (double*)malloc(m * sizeof(double));

    b = malloc(m * sizeof *b);
    x = malloc(m * sizeof *x) ;  
    f = malloc(m * sizeof *f) ;
    for (i = 0; i < n; i++) {
        for (j = 0; j < m; j++) { 
            printf("a[%d][%d]=", i, j); 
            scanf("%le", &a[i][j]); 
            if(fabs(a[i][i])<1.e-10) return 0 ; 
        }

        printf("\n") ;
    }
        
    printf("\n") ;
        
    for (i = 0; i < n; i++) {
        for (j = 0; j < m; j++) { 
            printf("a[%d][%d]=%le  ", i, j, a[i][j]); 
        }
         
        printf("\n") ;
    }
    
    for (j = 0; j < m; j++) { 
        printf("x[%d]=", j); 
        scanf("%le", &x[j]); 
    } //intial guess
    
    printf("\n") ;
    
    for (j = 0; j < m; j++) { 
        printf("b[%d]=", j); 
        scanf("%le", &b[j]); 
    }
    
    printf("\n")  ;

    while (1) {
        bool = 0;
        for (i = 0; i < n; i++) {
            c = 0.0;
            for (j = 0; j < m; j++) 
                if (j != i) 
                    c += a[i][j] * x[j];  
            f[i] = (b[i] - c) / a[i][i];
        }
       
        for (i = 0; i < m; i++)  
            if (fabs(f[i] - x[i]) > eps) 
                bool = 1;
       
        if (bool == 1) 
            for (i = 0; i < m; i++) 
                x[i] = f[i];
        else if (bool == 0) 
            break;
    }

    for (j = 0; j < m; j++) 
        printf("%le\n", f[j]);

    return 0;
}

The condition of stoping the loop is that previous approximation minus current approximation for all x is less than epsilon. It seems like i did everything according to algorithm,but the programm doesn't work. Where did i make a mistake?

Don't cast the value returned by `malloc`. https://stackoverflow.com/questions/605845/do-i-cast-the-result-of-malloc — William Pursell, Dec 01 '21 at 15:50
`a=(double**)(malloc(n*sizeof(double))) ;` is not correct. Perhaps `sizeof(double)` == `sizeof(double *)` and it's not a problem, but maybe not. You want `a = malloc(n * sizeof *a);` — William Pursell, Dec 01 '21 at 15:52
What does "doesn't work" mean? Certainly the division by `a[i][i]` is a concern, since there's no check that it is non-zero, but without a more descriptive explanation of the error it's not really worth looking further. — William Pursell, Dec 01 '21 at 16:01
@WilliamPursell There is a requirement that to implement this method diagonal elements must be nonzero. — Nickie, Dec 01 '21 at 16:05
If you have a requirement imposed on the input, it is incumbent on the program to verify that the requirement is met. This includes checking the values returned by every call to `scanf` and ensuring that all diagonal elements are non-zero. — William Pursell, Dec 01 '21 at 16:25
Note that the Jacobi and Gauss-Seidel methods are only guaranteed to work if the matrix is diagonally dominant. Check your test data for it, or better, implement a test for it after the input. — Lutz Lehmann, Dec 01 '21 at 16:35
@LutzLehmann found out that it's works when satisfies all the conditions.Thanks you very much.If you want,you can write it as an answer and i will accept it. — Nickie, Dec 01 '21 at 17:10
@YvesDaoust output was nan or inf.I didn't know about diagonal dominance condition.It works,when satisfes conditions. — Nickie, Dec 02 '21 at 10:14
@Nickie: this is a fixed-point method (`x' = Px + Q`), hence it must fulfill a convergence condition (i.e. be contracting). — , Dec 02 '21 at 10:28

score 3 · Accepted Answer · answered Dec 01 '21 at 18:45

3

While not the most strict condition, the usual condition requiered to guarantee convergence in the Jacobi and Gauss-Seidel methods is diagonal dominance,

abs(a[i][i]) > sum( abs(a[i][j]), j=0...n-1, j!=i)

This test is also easy to implement as a check to run before the iteration.

The larger the relative gap in all these inequalities, the faster the convergence of the method.

answered Dec 01 '21 at 18:45

Lutz Lehmann

25,219
2
22
51

In practice, if we want to solve `Ax = b`, we can transform it with `A^t Ax = A^t b`, which garanties a better convergence. In my understanding, Gauss himself used this trick. (No I was not born, I read it). – Damien Dec 02 '21 at 10:35
1

This in general is not a good idea, especially in large dimensions, as the product matrix has the squared condition number of `A`. What is done is to apply a pre-conditioner instead of `A^t`. That is, if available, an easy approximation of the inverse of the matrix `A`. As the resulting matrix is then closer to the identity matrix, diagonal dominance often results. For example for sparse matrices there exists the ILU, the incomplete LU factorization, that preserves or reduces the sparseness pattern. – Lutz Lehmann Dec 02 '21 at 10:41

Gauss-Jacobi iteration method

1 Answers1