c++实现高斯塞德尔迭代法的收敛判断问题

Question

最近在用c++写一些东西，其中要实现一个高斯塞德尔迭代函数，目前关于是否收敛的判断我个人不是很确定(就是下面程序的err判断部分)，我在网上找到了一些资料，但是感觉并不是很可信(目前下面的写法就是参考网上的例子)

其中我用到了eigen矩阵库。

std::vector<double> GaussSeidelSolution(SparseMatrix<int> A,std::vector <int> b,float tolerance,int max_iter, MatrixXf x0){
    
    int i,k,m,n,iterNumber;
    
    int length = (int)b.size();
    
    std::vector<double> Answer;
    std::vector<double> LastAnswer;
    
    double err;
    
    for(i=0;i<x0.rows();i++){
        Answer.push_back((double)x0(i,0));
    }
    
    for(iterNumber=0;iterNumber<max_iter;iterNumber++){
        LastAnswer=Answer;
        for(k=0;k<length;k++){
            double sum=0;
            for(m=0;m<length;m++){
                if(m==k)continue;
                sum+=A.coeff(k,m)*Answer.at(m);
            }
            Answer[k]=(b[k]-sum)/A.coeff(k,k);
        }
        //这一块判断收敛目前不是很确定
        err = 0;
        for(n=0;n<length;n++){
            err = fabs(Answer[n]-LastAnswer[n]) > err ? fabs(Answer[n]-LastAnswer[n]):err;
        }
        if(err < tolerance){
            std::cout<<"err < tolerance at "<<i<<std::endl;
            return Answer;
        }
    }
    
    return Answer;
    
}

希望有相关经验的认人士或者相关同学能够指点一二或者讨论一下，谢谢。

Answer 1

高斯塞德尔迭代的计算，最后是否收敛的判断依据，参考wiki的例程
Gauss–Seidel method

import numpy as np

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array([[10., -1., 2., 0.],
              [-1., 11., -1., 3.],
              [2., -1., 10., -1.],
              [0.0, 3., -1., 8.]])
# initialize the RHS vector
b = np.array([6., 25., -11., 15.])

# prints the system
print("System:")
for i in range(A.shape[0]):
    row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
    print(" + ".join(row), "=", b[i])
print()

x = np.zeros_like(b)
for it_count in range(ITERATION_LIMIT):
    print("Current solution:", x)
    x_new = np.zeros_like(x)

    for i in range(A.shape[0]):
        s1 = np.dot(A[i, :i], x_new[:i])
        s2 = np.dot(A[i, i + 1:], x[i + 1:])
        x_new[i] = (b[i] - s1 - s2) / A[i, i]

    if np.allclose(x, x_new, rtol=1e-8):   # 在这里判断收敛
        break

    x = x_new

print("Solution:")
print(x)
error = np.dot(A, x) - b
print("Error:")
print(error)

其中判断收敛的是这句if np.allclose(x, x_new, rtol=1e-8)
这里numpy.allclose函数的作用是【如果两个数组在元素级别在容差内相等，则返回True。】

题主的程序中，err部分的代码，把Answer与LastAnswer对应位置元素之差的最大值，找出来，然后与容差tolerance比较，小于容差则成功。

所以，题主程序关于收敛的判断是OK的～