Gauss Seidel iteration in matlab

Question

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?

Answer 1

Ok, let's start with the method.

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.