Solve Problems in Loops: Get Feedback & Robustness

  • Thread starter Thread starter mfarooq52003
  • Start date Start date
  • Tags Tags
    Loops
AI Thread Summary
The discussion revolves around two algorithms designed to solve a linear system of equations Ax=b using iterative methods. The user seeks feedback on potential errors, particularly within the loop structures of both algorithms. Key issues highlighted include concerns about the algorithms not executing the 'else' part of the loop when conditions are met, which may indicate a flaw in the termination condition. The user emphasizes the need for the algorithms to stop once the residual equals a specified epsilon value, suggesting a possible need to adjust the comparison operator from 'not equal' to 'greater than or equal to' or 'less than or equal to' for proper termination. Participants recommend using debugging tools to step through the code and identify any logical errors, and they encourage adding comments for clarity. Overall, the user aims to refine the algorithms for a comparative analysis of their performance.
mfarooq52003
Messages
21
Reaction score
0
Hi There

Can anyone check the following algorithm and tell me what I am doing wrong. If every thing is OK then please give me some feedback to make it robust. Please check both the algorithms

Algorithm-1

B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
r(:,1)=b(:,1)-A*x(:,1);
y(:,1)=y;
p=A*r(:,1);
p1=A*p;
c0=y'*r(:,1);
c1=y'*p;
c2=y'*p1;
c3=y'*A*p1;
delta=c1*c3-c2^2;
alpha=(c0*c3-c1*c2)/delta;
beta=(c0*c2-c1^2)/delta;
r(:,2)=r(:,1)-(c0*p)/c1;
x(:,2)=x(:,1)+(c0/c1)*r(:,1);
r(:,3)=r(:,1)-alpha*p+beta*p1;
x(:,3)=x(:,1)+alpha*r(:,1)-beta*p;
y(:,2)=A'*y(:,1);
y(:,3)=A'*y(:,2);
y(:,4)=A'*y(:,3);
eps=0.0001;
tic;
%Step-2
for k=3:m
y(:,k+2)=A'*y(:,k+1);
q1=A*r(:,k-1);
q2=A*q1;% q2=A*q1=A*(A*r(:,k-1))=A^2r(:,k-1)
q3=A*r(:,k-2);
a11=y(:,k-1)'*r(:,k-1);
a13=y(:,k-2)'*r(:,k-2);
a21=y(:,k)'*r(:,k-1);
a22=a11;
a23=y(:,k-1)'*r(:,k-2);
a31=y(:,k+1)'*r(:,k-1);
a32=a21;
a33=y(:,k)'*r(:,k-2);
s=y(:,k+2)'*r(:,k-1);
t=y(:,k+1)'*r(:,k-2);
F_{k+1}=-a11/a13;
b1=-a21-F_{k+1}*a23;
b2=-a31-F_{k+1}*a33;
b3=-s-F_{k+1}*t;
Delta_{k+1}=a11*(a22*a33-a32*a23)+a13*(a21*a32-a31*a22);
B_{k+1}=1/Delta_{k+1}*[b1*(a22*a33-a32*a23)+a13*(b2*a32-b3*a22)];
G_{k+1}=[b1-B_{k+1}*a11]/a13;
C_{k+1}=[b2-B_{k+1}*a21-G_{k+1}*a23]/a22;
A_{k+1}=1/[C_{k+1}+G_{k+1}];
r(:,k+1)=A_{k+1}*[q2+B_{k+1}*q1+C_{k+1}*r(:,k-1)+F_{k+1}*q3+G_{k+1}*r(:,k-2)];
x(:,k+1)=A_{k+1}*[C_{k+1}*x(:,k-1)+G_{k+1}*x(:,k-2)-(q1+B_{k+1}*r(:,k-1)+F_{k+1}*r(:,k-2))];
%Step-3
if (r(:,k+1))~=eps
continue;
else
x=x(:,k+1);
end
end
toc;
x1=A\b;
x2=x(:,k+1);
Res=[x1, x2];
E=norm(x1-x2)
Er=E/norm(x1)

