6-4 Compare Methods of Jacobi with Gauss-Seidel (50分)

Use Jacobi and Gauss-Seidel methods to solve a given n×n linear system A
​x​​ =​b​⃗​​ with an initial approximation ​x​⃗​​​(0)​​.
Note: When checking each a​ii​​ , first scan downward for the entry with maximum absolute value (a​ii​​ included). If that entry is non-zero, swap it to the diagonal. Otherwise if that entry is zero, scan upward for the entry with maximum absolute value. If that entry is non-zero, then add that row to the i-th row.

Format of functions:

int Jacobi( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );
int Gauss_Seidel( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );

where int n is the size of the matrix of which the entries are in the array double a[][MAX_SIZE] and MAX_SIZE is a constant defined by the judge; double b[] is for ​b​⃗​​ , double x[] passes in the initial approximation ​x​⃗​​(0)​​ and return the solution ​x​⃗​​ ; double TOL is the tolerance for ∥​x​⃗​​(k+1)​​ −​x​⃗​​(k)​​ ∥​∞​​ ; and finally int MAXN is the maximum number of iterations.
The function must return:
k if there is a solution found after k iterations;
0 if maximum number of iterations exceeded;
1 if the matrix has a zero column and hence no unique solution exists;
2 if there is no convergence, that is, there is an entry of
​x​⃗ ​(K)​​ that is out of the range [-bound, bound] where bound is a constant defined by the judge.
Sample program of judge:

#include <stdio.h>
#include <math.h>

#define MAX_SIZE 10
#define bound pow(2, 127)
#define ZERO 1e-9 /* X is considered to be 0 if |X|<ZERO */

enum bool { false = 0, true = 1 };
#define bool enum bool

int Jacobi( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );

int Gauss_Seidel( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );

int main()
{
    int n, MAXN, i, j, k;
    double a[MAX_SIZE][MAX_SIZE], b[MAX_SIZE], x[MAX_SIZE];
    double TOL;

    scanf("%d", &n);
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++)
            scanf("%lf", &a[i][j]);
        scanf("%lf", &b[i]);
    }
    scanf("%lf %d", &TOL, &MAXN);

    printf("Result of Jacobi method:\n");
    for ( i=0; i<n; i++ )
        x[i] = 0.0;
    k = Jacobi( n, a, b, x, TOL, MAXN );
    switch ( k ) {
        case -2:
            printf("No convergence.\n");
            break;
        case -1: 
            printf("Matrix has a zero column.  No unique solution exists.\n");
            break;
        case 0:
            printf("Maximum number of iterations exceeded.\n");
            break;
        default:
            printf("no_iteration = %d\n", k);
            for ( j=0; j<n; j++ )
                printf("%.8f\n", x[j]);
            break;
    }
    printf("Result of Gauss-Seidel method:\n");
    for ( i=0; i<n; i++ )
        x[i] = 0.0;
    k = Gauss_Seidel( n, a, b, x, TOL, MAXN );
    switch ( k ) {
        case -2:
            printf("No convergence.\n");
            break;
        case -1: 
            printf("Matrix has a zero column.  No unique solution exists.\n");
            break;
        case 0:
            printf("Maximum number of iterations exceeded.\n");
            break;
        default:
            printf("no_iteration = %d\n", k);
            for ( j=0; j<n; j++ )
                printf("%.8f\n", x[j]);
            break;
    }

    return 0;
}
/* Your function will be put here */

Sample Input 1:
3
2 -1 1 -1
2 2 2 4
-1 -1 2 -5
0.000000001 1000

Sample Output 1:
Result of Jacobi method:
No convergence.
Result of Gauss-Seidel method:
no_iteration = 37
1.00000000
2.00000000
-1.00000000

Sample Input 2:
3
1 2 -2 7
1 1 1 2
2 2 1 5
0.000000001 1000

Sample Output 2:
Result of Jacobi method:
no_iteration = 232
1.00000000
2.00000000
-1.00000000
Result of Gauss-Seidel method:
No convergence.

Sample Input 3:
5
2 1 0 0 0 1
1 2 1 0 0 1
0 1 2 1 0 1
0 0 1 2 1 1
0 0 0 1 2 1
0.000000001 100
Sample Output 3:
Result of Jacobi method:
Maximum number of iterations exceeded.
Result of Gauss-Seidel method:
no_iteration = 65
0.50000000
0.00000000
0.50000000
0.00000000
0.50000000

