Метод Гаусса-Зейделя в MATLAB

Question

Я пытаюсь реализовать метод Гаусса-Зейделя в MATLAB. Но в моем коде есть две основные ошибки, и я не мог их исправить:

1. Мой код очень хорошо сходится на маленьких матрицах, но никогда не сходится на больших матрицах.

2. Код делает избыточные итерации. Как я могу предотвратить избыточные итерации?

Метод Гаусса-Зайделя в Википедии.

N=5;
A=rand(N,N);
b=rand(N,1);
x = zeros(N,1);
sum = 0;
xold = x;
tic
for n_iter=1:1000
    for i = 1:N
        for j = 1:N
            if (j ~= i)
                sum = sum + (A(i,j)/A(i,i)) * xold(j);
            else
                continue;
            end
        end
        x(i) = -sum + b(i)/A(i,i);
        sum = 0;
    end
    if(abs(x(i)-xold(j))<0.001)
        break;
    end
    xold = x;
end
gs_time=toc;
prompt1='Gauss-Seidel Method Time';
prompt2='x Matrix';
disp(prompt2);
disp(x);
disp(prompt1);
disp(gs_time);

Answer 1

Во-первых, общность. В Гаусс-Зейделя и методы Якобов применяются только к диагонально доминирующим матрицам, а не родовой случайные. Итак, чтобы получить правильные тестовые примеры, вам необходимо конструктивно обеспечить выполнение этого условия, например, через

A = rand(N,N)+N*eye(N)

или похожие.

В противном случае метод будет расходиться к бесконечности в некоторых или всех компонентах.

---

Теперь о некоторых других странностях в вашей реализации. Что значит

if(abs(x(i)-xold(j))<0.001)

подлый? Обратите внимание, что эта инструкция находится вне циклов, гдеi а также jявляются переменными итерации, поэтому потенциально значения индекса не определены. По инерции они случайно оба будут иметь значениеN, так что этот критерий имеет хоть какой-то смысл.

То, что вы хотите проверить, - это некоторая норма разности векторов в целом, поэтому используйте sum(abs(x-xold))/N или max(abs(x-xold)). С правой стороны вы можете захотеть произвести умножение с той же конструкцией нормы, применяемой кx так что тест проводится на относительную ошибку с учетом масштаба проблемы.

---

Кроме того, вы можете сократить / упростить внутренний цикл

for i = 1:N
    sum = b(i);
    for j = 1:N
        if (j ~= i)
            sum = sum - A(i,j) * xold(j);
        end
    end
    x(i) = sum/A(i,i);
end

или даже короче, используя языковые особенности Matlab

for i = 1:N
    J = [1:(i-1) (i+1):N];
    x(i) = ( b(i) - A(i,J)*xold(J) )/A(i,i);
end   

Answer 2

%Gauss-seidal method for three equations
clc;
x1=0;
x2=0;
x3=0;
m=input('Enter number of iteration');
for i=1:1:m
    x1(i+1)=(-0.01-0.52*x2(i)-x3(i))/0.3
    x2(i+1)=0.67-1.9*x3(i)-0.5*x1(i+1)
    x3(i+1)=(0.44-0.1*x1(i+1)-0.3*x2(i+1))/0.5
    er1=abs((x1(i+1)-x1(i))/x1(i+1))*100
    er2=abs((x2(i+1)-x2(i))/x2(i+1))*100
    er3=abs((x3(i+1)-x3(i))/x3(i+1))*100

    if er1<=0.01
       er2<=0.01
       er3<=0.01
        break;
    end
end