Algorithm-2
B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
y(:,1)=y;
r(:,1)=b-A*x(:,1);
p(:,1)=r(:,1);
C=1;
A_{2}=(-1*y(:,1)'*r(:,1))/(y(:,1)'*A*r(:,1));
x(:,2)=x(:,1)-A_{2}*r(:,1);
r(:,2)=r(:,1)+A_{2}*A*r(:,1);
eps=0.0001;
tic;
%Step-2
for k=2:m
y(:,k)=A'*y(:,k-1);
D_{k+1}=(-1*y(:,k)'*r(:,k))/(C*y(:,k)'*p(:,k-1));
p(:,k)=r(:,k)+D_{k+1}*C*p(:,k-1);
A_{k+1}=(-1*y(:,k)'*r(:,k))/(y(:,k)'*A*p(:,k));
r(:,k+1)=r(:,k)+A_{k+1}*A*p(:,k);
x(:,k+1)=x(:,k)-A_{k+1}*p(:,k);
%Step-3
if r(:,k+1)~=eps,
C=C/A_{k};
fprintf('I am here')
else
x=x(:,k+1);
fprintf('You are here')
end
end
toc;
x2=x(:,k+1);
x1=A\b;
Res=[x1, x2];
E=norm(x1-x2)
Er=E/norm(x1)

Please let me know of any blunder if I have done there in the loops and I will appreciate if you can correct it and post it in correct form.

I will be grateful for your help and support.

Regards
 
Technology news on Phys.org
If you can only check the loop part of both the algorithm and see if there is any mistake, if you people don't have enough time for the whole algorithm.
 
mfarooq52003 said:
If you can only check the loop part of both the algorithm and see if there is any mistake, if you people don't have enough time for the whole algorithm.

It would sure help if you could:

1) Describe what you are trying to do with the two pieces of code.

2) Use "code" tags around your code, so that indentation is preserved. It's basically unreadable without indentations.

3) Tell us why you think there may be a problem. What are the symptoms that make you think some loop somewhere has a problem?
 
Thank you very much for such a prompt response. Actually I think my programme is not going through the else part of the algorithm even when the condition is satisfied, I guess. I want to get result from these two algorithms and then do the comparison that which one is better.
B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
r(:,1)=b(:,1)-A*x(:,1);
y(:,1)=y;
p=A*r(:,1);
p1=A*p;
c0=y'*r(:,1);
c1=y'*p;
c2=y'*p1;
c3=y'*A*p1;
delta=c1*c3-c2^2;
alpha=(c0*c3-c1*c2)/delta;
beta=(c0*c2-c1^2)/delta;
r(:,2)=r(:,1)-(c0*p)/c1;
x(:,2)=x(:,1)+(c0/c1)*r(:,1);
r(:,3)=r(:,1)-alpha*p+beta*p1;
x(:,3)=x(:,1)+alpha*r(:,1)-beta*p;
y(:,2)=A'*y(:,1);
y(:,3)=A'*y(:,2);
y(:,4)=A'*y(:,3);
eps=0.0001;
tic;
%Step-2
for k=3:m
y(:,k+2)=A'*y(:,k+1);
q1=A*r(:,k-1);
q2=A*q1;% q2=A*q1=A*(A*r(:,k-1))=A^2r(:,k-1)
q3=A*r(:,k-2);
a11=y(:,k-1)'*r(:,k-1);
a13=y(:,k-2)'*r(:,k-2);
a21=y(:,k)'*r(:,k-1);
a22=a11;
a23=y(:,k-1)'*r(:,k-2);
a31=y(:,k+1)'*r(:,k-1);
a32=a21;
a33=y(:,k)'*r(:,k-2);
s=y(:,k+2)'*r(:,k-1);
t=y(:,k+1)'*r(:,k-2);
F_{k+1}=-a11/a13;
b1=-a21-F_{k+1}*a23;
b2=-a31-F_{k+1}*a33;
b3=-s-F_{k+1}*t;
Delta_{k+1}=a11*(a22*a33-a32*a23)+a13*(a21*a32-a31*a22);
B_{k+1}=1/Delta_{k+1}*[b1*(a22*a33-a32*a23)+a13*(b2*a32-b3*a22)];
G_{k+1}=[b1-B_{k+1}*a11]/a13;
C_{k+1}=[b2-B_{k+1}*a21-G_{k+1}*a23]/a22;
A_{k+1}=1/[C_{k+1}+G_{k+1}];
r(:,k+1)=A_{k+1}*[q2+B_{k+1}*q1+C_{k+1}*r(:,k-1)+F_{k+1}*q3+G_{k+1}*r(:,k-2)];
x(:,k+1)=A_{k+1}*[C_{k+1}*x(:,k-1)+G_{k+1}*x(:,k-2)-(q1+B_{k+1}*r(:,k-1)+F_{k+1}*r(:,k-2))];
%Step-3
if (r(:,k+1))~=eps
continue;
else
x=x(:,k+1);
end
end
toc;
 
