For backward and forward elimination I used
% Now use a vector y to solve 'Ly=b'
for j=1:z
for k=1:j-1
b(j)=b(j)-L(j,k)*b(k);
end;
b(j) =b(j)/L(j,j);
end;
% Now we use this y to solve Rx = y
x = zeros( z, 1 );
for i=z:-1:1
x(i) = ( b(i) - R(i, :)*x )/R(i, i);
end
I put that into your code, and it works ;)
For helping you with pivot strategies it would be helpful to know what strategie you want/have to use as there are various versions. One simple version would be to just swap rows such that the diagonal element $a_{ii} \neq 0$ for all $i = 1, \dots, z$.
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.
Best Answer
Crout: