Метод простой итерации для решения СЛАУ в MATLAB

[sergey zhukov]

Попробовал запрограммировать метод Гаусса-Зейделя для итеративного решения СЛАУ в матлабе. В самой математике ничего сложного нет, все в общем работает. Но мне кажется, что я как-то коряво это сделал в плане работы с матрицами. Может, число циклов как-то можно сократить и использовать какие-то матричные операции? Что тут можно улучшить? Опыт программирования у меня так себе:

Используется синтаксис Matlab M

clear
clc
% Метод Гаусса-Зейделя для решения СЛАУ. Метод Гаусса-Зейделя является модификацией метода Якоби.
% Метод Якоби работает с вектором решения СЛАУ X, как с единым элементом.
% Из значения X(i) на i-ом шаге получается следующее значение X(i+1) простой итерацией, т.е. X(i+1)=C*X(i)+F
% На самом деле это можно улучшить, т.к. внутри каждого i-ого шага вычисление
% компонент вектора X происходит последовательно, при этом разумно уже
% вычисленные компоненты на i+1-ом шаге использовать для вычисления
% оставшихся компонент на этом же i+1-ом шаге.

    % Решим каку-нибудь случайную СЛАУ Ax=B:
    % Задаем размер матрицы:
N = 100;
    % Создаем единичную диагональную матрицу:
A1 = eye(N);
    % Создаем случайную матрицу с малыми компонентами. Сходимость не гарантируется, если эти коэффициенты недостаточно малы:
A2 = 1/12*randn(N);
    % Итоговая матрица СЛАУ:
A = A1+A2;
    % Создаем случайный вектор-столбец правой части:
B = randn(N,1);
    % Случайный вектор начального приближения:
X=10*randn(N,1);
    % Требуемая точность решения (порог разности длин вектора X двух последовательных итераций):
err=0.01;
    % Максимальное число итераций:
iter=100;

    % Для того, чтобы размер матрицы/вектора не менялся в цикле расчета, ее полезно сразу определить постоянного максимально возможного размера:
    % Матрица эволюции решения:
Z = zeros(iter,N);
    % Вектор эволюции точности решения:
V(1:iter) = 1;
    % Матрица С:
C = zeros(N,N);
    % Матрица F:
F = zeros(N);

    % Расчет матриц для итеративной схемы X(i+1)=C*X(i)+F:
for i=1:N
    
    
    
        % Расчет матрицы Сij=-Aij/Aii:
    for j=1:N
        if i~=j
            C(i,j)=-A(i,j)/A(i,i);
        end
    end
    
    
    
        % Расчет матрицы F(i)=Bi/Aii:
    F(i,1)=(B(i,1)/A(i,i));
end

    % Итеративный цикл, который останавливается либо по достижению предела числа итераций, либо по достижении заданной точности:
i=1;
while and(iter~=0,V(i)>err)
    
    iter=iter-1;
    
    i=i+1;
    
    
            % цикл последовательного вычисления вектора переменных на i-ом шаге:
        for j=1:N
            X(j,1)=C(j,:)*X+F(j,1);
        end
        
        
            % Матрица Z сохраняет историю итераций X (для визуализации):
        Z(i,:)=X;
        
        
            % Длина вектора ошибки (величина отклонения предыдущего решения
            % от последующего):
        if i~=1
            V(i)=norm(Z(i,:)-Z(i-1,:));
        end
end

    % Нарисуем ход решения:
subplot(2,1,1);
plot(Z);
subplot(2,1,2);
plot(V);
    % Вывод числа итераций:
disp('Число итераций до получения решения')
disp(i)
 

[GAA]

Пусть имеется система уравнений $Ax = b$. Eё приводят к виду
$x =Cx +f$, где $c_{ij} = -a_{ij}/a_{ii}$, если $i \ne j$ и $c_{ii} = 0$; $f_i = b_i/a_{ii}$.
В методе простых итераций $x^{(k+1)} = Cx^{(k)} +f$.
В методе Зейделя $x^{(k+1)} = C_1x^{(k+1)}+ C_2x^{k}+f$, где $C_1$ — нижняя треугольная часть матрицы $C$, $C_1$ — верхняя треугольная. [В файле по ссылке опечатка в матричной записи метода Зейделя]

Для реализации можно:
1. Получить вектор диагональных элементов исходной матрицы
с = diag(A);
2. Найти матрицу С
A = diag(c) – A;
А затем поделить строки матрицы на соответствующие диагональные элементы
for i=1:3
A(i, :) = A(i, :)/c(i);
end
3. Найти вектор f при помощи поэлементного деления
b = b./c

Дальше в методе простых итераций в цикле выполняем итерации.
x = x_new;
x_new = A*x + b;

В методе Зйделя реализовать итерации проще при помощи однократного цикла
x_old = x;
for i=1:n
x(i)= A(i, :)*x +b
end
Пример решения системы из книги Демидович Б.П., Марон И.А. Основы вычислительной математики, 1966.

Используется синтаксис Matlab M

n = 3;
eps = 0.001;
A = [10, 1, 1; 2, 10, 1; 2, 2, 10];
b = [12, 13, 14]';
x = [1.2, 0, 0]';
aa = diag(A);
A = diag(aa) - A;
b = b./aa;
for i=1:n
 A(i,:) = A(i, :)/aa(i);   
end
 
x_old = inf*ones(n, 1);
while max(abs(x - x_old)) > eps
 x_old = x;   
 for i=1:n
  x(i) = A(i, :)*x + b(i);   
 end    
end 
disp(x);
 

Результат
1.0000
1.0000
1.0000

Т.е. если в методе простой итерации легко реализовать итерации через матричное умножение, то в методе Зейделя использование операций над массивами позволяет во время итераций уйти от двойного цикла, но однократный цикл остаётся.

В общем, не видно смысла уходить от цикла в итерациях методом Зейделя. Так что у Вас всё по смыслу. Просто для скорости такие вещи надо писать на низком уровне и вызывать уже из Matlab "внешнюю" функцию.

[Ende]

Тема перемещена из форума «Программирование» в форум «Околонаучный софт»
Причина переноса: Матлаб сюда.