Dear Berkeman
Can you read it now or should I email you in attachment as M-file.

Cheers
Farooq
 
What language is this? Do you have single-step and watch capabilities in a debugger for it? You can watch to see why the condition seems to always be met.
 
mfarooq52003 said:
Dear Berkeman
Can you read it now or should I email you in attachment as M-file.

Cheers
Farooq

It's still pretty unreadable. Use [ code ] [ /code ] tags (without the spaces) around your code. Your source code is indented, right?
 
Code: 
B=gallery('tridiag',10,-6,4,4);n = 10;
BL = kron(eye(n),B);
[p1,q1] = size(B);
BLI = BL - diag(ones(1,n*q1-p1),-p1) - diag(1*ones(1,n*q1-p1),p1);
A=BLI;
X(:,1)=[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
b(:,1)=A*X(:,1);
x(:,1)=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
y=[-6;-5;-4;-3;-2;-1;0;1;2;3;-4;-3;-2;-6;-5;-1;0;1;2;3;3;-1;0;8;3;-1;0;8;0;0;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8;-5;-4;-3;-2;-1;0;1;2;3;8];
m=100;
%Step-1
r(:,1)=b(:,1)-A*x(:,1);
y(:,1)=y;
p=A*r(:,1);
p1=A*p;
c0=y'*r(:,1);
c1=y'*p;
c2=y'*p1;
c3=y'*A*p1;
delta=c1*c3-c2^2;
alpha=(c0*c3-c1*c2)/delta;
beta=(c0*c2-c1^2)/delta;
r(:,2)=r(:,1)-(c0*p)/c1;
x(:,2)=x(:,1)+(c0/c1)*r(:,1);
r(:,3)=r(:,1)-alpha*p+beta*p1;
x(:,3)=x(:,1)+alpha*r(:,1)-beta*p;
y(:,2)=A'*y(:,1);
y(:,3)=A'*y(:,2);
y(:,4)=A'*y(:,3);
eps=0.0001;
tic;
%Step-2
for k=3:m
    y(:,k+2)=A'*y(:,k+1);
    q1=A*r(:,k-1);
    q2=A*q1;% q2=A*q1=A*(A*r(:,k-1))=A^2r(:,k-1)
    q3=A*r(:,k-2);
    a11=y(:,k-1)'*r(:,k-1);
    a13=y(:,k-2)'*r(:,k-2);
    a21=y(:,k)'*r(:,k-1);
    a22=a11;
    a23=y(:,k-1)'*r(:,k-2);
    a31=y(:,k+1)'*r(:,k-1);
    a32=a21;
    a33=y(:,k)'*r(:,k-2);
    s=y(:,k+2)'*r(:,k-1);
    t=y(:,k+1)'*r(:,k-2);
    F_{k+1}=-a11/a13;
    b1=-a21-F_{k+1}*a23;
    b2=-a31-F_{k+1}*a33;
    b3=-s-F_{k+1}*t;
    Delta_{k+1}=a11*(a22*a33-a32*a23)+a13*(a21*a32-a31*a22);
    B_{k+1}=1/Delta_{k+1}*[b1*(a22*a33-a32*a23)+a13*(b2*a32-b3*a22)];
    %a=a22*a33-a32*a23;
    %Delta_{k+1}=a11*a+a13*(a21*a32-a31*a22);
    %B_{k+1}=1/Delta_{k+1}*(b1*a+a13*(b2*a32-b3*a22));
    G_{k+1}=[b1-B_{k+1}*a11]/a13;
    C_{k+1}=[b2-B_{k+1}*a21-G_{k+1}*a23]/a22;
    A_{k+1}=1/[C_{k+1}+G_{k+1}];
    r(:,k+1)=A_{k+1}*[q2+B_{k+1}*q1+C_{k+1}*r(:,k-1)+F_{k+1}*q3+G_{k+1}*r(:,k-2)];
    x(:,k+1)=A_{k+1}*[C_{k+1}*x(:,k-1)+G_{k+1}*x(:,k-2)-(q1+B_{k+1}*r(:,k-1)+F_{k+1}*r(:,k-2))];
    %Step-3
    if (r(:,k+1))~=eps
        continue;
    else
        x=x(:,k+1);
    end
end
    toc;
    x1=A\b;
    x2=x(:,k+1);
    Res=[x1, x2];
    E=norm(x1-x2)
    Er=E/norm(x1)
 
Thank you very much for your guidance, support and help berkeman
 
  • #10
That's much more readable. You said in your PM that this is for MATLAB. Do you have a debugger with single-step and watch capability? I'm not that familiar with MATLAB. You might PM MATLABdude and point him to this thread -- as the user name implies, he's pretty conversant with that package.
 
  • #11
Also, are you sure you want ~= not equal for that comparison with epsilon, and not <= ?
 
  • #12
Well, I want it as not equal but we can do it other way as well but then it will be >=. Can you check the step-2, the loop bit only please, if there is any mistake. It runs but I want to share it with experts so that I can go for comparison and analysis after their expert opinion, if it is OK.
 
  • #13
mfarooq52003 said:
Well, I want it as not equal but we can do it other way as well but then it will be >=. Can you check the step-2, the loop bit only please, if there is any mistake. It runs but I want to share it with experts so that I can go for comparison and analysis after their expert opinion, if it is OK.

Please describe in words what your algorithm is supposed to be doing. What is the overall problem, and how are you approaching it? Most code contains comments, after all. At least my code always does...
 
  • #14
My algorithms solve a linear system of equation Ax=b using iterative approach. I just want some one, to check the loop bit and see if I have done any thing wrong, if not just tell me. The algorithm has to stop once the residual becomes equal to the epsilon and then we get the solution, otherwise we have to start another iteration.

Thank you very much for your help and prompt responses berkeman.


Farooq
 
  • #15
mfarooq52003 said:
My algorithms solve a linear system of equation Ax=b using iterative approach. I just want some one, to check the loop bit and see if I have done any thing wrong, if not just tell me. The algorithm has to stop once the residual becomes equal to the epsilon and then we get the solution, otherwise we have to start another iteration.

Well, without understanding in detail that comparison, I definitely think it should not be ~=, if you are trying for a termination condition. Use the >= or <= comparison, whichever is appropriate.

Beyond that, I'd recommend that you PM MATLABdude and ask him to take a look starting at post #8 where the code is cleaned up. He might be able so spot some other problem with your loop and termination condition.
 

Similar threads

Debugging a Backprop Neural Network: Issues & Solutions
Replies
3
Views
2K
Wrong solution order using Runge Kutta 4
Replies
15
Views
3K
C/C++ C++ matrix boundary condition problems
Replies
3
Views
2K
Python The Wimol-Banlue attractor via the Taylor Series Method
Replies
6
Views
2K
Python Lorenz ODE solver in 35 lines of pure Python
Replies
17
Views
3K
Python Thermal lattice Boltzmann model ignoring source term -- python code help please
Replies
1
Views
2K
Loop integral computation
Replies
2
Views
1K
Find the equation of the tangent plane and normal to a surface
Replies
1
Views
1K
I Understanding Legendre differential equation
Replies
2
Views
5K
B How to multiply matrix with row vector?
Replies
4
Views
3K

Hot Threads

Recent Insights

Back
Top