1. Linear Equations and Matrices

>	A := RandomMatrix(3,3);

>	b := convert(RandomMatrix(3,1),Vector);

We will define the linear system A*x = b.  First we pick variables:

>	X := [seq(x[i],i=1..3)];

>	eqs := GenerateEquations(A,X,b);

This is the symbolic view of a linear system.  To work with a linear system, we typically work with the coefficient matrix A and the right hand side vector b.

>	AA,bb := GenerateMatrix(eqs,X);

>

AA; bb;

2. Vector and Matrix Arithmethic

We noticed already the distinction between the type Vector and the type Matrix:

>	whattype(bb); whattype(AA);

>	mb := convert(bb,Matrix); whattype(mb);

But still, we cannot add AA to mb:

>

AA+mb;

Error, (in rtable/Sum) invalid arguments

>	MatrixAdd(AA,mb);

Error, (in LinearAlgebra:-MatrixAdd) matrix dimensions don't match: 3 x 3 vs 3 x 1

There are two errors, first the overloading of the "+" did not work.  The MatrixAdd complains about the wrong dimensions.  

>	B := RandomMatrix(3);

>	C := MatrixAdd(A,B);

>	MatrixMatrixMultiply(A,mb);

>	MatrixVectorMultiply(A,bb);

Notice the dimensions.  If we multiply a 3-by-3 matrix with a 3-by-1 matrix, we obtain a 3-by-1 result, this is a matrix of 3 rows and 1 column.  Notice that a column vector is not automatically converted into a matrix.

3. Basic Matrix Functions

To determine whether a matrix is singular or not, to find there exists a unique solution to the linear system with coefficient matrix A, we compute the determinant:

>	Determinant(A);

To solve the system A*x = b, we can apply LU decomposition:

>	P,L,U := LUDecomposition(A);

>	MatrixMatrixMultiply(L,U) - A;

With this LU decomposition, the system A*x = b becomes (L*U)*x = b, or L*(U*x) = b.  We can first solve the lower triangular system L*y = b, where y is auxiliary, and then solve U*x = y.

>	y := ForwardSubstitute(L,b); # solve L*y = b

>	sol := BackwardSubstitute(U,y); # solve U*x = y

>	MatrixVectorMultiply(A,sol); b; # verify if sol is a solution

>	fA := convert(A,float);

>	P,L,U := LUDecomposition(fA);

With floating-point numbers on input, the result looks quite different.  Now, for numerical stability, Maple has permuted the first and third row.

>	y := ForwardSubstitute(L,MatrixVectorMultiply(P,b));

>	nsol := BackwardSubstitute(U,y);

>	MatrixVectorMultiply(A,nsol); b;

Up to permutations, we get the solution...

4. Vector Operations

We can apply the Gram-Schmidt orthogonalization procedure:

>	GramSchmidt(A);

Error, (in LinearAlgebra:-GramSchmidt) invalid input: LinearAlgebra:-GramSchmidt expects its 1st argument, V, to be of type {list(Vector), set(Vector)} but received Matrix(3, 3, [[...],[...],[...]], datatype = anything)

>	v := seq(Column(A,i),i=1..ColumnDimension(A));

>	basis := GramSchmidt({v});

>	Matrix(3,3,(i,j) -> DotProduct(basis[i],basis[j]));

We just verified that the three vectors returned by GramSchmidt are perpendicular to each other.

5. Smith Normal Form and the SVD