用MATLAB实现Gauss-Seidel迭代
用Gauss-Seidel迭代计算如下方程组
迭代推导
设计步骤
在求解线性方程组时, 程序中输入的是线性方程组的增广矩阵,并设定送
代结果的精度要求, 在开始迭代之前则可以进行以下判断:
- 检测该用户输入矩阵是否有误,并判断该线性方程组是否有唯一解;
- 判断系数矩阵的对角元是否有零元素。
由于迭代过程的本质是矩阵乘法,因此程序中直接调用 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函数和主函数是同一文件夹路径下的文件