> with(linalg);
[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp,
Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub,
band, basis, bezout, blockmatrix, charmat, charpoly, cholesky,
col, coldim, colspace, colspan, companion, concat, cond, copyinto,
crossprod, curl, definite, delcols, delrows, det, diag, diverge,
dotprod, eigenvals, eigenvalues, eigenvectors, eigenvects,
entermatrix, equal, exponential, extend, ffgausselim, fibonacci,
forwardsub, frobenius, gausselim, gaussjord, geneqns, genmatrix,
grad, hadamard, hermite, hessian, hilbert, htranspose, ihermite,
indexfunc, innerprod, intbasis, inverse, ismith, issimilar,
iszero, jacobian, jordan, kernel, laplacian, leastsqrs, linsolve,
matadd, matrix, minor, minpoly, mulcol, mulrow, multiply, norm,
normalize, nullspace, orthog, permanent, pivot, potential,
randmatrix, randvector, rank, ratform, row, rowdim, rowspace,
rowspan, rref, scalarmul, singularvals, smith, stack, submatrix,
subvector, sumbasis, swapcol, swaprow, sylvester, toeplitz, trace,
transpose, vandermonde, vecpotent, vectdim, vector, wronskian]
# Matrix Operations
#
# Here is how to input a matrices A and B.
#
> A:=matrix(3,3,[1,2,1,0,-1,2,4,7,-6]);
[1 2 1]
[ ]
A := [0 -1 2]
[ ]
[4 7 -6]
> B:=matrix(3,3,[4,2,1,-2,0,2,3,7,2]);
[ 4 2 1]
[ ]
B := [-2 0 2]
[ ]
[ 3 7 2]
#
# Scalar multiply A by -2
> evalm(2*A);
[2 4 2]
[ ]
[0 -2 4]
[ ]
[8 14 -12]
# Note that without evalm, all you would get is 2A
> 2*A;
2 A
# Scalar division works too.
> evalm(A/2);
[1/2 1 1/2]
[ ]
[ 0 -1/2 1 ]
[ ]
[ 2 7/2 -3 ]
# Add two matrices: A+B
> evalm(A+B);
[ 5 4 2]
[ ]
[-2 -1 4]
[ ]
[ 7 14 -4]
# Multiply two matrices: A*B
> evalm(A&*B);
[ 3 9 7]
[ ]
[ 8 14 2]
[ ]
[-16 -34 6]
# The matrix A and then its transpose:
> evalm(A); transpose(A);
[1 2 1]
[ ]
[0 -1 2]
[ ]
[4 7 -6]
[1 0 4]
[ ]
[2 -1 7]
[ ]
[1 2 -6]
# The inverse of A
> inverse(A);
[ 19 ]
[-2/3 -- 5/12]
[ 12 ]
[ ]
[2/3 -5/6 -1/6]
[ ]
[ -1 ]
[1/3 1/12 -- ]
[ 12 ]
# The negative power also works to find the inverse
> evalm(A^(-1));
[ 19 ]
[-2/3 -- 5/12]
[ 12 ]
[ ]
[2/3 -5/6 -1/6]
[ ]
[ -1 ]
[1/3 1/12 -- ]
[ 12 ]
# Solving Matrix Equations
#
# Let's solve the equation Ax=b where b is the column vector (1,4,-4)
# (transposed).
# First, define A and b
> A:=matrix(3,3,[1,2,1,0,-1,2,4,7,-6]);
[1 2 1]
[ ]
A := [0 -1 2]
[ ]
[4 7 -6]
> b:=matrix(3,1,[1,4,-4]);
[ 1]
[ ]
b := [ 4]
[ ]
[-4]
# Then augment A with b (and call the result AA)
> AA:=augment(A,b);
[1 2 1 1]
[ ]
AA := [0 -1 2 4]
[ ]
[4 7 -6 -4]
# Then apply Gaussian elimination
> G:=gaussjord(AA);
[1 0 0 4]
[ ]
G := [0 1 0 -2]
[ ]
[0 0 1 1]
# The answer is x=(4,-2,1) (transposed) - the last column of G. Let's
# check this.
> x:=matrix(3,1,[4,-2,1]);
[ 4]
[ ]
x := [-2]
[ ]
[ 1]
> evalm(A&*x);
[ 1]
[ ]
[ 4]
[ ]
[-4]
# which equals b.
#
# Maple can extract the fourth column of G with the col command
# so that you won't have to write out the solution by hand.
> x:=col(G,4);
x := [4, -2, 1]
#
# Another way to solve Ax=b is to solve for x as
# x=inverse(A)*b. Note that the order is important
# (b*inverse(A) would not make sense)
> x:=evalm(inverse(A)&*b);
[ 4]
[ ]
x := [-2]
[ ]
[ 1]
# Here is another example
> A:=matrix(3,3,[1,2,-1,0,1,-3,2,3,1]);
[1 2 -1]
[ ]
A := [0 1 -3]
[ ]
[2 3 1]
> b:=matrix(3,1,[-10,-13,-7]);
[-10]
[ ]
b := [-13]
[ ]
[ -7]
> AA:=augment(A,b);
[1 2 -1 -10]
[ ]
AA := [0 1 -3 -13]
[ ]
[2 3 1 -7]
> gaussjord(AA);
[1 0 5 16]
[ ]
[0 1 -3 -13]
[ ]
[0 0 0 0]
# There are an infinite number of solutions.
# The third component can be arbitrary, x_3=c.
# then x_2=3c-13 and x_1=16-5c.
#
# On the other hand, take the following b
> b:=matrix(3,1,[2,4,-5]);
[ 2]
[ ]
b := [ 4]
[ ]
[-5]
> AA:=augment(A,b);
[1 2 -1 2]
[ ]
AA := [0 1 -3 4]
[ ]
[2 3 1 -5]
> gaussjord(AA);
[1 0 5 0]
[ ]
[0 1 -3 0]
[ ]
[0 0 0 1]
# Here, there is no solution.
>
> C:=matrix(3,3,[1,2,-1,0,0,0,2,3,1]);
[1 2 -1]
[ ]
C := [0 0 0]
[ ]
[2 3 1]
> E:=matrix(3,3,[1,0,0,2,1,-1,0,0,1]);
[1 0 0]
[ ]
E := [2 1 -1]
[ ]
[0 0 1]
> evalm(E&*C);
[1 2 -1]
[ ]
[0 1 -3]
[ ]
[2 3 1]
>
# Determinants and Other Operations
#
# Determinants. Here is a matrix
> A:=matrix(3,3,[1,2,1,0,-1,2,4,7,-6]);
[1 2 1]
[ ]
A := [0 -1 2]
[ ]
[4 7 -6]
# Here is its determinant
> det(A);
12
# Here is a singular matrix
> A:=matrix(3,3,[1,2,-1,0,1,-3,2,3,1]);
[1 2 -1]
[ ]
A := [0 1 -3]
[ ]
[2 3 1]
> det(A);
0
# Note that its inverse does not exist
> inverse(A);
Error, (in inverse) singular matrix
# The adjoint of A
> A:=matrix(3,3,[1,2,1,0,-1,2,4,7,-6]);
[1 2 1]
[ ]
A := [0 -1 2]
[ ]
[4 7 -6]
> adjoint(A);
[-8 19 5]
[ ]
[ 8 -10 -2]
[ ]
[ 4 1 -1]
# Note that adjoint(A)/det(A) is the same as the inverse of A
> evalm(adjoint(A)/det(A));
[ 19 ]
[-2/3 -- 5/12]
[ 12 ]
[ ]
[2/3 -5/6 -1/6]
[ ]
[ -1 ]
[1/3 1/12 -- ]
[ 12 ]
> inverse(A);
[ 19 ]
[-2/3 -- 5/12]
[ 12 ]
[ ]
[2/3 -5/6 -1/6]
[ ]
[ -1 ]
[1/3 1/12 -- ]
[ 12 ]
>
# Kernel, Rowspace and Column Space
#
> A:=matrix(3,5,[1,-1,0,2,3,2,0,2,2,8,1,1,2,0,5]);
[1 -1 0 2 3]
[ ]
A := [2 0 2 2 8]
[ ]
[1 1 2 0 5]
> gaussjord(A);
[1 0 1 1 4]
[ ]
[0 1 1 -1 1]
[ ]
[0 0 0 0 0]
> NA:=kernel(A);
NA := {[-1, -1, 1, 0, 0], [-4, -1, 0, 0, 1], [-1, 1, 0, 1, 0]}
> v:=matrix(5,1,op(2,NA));
[-4]
[ ]
[-1]
[ ]
v := [ 0]
[ ]
[ 0]
[ ]
[ 1]
> evalm(A&*v);
[0]
[ ]
[0]
[ ]
[0]
> gaussjord(A);
[1 0 1 1 4]
[ ]
[0 1 1 -1 1]
[ ]
[0 0 0 0 0]
>
> rowspace(A);
{[1, 0, 1, 1, 4], [0, 1, 1, -1, 1]}
> colspace(A);
{[0, 1, 1], [1, 0, -1]}
>
>