Simple Jacobi Iteration (help please...)

[iahn13]

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;
}

[JohnW@Wessex]

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?

[JohnW@Wessex]

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

Your arrays go from 0 to n - 1

[iahn13]

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...

#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...

[JohnW@Wessex]

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.

[monarch_dodra]

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

[D_Drmmr]

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.

[Paul McKenzie]

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

[ninja9578]

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:

'this is Visual Basic
dim myarray(1 to 15) as int

But almost every language starts at zero.


Also

#include <iostream.h>
#include <math.h>

should be

#include <iostream>
#include <cmath>

[iahn13]

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..

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...

[monarch_dodra]

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.

[ninja9578]

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.