NumericalMethods[MatrixEliminate] - applies elimination to a matrix

Calling Sequence
    MatrixEliminate(A,options)
    MatrixEliminate(<A|b>,options)
    MatrixEliminate(<A|B>,options)

Parameters
    A - Matrix; a numerical matrix
    b - Vector; typically right-hand side of the system
    B - Matrix; a numerical matrix
    options - (optional) parameters of the form keyword = value, where keyword is one of output, pivoting.

Description
The MatrixEliminate command applies elimination to a matrix A. Several variations of basic elimination are available, controlled by the choice of pivoting strategy and desired output. 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 interesting options, not performace.

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 MatrixEliminate(..) only after executing the command with(NumericalMethods). However, it can always be accessed through the long form of the command by using NumericalMethods[MatrixEliminate](..).

Options

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

output = gem, gaussjordan, lu, or lup
Specifies the return of the procedure.
- gem returns a matrix in a row echelon form. The standard Gaussian elimination is used. This is the forced output when Euclidean elimination is used.
- gaussjordan returns a matrix in a reduced row echelon form. The standard Gaussian elimination is used.
- lu returns matrices L and U forming the LU decomposition A = LU, if it exists. With Euclidean elimination this option is not available.
- lup returns matrices L, U, and P forming the LUP decomposition, that is, PA = LU. With Euclidean elimination this option is not available.
The default value is output = gem.

pivoting = partial, minimal, smallpivot, or euclidean
Specifies the pivoting strategy applied by the elimination, in one particular case it influences the whole elimination process.
- partial is the standard partial pivoting elimination: From the active column, starting from the pivoting position down, the largest number (in absolute value) is used. This is the Maple's GaussianElimination default for floating point matrices.
- minimal is the minimal effort pivoting: If there is a zero at the pivot place, we choose the first non-zero candidate from the active column, starting from the pivoting position down. This is the Maple's GaussianElimination default for integer and fractional matrices.
- smallpivot: From the active column, starting from the pivoting position down, the smallest (in absolute value) non-zero number is used.
- euclidean switches to Euclidean elimination. For the active column, starting from the pivoting position down, the Euclidean algorithm is applied, finding the greatest common divisor of these numbers, leaving it at the pivoting position. This elimination can be used to solve systems of linear Diophantic equations.
The default value is pivoting = partial.

Examples

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

First we try to obtain the solution [1,-1,2] using Gaussian elimination and backward substitution.
> Ael:=MatrixEliminate(<A|b>);
  BackSubstitute(Ael);

Next we try Gauss-Jordan elimination.
> Ael:=MatrixEliminate(<A|b>, output=gaussjordan);

When eliminating by hand, we usually prefer smaller pivots.
> Ael:=MatrixEliminate(<A|b>, pivoting=smallpivot);

The LU decomposition of A can also be found easily. We check that it works.
> (L,U,P):=MatrixEliminate(A, output=lup);
  L.U-P.A;

> A:=<<0,2,4>|<1,-1,1>|<2,2,0>>;

> MatrixEliminate(A, output=lu);
Error, (in pswitchLU) LU decomposition not possible, try LUP.

Now we compare various pivoting strategies.
> A:=<<18,-12,24>|<1,-1,1>|<2,2,0>>;

> MatrixEliminate(A);
  MatrixEliminate(A, pivoting=smallpivot);
  MatrixEliminate(A, pivoting=euclidean);

Partial pivoting is used to help alleviate numerical errors in floating point elimination. We try a random experiment.
> n:=100:
  A:=RandomMatrix(n, generator=-9...9.):
  while Rank(A)<n do A:=RandomMatrix(n,m, generator=-9...9.): end do:
  xsol:=RandomVector(n, generator=-5..5): b:=A.xsol:
  gjsol:=MatrixEliminate(<A|b>, pivoting=minimal, output=gaussjordan):
  errorPminimal:=add(abs(xsol[i]-gjsol[i,n+1]),i=1..n);
  gjsol:=MatrixEliminate(<A|b>, pivoting=partial, output=gaussjordan):
  errorPpartial:=add(abs(xsol[i]-gjsol[i,n+1]),i=1..n);
  errorPminimal/errorPpartial;