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

релаксации. Как и следовало ожидать, самым эффективным из рассмотрен-
ных методов оказался метод сопряженных градиентов, потребовавший всего
6 итераций.
    Программа sam3_02.m демонстрирует решение СЛАУ с матрицей, полу-
ченной в результате конечно-разностной аппроксимации оператора Лапласа
в квадрате, методом сопряженных градиентов (функция con_grad). Взята
область размером n = 22 узловых точек. Размерность матрицы в этом случае

равна ( n − 2 ) = 400 . Причем, с помощью функции sparse матрица преобра-
              2



зована в разреженный формат. Вектор правой части состоит из единиц.
% Пример решения СЛАУ с разреженной матрицей СЛАУ с матрицей, полученной в
% результате конечно-разностной аппроксимации оператора Лапласа в квадрате,
% явным методом сопряженных градиентов (двучленная формула).
% Вектор правой части состоит из единиц.
clear % очистка рабочего пространства
clc % очистка командного окна
close all % Закрытие всех окон
n=22; % Число узловых точек
% матрица, полученная в результате конечно-разностной аппроксимации оператора
Лапласа в квадрате
A=delsq(numgrid('S',n));
AS=sparse(A);% AS - матрица в разреженном формате хранения
b=ones(length(A),1); % Правая часть из единиц
x0=zeros(length(A),1); % Нулевое начальное приближение
tol = 1e-6; % Точность
maxit = 1000; % Максимальное число итераций
[x,flag,relres,iter,resvec,alpha,beta] = con_grad(A,b,tol,maxit,x0);
if flag == 0 % Решение не достигнуто
    disp ('Решение не достигнуто')
    fprintf(1, '\nДостигнута относительная норма невязки %e за %u итера-
ций\n', relres, iter)
end
if flag == 1 % Решение достигнуто
    disp ('Решение достигнуто')
    fprintf(1, '\nДостигнута относительная норма невязки %e за %u итера-
ций\n', relres, iter)
    figure;
    subplot(3,1,1);
    semilogy (0:iter, resvec / norm (b,inf)) % График изменения относительной
нормы невязки
    grid on
    title('Метод сопряженных градиентов')% заголовок окна
    xlabel('Номер итерации') % надпись по оси X
    ylabel('Относ. норма невязки')% надпись по оси Y
    % График alpha
    subplot(3,1,2);
    plot(1:iter,alpha);
    grid on
    xlabel('Номер итерации') % надпись по оси X
    ylabel('alpha')% надпись по оси Y
    % График beta
    subplot(3,1,3);
    plot(1:iter-1,beta);
    xlabel('Номер итерации') % надпись по оси X

                                                                           71