Solve Problems in Loops: Get Feedback & Robustness

  • Thread starter Thread starter mfarooq52003
  • Start date Start date
  • Tags Tags
    Loops
Click For Summary

Discussion Overview

The discussion revolves around two algorithms intended for solving problems involving loops in numerical computations. Participants seek feedback on the correctness and robustness of the algorithms, particularly focusing on the loop structures and any potential issues that may arise during execution.

Discussion Character

  • Technical explanation
  • Debate/contested
  • Mathematical reasoning

Main Points Raised

  • One participant requests a review of both algorithms, specifically looking for mistakes and suggestions for improving robustness.
  • Another participant suggests focusing solely on the loop parts of the algorithms if time is limited for a full review.
  • A participant emphasizes the need for clarity in the code presentation, suggesting the use of "code" tags for better readability and asking for a description of the intended functionality of the algorithms.
  • One participant expresses concern that their program is not executing the expected "else" part of the algorithm, indicating a potential logical issue in the loop condition.
  • Another participant inquires about the programming language used and suggests utilizing debugging tools to investigate why the loop condition appears to always be satisfied.

Areas of Agreement / Disagreement

There is no consensus on the correctness of the algorithms or the specific issues within the loops, as participants have raised various concerns and suggestions without resolving them.

Contextual Notes

Participants have noted potential issues with the loop conditions and the overall structure of the algorithms, but specific assumptions or definitions that may affect the outcomes have not been fully explored.

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 ·
Replies
3
Views
2K
Wrong solution order using Runge Kutta 4
  • · Replies 15 ·
Replies
15
Views
3K
C/C++ C++ matrix boundary condition problems
  • · Replies 3 ·
Replies
3
Views
2K
Python The Wimol-Banlue attractor via the Taylor Series Method
  • · Replies 6 ·
Replies
6
Views
2K
Python Lorenz ODE solver in 35 lines of pure Python
  • · Replies 17 ·
Replies
17
Views
3K
Python Thermal lattice Boltzmann model ignoring source term -- python code help please
  • · Replies 1 ·
Replies
1
Views
3K
Loop integral computation
  • · Replies 2 ·
Replies
2
Views
2K
Find the equation of the tangent plane and normal to a surface
  • · Replies 1 ·
Replies
1
Views
2K
Undergrad Understanding Legendre differential equation
  • · Replies 2 ·
Replies
2
Views
6K
Graduate Have I solved this iterative equation correctly?
  • · Replies 2 ·
Replies
2
Views
1K