Hi, I am trying to implement Gauss-Seidel in Matlab without using too many built in functions.

The way I am approaching it is to rewrite the system $\displaystyle Ax=b$ as

$\displaystyle x = D^{-1}(b-A_{off}\cdot x)$ where $\displaystyle A_{off}$ is the matrix $\displaystyle A$ with it's diagonal zeroed out.

Since I assume that the diagonal element of $\displaystyle A$ are non-zero, the elements in $\displaystyle D^{-1}$ are simply the reciprocals of the diagonal elements of $\displaystyle A$.

The first part of my code is just creating a strictly diagonally dominant matrix.

The loop converges to a solution, but it is not the right one, $\displaystyle Ax \neq b$. I have stared at it for quite some time now, and do not see what I have done wrong.

Code:

format long;
%-------------------------------TESTMATRIX------------------------------
n = 3;
A = zeros(n);
% The diagonal of A.
for i = 1:n
A(i,i) = 2*i;
end
% Superdiagonal
for i = 1:n-1
A(i,i+1)=i-1;
end
% Subdiagonal
for i = 2:n
A(i,i-1) = -i+1;
end
% Guessing a vector b in the system Ax=b.
b = zeros(1,n);
for i = 1:n
b(i) = i/2;
end
%-----------------------------TESTMATRIX END----------------------------
d = diag(A);
Aoff = A - diag(d);
x = zeros(1,n); % My first guess at a solution.
N = 10;
e = 10e-5;
while N >= 0
x_previous = x;
for i = 1:n
for j = 1:n
x(i) = x(i) + (Aoff(i,j)*x(j));
end
x(i) = (b(i)-x(i))/A(i,i);
end
if norm(x-x_previous,inf) < e
x
return
end
N = N - 1;
end