Solve Problems in Loops: Get Feedback & Robustness

  • Thread starter Thread starter mfarooq52003
  • Start date Start date
  • Tags Tags
    Loops
Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
Sign up Log in
14 replies · 2K views
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
 
Physics 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
 
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.
 
Also, are you sure you want ~= not equal for that comparison with epsilon, and not <= ?
 
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.
 
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...
 
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
 
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.