NumericalMethods[MatrixIterate] - approximates solution of a system of linear equations by iteration

Calling Sequence
    MatrixIterate(A, b, xinit=v, method=m, tolerance=t, options)
    MatrixIterate(A, method=m, output=matrix)

Parameters
    A - Matrix; a square numerical matrix of the system
    b - Vector; right-hand side of the system
    v - Vector; initial approximate eigenvector of appropriate dimension
    m - name; method to be used
    t - positive; tolerance for the stopping condition
    options - (optional) parameters of the form keyword = value, where keyword is one of stoppingcriterion, output, maxiterations, maxnumber, digits.

Description
The MatrixIterate command approximates the solution of a system of linear equations Ax = b, where diagonal terms of A are not zero, using iterative methods: Jacobi, Gauss-Seidel or SOR (relaxation). They all fit the general pattern x_k+1 = Bx_k + c for a certain matrix B and a shift vector c, and the procedure MatrixIterate can also provide B and c for these methods. Its main purpose is to exhibit behaviour of basic methods and be accessible to students not versed in Maple. Acordingly, the main focus is on ease of use and informative output, not performace. Controls are similar to other iterative procedures in the package.

There is no default method or default tolerance, they must be user specified.

It is assumed that the package LinearAlgebra has ben loaded before this procedure is read.

This command is part of the NumericalMethods package, so it can be used in the form MatrixIterate(..) only after executing the command with(NumericalMethods). However, it can always be accessed through the long form of the command by using NumericalMethods[MatrixIterate](..).

Options

method = jacobi, gaussseidel, or [sor,numeric]
- jacobi: the Jacobi iteration method (JIM):
(x_k+1)_i = b_i /a_i,i - (1/a_i,i)[a_i,1(x_k)₁ + ... + a_{i, j-1}(x_k)_j-1 + a_{i, j+1}(x_k)_j+1 + ... + a_i,n(x_k)_n ].
- gaussseidel: the Gauss-Seidel iteration method (GSM):
(x_k+1)_i = b_i /a_i,i - (1/a_i,i)[a_i,1(x_k+1)₁ + ... + a_{i, j-1}(x_k+1)_j-1 + a_{i, j+1}(x_k)_j+1 + ... + a_i,n(x_k)_n ].
- [sor,omega], where omega = numeric: the Successive over-relaxation method (SOR):
x_k+1 = omega*GS(x_k ) + (1 - omega)*x_k , where by GS(x_k ) we denote the outcome of Gauss-Seidel iteration applied to vector x_k.

The options argument can be one or more of the following equations.

stoppingcriterion = relative, absolute, or residual
The criterion that the iteration must meet to stop (max-norm is used).
- relative: |x_k - x_k-1| / |x_k| < tolerance
- absolute: |x_k - x_k-1| < tolerance
- residual: |b - Ax_k| < tolerance
The default value is stoppingcriterion = relative.

output = information, iterates, value, matrix, or matrices
Specifies what information is given by the procedure. In all cases, the iteration index k is printed at every iteration step, so that the user can see the progress and knows how many iterations it took to produce the final approximation.
- information returns the final numerical approximation of the solution. During iterations, data are shown at each step: the generated solution approximation, the residual of the current approximation and value test used in stoppingcriterion to determine whether the algorithm should stop.
- iterates returns the final numerical approximation of the solution. During iterations, generated solution approximations are shown.
- value returns the final numerical approximation of the solution.
- matrix: returns the matrix B corresponding to the iterative formula  x_k+1 = Bx_k + c of the chosen method. For this output option, a shortened form of a procedure call is provided so that irrelevant data need not be inputed.
- matrices: returns the matrix B and vector c corresponding to the iterative formula  x_k+1 = Bx_k + c of the chosen method. There is no shorter form provided.
The default value is output = information if matrix dimension is less than 5 for real valued approximations and less than 3 for complex valued approximations, and output = iterates if matrix dimension is less than 7 for real valued approximations and less than 4 for complex valued approximations. Otherwise the default is output = value.

maxiterations = posint
The maximum number of iterations to perform. If it is reached, the algorithm is interrupted and a message is displayed. The return as defined by output is still produced.
The default value is maxiterations = 50.

maxnumber = positive
The maximal number that is not to be exceeded by vector entries during iterations. If it does happen, the algorithm is interrupted and a message is displayed. The return as defined by output is still produced.
The default value is maxnumber = 1000000.

digits = posint
Specifies how many decimal places should be shown when output = information or output = iterates is specified.
If digits is set to a value greater than Digits, the latter is set to digits for the duration of this procedure.
The default value is digits = Digits.

Examples

> with(LinearAlgebra):
  A:=<<4,-1,1>|<0,5,-1>|<2,1,3>>;
  b:=A.<1,-1,2>;

