用MATLAB实现Gauss-Seidel迭代

用MATLAB实现Gauss-Seidel迭代

用Gauss-Seidel迭代计算如下方程组

迭代推导

设计步骤

在求解线性方程组时， 程序中输入的是线性方程组的增广矩阵，并设定送
代结果的精度要求， 在开始迭代之前则可以进行以下判断:

1. 检测该用户输入矩阵是否有误，并判断该线性方程组是否有唯一解;
2. 判断系数矩阵的对角元是否有零元素。
   由于迭代过程的本质是矩阵乘法，因此程序中直接调用 MATLAB 的矩阵乘法来完成送代，选代均以零向量为初始向量，并且在送代之前判断该送代过程是否收敛。 倘若收敛，则输出结果并打印迭代次数。

具体实现

function函数

function solution = GaussSeidel(Ab,epsilon)
%% 输入参数检查
if nargin ==1
    disp('请输入精度要求epsilon')
    return 
end

row = size(Ab,1);
col = size(Ab,2);
if ndims(Ab)~= 2 | col - row ~= 1
    disp('矩阵的大小有误，不能采用Gauss-Seidel迭代法')
    return
end

A = Ab(:,1:row);
b = Ab(:,col);

ddet = abs(det(A));
ddiag = abs(det(diag(diag(A))));
if ddet < eps | ddiag < eps
    disp('该方程的系数矩阵行列式不为0，方程组无解或有无穷解，或系数矩阵的对角线有零元，不能采用Gauss-Seidel迭代法')
    return
end

%%提取上下三角矩阵及对角矩阵
U = -triu(A,1);   %提取对角元素为0的上三角矩阵
L = -tril(A,-1);  %提取对角元素为0的下三角矩阵
D = diag(diag(A));%提取对角阵
G = (D - L)^(-1) * U;
if max(abs(eig(G))) > 1
    disp('迭代法不收敛！！！')
    return
end

%% 迭代过程
error = 10;
n1 = 0;
n = row;
start = zeros(row,1);
xk = start;
xknext = start;
while error > epsilon
    xk = xknext;
    xknext(1) = 1/A(1,1)*(b(1) - sum(A(1,2:n).* xk(2:n)'));
    for i = 2:n - 1
        Ssum1 = sum(A(i,1:i - 1) .* xk(1:i - 1)');
        Ssum2 = sum(A(i,i + 1:n) .* xknext(i + 1:n)');
        xknext(i) = 1/A(i,i) * (b(i) - Ssum1 - Ssum2);
    end
        xknext(n) = 1/A(n,n)*(b(n) - sum(A(n,1:n-1).* xknext(1:n-1)'));
        error = norm(xk - xknext);
        n1 = n1 + 1;
end
    fprintf('GS迭代次数：%d',n1)
    solution = xknext;

主函数：

Ab = [10 -2 -1 3;-2 10 -1 15;-1 -2 5 10]
GaussSeidel(Ab,1e-8)

注：function函数和主函数是同一文件夹路径下的文件