Chapter 7. Iterative Techniques in Matrix Algebra Sample Program in the Tex #include /* Jacobi iterative method */ void matvec(double a[][], double x[], double temp[], int n ) main() { int i, j, iter, n = 4; double x[4], temp[4] ; double r0, rn, eps, numer1, numer2 ; /* A matrix and b vector */ double b[4] = {6, 25, -11, 15}; double a[4][4] = {10, -1, 2, 0, -1, 11, -1, 3, 2, -1, 10, -1, 0, 3, -1, 8 };
;
/* Initialize x to 0 */ for (i=0; i < n; i++) { x[i] =0. ; } /* Compute R0 = b - A x0 */ matvec(a, x, temp, n) ; for (i = 0; i < n ; i++) { temp[i] = b[i] - temp[i]; } /* Compute sqrt(R0) */ r0 = 0. ; for (i = 0 ; i < n ; i++) { r0 = r0 + temp[i]*temp[i]; } r0 = sqrt(r0); rn = r0 ; iter = 0; eps = 1.e-8 ; while ((rn/r0) > eps ) { for (i=0; i < n ; i++) * Compute sum(a(i,j)*x(j), j=1, i-1) * numer1 = 0.; j = 0 ; while (j < i ) { numer1 = numer1 + a[i][j]*x[j]; j++; }
1
/* Compute sum(a(i,j)*x(j), j=i+1, n) */ numer2 = 0.; j = i+1 ; while (j < n ) { numer2 = numer2 + a[i][j]*x[j]; j++; } temp[i] = (-numer1 - numer2 + b[i] )/a[i][i];
} for (i=0; i < n; i++) { x[i] = temp[i]; } /* Compute Rn = n - A Xn */ matvec(a, x, temp, n); rn = 0. ; for (i=0; i < n; i++) { rn = rn + (b[i]-temp[i])*(b[i]-temp[i]); } rn = sqrt(rn); iter++; printf("At iter = %d residual =%f %f %f\n", iter, r0, rn, rn/ printf(" X = %f %f %f %f\n", x[0], x[1], x[2], x[3]); } }
void matvec(double a[][4], double x[], double temp[], int n ) /* Compute temp = A x */ { int i, j ; for (i=0; i < n; i++) { temp[i] = 0. ; for (j=0; j < n; j++) {
2
temp[i] = temp[i] + a[i][j]*x[j] ; } } }
3
Outputs with eps = 1.e-8 At X At X At X At X At X At X At X At X At X At X At X At X At X At X At X At X At iter = 1 residual =31.733263 11.353749 0.357787 = 0.600000 2.272727 -1.100000 1.875000 iter = 2 residual =31.733263 4.990955 0.157278 = 1.047273 1.715909 -0.805227 0.885227 iter = 3 residual =31.733263 2.029878 0.063967 = 0.932636 2.053306 -1.049341 1.130881 iter = 4 residual =31.733263 0.891141 0.028082 = 1.015199 1.953696 -0.968109 0.973843 iter = 5 residual =31.733263 0.368628 0.011616 = 0.988991 2.011415 -1.010286 1.021351 iter = 6 residual =31.733263 0.160473 0.005057 = 1.003199 1.992241 -0.994522 0.994434 iter = 7 residual =31.733263 0.067106 0.002115 = 0.998128 2.002307 -1.001972 1.003594 iter = 8 residual =31.733263 0.029022 0.000915 = 1.000625 1.998670 -0.999036 0.998888 iter = 9 residual =31.733263 0.012221 0.000385 = 0.999674 2.000448 -1.000369 1.000619 iter = 10 residual =31.733263 0.005261 0.000166 = 1.000119 1.999768 -0.999828 0.999786 iter = 11 residual =31.733263 0.002225 0.000070 = 0.999942 2.000085 -1.000068 1.000109 iter = 12 residual =31.733263 0.000955 0.000030 = 1.000022 1.999959 -0.999969 0.999960 iter = 13 residual =31.733263 0.000405 0.000013 = 0.999990 2.000016 -1.000013 1.000019 iter = 14 residual =31.733263 0.000173 0.000005 = 1.000004 1.999993 -0.999994 0.999992 iter = 15 residual =31.733263 0.000074 0.000002 = 0.999998 2.000003 -1.000002 1.000003 iter = 16 residual =31.733263 0.000032 0.000001 = 1.000001 1.999999 -0.999999 0.999999 iter = 17 residual =31.733263 0.000013 0.000000
4
X At X At X At X At X At X
= 1.000000 2.000001 -1.000000 iter = 18 residual =31.733263 = 1.000000 2.000000 -1.000000 iter = 19 residual =31.733263 = 1.000000 2.000000 -1.000000 iter = 20 residual =31.733263 = 1.000000 2.000000 -1.000000 iter = 21 residual =31.733263 = 1.000000 2.000000 -1.000000 iter = 22 residual =31.733263 = 1.000000 2.000000 -1.000000
1.000001 0.000006 1.000000 0.000002 1.000000 0.000001 1.000000 0.000000 1.000000 0.000000 1.000000
0.000000 0.000000 0.000000 0.000000 0.000000
5
Gauss-Seidel Code #include /* Gauss-Seidel Iterative method */ void matvec(double a[][], double x[], double temp[], int n ) main() { int i, j, iter, n = 4; double x[4], temp[4] ; double r0, rn, eps, numer1, numer2 ; /* A matrix and b vector */ double b[4] = {6, 25, -11, 15}; double a[4][4] = {10, -1, 2, 0, -1, 11, -1, 3, 2, -1, 10, -1, 0, 3, -1, 8 }; /* Initialize x to 0 */ for (i=0; i < n; i++) { x[i] =0. ; } /* Compute R0 = b - A x0 */ matvec(a, x, temp, n) ; for (i = 0; i < n ; i++) { temp[i] = b[i] - temp[i]; } /* Compute sqrt(R0) */ r0 = 0. ; for (i = 0 ; i < n ; i++) { r0 = r0 + temp[i]*temp[i]; } r0 = sqrt(r0); rn = r0 ; iter = 0; eps = 1.e-8 ; while ((rn/r0) > eps ) { for (i=0; i < n ; i++) { numer1 = 0.; /* Compute sum(a(i,j)*x(j), j=1, i-1) */ j = 0 ; while (j < i ) { numer1 = numer1 + a[i][j]*x[j];
6
;
j++; } * Compute sum(a(i,j)*x(j), j=i+1, n) */ numer2 = 0.; j = i+1 ; while (j < n ) { numer2 = numer2 + a[i][j]*x[j]; j++; }
/* Here, temp is replaced by x. */ x[i] = (-numer1 - numer2 + b[i] )/a[i][i]; } /* Compute Rn = n - A Xn */ matvec(a, x, temp, n); rn = 0. ; for (i=0; i < n; i++) { rn = rn + (b[i]-temp[i])*(b[i]-temp[i]); } rn = sqrt(rn); iter++; printf("At iter = %d residual =%f %f %f\n", iter, r0, rn, rn/ printf(" X = %f %f %f %f\n", x[0], x[1], x[2], x[3]); } }
void matvec(double a[][4], double x[], double temp[], int n ) /* Compute temp = A x */ { int i, j ; for (i=0; i < n; i++) { temp[i] = 0. ; for (j=0; j < n; j++) { temp[i] = temp[i] + a[i][j]*x[j] ;
7
} } }
8
Outpus with eps = 10**(-8) At X At X At X At X At X At X At X At X At X iter = 1 residual =31.733263 5.693016 0.179402 = 0.600000 2.327273 -0.987273 0.878864 iter = 2 residual =31.733263 0.429975 0.013550 = 1.030182 2.036938 -1.014456 0.984341 iter = 3 residual =31.733263 0.066172 0.002085 = 1.006585 2.003555 -1.002527 0.998351 iter = 4 residual =31.733263 0.008165 0.000257 = 1.000861 2.000298 -1.000307 0.999850 iter = 5 residual =31.733263 0.000852 0.000027 = 1.000091 2.000021 -1.000031 0.999988 iter = 6 residual =31.733263 0.000078 0.000002 = 1.000008 2.000001 -1.000003 0.999999 iter = 7 residual =31.733263 0.000006 0.000000 = 1.000001 2.000000 -1.000000 1.000000 iter = 8 residual =31.733263 0.000000 0.000000 = 1.000000 2.000000 -1.000000 1.000000 iter = 9 residual =31.733263 0.000000 0.000000 = 1.000000 2.000000 -1.000000 1.000000
Jacobi required 22 iterations, while Gauss-Seidel needed only 9 iterations.
9