The matrix is strictly diagonally dominant, so iteration should yield the solution [1,-1,2]. We start with the Jacobi iteration method. Since for such small matrices the default output is information, we restrict the number of digits shown.
> MatrixIterate(A, b, xinit=<1.3,1.4,-2.3>, method=jacobi, tolerance=0.001, digits=5);
k=00  x=[  1.30000,  1.40000, -2.30000]
k=01  x=[  3.15000, -0.08000,  2.70000],   res= 10.00000   test= 1.58730,
k=02  x=[  0.65000, -0.71000,  1.59000],   res= 2.22000   test= 1.57233,
k=03  x=[  1.20500, -0.98800,  2.21333],   res= 1.24667   test= 0.28163,
k=04  x=[  0.89333, -1.00167,  1.93567],   res= 0.55533   test= 0.16101,
k=05  x=[  1.03217, -1.00847,  2.03500],   res= 0.19867   test= 0.06822,
k=06  x=[  0.98250, -1.00057,  1.98646],   res= 0.09709   test= 0.02500,
k=07  x=[  1.00677, -1.00079,  2.00564],   res= 0.03838   test= 0.01210,
k=08  x=[  0.99718, -0.99977,  1.99748],   res= 0.01633   test= 0.00480,
k=09  x=[  1.00126, -1.00006,  2.00102],   res= 0.00707   test= 0.00204,
k=10  x=[  0.99949, -0.99995,  1.99956],   res= 0.00291   test= 0.00088,

The Gauss-Seidel iteration should be faster. We restrict output to iterates. Note that the residual of the last approximation above is larger than tolerance, so we try to use it as stoppingcriterion.
> MatrixIterate(A, b, xinit=<1.3,1.4,-2.3>, method=gaussseidel, tolerance=0.001, stoppingcriterion=residual, output=iterates);
k=00  x=[  1.3000000000,  1.4000000000, -2.3000000000]
k=01  x=[  3.1500000000,  0.2900000000,  1.7133333330]
k=02  x=[  1.1433333340, -0.9139999998,  1.9808888890]
k=03  x=[  1.0095555560, -0.9942666666,  1.9987259260]
k=04  x=[  1.0006370370, -0.9996177778,  1.9999150620]
k=05  x=[  1.0000424690, -0.9999745186,  1.9999943370]

Of course, for badly conditioned matrices, residual is not a good guide.

> A:=<<1,-1,1,0,1>|<1,-1,3,-1,1>|<-1,1,1,2,1>|<1,-1,0,1,-1>|<1,5,-1,1,3>>:
  b:=A.<1,-1,0,-1,1>:
For systems with dimension 5 and 6, output = iterates is the default.
> MatrixIterate(A, b, xinit=<1.,-1.3,1.4,1.3,-1.4>, method=gaussseidel, tolerance=0.001);
k=00  x=[  1.0000000000, -1.3000000000,  1.4000000000,  1.3000000000, -1.4000000000]
k=01  x=[  2.8000000000,-15.7000000000, 39.9000000000,-93.1000000000,-38.7000000000]
k=02  x=[ 187.4000000000,-253.9000000000, 532.6000000000,-1279.4000000000,-580.5000000000]
k=03  x=[ 2646.4000000000,-3742.9000000000, 7998.8000000000,-19159.0000000000,-8685.7666670000]
k=04  x=[ 39586.4666700000,-55863.5000100000, 119315.2666000000,-285807.2665000000,-129613.8333000000]
k=05 Numbers in iteration too large.

This does not look good, let's check on this iteration closer. We add output parameter to get the iteration matrix.
> B:=MatrixIterate(A, b, xinit=<1.,-1.3,1.4,1.3,-1.4>, method=gaussseidel, tolerance=0.001, output=matrix);
  SpectralRadius(B);

The spectral radius  larger than 1, which explains the divergence. The shorter form would work as well.
> B:=MatrixIterate(A, method=gaussseidel, output=matrix);

Let's go back to the first matrix.
> A:=<<4,-1,1>|<0,5,-1>|<2,1,3>>:
  b:=A.<1,-1,2>:
The Gauss-Seidel method required five iterations. A bit of experimentation shows that SOR can improve on this. We try another output option.
> MatrixIterate(A, b, xinit=<1.3,1.4,-2.3>, method=[sor,21/20], tolerance=0.001, output=value);
k=00 k=01 k=02 k=03 k=04

Now we will try it manually.
> (B,c):=MatrixIterate(A, b, xinit=<1.3,1.4,-2.3>, method=[sor,21/20], tolerance=0.001, output=matrices);

> y:=<1.3,1.4,-2.3>:
  for i from 1 to 4 do
   y:=evalf(B.y+c):
  end do:
  y;

This matches our previous result quite well. Note that matrix calculations are automatically done with higher precision than the one specified by Digits when the LinearAlgebra package is used.

The last example will show that this procedure can handle complex iterations as well. For 3x3 matrices and complex iteration, the default output is iterates. We will ask for a smaller tolerance and adjust digits to it.
> MatrixIterate(A, b, xinit=<I,-1.3,1.4-I>, method=gaussseidel, tolerance=0.00001, digits=5);
k=00  x=[  0.00000 +1.00000I, -1.30000 +0.00000I,  1.40000 -1.00000I]
k=01  x=[  1.30000 +0.50000I, -0.82000 +0.30000I,  1.96000 -0.06667I]
k=02  x=[  1.02000 +0.03333I, -0.98800 +0.02000I,  1.99733 -0.00444I]
k=03  x=[  1.00133 +0.00222I, -0.99920 +0.00133I,  1.99982 -0.00030I]
k=04  x=[  1.00009 +0.00015I, -0.99995 +0.00009I,  1.99999 -0.00002I]
k=05  x=[  1.00001 +0.00001I, -1.00000 +0.00001I,  2.00000 -0.00000I]
k=06  x=[  1.00000 +0.00000I, -1.00000 +0.00000I,  2.00000 -0.00000I]