C++

Forum

 
Jacobi Iterative Method Issue

hans4354 (2)
Hi all,

Attempting to create a program that uses the Jacobi Iterative Method to solve an 'n'-dimensional A.x=b system (which I can then base Gauss-Seidel program on). I wish to use user input to determine not only the coefficient matrix and constant vector, but also the size of the system; I also would like to use the two norm of the difference between X[m] and X[m-1] (iterate m and m-1 for vector x that solves A.x=b) to decide when to exit the algorithm.

I have it written, and it almost works except that it does not return the correct X vector that solves the system -- it seems that the issue is the program is ignoring the if statements controlling the dummy variable 'sum[j]'.

I was hoping someone could help me troubleshoot as I have been trying for several hours and just can't find the issue...

code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
 
// Created by Kyle James Hansen on 9 May 2016
// Copyright © 2016 Kyle James Hansen. All rights reserved.

#include <iostream>
#include <string>
#include <stdio.h>
#include <math.h>
#include <iomanip>
#include <stdlib.h>
using namespace std;

double abs(double ABS)
{
    return sqrt(ABS*ABS);
}

double chop(double CHOP)
{
    if (abs(CHOP) < pow(10,-15))
    {
        CHOP = 0;
    }
    return CHOP;
}

int main()
{
    int i,j,m,N,n,w=15;
    
    // Define size of system
    cout << "Solve system [A.x = b] using Jacobi Iteration"<< endl << endl << "Enter 'N', the size of the system:" << endl;
    cin >> N;
    n = N - 1;
    
    double A[N][N],b[N],x[N][1000],xnorm,errtol=pow(10,-10),sum[N];
    
    // Define matrix A (coefficient matrix)
    cout << endl << "Enter coefficients of " << N << "*" << N << " system:" << endl;
    for(j=0;j<=n;j=j+1)
    {
        for (i=0;i<=n;i=i+1)
        {
            cin >> A[j][i];
        }
    }
    
    // Define constant vector b
    cout << endl << "Enter " << N << " constants of RHS of system:" << endl;
    for(j=0;j<=n;j=j+1)
    {
        cin >> b[j];
    }
    
    // Define initial guess for vector x that solves [A.x = b]
    m=0;
    for(j=0;j<=n;j=j+1)
    {
        x[j][m] = j + 1;
    }
    
    // Print system as check
    cout << endl << "System:" << endl;
    for(j=0;j<=n;j=j+1)
    {
        for(i=0;i<=n;i=i+1)
        {
            cout << A[j][i];
            if (i<n)
            {
                cout << setw(5);
            }
        }
        cout << " | " << b[j] << endl;
    }
    
    // Start output table
    cout << endl << "m";
    for(j=0;j<=n;j=j+1)
    {
        cout << setw(w-1) << "x" << j;
    }
    cout << setw(w) << "||X||2" << endl;
    cout << m;
    for(j=0;j<=n;j=j+1)
    {
        cout << setw(w) << x[j][m];
    }
    cout << setw(w) << "-" << endl;
    
    // Jacobi algorithm
    do
    {
        m = m + 1;
        
        // Iterate one row at a time (x0,x1,...,xn)
        for (j = 0; j <= n; j = j + 1)
        {
            sum[j] = 0;
            
            // for x0
            if ((j = 0))
            {
                for (i = 1; i <= n; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m - 1]);
                }
                x[j][m] = (1/A[j][j])*(b[j] - sum[j]);
            }
            
            // for x1,...,xn-1
            else if ((j > 0) && (j < n))
            {
                for (i = 0; i < j; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m - 1]);
                }
                for (i = j + 1; i <= n; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m - 1]);
                }
                x[j][m] = (1/A[j][j])*(b[j] - sum[j]);
            }
            
            // for xn
            else if ((j = n))
            {
                for (i = 0; i < n; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m - 1]);
                }
                x[j][m] = (1/A[j][j])*(b[j] - sum[j]);
            }
        }
        
        // Calculate two norm of error vector (X[m]-X[m-1])
        for(xnorm=0,j=0;j<=n;j=j+1)
        {
            xnorm = xnorm + pow(x[j][m]-x[j][m-1],2);
        }
        xnorm = sqrt(xnorm);
        
        // Print results of current iterate
        cout << m;
        for(j=0;j<=n;j=j+1)
        {
            cout << setw(w) << x[j][m];
        }
        cout << setw(w) << xnorm << endl;
    } while ((xnorm > errtol) && (m<1000));
    
    // Format final results
    cout << endl << "x[" << m << "] that solves [A.x = b] = {";
    for (j=0;j<=n;j=j+1)
    {
        cout << x[j][m];
        if (j<n)
        {
            cout << ", ";
        }
    }
    cout << "}T" << endl;
    
    // Print info
    cout << endl << "Copyright © 2016 Kyle James Hansen. All rights reserved. ";
}


