有限差分的 Jacobi 迭代:何时停止?

计算科学 有限差分 迭代法 C 泊松
2021-12-24 21:16:00

我实现了一个有限差分方案来求解 C 中二维网格中的泊松方程。我使用 Jacobi 迭代来求解系统。一切正常,直到我使用 while 循环检查是否该停止迭代(使用 for 循环很容易)。在我关注的笔记上写着我必须计算以下内容:

δ=max||uNewu||,

为了1<=i,j<=n. uNew作为当前的解决方案和u上一次迭代。显然它们是二维数组。

我尝试执行以下操作:

// Iterations
    double delta=(tau+1), temp_delta;
    #ifdef _OPENMP
            wt1=omp_get_wtime();
    #endif
    do {

            for(i=1;i<y-1;i++) {
                    for(j=1;j<x-1;j++) {
                            uNew[i*x+j] = 0.25 * (u[(i-1)*x+j] + u[i*x+(j+1)] + u[i*x+(j-1)] + u[(i+1)*x+j] - dx*dy*func(i,j,dx,dy));
                    }
            }

            // Boundary conditions using g(x,y)
            for(j=0;j<x;j++) {
                    uNew[j] = gunc(0,j,dx,dy);
                    uNew[(y-1)*x+j] = gunc(y-1,j,dx,dy);
            }
            for(i=0;i<y;i++) {
                    uNew[i*x] = gunc(i,0,dx,dy);
                    uNew[i*x+(x-1)] = gunc(i,x-1,dx,dy);
            }

            // Check if to terminate Jacobi iteration
            for(i=1;i<y-1;i++) {
                    for(j=1;j<x-1;j++) {
                            temp_delta = abs(uNew[i*x+j]-u[i*x+j]);
                            printf("%f ", temp_delta);
                            if (delta <= temp_delta) {
                                    delta = temp_delta;
                            }
                    }
            }

            // Update solution
            for(i=0;i<y;i++) {
                    for(j=0;j<x;j++) {
                            u[i*x+j] = uNew[i*x+j];
                    }
            }
    } while(delta > tau);

在哪里τ是公差,δ是上式的结果,temp_delta用于找到最大值和uNewu只是包含网格点处解的矩阵。问题是它不起作用。有人可以给我一个提示吗?

1个回答

对于求解系统,您不应该比较每次迭代之间的结果,而是计算残差。

如果您认为矩阵表示为标准形式:

Ay=b

然后您可以在某个迭代中定义残差(yi)

r=bAyi

然后你“只”需要根据残差向量确定停止方法,一些选择是

δ<max(|r|)

δ<sum(|r|)

关于残差的维基百科文章