主页 / 问题 / 1249613
Accepted
kontsev_
Asked:2022-02-28 12:55:02 +0000 UTC2022-02-28 12:55:02 +0000 UTC

使用 Jacobi 方法求解 SLE

  • 772

任务是使用 Jacobi 方法找到一个近似解，但是程序显示为零并且不继续迭代，我不明白如何修复此错误。

#include <stdio.h>
double diff(int n, double x1[n], double x2[n]) {
    double s = 0;
    int i;
        for( i = 0; i < n; i++) {
            s += (x2[i]-x1[i]) * (x2[i]-x1[i]);
        }
        return s;
}
double jacobi(int n, double x1[n], double x2[n]) {
    int i, j;
    double s = 0, a[n][n], f[n];
        for(i = 0; i < n; i++) {
            for(j = 0; j < i-1; j++) {
                s += x1[j]*a[i][j];
            }
            for(j = i+1; j < n; j++) {
                s += x1[j]*a[i][j];
            }
            x2[i] = (f[i] - s) / a[i][i];
        }
}
int main () {
    int n, i, j, k;
    double eps;
    printf ("Add N:\n");
    scanf ("%d", &n);
    printf ("Add epsilon:\n");
    scanf("%lf", &eps);
    double a[n][n], f[n], x1[n], x2[n];
        for(i = 0; i < n; i++) {
            for(j = 0; j < n; j++) {
                scanf("%lf", &a[i][j]);
            }
        }
        for(i = 0; i < n; i++) {
            scanf ("%lf", &f[i]);
        }
        for(i = 0; i < n; i++) {
            x2[i] = 0;
        }
        while(diff(n, x1, x2) > eps*eps) {
            for(i = 0; i < n; i++) {
                x1[i] = x2[i];
            }
            jacobi(n, x1, x2);
        }
        for(i = 0; i < n; i++)
            printf ("%lf", x2[i]);
}
c

1 个回答

  1. Best Answer

    您正在将矩阵注入到a[n][n]函数中的本地数组中main()。然后你调用你的函数jacobi()，它定义了另一个名为 a[n][n] 的数组，你甚至不用初始化它就可以使用它。因此，你得到零（任何垃圾都可以在最后一个矩阵中）。因此，即使您将工作矩阵设为全局数组，该函数jacobi()中同名的本地数组也会破坏一切。您需要删除最后一个数组，并且需要将工作矩阵作为参数从main(). 也就是说，您需要更改：1）a[n][n]从函数中删除数组jacobi()，2）通过将数组添加到参数来更改此函数的标头a[n][n]，以及 3）更正其调用。

    UPD：该程序不起作用，因为向量 f[n] 和另一个错误出现了相同的错误 - 您没有重置s矩阵每一行的总和（对于每个新的 x2[i]）。

    并且你还需要在main()( while(diff(n, x1, x2) > eps*eps))中的主循环加上后置条件do ... while()，否则你不会进入它。

    纠正错误后，使用 Jacobi 方法的函数应如下所示：

    void jacobi(int n, double x1[n], double x2[n], double g[n], double m[n][n]) {
    double s = 0;
        for(int i = 0; i < n; i++) {
            s=0;
            for(int j = 0; j < i-1; j++) {
                s += x1[j]*m[i][j];
            }
            for(int j = i+1; j < n; j++) {
                s += x1[j]*m[i][j];
            }
            x2[i] = (g[i] - s) / m[i][i];
        } 
}

    我会让这个功能更紧凑：

    void jacobi(int n, double x1[n], double x2[n], double g[n], double m[n][n]) 
{
    double row[n];
    
    for(int i = 0; i < n; i++)
    {
        row[i] = g[i];
        for(int j = 0; j < n; j++)
            row[i] -= (i != j) ? m[i][j] * x1[j] : 0.0;             
        x2[i] = row[i] / m[i][i];
    }
}

    检查，这两个选项似乎都可以正常工作。并且不要忘记 Jacobi 方法不会对每个矩阵都收敛。

    • 1

相关问题