Basic iterative schemes in matlab

Hello,

I am becoming acquainted with solving basic linear problems in matlab iteratively. So, I apologize for this very basic question.

I am writing a program to solve Ax=f, for a randomly generated matrix (say, of size 6).

I want to extract to break A into M and N, where A=M-N.

I have written a code for this, and I would like to determine the number of iterations that are required before the solution converges. However, my solution does not converge! It blows up to infinity. I have checked to make sure the matrix is non-singular, and I think I have the iterative scheme correct. I would be very grateful if someone could please check my work and let me know where I might be going wrong. Here is my code:

scale=6; 

%A is a random matrix 
A = randn(scale,scale);
if det(A)==0
    fprintf('Matrix is singular');
    %break
else
    fprintf('Matrix is non singular');
end

n=size(A,1);
%need to construct M and subtract N from it

%M contains diagonal elements, here a tridiagonal;

%N contains off diagonal elements of A

M=tril(triu(A,-1),1); %tri diagonal
%tril(triu(A,-2),2) %penta diagonal
N=-(M-A);

%B is a random name for a check to ensure that A=M-N
B=M-N; %Check that we infact get back matrix A

%Solve A*x=f
f=randn(length(M),1);
x=zeros(length(M),1); %starting vector
k=1; %iteration parameter
r=f-A*x(:,k); 
tol=1e-3; %tolernace of the convergence, dictated by the 2-norm of r

while norm(r)>=tol
    x(:,k+1)=M\-N*(x(:,k) + f)
    r(:,k)=f-A*x(:,k);
    norm(r)
    k=k+1;
end

Edit:

Added a condition:

if abs(max(E))>=1
    fprintf('Conditions do not hold')
    break
end

However, there are times when the random matrix A satisfied this condition, and runs through the while look, and still blows up! Any ideas?

[John D'Errico]

Just wondering why you would define a variable named gamma, which overloads the very useful FUNCTION with the same name, then not even bother to use the result?

I might also point out that your code would be FAR more easy to read if you used variable names that actually meant something. Last I checked, MATLAB does not charge extra for longer names than one character. And since your question is asking for us to read and follow your code, you might have made it easy to do so.

[John D'Errico]

Oh, by the way, don't bother testing to see if the determinant of a randomly generated matrix is EXACTLY zero. That will essentially NEVER happen, even if the matrix was numerically singular.

Use rank instead.

[A]

Thank you! I'll implement the rank comment. And will clean up the naming of my parameters momentarily.

[John D'Errico]

And finally, while I usually don't bother to up-vote questions (not sure why not) I did do so here. That is because I did like that you included your code, that despite the terrible choices for variable names and use of det, it was moderately easy to read. And finally, you wrote some code, then decided to think about what it was doing, and why it was doing that. One who thinks about what it is they are doing gains my respect.

Answer 1

[John D'Errico]

(Note: I've not checked to see if your iterative code actually represents the scheme you seem to want to use. Really, I had no need to do so, because the answer is, it WILL diverge almost always for such a random mnatrix.)

Your question comes down to, under what circumstances would such a scheme diverge? Perhaps you want to consider if the Jacobi method always converges, for any matrix. (No, it will not in general converge for such a random matrix.)

As the wiki page points out,

scale=6; 
A = randn(scale,scale);

M=tril(triu(A,-1),1); %tri diagonal
N=-(M-A);

Now, what does that page point out as a general requirement for convergence for such a scheme?

max(abs(eig(M\N)))
ans =
   3.3047

Significantly greater than 1. Not gonna converge.