Simple Jacobi Iteration (help please...)
CodeGuru Home VC++ / MFC / C++ .NET / C# Visual Basic VB Forums Developer.com
Results 1 to 12 of 12

Thread: Simple Jacobi Iteration (help please...)

  1. #1
    Join Date
    Jul 2010
    Posts
    3

    Simple Jacobi Iteration (help please...)

    I'm very new to c++ programming, and my first Homework is to implement some iterations for linear system.. is it too much to ask if you could check what's wrong with my code?

    here's the code:

    my target is to solve a predefined 3500x3500 matrix using jacobi iteration... the code works perfectly for a simple 3x3 matrix (correct answer up to 5 decimal places) however, problem arises when the matrix gets larger (up to 120x120, the matrix still outputs an answer)... iteration always end up with "-1.#IND" error...please please...

    any comment or help would be greatly appreciated..many thanks

    #include <iostream.h>
    #include <math.h>
    #define n 121 //number of equations (ok up to 120) [target n=3500]

    float a[n][n];
    int main()
    {
    int i, j;
    float b[n], h = 1.0/(float)(n+1);
    //left hand side matrix of the equation ax=b
    for (i=1;i<=n-1;i++)
    for (j=1;j<=n-1;j++)
    a[i][j] = 0.0;
    a[1][1] = -(2.0 + 4.0*h*h);
    a[1][2] = 1.0 ;
    for (i=2;i<=n-1;i++)
    {
    a[i][i-1] = 1.0 ;
    a[i][i] = -(2.0 + 4.0*h*h);
    a[i][i+1] = 1.0 ;
    }
    a[n][n-1] = 1.0;
    a[n][n] = -(2.0 + 4.0*h*h);

    //right hand side of the vector of the equation ax=b
    for(i=1;i<=n-1;i++)
    b[i] = h*h;
    b[n] = h*h - 1.0;

    //calculation using Jacobi Iteration
    double pold[n], pnew[n], tol = 0.000001; //tol = set as tolerance
    //pold[]=vector of old values
    //pnew[]=vector of new results
    int counter = 0, s;
    s = 0;
    //initialize
    for (i=1;i<=n;i++)
    pnew[i] = pold[i] = 0.0;
    //computation loop
    while (s!=1)
    {
    for(i=1;i<=n;i++)
    {
    double sum = b[i];
    for(j=1;j<=i-1;j++)
    sum -= a[i][j]*pold[j];
    for(j=i+1;j<=n;j++)
    sum -= a[i][j]*pold[j];
    pnew[i] = sum/(a[i][i]);
    if (fabs(pnew[i]-pold[i])<=tol)
    s = 1;
    }
    counter++;
    if (counter%1000 == 0) //indication if program still running
    cout << "|";
    for (i=1;i<=n;i++)
    pold[i] = pnew[i];
    }
    //print output
    cout<<"\n";
    for (i=1;i<=n;i++)
    cout<<pnew[i]<<endl;
    cout<<"\nno. of iterations = "<<counter<<endl;
    cin>>i; //system 'PAUSE'
    return 0;
    }

  2. #2
    Join Date
    Jul 2002
    Location
    Portsmouth. United Kingdom
    Posts
    2,712

    Re: Simple Jacobi Iteration (help please...)

    Code:
    for (i = 1; i <= n; i++)
    
    should be
    
    for (i = 0; i < n; i++)
    <iostream.h> is an old and now non-standard header.
    The correct header is <iostream>

    Remember to use code tags when posting code

    What compiler are you using?
    "It doesn't matter how beautiful your theory is, it doesn't matter how smart you are. If it doesn't agree with experiment, it's wrong."
    Richard P. Feynman

  3. #3
    Join Date
    Jul 2002
    Location
    Portsmouth. United Kingdom
    Posts
    2,712

    Re: Simple Jacobi Iteration (help please...)

    Also a[n][n] is illegal as it is beyond the end of the array.

    Your arrays go from 0 to n - 1
    "It doesn't matter how beautiful your theory is, it doesn't matter how smart you are. If it doesn't agree with experiment, it's wrong."
    Richard P. Feynman

  4. #4
    Join Date
    Jul 2010
    Posts
    3

    Re: Simple Jacobi Iteration (help please...)

    Thanks for replying..

    i'm using the freeware Dev-C++ 4.9.9.2 for my compiler

    as for the array.. yeah, I know that they go from 0 to n-1, however my Professor insists on using arrays from 1 to n (we are allowed to extend the array to any number but we have to comply with the 3500x3500 [i.e. "n" could be set to 3505])

    and...I did change my library to iostream...still works the same
    maybe, the problem lies with the arrays...i'll try to change the arrays, i'll post the result later...

    btw..sorry for the code...
    here's the repost...

    Code:
    #include <iostream>
    #include <math.h>
    #define n 3 //number of equations (ok up to 120) [target n=3500]
    using namespace std;
    float a[n][n];
    int main()
    {
        int i, j;
        float b[n], h = 1.0/(float)(n+1);
        //left hand side matrix of the equation ax=b
        for (i=1;i<=n-1;i++)
            for (j=1;j<=n-1;j++)
                a[i][j] = 0.0;
        a[1][1] = -(2.0 + 4.0*h*h);
        a[1][2] = 1.0 ;
        for (i=2;i<=n-1;i++)
                {
                    a[i][i-1] = 1.0 ;
                    a[i][i] = -(2.0 + 4.0*h*h);
                    a[i][i+1] = 1.0 ;
                }
        a[n][n-1] = 1.0;
        a[n][n] = -(2.0 + 4.0*h*h);
        
        //right hand side of the vector of the equation ax=b
        for(i=1;i<=n-1;i++)
            b[i] = h*h;
        b[n] = h*h - 1.0;    
        
        //calculation using Jacobi Iteration
        double pold[n], pnew[n], tol = 0.000001; //tol = set as tolerance
        //pold[]=vector of old values
        //pnew[]=vector of new results    
        int counter = 0, s;
        s = 0;
        //initialize
        for (i=1;i<=n;i++)
            pnew[i] = pold[i] = 0.0;
        //computation loop
        while (s!=1)
        {
            for(i=1;i<=n;i++)
            {
              double sum = b[i];
              for(j=1;j<=i-1;j++)
                 sum -= a[i][j]*pold[j];
              for(j=i+1;j<=n;j++) 
                 sum -= a[i][j]*pold[j];
              pnew[i] = sum/(a[i][i]);
              if (fabs(pnew[i]-pold[i])<=tol)
                 s = 1;
            }
            counter++;
            if (counter&#37;1000 == 0)           //indication if program still running
                    cout << "|";
            for (i=1;i<=n;i++) 
                pold[i] = pnew[i];
        }
        //print output
        cout<<"\n";
        for (i=1;i<=n;i++) 
            cout<<pnew[i]<<endl;
        cout<<"\nno. of iterations = "<<counter<<endl;
        cin>>i;      //system 'PAUSE'
        return 0;
    }
    i'll get back after I finish the replacement of the array...

  5. #5
    Join Date
    Jul 2002
    Location
    Portsmouth. United Kingdom
    Posts
    2,712

    Re: Simple Jacobi Iteration (help please...)

    Quote Originally Posted by iahn13 View Post
    i'm using the freeware Dev-C++ 4.9.9.2 for my compiler
    That IDE is no longer under development.

    If you want to learn modern C++ then you should use something like Code:Blocks or Visual Studio Express, both of which are free.
    and...I did change my library to iostream...still works the same
    Yes it would, it's just that iostream.h is no longer a compliant or supported header. Dev-C++ supports it because it is also old and non-compliant.
    maybe, the problem lies with the arrays..
    Yes, you are stepping over the boundaries of the array when you use 'n' as an index.
    "It doesn't matter how beautiful your theory is, it doesn't matter how smart you are. If it doesn't agree with experiment, it's wrong."
    Richard P. Feynman

  6. #6
    Join Date
    Jun 2009
    Location
    France
    Posts
    2,290

    Re: Simple Jacobi Iteration (help please...)

    Quote Originally Posted by iahn13 View Post
    however my Professor insists on using arrays from 1 to n
    He's right to do so. He wants you to calculate a complicated mathematical function, for which all formulas are conventionally calculated with matrixes who's indexes go from 1 to n (inlcuded).

    Changing the indexes from 0 to n-1 would just make the code hell, and a grader's nightmare.

    Simply declare your array as int[n+1][n+1], and you'll be able to index from 1 to n included.

    If you feel you are wasting space (even ever so little), you can create a wrapper object, which translates your indexes from [1, n] to [0, n-1]. But try to stick with an interface from 1 to n
    Is your question related to IO?
    Read this C++ FAQ LITE article at parashift by Marshall Cline. In particular points 1-6.
    It will explain how to correctly deal with IO, how to validate input, and why you shouldn't count on "while(!in.eof())". And it always makes for excellent reading.

  7. #7
    Join Date
    Jul 2005
    Location
    Netherlands
    Posts
    2,012

    Re: Simple Jacobi Iteration (help please...)

    Quote Originally Posted by iahn13 View Post
    as for the array.. yeah, I know that they go from 0 to n-1, however my Professor insists on using arrays from 1 to n
    Well, your professor is wrong. In C++, indices are 0-based... always. There is nothing to be gained from these types of tricks to make indices 1-based. At best, you (and all other students) will spend 3 times as long learning how to use arrays. First, you learn to use them with 1-based indices. Then you have to unlearn using 1-based indices and, finally, you have to learn using 0-based indices. That's a total waste of effort.
    Cheers, D Drmmr

    Please put [code][/code] tags around your code to preserve indentation and make it more readable.

    As long as man ascribes to himself what is merely a posibility, he will not work for the attainment of it. - P. D. Ouspensky

  8. #8
    Join Date
    Apr 1999
    Posts
    27,427

    Re: Simple Jacobi Iteration (help please...)

    Quote Originally Posted by D_Drmmr View Post
    Well, your professor is wrong. In C++, indices are 0-based... always. There is nothing to be gained from these types of tricks to make indices 1-based.
    I can vouch for not trying to make C++ act like math arrays, i.e. making C++ arrays act like 1-based arrays.

    There is too much potential for "off-by-one" errors when you setup a program this way. In my experience, for all programs that I've seen students write that try to fake a 1-based array system in C++, more than half of them have one-off bugs -- not general off-by-one bugs that we all encounter, but explicitly due to the 1-based faking out that's going on.

    Regards,

    Paul McKenzie

  9. #9
    Join Date
    Jan 2009
    Posts
    1,689

    Re: Simple Jacobi Iteration (help please...)

    Quote Originally Posted by iahn13 View Post
    my Professor insists on using arrays from 1 to n
    O.o I'm scared

    The only language that I can think of that starts at 1 is Pascal. A few languages let you define a range:
    Code:
    'this is Visual Basic
    dim myarray(1 to 15) as int
    But almost every language starts at zero.


    Also
    Code:
    #include <iostream.h>
    #include <math.h>
    should be
    Code:
    #include <iostream>
    #include <cmath>
    Last edited by ninja9578; July 12th, 2010 at 03:04 PM.

  10. #10
    Join Date
    Jul 2010
    Posts
    3

    Smile Re: Simple Jacobi Iteration (help please...)

    thank you guys for the reply...

    here's the thing.. remember that we are allowed to adjust the matrix in any size but we have to keep in mind the inclusion of the required n=3500..

    Quote Originally Posted by monarch_dodra

    Simply declare your array as int[n+1][n+1], and you'll be able to index from 1 to n included.
    my assumption is that the error came from the fact that the matrix is having trouble at the end of the matrix (well, at least based on inspecting/printing the code per line).. I remember now that we were advised to note that problem... anyway, I just declared my a[n][n] to a[n+5][n+5]...

    the program's now working...thank you all for the things I learned here.. all of the comments were very helpful...thanks a lot...

    btw..the reason my Prof wants the arrays go from 1 to n, is that (according to him) zeros can't be counted using fingers...really weird...

  11. #11
    Join Date
    Jun 2009
    Location
    France
    Posts
    2,290

    Re: Simple Jacobi Iteration (help please...)

    Quote Originally Posted by iahn13 View Post
    the reason my Prof wants the arrays go from 1 to n, is that (according to him) zeros can't be counted using fingers...really weird...
    Well, that's not really a good reason actually...

    I was saying that if you have to transform a really complicated mathematical formulal, the type that's a double sum of 1 to n, it might be more dangerous to try to shift the indexes of the formula by 1, then just using a "start at 1 array".

    Yes, it is less than optimal, but it is impossible to debug a program where the code is right, but the formula wrong.

    Its always a question of optimal choice, there is never a solution that fit alls.

    Where I work, we have a class that wraps a vector to create a new vector class that starts at 1. We only use it like 1% of the time, but for those 1%, they are a life saver.
    Is your question related to IO?
    Read this C++ FAQ LITE article at parashift by Marshall Cline. In particular points 1-6.
    It will explain how to correctly deal with IO, how to validate input, and why you shouldn't count on "while(!in.eof())". And it always makes for excellent reading.

  12. #12
    Join Date
    Jan 2009
    Posts
    1,689

    Re: Simple Jacobi Iteration (help please...)

    Quote Originally Posted by iahn13 View Post
    btw..the reason my Prof wants the arrays go from 1 to n, is that (according to him) zeros can't be counted using fingers...really weird...
    Uh... this is how you count to zero using fingers.


    Also, I've never seen a programmer count on his fingers. If you can't do math to 10 in your head, you shouldn't be an engineer. You also can't count to 3500 on your fingers. :P


    Disclaimer: Don't ever search for "fist" in google at work without safe-search on. Ewww.

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  


Azure Activities Information Page

Windows Mobile Development Center


Click Here to Expand Forum to Full Width

This is a CodeGuru survey question.


Featured


HTML5 Development Center