c - C 中的高斯赛德尔(特定方程求解器)

Question

很抱歉在此尝试之前发布了一个格式错误的问题。

我正在尝试让 Gauss Seidel 方法在 C 中工作，以检查它比高级解释语言(即 python)快多少，但我对获得的结果有一些疑问.

我的输入矩阵是

• 对称正定式
• & 对角优势

所以我相信它应该收敛。

问题试图解决 "Ax=b",

(其中“A”=“a[][]”，“b”=“b[]”，而“x”=“x[]”)

最终数组 'check [ ]' 是通过 'a' 和 'x' 之间的点积获得的，以查看它是否返回接近 'b' 的值。

下面的代码是完全可执行的。

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


int main(void) 
{   
int i=0,j=0;
int num=0;

double h = 1.0/(3+1);
double h2 = pow(h,2);
double w=1.5, sum=0.0;

long double x[9],y[9], check[9];
long double tol = pow(10, -10);
long double  a[9][9] = {{-4, 1, 0, 1, 0, 0, 0, 0, 0} ,
                        {1, -4, 1, 0, 1, 0, 0, 0, 0} ,
                        {0, 1, -4, 0, 0, 1, 0, 0, 0} ,
                        {1, 0, 0, -4, 1, 0, 1, 0, 0} ,
                        {0, 1, 0, 1, -4, 1, 0, 1, 0} ,
                        {0, 0, 1, 0, 1, -4, 0, 0, 1} ,
                        {0, 0, 0, 1, 0, 0, -4, 1, 0} ,
                        {0, 0, 0, 0, 1, 0, 1, -4, 1} ,
                        {0, 0, 0, 0, 0, 1, 0, 1, -4}};


long double b[9] =   {0.000000,
                      0.000000,
                      0.000000,
                      0.000000,
                      0.125000,
                      0.000000,
                      0.000000,
                      0.000000,
                      0.000000 };



for(i=0;i<9;i++){               // initialise the arrays to 0
    x[i]=0;
    y[i]=2*tol;                 
}   


for(i=0;i<9;i++){               // prints 'a' matrix, to check if its right
    for(j=0;j<9;j++){
        printf("a[%d][%d] = %LF ",i,j,a[i][j]);
    }
    printf("\n" );
}

printf("\n\n");

for(i=0;i<9;i++){               // prints b matrix 
    printf("b[%d] = %LF \n",i,b[i]);
}

do{                             // the Gauss seidel Solver
    for(i =0;i<9;i++){
        for(j=0; j<9; j++){
            if(i!=j){
                sum += a[i][j]*x[j];
            }
            x[i] = (w/a[i][i])* (b[i] - sum + a[i][i]*x[i]) + (1-w)*x[i];   
        }
    }
num=num+1;
}while (fabs(y[i]-x[i])>tol);

printf("\n\n\n"); 

for(i=0;i<9;i++){               // Prints the solution X
    printf("x[%d] = %LF \n",i,x[i]);
}

printf("\n\n");

for(i=0;i<9;i++){               // Ititialises matrix 
        check[i]=0;
}

for (i = 0; i < 9; i++){        // performs a dot product of
                                // 'a' and 'x', to see if we get
                                // 'b' as the result 
    for(j = 0; j< 9; j++){
        check[i]+= a[i][j] * x[j];
    }   
    check[i] = check[i]/h2;     // the 'b' matrix was multiplied by h2, 
                                // hence to get it back, a division is necessary
    printf("check[%d] = %LF\n",i,check[i]);
}

printf("num=%d\n",num );
return 0;

}

我得到的输出即“x”是:

x[0] = -0.000000 
x[1] = -0.000000 
x[2] = -0.000000 
x[3] = -0.000000 
x[4] = -0.421875 
x[5] = -0.791016 
x[6] = -1.423828 
x[7] = -3.816650 
x[8] = -11.702087 

我得到的“检查”输出是:

check[0] = 0.000000
check[1] = -4.500000
check[2] = -5.625000
check[3] = -14.625000
check[4] = -10.968750
check[5] = -42.328125
check[6] = 17.156250
check[7] = 18.421875
check[8] = 212.343750

理想情况下，如果一切正常，check[4] 应该输出 2(原因在代码中输出 'check' 时的注释中给出)，并且 check 的每个其他元素都应该为 0。

有什么建议吗？

Answer 1

在开始下一行之前，

sum 应该在 for 循环内重新初始化为 0，等式不正确。您在 python 实现中使用的等式假设添加了 a[i][i]*x[i] 来生成完整的点积，他们使用了 numpy获取点积而不是循环，因此他们没有机会执行 i != j。此外，我不确定该实现中的方程是高斯赛德尔方法，由于 w 和 (1 - w)，它看起来更像是连续松弛。无论如何，这是我修改后的解决方案。我使用错误 |Ax - b| 检查收敛性< 所有条目的 tol。 for循环被分成两部分作为一个小的优化。 a[i][i] * x[i] 被添加到 sum 以获得错误中 (Ax)_i 的当前值查看。

int converged;                                                                  
do {                                                                
  converged = 1;                                                                    
  for (i = 0; i < 9; i++) {                                                         
    sum = 0;                                                                       
    for (j = 0; j < i; j++) {
      sum += a[i][j] * x[j];
    }                                 
    for (j = i + 1; j < 9; j++) {
      sum += a[i][j] * x[j];
    }                             
    x[i] += w * ((b[i] - sum) / a[i][i] - x[i]);                                   
    if (fabs(sum + a[i][i] * x[i] - b[i]) > tol) {
      converged = 0;
    }                    
  }                                                                      
} while (!converged); 

给出输出:

x[0] = -0.007812 
x[1] = -0.015625 
x[2] = -0.007812 
x[3] = -0.015625 
x[4] = -0.046875 
x[5] = -0.015625 
x[6] = -0.007812 
x[7] = -0.015625 
x[8] = -0.007813 


check[0] = 0.000000
check[1] = -0.000000
check[2] = -0.000000
check[3] = -0.000000
check[4] = 2.000000
check[5] = 0.000000
check[6] = -0.000000
check[7] = 0.000000
check[8] = 0.000000
num=31