Solve Problems in Loops: Get Feedback & Robustness

  • Thread starter mfarooq52003
  • Start date
  • Tags
    Loops
In summary, to efficiently solve problems using loops, it is important to break the problem down into smaller parts and use the appropriate loop structure for each part. Regularly testing and debugging can improve the robustness of loop-based solutions. Some common mistakes to avoid include forgetting to initialize or update loop variables and creating infinite loops. To make loop-based solutions more efficient, consider using more efficient loop structures and minimizing the number of times variables are accessed or calculations are performed. While loops can be used to solve a variety of problems, it is important to consider other methods or algorithms that may be more efficient or appropriate for certain problems.
  • #1
mfarooq52003
22
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
  • #2
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.
 
  • #3
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?
 
  • #4
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;
 
  • #5
Dear Berkeman
Can you read it now or should I email you in attachment as M-file.

Cheers
Farooq
 
  • #6
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.
 
  • #7
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?
 
  • #8
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)
 
  • #9
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.
 

1. How can I efficiently solve problems using loops?

To efficiently solve problems using loops, it is important to break the problem down into smaller, more manageable parts and use the appropriate loop structure for each part. Additionally, regularly testing and debugging your code can help identify and fix any errors or inefficiencies.

2. How do I get feedback and improve the robustness of my loop-based solutions?

To get feedback and improve the robustness of your loop-based solutions, it is helpful to have a clear understanding of the problem you are trying to solve and the expected output. You can also seek feedback from peers or use test cases to ensure your solution works for different scenarios.

3. What are some common mistakes to avoid when using loops to solve problems?

Some common mistakes to avoid when using loops to solve problems include forgetting to initialize or update loop variables, creating infinite loops, and not considering edge cases or unexpected inputs. It is important to carefully plan and test your code to avoid these mistakes.

4. How can I make my loop-based solutions more efficient?

To make your loop-based solutions more efficient, you can consider using more efficient loop structures such as while loops instead of for loops, minimizing the number of times you access variables or perform calculations, and using built-in functions or methods when applicable.

5. Can loops be used to solve any type of problem?

While loops can be used to solve a wide range of problems, they may not always be the most efficient or appropriate solution. It is important to consider the problem at hand and determine if a loop is the best approach or if there may be other methods or algorithms that can be used.

Similar threads

Debugging a Backprop Neural Network: Issues & Solutions
    • Programming and Computer Science
Replies
3
Views
1K
Wrong solution order using Runge Kutta 4
    • Programming and Computer Science
Replies
15
Views
2K
Thermal lattice Boltzmann model ignoring source term -- python code help please
    • Programming and Computer Science
Replies
1
Views
1K
Python The Wimol-Banlue attractor via the Taylor Series Method
    • Programming and Computer Science
Replies
6
Views
1K
Lorenz ODE solver in 35 lines of pure Python
    • Programming and Computer Science
Replies
17
Views
2K
B How to multiply matrix with row vector?
    • Linear and Abstract Algebra
Replies
4
Views
1K
C/C++ C++ matrix boundary condition problems
    • Programming and Computer Science
Replies
3
Views
1K
Condition for coplanarity of two lines
    • Precalculus Mathematics Homework Help
Replies
14
Views
272
Python "ImportError: bad magic number ...."
    • Programming and Computer Science
Replies
1
Views
2K
C/C++ How to Implement Polynomial Operations in C++?
    • Programming and Computer Science
Replies
5
Views
2K

Hot Threads

Recent Insights

Back
Top