Next: Numerical Computation in Up: Programming in Maple Previous: Programming in Maple

Matrix and Vector Computation in Maple

A vector in Maple is represented by a one-dimensional array indexed from 1, and a matrix is represented by a two-dimensional array, with row and column indices starting from 1. Here is one way to create a 5 by 5 Hilbert matrix. Recall that a Hilbert matrix is a symmetric matrix whose (i,j)'th entry is .


> H := array(1..5,1..5):
> for i to 5 do for j to 5 do H[i,j] := 1/(i+j-1) od od;
> H;

                                       H

> eval(H);
                          [  1   1/2  1/3  1/4  1/5 ]
                          [                         ]
                          [ 1/2  1/3  1/4  1/5  1/6 ]
                          [                         ]
                          [ 1/3  1/4  1/5  1/6  1/7 ]
                          [                         ]
                          [ 1/4  1/5  1/6  1/7  1/8 ]
                          [                         ]
                          [ 1/5  1/6  1/7  1/8  1/9 ]

The linalg package contains many functions for computing with vectors and matrices in Maple. This matrix could also have been created in the following way using the matrix command in the linear algebra package.


linalg[matrix](5,5,(i,j) -> 1/(i+j-1));

Maple can compute the determinant and inverse of the matrix H and many other matrix operations. See ?linalg for a list of the operations.

Here we show a program for doing a so called row reduction on a matrix using Gaussian elimination. GaussianElimination(A,'r') computes the reduced matrix, an upper triangular matrix. It also has an optional 2nd argument which, if given, is assigned the rank of the matrix.


GaussianElimination := proc(A,rank)
local m,n,i,j,B,r,c;

    if not type(A,'matrix(rational)') then
        ERROR(`expecting a matrix of rational numbers`,A) fi;

    m := linalg[rowdim](A); # the number of rows of the matrix
    n := linalg[coldim](A); # the number of columns of the matrix

    B := array(1..m,1..n);
    for i to m do for j to n do B[i,j] := A[i,j] od od;

    r := 1; # r and c are row and column indices
    for c to m while r <= n do

        for i from r to n while B[i,c] = 0 do od;  # search for a pivot
        if i <= n then
            if i <> r then # interchange row i with row r
                for j from c to m do
                    t := B[i,j]; B[i,j] := B[r,j]; B[r,j] := t
                od
            fi;
            for i from r+1 to n do
                if B[i,c] <> 0 then
                    t := B[i,c]/B[r,c];
                    for j from c+1 to m do B[i,j] := B[i,j]-t*B[r,j] od;
                    B[i,c] := 0
                fi
            od;
            r := r + 1                # go to next row
        fi
    od;                       # go to next column

    if nargs>1 then rank := r-1 fi;
    eval(B)

end:

The type matrix(rational) specifies a matrix (a 2-dimensional array in Maple) whose entries are all rational numbers.


bondaren@thsun1.jinr.ru