# Gauss Seidel iteration in matlab

I've posted this question before for crout factorization. Now, I need help with Gauss-Seidel iteration.

Write a program that takes a value for n and solves for x using the following method:

Gauss-Seidel iteration starting with $x_{0} = 0$ and terminating when the residual is less than $10^8$ in $\infty$ norm. The program should output the $\infty$ norm of the residual of your computed solution and the number of iterations used.

I am not posting any codes because I don't have one. I need help with figuring out the codes. Thanks.

And, if I need p (permutation matrix), would it be the same as identity matrix in this case?

Gauss-Seidel iteration uses a decomposition of a matrix into a lower triangular matrix $L$ and a strictly upper triangular matrix $U$, $A = L + U$. Let's begin by defining our function and forming $A$, $L$, and $U$ based on a user's choice of $n$

function X = MyGaussSeidelExample(n)
A = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1);  % We don't actually need to form this explicitly!
L = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1);
U = -diag(ones(n-1,1),1)

b = 1/(n+1)^4*(1:n)'.^2; b(1) = b(1) + 1; b(n) = b(n) + 6; % Form b according to problem statement.
% (1:n) gives the vector [1 2 3 ... n]. (1:n)' makes it a column vector.
% (1:n)'.^2 squares each element. Pre-multiplying by 1/(n+1)^4 gives us the desired form.
% Manually add 1 to the first entry and 6 to the last.
x = zeros(n,1);

X = x;
end


This sets up your problem.

Now, the iteration rule is $Lx_{k+1} = b-Ux_k$. So, naively, we can solve $x_{k+1} = L^{-1}(b-Ux_k)$ at each step. Note, this won't end up being the best way to do it, but let's start with it.

function X = MyGaussSeidelExample(n)
A = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1);  % We don't actually need to form this explicitly!
L = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1);
U = -diag(ones(n-1,1),1);

b = 1/(n+1)^4*(1:n)'.^2; b(1) = b(1) + 1; b(n) = b(n) + 6; % Form b according to problem statement.
% (1:n) gives the vector [1 2 3 ... n]. (1:n)' makes it a column vector.
% (1:n)'.^2 squares each element. Pre-multiplying by 1/(n+1)^4 gives us the desired form.
% Manually add 1 to the first entry and 6 to the last.
x = zeros(n,1);
X = ones(n,1);

tol = 10e-8; % I assume you want an error tolerance to be small.
epsilon = 10*tol;
while (epsilon > tol)
x = inv(L)*(b-U*x);
epsilon = max(abs(X-x));
X = x;
end

total_error = max(abs(A\b - X));  % Compare the error to MATLAB's built in solver

end


Computing this in this manner using the inverse of $L$ is still very crappy. However, if you read the commentary here, it should be easy enough to directly compute each element of $x_{k+1}$ at each step without forming $L^{-1}$ directly. I'll leave that as an exercise for you. I have created the framework you need.

• Thank you so much and since you answered, I think I can ask you another question if you don't mind. If it were a Crout factorization, what does the input of this sentence "The program should output the $\infty$ norm of the residual of your computed solution and the number of iterations used" would be? Thanks. – user136422 Apr 23 '14 at 0:15
• I mean, just look up the definitions of those terms. It's not hard. It's actually really easy. If you can't figure out how to count iterations or find maxima, you need to take a more elementary programming course. – Emily Apr 23 '14 at 4:46
• Yes. I already figured it out. Thanks again for your help. – user136422 Apr 23 '14 at 4:47