Gauss-Jacobi 迭代法

Posted 2023-02-18

【中文标题】Gauss-Jacobi 迭代法

【英文标题】：Gauss-Jacobi iteration method

【发布时间】：2022-01-08 05:07:24

【问题描述】：

我正在尝试编写一个使用 求解方程组 Ax=B 的程序。

#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;

停止循环的条件是所有 x 的先前近似值减去当前近似值小于 epsilon。 似乎我按照算法做了所有事情，但是程序不起作用。 我哪里做错了？

【问题讨论】：

不要转换malloc返回的值。 ***.com/questions/605845/…

a=(double**)(malloc(n*sizeof(double))) ; 不正确。也许sizeof(double) == sizeof(double *) 这不是问题，但也许不是。你想要a = malloc(n * sizeof *a);

“不起作用”是什么意思？当然，a[i][i] 的除法是一个问题，因为没有检查它是否非零，但是如果没有对错误的更描述性解释，它真的不值得进一步研究。

谢谢。从未听说过以这种方式使用 malloc。

@WilliamPursell 要求实现此方法的对角线元素必须非零。

【参考方案1】：

虽然不是最严格的条件，但在 Jacobi 和 Gauss-Seidel 方法中保证收敛所需的通常条件是对角优势，

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

这个测试也很容易实现，作为在迭代之前运行的检查。

所有这些不等式的相对差距越大，该方法的收敛速度越快。

【讨论】：

在实践中，如果我们想解决Ax = b，我们可以将其转换为A^t Ax = A^t b，这样可以保证更好的收敛性。据我了解，高斯本人就是使用了这个技巧。 （不，我不是出生的，我读过它）。

这通常不是一个好主意，尤其是在大维度上，因为乘积矩阵的平方条件数为A。所做的是应用预处理器而不是A^t。也就是说，如果可用，矩阵的逆矩阵A 的简单近似。由于得到的矩阵更接近单位矩阵，因此通常会导致对角线优势。例如，对于稀疏矩阵，存在 ILU，即不完全 LU 分解，它保留或减少了稀疏模式。