Вычислительные методы линейной алгебры. Горбаченко В.И - 65 стр.

%   x0 - Начальное приближение
%ВОЗВРАЩАЕМЫЕ ЗНАЧЕНИЯ
%   x - вектор решения
%   flag - признак достижения решения (flag = 1, если достигнута точность tol
%          в пределах числа итераций maxit; flag = 0, если за число итераций
%          maxit точность tol не достигнута)
%   relres - достигнутая относительная норма невязки
%   iter - Число проделанных итераций
%   resvec - Вектор полученных невязок
x = x0;
flag = 0;
N=length(b);
nb=norm(b,inf); %Норма вектора правой части
r = b - A * x; % Расчет невязки начального приближения
nr = norm (r,inf); % Норма невязки
resvec (1) = nr; % результирующий вектор невязок
for iter = 1 : maxit
    % Уточнение решения
     xk=x;
     for j=1:N
         x(j)=(b(j)-A(j,[1:j-1,j+1:N])*xk([1:j-1,j+1:N]))/A(j,j);
     end
     r = b - A * x; % Расчет невязки нового приближения
     nr = norm (r,inf); % Норма невязки
     resvec (iter+1) = nr; % результирующий вектор невязок
     relres = nr / nb; % Достигнутая относительная норма невязки
     if relres <= tol % Проверка условия окончания процесса
         flag = 1;
         break;
     end
end

    Функция seidel, реализующая метод Зейделя, имеет вид:
function [x,flag,relres,iter,resvec] = seidel(A,b,tol,maxit,x0)
% Метод Зейделя
%ПЕРЕДАВАЕМЫЕ ЗНАЧЕНИЯ
%   A,b – Матрица и вектор правой части (соответственно)
%   tol - требуемая точность (относительно нормы невязки)
%   maxit - Максимальное число итераций
%   x0 - Начальное приближение
%ВОЗВРАЩАЕМЫЕ ЗНАЧЕНИЯ
%   x - вектор решения
%   flag - признак достижения решения (flag = 1, если достигнута точность tol
%          в пределах числа итераций maxit; flag = 0, если за число итераций
%          maxit точность tol не достигнута)
%   relres - достигнутая относительная норма невязки
%   iter - Число проделанных итераций
%   resvec - Вектор полученных невязок
x = x0;
flag = 0;
N=length(b);
nb=norm(b,inf); % Норма вектора правой части
r = b - A * x; %Расчет невязки начального приближения
nr = norm (r,inf); %Норма невязки
resvec (1) = nr; % результирующий вектор невязок
for iter = 1 : maxit
    % Уточнение решения
    for j=1:N
        x(j)=(b(j)-A(j,[1:j-1,j+1:N])*x([1:j-1,j+1:N]))/A(j,j);
    end
    r = b - A * x; %Расчет невязки нового приближения
    nr = norm (r,inf); % Норма невязки
    resvec (iter+1) = nr; % результирующий вектор невязок

                                                                           65