解法：

1.先对输入的矩阵进行变换，从第一行开始，第i行找$\underset{i\le j\le n}{\max }(a_{j,i})$，如果不是0，则换到第i行。
这里注意b也要跟着换。

2.如果是0，则找$\underset{1\le j\le i-1}{\max }(a_{j,i})$,如果不是0则加到第i行，否则返回-1

3.算法按书上的伪代码写，不论是雅可比法还是高斯赛德迭代法，都应该另设一个数组x2，因为需要比较两次答案之间的误差。

4.如果达到的条件是最大的误差小于TOL，必须确保当时的x是已经迭代过的新值（x是传进来的，并不需要返回，但是被主函数直接引用）

测试点：

倒数第二个返回0 0
最后一个返回-1 -1（n=15）

总结

总而言之并不是很难，但是得仔细，小心。

上代码

int _n;
double _a[MAX_SIZE][MAX_SIZE];
double _b[MAX_SIZE];

int changeline(){
	for (int i = 0; i < _n; i++){
        double max = fabs(_a[i][i]);
        int index = i;
        for (int j = i; j < _n; j++){
            if (fabs(_a[j][i]) > max){
                max = fabs(_a[j][i]);
                index = j;
            }
        }
        if(max == 0){
            for (int j = 0; j < i; j++){
                if (fabs(_a[j][i]) > max){
                    max = fabs(_a[j][i]);
                    index = j;
                }
            }
            if(max==0) return -1;//可能是这里
            for(int j = 0;j < _n; j++) 
                _a[i][j]+=-_a[index][j];
            _b[i] += _b[index];
        }else{
            if(index == i) continue;//如果aii已经是i列最大元素，继续 
		    for (int k = 0; k < _n; k++){
        	    double temp = _a[i][k]; _a[i][k] = _a[index][k]; _a[index][k] = temp;
    	    }
    	    double temp = _b[i] ;_b[i] = _b[index]; _b[index] = temp;
	       }
        }
}

int Jacobi(int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN){
	_n = n;
	for(int i = 0; i < n; i++){
		for(int j = 0; j < n; j++)
			_a[i][j] = a[i][j];
		_b[i] = b[i];
	}
	if(changeline()==-1) return -1;
    int k = 1;
    double x2[MAX_SIZE];
    while (k <= MAXN){//迭代 
        for (int i = 0; i < n; i++){//行 
            double sum = 0;
            for(int j = 0; j < n; j++){//每行的所有元素 
                if (j == i) continue;
                sum += x[j] * _a[i][j];
            }
            x2[i] = (_b[i] - sum) / _a[i][i];
		}
		double maxerror = fabs(x2[0] - x[0]);
        for (int j = 0; j < n; j++){
            if (fabs(x2[j] - x[j]) > maxerror) 
				maxerror = fabs(x2[j]-x[j]);
			if (fabs(x[j]) > bound)  return -2;
			x[j] = x2[j];  				
        }
        if(maxerror < TOL) return k; 
        k++;
    }
	return 0;
}

int Gauss_Seidel(int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN){
	_n = n;
	for(int i=0;i<n;i++){
		for(int j=0;j<n;j++)
			_a[i][j] = a[i][j];
		_b[i] = b[i];
	}
    
	if(changeline()==-1) return -1;//这里出去。
    int k = 1;
    double x2[MAX_SIZE];
    while ( k <= MAXN ){//迭代 
        for (int i = 0; i < _n; i++){//行 
            double sum = 0;
            for(int j = 0; j < i; j++)//每行的所有元素 
                sum += x2[j] * _a[i][j];
            for(int j=i+1;j<n;j++)
            	sum += x[j] * _a[i][j];
            x2[i] = (_b[i] - sum) / _a[i][i];
		}
		double maxerror = fabs(x2[0] - x[0]);
        for (int j = 0; j < _n; j++){
            if (fabs(x2[j] - x[j]) > maxerror) 
				maxerror = fabs(x2[j] - x[j]);
			if (fabs(x2[j] > bound)) return -2;
			x[j] = x2[j];
        }
        if(maxerror < TOL) return k;
        k++;
    }
    return 0;
}