MATLAB Answers

Basic iterative schemes in matlab

9 views (last 30 days)
A
A on 5 Mar 2015
Commented: John D'Errico on 5 Mar 2015
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?
  4 Comments
John D'Errico
John D'Errico on 5 Mar 2015
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.

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 5 Mar 2015
Edited: John D'Errico on 5 Mar 2015
(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.

More Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by