STHB STHB - 2 months ago 9
C++ Question

Iterations in a Jacobian method function

I'm making a function for doing the Jacobian Method. My problem is that my iterations are supposed to let the x vector elements converge to 1s, however instead of 1s the x vector elements in my functions just become 'nan's ?

void Jacobi(int n, double (* A)[20], double * b, double * x)
{
double sum[20];
double xnext[20];

cout << "sum" << endl;
for (int i = 0 ; i < 20 ; i++)
{
sum[i]=0;
cout << sum[i]<<endl;
}

cout << endl;
cout << endl;

cout << "xnext" << endl;
for (int i = 0 ; i < 20 ; i++)
{
xnext[i]=0;
cout << xnext[i]<< endl;
}

//iterations loop
for (int iter = 0 ; iter < n ; iter++)
{
for (int i = 0; i < 20 ; i++)
{
sum[i] = 0;

for (int j = 0 ; j < 20 ; j++)
{
if (i != j)
sum[i] = sum[i] + (A[i][j] * x[j]) ;
}

xnext[i] = (1/A[i][i])*(b[i]-sum[i]);

}

cout << endl;

cout << "xnext "<< iter <<endl;
for(int i= 0; i < 20 ; i++)
{
x[i]=xnext[i];
cout << x[i]<< endl;
}

}

cout << endl;

cout << "New x vector" << endl;
for (int i = 0 ; i < 20 ; i++)
{
cout << x[i] << endl;
}

return;
}

Answer

It is almost right, but there are some problems.

The matrix A is not initialized correctly. The main diagonal is filled in with 2.0. However, the code that fills the next closest diagonals with -1.0 goes out of bounds when i is 19. Furthermore, none of the other elements of A are initialized. This leads to the following changes:

// initialize everything to zero
memset(A[0], 0, 20*20*sizeof(double));
// set elements for main three diagonals
for (int i = 0 ; i< 19; ++i)
{
    A[i][i] = 2.0; 
    A[i][i+1] = -1.0; 
    A[i+1][i] = -1.0; 
}
// set last element
A[19][19] = 2.0;

The next change is to the Jacobi method itself. xnext[i] = (1.0/A[i][i])*(b[i] - sum[i]); is the correct formula for the next x. However, the next line replaces xnext with the delta between the iterates. After increasing the number of iterations to 600, the resulting vector ended up being pretty close to what Octave came up for a solution.

Comments