sample output:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 
Solve system [A.x = b] using Jacobi Iteration

Enter 'N', the size of the system:
2

Enter coefficients of 2*2 system:
1 2
3 4

Enter 2 constants of RHS of system:
3
5

System:
1    2 | 3
3    4 | 5

m             x0             x1         ||X||2
0              1              2              -
1              0            0.5        1.80278
2              0            0.5              0

x[2] that solves [A.x = b] = {0, 0.5}T

Copyright © 2016 Kyle James Hansen. All rights reserved. Program ended with exit code: 0


-> this should return x={-1,2}T, which it does if I replace all the if statements with manual formulas for x0 and x1...

Thanks in advance for any help!
Last edited on
Chervil (7320)
There's a whole lot of stuff there which I don't understand, so I will only comment on the most obvious.

1
2
 
// for x0
if ((j = 0))


1
2
 
// for xn
else if ((j = n))


Here = is the assignment operator. But it looks like it should be the == operator which tests for equality.

There are a few other things in the code which look a little odd, such as the abs() function (which is never used), and the use of variable-length arrays which are not valid in standard C++, I assume you use a compiler feature to permit that.

Also, please use code tags for legibility:
http://www.cplusplus.com/articles/jEywvCM9/
hans4354 (2)
Reformatted per the link (did try when I posted but the button didn't do anything, also told me my original post was zero characters long).

I was in a programming-oriented class this semester and used the abs and chop fcn's so often I just always define them now.

Will try that fix in a bit and get back to you, thanks!

Edit:
@Chervil 's fix and changing it to a Gauss-Seidel iterative method instead of the original Jacobi worked! After the 
1
2
 
// for xn
else if ((j==n))
fix the program converged on correct values for x0 and x1 but not on the same iteration so the two norm did not have a limit of zero. Changed it Gauss-Seidel using the following:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 
// Iterate one row at a time (x0,x1,...,xn)
        for (j = 0; j <= n; j = j + 1)
        {
            sum[j] = 0;
            
            // for x0
            if ((j == 0))
            {
                for (i = 1; i <= n; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m - 1]);
                }
                x[j][m] = (1/A[j][j])*(b[j] - sum[j]);
            }
            
            // for x1,...,xn-1
            else if ((j > 0) && (j < n))
            {
                for (i = 0; i < j; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m]);
                }
                for (i = j + 1; i <= n; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m - 1]);
                }
                x[j][m] = (1/A[j][j])*(b[j] - sum[j]);
            }
            
            // for xn
            else if ((j == n))
            {
                for (i = 0; i < n; i++)
                {
                    sum[j] = sum[j] + (A[j][i]*x[i][m]);
                }
                x[j][m] = (1/A[j][j])*(b[j] - sum[j]);
            }
        }

which generally converges faster anyways.

Thanks for the help!
Last edited on
Chervil (7320)
I'm glad that worked.

I can now comment on something which was not directly related, but may be useful for future reference.

This function is heavy on the amount of calculation done:
1
2
3
4
 
double abs(double ABS)
{
    return sqrt(ABS*ABS);
}

all that is required is to reverse the sign if the value is negative:
1
2
3
4
 
double abs(double ABS)
{
    return (ABS < 0.0) ? -ABS : ABS;
}

Above, the ternary operator is just a concise way to write an if-else statement.

Also calculation-intensive, pow(10,-15) can be written directly in scientific notation as 1e-15, same with double errtol = 1e-10;


doug4 (1538)
Better than defining the double abs(double ABS) function would be to use the standard library function double fabs(double x) found in the header <cmath>
Topic archived. No new replies allowed.