CODEWITHSUNDEEP
c++mathlinear-algebraphysicsscientific-computing
Mohsen Joodaki
Mohsen Joodaki

Reputation: 19

jacobi iterative method has wrong answer in c++

I'm writing Jacobi iterative method to solve any linear system of equations. this program works for some examples but doesn't work for others. for example for

A=     and B=
7  3       5
2  3       4

this will works and answers are true but for

A=     and B=
1  2       3
3  4       7

the answers are wrong and huge numbers. I really don't know what should I do to make a correct calculation. I used some other codes but still I have this issue with codes.

#include <iostream>
using namespace std;

int main(){
    double A[10][10], alpha[10][10], B[10], betha[10], x[10][100], sum[10];
    int i, j, n, k, kmax;
    cout << "insert number of equations \n";
    cin >> n;
    cout << "insert LHS of equations (a11,a12,...,ann)\n";
    for (i = 1; i <= n; i++){
        for (j = 1; j <= n; j++){
            cin >> A[i][j];
        }
    }
    cout << "A=\n";
    for (i = 1; i <= n; i++){
        for (j = 1; j <= n; j++){
            cout << A[i][j] << "\t\t";
        }
        cout << "\n\n";
    }
    cout << "alpha=\n";
    for (i = 1; i <= n; i++){
        for (j = 1; j <= n; j++){
            if (i == j){
                alpha[i][j] = 0;
            }
            else{
                alpha[i][j] = -A[i][j] / A[i][i];
            }
        }
    }
    for (i = 1; i <= n; i++){
        for (j = 1; j <= n; j++){
            cout << alpha[i][j] << "\t\t";
        }
        cout << "\n\n";
    }
    cout << "insert RHS of equations";
    for (i = 1; i <= n; i++){
        cin >> B[i];
    }
    cout << "\nbetha=\n";
    for (i = 1; i <= n; i++){
        betha[i] = B[i] / A[i][i];
        cout << betha[i] << endl;
    }
    cout << "Enter the number of repetitions." << endl;
    cin >> kmax;
    k = 0;
    for (i = 1; i <= n; i++){
        sum[i] = 0;
        x[i][k] = betha[i];    //initial values 
    }
    for (k = 0; k <= kmax; k++){
        for (i = 1; i <= n; i++){
            for (j = 1; j <= n; j++){
                sum[i] += alpha[i][j] * x[j][k];
            }
            x[i][k] = betha[i] + sum[i];
            sum[i] = 0; 
        }
    }
    cout << "answers:\n\n";
    for (i = 1; i <= n; i++){
        cout << x[i][kmax] << endl;
    }
    return 0;
}

Upvotes: 0

Views: 1433

Answers (2)

duffymo
duffymo

Reputation: 308763

Your matrix, in row order, is: [{1, 2} {3, 4}]

It has determinant equal to -2; clearly it's not singular.

It has inverse: [{4, -2}, {-3, 1}]/(-2)

The correct solution is: {1, 1}

You can verify this by substituting back into the original equation and checking to make sure you have an identity: [{1, 2} {3, 4}]{1, 1} = {3, 7}

Iterative methods can be sensitive to initial conditions.

The point about diagonally dominant is correct. Perhaps a more judicious choice of initial condition, closer to the right answer, would allow you to converge.

Update:

Jacobi iteration decomposes the matrix into diagonal elements D and off-diagonal elements R:

enter image description here

Jacobi will converge if:

Convergence criteria

Since this is not the case for the first row of your sample matrix, you might have a problem.

You still get there in a single step if you use the correct answer as your initial guess. This says that even Jacobi will work with a judicious choice.

If I start with {1, 1} I converge to the correct answer in a single iteration.

Upvotes: 0

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

You should again check the condition for convergence. There you will find that usually the method only converges for diagonally dominant matrices. The first example satisfies that condition, while the second violates it clearly.

If convergence is not guaranteed, divergence might happen, as you found.

More specifically, the Jacobi iteration in the second example computes

xnew[0] = (3 - 2*x[1])/1;
xnew[1] = (7 - 3*x[0])/4;

Over two iterations the composition of steps gives

xtwo[0] = (3 - 2*xnew[1])/1 = -0.5 + 1.5*x[0];
xtwo[1] = (7 - 3*xnew[0])/4 = -0.5 + 1.5*x[1];

which is clearly expanding the initial errors with factor 1.5.

Upvotes: 1

Related Questions