
July 12th, 2010, 03:49 AM
#1
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<=n1;i++)
for (j=1;j<=n1;j++)
a[i][j] = 0.0;
a[1][1] = (2.0 + 4.0*h*h);
a[1][2] = 1.0 ;
for (i=2;i<=n1;i++)
{
a[i][i1] = 1.0 ;
a[i][i] = (2.0 + 4.0*h*h);
a[i][i+1] = 1.0 ;
}
a[n][n1] = 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<=n1;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<=i1;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;
}

July 12th, 2010, 04:23 AM
#2
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 nonstandard 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

July 12th, 2010, 04:30 AM
#3
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

July 12th, 2010, 05:11 AM
#4
Re: Simple Jacobi Iteration (help please...)
Thanks for replying..
i'm using the freeware DevC++ 4.9.9.2 for my compiler
as for the array.. yeah, I know that they go from 0 to n1, 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<=n1;i++)
for (j=1;j<=n1;j++)
a[i][j] = 0.0;
a[1][1] = (2.0 + 4.0*h*h);
a[1][2] = 1.0 ;
for (i=2;i<=n1;i++)
{
a[i][i1] = 1.0 ;
a[i][i] = (2.0 + 4.0*h*h);
a[i][i+1] = 1.0 ;
}
a[n][n1] = 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<=n1;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<=i1;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;
}
i'll get back after I finish the replacement of the array...

July 12th, 2010, 05:37 AM
#5
Re: Simple Jacobi Iteration (help please...)
Originally Posted by iahn13
i'm using the freeware DevC++ 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. DevC++ supports it because it is also old and noncompliant.
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

July 12th, 2010, 06:57 AM
#6
Re: Simple Jacobi Iteration (help please...)
Originally Posted by iahn13
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 n1 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, n1]. 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 16.
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.

July 12th, 2010, 07:04 AM
#7
Re: Simple Jacobi Iteration (help please...)
Originally Posted by iahn13
as for the array.. yeah, I know that they go from 0 to n1, however my Professor insists on using arrays from 1 to n
Well, your professor is wrong. In C++, indices are 0based... always. There is nothing to be gained from these types of tricks to make indices 1based. 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 1based indices. Then you have to unlearn using 1based indices and, finally, you have to learn using 0based 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

July 12th, 2010, 02:34 PM
#8
Re: Simple Jacobi Iteration (help please...)
Originally Posted by D_Drmmr
Well, your professor is wrong. In C++, indices are 0based... always. There is nothing to be gained from these types of tricks to make indices 1based.
I can vouch for not trying to make C++ act like math arrays, i.e. making C++ arrays act like 1based arrays.
There is too much potential for "offbyone" 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 1based array system in C++, more than half of them have oneoff bugs  not general offbyone bugs that we all encounter, but explicitly due to the 1based faking out that's going on.
Regards,
Paul McKenzie

July 12th, 2010, 02:56 PM
#9
Re: Simple Jacobi Iteration (help please...)
Originally Posted by iahn13
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.

July 12th, 2010, 10:06 PM
#10
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..
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...

July 13th, 2010, 01:05 AM
#11
Re: Simple Jacobi Iteration (help please...)
Originally Posted by iahn13
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 16.
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.

July 13th, 2010, 08:05 AM
#12
Re: Simple Jacobi Iteration (help please...)
Originally Posted by iahn13
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 safesearch 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

Forum Rules

Click Here to Expand Forum to Full Width
