Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Problem in loops

  1. Aug 4, 2009 #1
    Hi There

    Can any one 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
     
  2. jcsd
  3. Aug 4, 2009 #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.
     
  4. Aug 4, 2009 #3

    berkeman

    User Avatar

    Staff: Mentor

    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?
     
  5. Aug 4, 2009 #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.
     
  6. Aug 4, 2009 #5
    Dear Berkeman
    Can you read it now or should I email you in attachment as M-file.

    Cheers
    Farooq
     
  7. Aug 4, 2009 #6

    berkeman

    User Avatar

    Staff: Mentor

    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.
     
  8. Aug 4, 2009 #7

    berkeman

    User Avatar

    Staff: Mentor

    It's still pretty unreadable. Use [ code ] [ /code ] tags (without the spaces) around your code. Your source code is indented, right?
     
  9. Aug 4, 2009 #8
    Code (Text):
    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)
     
  10. Aug 4, 2009 #9
    Thank you very much for your guidance, support and help berkeman
     
  11. Aug 4, 2009 #10

    berkeman

    User Avatar

    Staff: Mentor

    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.
     
  12. Aug 4, 2009 #11

    berkeman

    User Avatar

    Staff: Mentor

    Also, are you sure you want ~= not equal for that comparison with epsilon, and not <= ?
     
  13. Aug 4, 2009 #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.
     
  14. Aug 4, 2009 #13

    berkeman

    User Avatar

    Staff: Mentor

    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...
     
  15. Aug 4, 2009 #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
     
  16. Aug 4, 2009 #15

    berkeman

    User Avatar

    Staff: Mentor

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Problem in loops
Loading...