mygs.m

解线性方程组的迭代法 C语言

function [x,iter] = mygs(A, b, x0, epsilon, iter_max)
%  Gauss-Seidel 迭代解线性方程组
%  输入：A 为系数矩阵， b为右端项，x0为初值，
%        epsilon 为相对误差精度，iter_max 为最大迭代步数
%  输出：x 为解，iter 为迭代步数

if prod(diag(A))==0
    disp('Error: diag(A) is not invertible!')
    return
end
n = length(b);
x_old = reshape(x0,n,1);
x_new = zeros(n,1);
for k = 1 : iter_max
    for i = 1 : n
        tmp = b(i) - A(i,[1:i-1])*x_new([1:i-1]) - A(i,[i+1:n])*x_old([i+1:n]);
        x_new(i) = tmp/A(i,i); 
    end
    error1 = norm(x_new - x_old);
    if error1 < epsilon, break, end
    x_old = x_new;
end
iter = k;    
x = x_new;

% END