In JACAL, a matrix is just a bunch
of equal length bunch
s,
and this is the structure that the matrix operations currently supported
by JACAL (ncmult(), ^^, transpose(), etc.) expect.
[elt_1, elt_2, ...]
To collect any number of Jacal objects into a bunch, simply enclose them
in square brackets. For example, to make the bunch whose elements are
1
, 2
, 4
, type [1, 2, 4]
. One can also nest
bunches, for example, [1, [[1, 3], [2, 5]], [1, 4]]
. Note
however that the bunch whose only element is [1, 2, 3]
is
[1 2 3]
. It is importance to notice that one has commas and the
other doesn't.
e3 : a:bunch(1, 2, 3); e3: [1, 2, 3] e4 : b:[a]; e4: [1 2 3] e5 : c:[b]; e5: [[1, 2, 3]] e6 : [[[1, 2, 3]]]; e6: [[1, 2, 3]]
ident
takes as argument a positive integer n
and returns an nxn identity matrix. This is sometimes more
convenient than obtaining this same matrix using the command
scalarmatrix
.
e6 : ident(4); [1 0 0 0] [ ] [0 1 0 0] e6: [ ] [0 0 1 0] [ ] [0 0 0 1]
scalarmatrix
takes as inputs a positive integer
size and an algebraic expression entry and returns an
n * n
diagonal matrix whose diagonal entries are all
equal to entry, where n = size
.
e1 : scalarmatrix(3, 6); [6 0 0] [ ] e1: [0 6 0] [ ] [0 0 6]
diagmatrix
takes as input a list of objects and
returns the diagonal matrix having those objects as diagonal entries. In
case one wants all of the diagonal entries to be equal, it is more
convenient to use the command scalarmatrix
.
e3 : diagmatrix(12,3,a,s^2); [12 0 0 0 ] [ ] [0 3 0 0 ] e3: [ ] [0 0 a 0 ] [ ] [0 0 0 2] [ s ] e4 : diagmatrix([1,2],2); [[1, 2] 0] e4: [ ] [ 0 2]
sylvester
returns the matrix introduced
by Sylvester (A Method of Determining By Mere Inspection the
Derivatives from Two Equations of Any Degree, Phil.Mag. 16
(1840)
pp. 132-135, Mathematical Papers, vol. I, pp. 54-57) for computing the
resultant of the two polynomials poly_1 and poly_2 with
respect to the variable var. If one wants to compute the resultant
itself, one can simply use the command resultant
with the same
syntax.
e5 : sylvester(a0 + a1*x + a2*x^2 + a3*x^3, b0 + b1*x + b2*x^2, x); [a3 a2 a1 a0 0 ] [ ] [0 a3 a2 a1 a0] [ ] e5: [b2 b1 b0 0 0 ] [ ] [0 b2 b1 b0 0 ] [ ] [0 0 b2 b1 b0]
genmatrix
takes as arguments a function
function of two variables and two positive integers, rows
and cols. It returns a matrix with the indicated numbers of rows
and columns in which the $(i,j)$th entry is obtained by evaluating
function at $(i,j)$. The function may be defined in any of the
ways available in Jacal, i.e previously by an explicit algebraic
definition, by an explicit lambda expression or by an implicit lambda
expression.
e4 : @1^2+@2^2; 2 2 e4: lambda([@1, @2], @1 + @2 ) e5 : genmatrix(e4,3,5); [2 5 10 17 26] [ ] e5: [5 8 13 20 29] [ ] [10 13 18 25 34]
row
returns the ith row of the matrix
matrix, where i = int
. If int is larger than
the number of rows of matrix, then Jacal prints an error message.
The corresponding command for columns of a matrix is col
.
e3 : u:[[1, 2, 3], [1, 5, 3]]; [1 2 3] e3: [ ] [1 5 3] e4 : row(u, 2); e4: [1, 5, 3]
col
is used to extract a column of a matrix. Here,
matrix is a matrix and integer is a positive integer. If
that integer exceeds the number of columns, an error message such
as
ERROR: list-ref: Wrong type in arg1 ()
appears. Here is an example of correct use of the command
col
:
e19 : a:[[1,2,4],[2,5,6]]; [1 2 4] e19: [ ] [2 5 6] e20 : col(a,2); [2] e20: [ ] [5]
minor
returns the submatrix of matrix
obtained by deleting the ith row and the jth column, where
i = row and j = col.
e21 : b:[[1,2,3],[3,1,5],[5,2,7]]; [1 2 3] [ ] e21: [3 1 5] [ ] [5 2 7] e22 : minor(b,3,1); [2 3] e22: [ ] [1 5]
rapply
is used to access elements of bunches. It
can also access elements nested at lower levels in a bunch. In
particular, it can also access matrix elements. In the above syntax,
bunch is the bunch whose parts one wishes to access, and n,
int_1, int_2, ..., int_n are positive integers.
It returns the int_n-th element of the int_{n-1}-th element
of ... of the int_2-th element of the int_1-th element
of bunch. One can have n = 0
. In that case, rapply
simply returns the bunch.
e2 : rapply([[1,2,3],[1,4,6],3],2,3); e2: 6 e6 : rapply([a,b],2); e6: b e7 : rapply([a,b]); e7: [a, b]
Matrix multiplication.
e1 : a:[[1, 2, 3], [5, 2, 7]]; [1 2 3] e1: [ ] [5 2 7] e2 : b:[[3, 2], [6, 4]]; [3 2] e2: [ ] [6 4] e3 : b . a; [13 10 23] e3: [ ] [26 20 46]
The infix operator ^^
is used for raising a square matrix to an
integral power. It can also be used on nonsquare matrices, but the
results are not documented.
e8 : a:[[1, 0], [-1, 1]]; [1 0] e8: [ ] [-1 1] e9 : a^^3; [1 0] e9: [ ] [-3 1]
dotproduct
returns the dot product of two
row vectors of the same length. It will also give the dot product of
two matrices of the same size by computing the sum of the dot products
of the corresponding rows or, what is the same, the trace of one
matrix times the transpose of the other one.
e28 : a:[1,2,3]; b:[3,1,5]; e28: [1, 2, 3] e29 : e29: [3, 1, 5] e30 : dotproduct(a,b); e30: 20
crossproduct
computes the cross product of two
vectors. By definition, the two vectors must each have three
components.
e24: [2 x, y + 2 x y ] e25 : crossproduct([1,2,3],[4,2,5]); e25: [4, 7, -6]
determinant
computes the determinant of a
square matrix. Attempting to take the determinant of a non-square
matrix will produce an error message.
e1 : a:[[1,2],[6,7]]; [1 2] e1: [ ] [6 7] e2 : determinant(a); e2: -5
(matrix)
.
The tensors
supported by JACAL are an extension of the matrix
structure (i.e., a bunch of bunches of bunches ...) with the added
stipulation that all dimensions
of the tensor be the same length
(e.g., 4x4x4). The number of dimensions (indices) in a tensor is its
rank: A scalar is a tensor of rank 0; a vector is a rank 1 tensor; a
matrix has rank 2; and so on.
Further, just as matrix binary operations place restrictions on the matrices involved (e.g., the row/column length requirement for matrix multiplication), the tensor binary operations require that the dimensions of each tensor be of the same length. For example, you could not multiply a 3x3 tensor and a 4x4x4 tensor.
JACAL's tensors do not support the construct of contravariant and covariant indices. Users must keep track of this information themselves, and perform the necessary operations with an appropriate metric so that the "index gymnastics" is performed correctly.
Before using any of JACAL's tensor operations, execute the following command from the JACAL prompt:
require("tensor");
This loads the file `tensor.scm' into JACAL, and makes the tensor operations available for use.
JACAL currently supports four tensor operations: tmult
,
contract
, indexshift
, and indexswap
. Each of these is
described in detail below.
tmult
takes a minimum of two arguments which are the tensors on
which the multiplication operation is to be performed.
With no additional arguments, tmult
will produce the outer
product of the two input tensors. The rank of the resulting tensor is
the sum of the inputs' ranks, and the components of the result are
formed from the pair-wise products of components of the inputs. For
example, for the input tensors x[a,b]
and y[c]
@result{z:tmult(x,y);} z[a,b,c] = x[a,b]*y[c]
With an additional argument, tmult
will produce the inner product
of the two tensors on the specified index. For example, given
x[i,j]
and y[k,l,m]
z:tmult(x,y,3); => length ----- \ z[a,b,c] = > x[a,q] * y[b,c,q] / ----- q = 1
Note that in this case x only has 2 indices. All of JACAL's tensor
operations modify index inputs to be between 1 and the rank of the
tensor. Thus, in this example, the 3 is modified to 2 in the case of x.
As another example, with x[i,j,k]
and y[l,m,n]
z:tmult(x,y,2); => length ----- \ z[a,b,c,d] = > x[a,q,b] * y[c,q,d] / ----- q = 1
With four arguments, tmult
produces an inner product of the two
tensors on the specified indices. For example, for x[i,j]
and
y[k,l,m]
z:tmult(x,y,1,3); => length ----- \ z[a,b,c] = > x[q,a] * y[b,c,q] / ----- q = 1
Note that matrix multiplication is the special case of an inner product (of
two "two dimensional matrices") on the second and first indices,
respectively: tmult(x,y,2,1) == ncmult(x,y)
Finally, tmult handles the case of a scalar times a tensor, in which case each component of the tensor is multiplied by the scalar.
The contraction operation produces a tensor of rank two less than a given tensor. It does this by performing a summation over two of the indices of the given tensor, as clarified in the examples below.
contract
takes at least one argument which is the tensor on which the
contraction operation is to be performed. One or two additional arguments
may be provided to specify the indices to be used in the summation. If no
additional arguments are provided, the summation is performed over the
first and second indices. With one additional argument, the summation is
over the specified index and the one following it (e.g., if 3 is specified,
the third and fourth indices are used). With two additional arguments, the
summation is performed over the indices specified. The actual indices used
will be constrained to be between 1 and the rank of the tensor.
Examples:
1) For a square matrix (tensor of rank 2), contract
returns a scalar that
is the sum of the diagonal elements of the matrix.
2) Given x[i,j,k,l]
, the command
y:contract(x,2,4);
produces:
length ----- \ y[a,b] = > x[a,q,b,q] / ----- q = 1
Special cases: If contract
is given a scalar (rank 0 tensor) as input,
it just returns the scalar. For a vector (tensor of rank 1), contract
returns a scalar that is the sum of the elements of the vector.
indexshift
rearranges the indices of a tensor. It is one of two
generalizations of the matrix transpose operation (cf. indexswap
).
indexshift
takes at least one argument which is the tensor on which the
index shifting is to be performed. One or two additional arguments may be
provided to specify the index and the position to which it is to be
shifted. If no additional arguments are provided, the first index of the
tensor is shifted to the second position (equivalent the matrix transpose
operation). If one additional argument is provided, it specifies the index
to be shifted, and that index will be shifted "to the right" one position
(e.g., if 3 is specified, the third index will be shifted to the forth
position). If two additional arguments are provided, the first specifies
the index and the second specifies the position to which it is to be
shifted. The actual index shifted and its shifted position will be
constrained to be between 1 and the rank of the tensor.
For example, given x[a,b,c,d]
, the command
y:indexshift(x,1,3);
produces a tensor y such that
y[a,b,c,d] == x[b,c,a,d]
. In this example, the element
that was in position [a,b,c,d]
in x will be in position
[b,c,a,d]
in y.
Special cases: If indexshift
is given a scalar (rank 0 tensor) as input,
it just returns the scalar. For a vector (tensor of rank 1), indexshift
transposes the 1-by-n matrix (row vector) to an n-by-1 matrix (column
vector).
indexswap
rearranges the indices of a tensor. It is one of two
generalizations of the matrix transpose operation (cf. indexshift
).
indexswap
takes at least one argument which is the tensor on which index
swapping is to be performed. One or two additional arguments may be
provided to specify the indices to be swapped. If no additional arguments
are provided, the first and second indices of the tensor are swapped
(equivalent the matrix transpose operation). With one additional
argument, the specified index is swapped with the one following it (e.g.,
if 2 is specified, the second and third indices will be swapped). If two
additional arguments are provided, they specify the indices to be swapped.
The actual indices used will be constrained to be between 1 and the rank of
the tensor.
For example, given x[a,b,c,d], the command y:indexswap(x,2,4);
produces a tensor y such that y[a,b,c,d] = x[a,d,c,b]. In this
example, the element that was in position [a,b,c,d] in x will be
in position [a,d,c,b] in y.
Special cases: If indexswap
is given a scalar (rank 0 tensor) as input,
it just returns the scalar. For a vector (tensor of rank 1), indexswap
transposes the 1-by-n matrix (row vector) to an n-by-1 matrix (column
vector).
Go to the first, previous, next, last section, table of contents.