Topics

• Background
• Methods
• Results
• Extra Credit search for a convergent matrix
• Movie

Background

Systems of linear equations appear in a wide variety of fields -- from science and engineering to social sciences. From determining electrical current through a system to figuring reproduction rates of bacteria, linear equations provide useful tools. Various techniques have arisen to solve for the unknowns of these systems. For a system of equations:

E_1: a_11(x_1) + a_12(x_2) + ... + a_1n(x_n) = b_1.

E_2: a_21(x_1) + a_22(x_2) + ... + a_2n(x_n) = b_2.

  .     .          .            .        .

  .     .          .            .        .

  .     .          .            .        .

when m<=n, the system can be solved directly. Pehaps the most common method for this is Gaussian elimination with backward subsitution. In this method, augmented matrix[A,b] is constructed, such that

                       a_11     a_12    ...     a_1n    :    b_1
                       a_21     a_22    ...     a_2n    :    b_2
              [A,b] =   .        .                .           .
                        .        .                .           .
                        .        .                .           .
                       a_m1     a_m2    ...     a_mn    :    b_n

While the direct method for solving matrices guarantees solutions, provided all the calculations are performed correctly, this method is not aways feasible. For instance in 1000x 1000 matrices, the chances for making an error are high, and many computers cannot store matrices of such size, making computer manipulation impossible. However, problems in the real world often produce such large matrices. For this reason, various iterative methods have been developed. In these methods, initial values are estimated, and successive iterations of the method produce improved results. Or, such is the hope.

In this project, we looked at the Jacobi iterative method. This method splits A into three matrices: the diagonal (D), an upper triangular (U), and a lower triangular (L), such that Diis the same as the diagonal of A, -U is equal to the upper triangular part of A, and -L is equal to the lower triangular part of A. Therefore,

     A = D - L - U.

     Ax = b

     (D - L - U)x = b.

     Dx = (L + U)x + b

     x = d^-1 (L + U)x + D^-1b,

     x^(k) = D^-1 (L + U) x^(k - 1) + D^-1b,    k = 1, 2,...

Methods: In an attempt to solve the given matrix by the Jacobi method, we used the following two programs:

function y = jacobi(S,b,N)
    %This function performs the Jacobi iterative on the (sparse) matrix S, to solve the
    system Sx = b, with N iterations.  Initial vector is X_0.
    n = size(S,1);
    for m = 1:n
     X_0(m) = 0;
    end
    iter = 1;
    while iter <=N
     sum = 0;
     for m = 1:n
      for k = 1:n
       if k~=m
        sum = sum + (S(m,k)*X_0));
      end
     end
      X(m) = (-sum + b(m))/S(m,m);
     end
      iter = iter + 1;
     for m = 1:n
      X_0(m) = X(m);
     end
    end
    y = X_0;

clear all
    A(1,:) = [10 -1 2 0];
    A(2,:) = [-1 11 -1 3];
    A(3,:) = [2 -1 10 -1];
    A(4,:) = [0 3 -1 8];
    b = [6 25 -11 15]';

    x = [0,0,0,0]';

    d = diag(diag(A));
    up = triu(A) - d;
    lo = tril(A) - d;

    for i = 1:4
     x = (-(up + lo)/d)'*x + 1./diag(A).*b
     r = b - A * x
    end

Results

From our data, we have seen that while the Jacobi method may successfully solve strictly diagonally dominant matrices, the rate at which it does so may make the method virtually ineffective in some cases. As shown, when we ran the nos_6 matrix with our program over night, the x values were hardly scratched. The rate of convergence is so slow that between the tenth and the hundred thousandth iteration, the norms of the vectors are close enough that Matlab shows their difference to be zero. Matlab can solve the same matrix in only a matter of moments when allowed to use its own devices. This problem, therefore provides a good example of a problem which is possible solve in the theoretical sense, but it is difficult to solve practically.

The output of the function jaco1.m is the infinity norm of x at iteration k. We denote this by [x]_k. Our results are as follows:

[x]_10 = jaco1(S,b,10) = 9.996002098750798e01

[x]_100 = jaco1(S,b,100) = 9.996002098750798e01

[x]_1000 = jaco1(S,b,1000) = 9.996002098750798e01

[x]_100000 = jaco1(S,b,100000) = 9.996002098750798e01

The difference between [x]_10 and [x]_100000 is nearly zero.

Why is the rate of convergence so slow for this matrix?

By corollary 7.20 from the text, if [T]<1 for any natural matrix norm and c is a given vector, then the sequence ({x(k)}_x)^inf = 0 defined by x^(k) = Tx(k-1) + c converges, for any x^(0) in the set of reals of dimension n, to a vector x also in the set of reals of dimension n, and the following error bound holds:

[-]<=°

T is the matrix given by inv(D)*(L + U), and [ ] is a natural matrix norm ( for us, the infinity norm). In the case of the nos6_mtx matrix, the quantity on the right equals:

(0.9999997500000626)^k*(3992002.55776)

As the right hand side of the equation shows, the error bound for the kth iteration is extremely large. Since the norm of the T matrix is so close to 1, this quantity shrinks to zero very slowly (but it does go to zero. Thus, the difference between the actual answer and the kth iteration will approach each other very slowly.