Study Guide
SCMUTILS Reference Manual
-
University:
Massachusetts Institute of Technology -
Course:
12.620J | Classical Mechanics: A Computational Approach Academic year:
2020
-
Views:
278
Pages:
47
Author:
Mara Mathis
1
SCMUTILS Reference Manual
This is a description of the Scmutils system, an integrated library of
procedures, embedded in the programming language Scheme, and intended
to support teaching and research in mathematical physics and electrical
engineering. Scmutils and Scheme are particularly effective in work
where the almost-functional nature of Scheme is advantageous, such as
classical mechanics, where many of the procedures are most easily
formulated as quite abstract manipulations of functions.
Many people contributed to the development of the Scmutils library
over many years, so we may miss some of them. The principal
contributors have been:
Gerald Jay Sussman, Harold Abelson, Jack Wisdom, Jacob Katzenelson,
Hardy Mayer, Christopher P. Hanson, Matthew Halfant, Bill Siebert,
Guillermo Juan Rozas, Panayotis Skordos, Kleanthes Koniaris, Kevin
Lin, Dan Zuras
Scheme and Functional Programming
Scheme is a simple computer language in the Lisp family of languages,
with important structural features derived from Algol-60. We will not
attempt to document Scheme here, as it is adequately documented in the
IEEE standard (IEEE-1178) and in numerous articles and books. The
R^4RS is a terse document describing Scheme in detail. It is included
with this document. We assume that a reader of this document is
familiar with Scheme and has read a book such as
Harold Abelson, Gerald Jay Sussman and Julie Sussman
Structure and Interpretation of Computer Programs
MIT Press and McGraw-Hill (1985, 1996)
As a reminder, Scheme is an expression-oriented language. Expressions
have values. The value of an expression is constructed from the
values of the constituent parts of the expression and the way the
expression is constructed. Thus, the value of
(+ (* :pi (square r)) 1)
is constructed from the values of the symbols
+, *, :pi, square, r, 1
and by the parenthetical organization. 2
In any Lisp-based language, an expression is constructed from parts
with parentheses. The first subexpresson always denotes a procedure
and the rest of the subexpressions denote arguments. So in the case
of the expression above, "square" is a symbol whose value is a
procedure that is applied to a thing (probably a number) denoted by
"r". That value of the application of square to r is then combined
with the number denoted by ":pi" using the procedure denoted by "*"
to make an object (again probably a number) that is combined with the
number denoted by "1" using the procedure denoted by "+". Indeed, if
the symbols have the values we expect from their (hopefully mnemonic)
names,
+
*
:pi
square
1
r
= a means of addition
= a means of multiplication
= a number approximately equal to 3.14159265358979
= a means of squaring
= the multiplicative identity in the reals
some number, say for example 4
then this expression would have the approximate value denoted by the
numeral 51.26548245743669.
We can write expressions denoting procedures.
procedure for squaring can be written
For example the
(lambda (x) (* x x)) ,
which may be read
"the procedure of one argument x that multiplies x by x" .
We may bind a symbol to a value by definition
(define square
(lambda (x) (* x x))) ,
or equivalently, by
(define (square x) (* x x)) .
The application of a defined procedure to operands will bind the
symbols that are the formal parameters to the actual arguments that
are the values of the operands:
(+ (square 3) (square 4)) => 25
This concludes the reminders about Scheme.
alternate source for more information.
You must consult an
One caveat: unlike the Scheme standard the Scmutils system is case
sensitive. 3
Generic Arithmetic
In the Scmutils library arithmetic operators are generic over a
variety of mathematical datatypes. For example, addition makes sense
when applied to such diverse data as numbers, vectors, matrices, and
functions. Additionally, many operations can be given a meaning when
applied to different datatypes. For example, multiplication makes
sense when applied to a number and a vector.
The traditional operator symbols, such as "+" and "*" are bound to
procedures that implement the generic operations. The details of
which operations are defined for which datatypes is described below,
organized by the datatype.
In addition to the standard operations, every piece of mathematical
data, x, can give answers to the following questions:
(type x)
Returns a symbol describing the type of x.
(type 3.14)
(type (vector 1 2 3))
For example,
=> *number*
=> *vector*
(type-predicate x)
Returns a predicate that is true on objects that
are the same type as x
(arity p)
Returns a description of the number of arguments that p,
interpreted as a procedure, accepts, compatible with the MIT
Scheme procedure-arity procedure, except that it is extended for
datatypes that are not usually interpreted as procedures. A
structured object, like a vector, may be applied as a vector of
procedures, and its arity is the intersection of the arities of
the components.
An arity is a newly allocated pair whose car field is the minimum
number of arguments, and whose cdr field is the maximum number of
arguments. The minimum is an exact non-negative integer. The
maximum is either an exact non-negative integer, or ‘#f’ meaning
that the procedure has no maximum number of arguments. In our
version of Scheme #f is the same as the empty list, and a pair
with the empty list in the cdr field is a singleton list, so the
arity will print as shown in the second column.
(arity
(arity
(arity
(arity
(arity
(arity
(arity
(lambda
(lambda
car)
(lambda
(lambda
(lambda
(vector
() 3))
(x) x))
=>
=>
=>
x x))
=>
(x . y) x))
=>
(x #!optional y) x)) =>
cos sin))
=>
(0
(1
(1
(0
(1
(1
(1
.
.
.
.
.
.
.
0)
1)
1)
#f)
#f)
2)
1)
=
=
=
=
=
=
=
(0 .
(1 .
(1 .
(0)
(1)
(1 .
(1 .
0)
1)
1)
2)
1) 4
We will now describe each of the generic operations. These operations
are defined for many but not all of the mathematical datatypes. For
particular datatypes we will list and discuss the operations that only
make sense for them.
(inexact? x)
This procedure is a predicate -- a boolean-valued procedure.
See the R^4RS for an explanation of exactness of numbers.
A compound object, such as a vector or a matrix, is
inexact if it has inexact components.
(zero-like x)
In general, this procedure returns the additive identity of the
type of its argument, if it exists. For numbers this is 0.
(one-like x)
In general, this procedure returns the multiplicative identity of
the type of its argument, if it exists. For numbers this is 1.
(zero? x)
Is true if x is an additive identity.
(one? x)
Is true if x is a multiplicative identity.
(negate x)
=
(- (zero-like x) x)
Gives an object that when added to x yields zero.
(invert x)
=
(/ (one-like x) x)
Gives an object that when multiplied by x yields one. 5
Most of the numerical functions have been generalized to many of the
datatypes, but the meaning may depend upon the particular datatype.
Some are defined for numerical data only.
(=
(+
(*
((/
x1
x1
x1
x1
x1
x2
x2
x2
x2
x2
...)
...)
...)
...)
...)
==>
(expt x1 x2)
(gcd n1 n2 ...)
(sqrt x)
Gives a square root of x, or an approximation to it.
(exp x)
(exp10 x)
(exp2 x)
=
=
=
:e^x
10^x
2^x
(log x)
(log10 x)
(log2 x)
=
=
(/ (log x) (log 10))
(/ (log x) (log 2))
(sin x), (cos x), (tan x)
(sec x), (csc x)
(asin x), (acos x), (atan x)
(atan x1 x2) = (atan (/ x1 x2)) but retains quadrant information
(sinh x), (cosh x), (tanh x)
(sech x), (csch x)
(make-rectangular a1 a2)
(make-polar a1 a2)
(real-part z)
(imag-part z)
(magnitude z)
(angle z)
=
=
a1+ia2
a1*:e^(* +i a2)
(conjugate z)
If M is a quantity that can be interpreted as a square matrix,
(determinant M)
(trace M) 6
Scheme Numbers
Operations on the Scheme Number datatype that are part of standard
Scheme are listed here without comment; those that are not part of
standard Scheme are described. In the following is (any
expression that denotes) an integer. is any real number, is
any complex number, and and are any kind of number.
(type )
(inexact? )
(zero-like )
(one-like )
(zero? )
(one? )
= *number*
==>
= 0
= 1
==>
==>
(negate ), (invert ), (sqrt )
(exp ), (exp10 ), (exp2 )
(log ), (log10 ), (log2 )
(sin ), (cos ), (tan ), (sec ), (csc )
(asin ), (acos ), (atan )
(atan )
(sinh ), (cosh ), (tanh ), (sech ), (csch )
(=
(+
(*
((/
...)
...)
...)
...)
...)
==>
(expt )
(gcd ...)
(make-rectangular )
(make-polar )
(real-part )
(imag-part )
(magnitude )
(angle )
(conjugate )
=
=
+i
*:e^(* +i ) 7
Structured Objects
Scmutils supports a variety of structured object types, such as
vectors, up and down tuples, matrices, and power series.
The explicit constructor for a structured object is a procedure whose
name is what we call objects of that type. For example, we make
explicit vectors with the procedure named "vector", and explicit lists
with the procedure named "list". For example
(list
1 2 3 4 5)
(vector 1 2 3 4 5)
(down 10 3 4)
a list of the first five positive integers
a vector of the first five positive integers
a down tuple with three components
There is no natural way to notate a matrix, except by giving its rows
(or columns). To make a matrix with three rows and five columns:
(define M
(matrix-by-rows (list 1 2 3 4 5)
(list 6 7 8 9 10)
(list 11 12 13 14 15)))
A power series may be constructed from an explicit set of coefficients
(series 1 2 3 4 5)
is the power series whose first five coefficients are the first five
positive integers and all of the rest of the coefficients are zero.
Although each datatype has its own specialized procedures, there are a
variety of generic procedures for selecting the components from
structured objects. To get the n-th component from a linear data
structure, v, such as a vector or a list, one may in general use the
generic selector, "ref":
(ref x n)
All structured objects are accessed by zero-based indexing, as is the
custom in Scheme programs and in relativity. For example, to get the
third element (index = 2) of a vector or a list we can use
(ref (vector 1 2 3 4 5) 2) = 3
(ref (list
1 2 3 4 5) 2) = 3
If M is a matrix, then the component in the i-th row and j-th column
can be obtained by (ref M i j). For the matrix given above
(ref M 1 3) = 9
Other structured objects are more magical
(ref cos-series 6)
= -1/720 8
The number of components of a structured object can be found with the
"size" generic operator:
(size (vector 1 2 3 4 5)) = 5
Besides the extensional constructors, most structured-object datatypes
can be intentionally constructed by giving a procedure whose values
are the components of the object. These "generate" procedures are
(list:generate
(vector:generate
(matrix:generate
(series:generate
n
proc)
n
proc)
m n proc)
proc)
For example, one may make a 6 component vector each of whose
components is :pi times the index of that component, as follows:
(vector:generate 6 (lambda (i) (* :pi i)))
Or a 3X5 matrix whose components are the sum of :pi times the row
number and the speed of light times the column number:
(matrix:generate 3 5 (lambda (i j) (+ (* :pi i) (* :c j))))
Also, it is commonly useful to deal with a structured object in an
elementwise fashion. We provide special combinators for many
structured datatypes that allow one to make a new structure, of the
same type and size of the given ones, where the components of the new
structure are the result of applying the given procedure to the
corresponding components of the given structures.
((list:elementwise proc) ... )
((vector:elementwise proc) ... )
((structure:elementwise proc) ... )
((matrix:elementwise proc) ... )
((series:elementwise proc) ... )
Thus, vector addition is equivalent to (vector:elementwise +). 9
Scheme Vectors
We identify the Scheme vector data type with mathematical
n-dimensional vectors. These are interpreted as up tuples when a
distinction between up tuples and down tuples is made.
We inherit from Scheme the constructors VECTOR and MAKE-VECTOR, the
selectors VECTOR-LENGTH and VECTOR-REF, and zero-based indexing. We
also get the iterator MAKE-INITIALIZED-VECTOR, and the type predicate
VECTOR? In the documentation that follows, will stand for a
vector-valued expression.
(vector? )
==>
(type )
==>
(inexact? )
==>
Is true if any component
*vector*
of is inexact, otherwise it is false.
(vector-length )
==> <+integer>
gets the number of components of
(vector-ref )
gets the th (zero-based) component of vector
(make-initialized-vector )
this is also called (v:generate )
and (vector:generate )
generates an -dimensional vector whose th component is the
result of the application of the to the number .
(zero-like )
==>
Gives the zero vector of the dimension of vector .
(zero? )
(negate )
==>
==>
(conjugate )
==>
Elementwise complex-conjugate of
Simple arithmetic on vectors is componentwise
(= ...)
(+ ...)
(- ...)
==>
==>
==> 10
There are a variety of products defined on vectors.
(dot-product
)
==>
(cross-product )
Cross product only makes sense for 3-dimensional vectors.
(* )
(* )
=
=
(scalar*vector )
(vector*scalar )
==>
==>
(/ )
=
(vector*scalar (/ 1 )) ==>
The product of two vectors makes an outer product structure.
(* )
=
(outer-product )
==>
(euclidean-norm ) = (sqrt (dot-product ))
(abs )
= (euclidean-norm )
(v:inner-product ) = (dot-product (conjugate ) )
(complex-norm )
(magnitude )
= (sqrt (v:inner-product ))
= (complex-norm )
(maxnorm )
gives the maximum of the magnitudes of the components of
(v:make-unit )
=
(/ (euclidean-norm ))
(v:unit? )
=
(one? (euclidean-norm ))
(v:make-basis-unit )
Makes the n-dimensional basis unit vector with zero in all
components except for the ith component, which is one.
(v:basis-unit? )
Is true if and only if is a basis unit vector. 11
Up Tuples and Down Tuples
Sometimes it is advantageous to distinguish down tuples and up tuples.
If the elements of up tuples are interpreted to be the components of
vectors in a particular coordinate system, the elements of the down
tuples may be thought of as the components of the dual vectors in that
coordinate system. The union of the up tuple and the down tuple data
types is the data type we call "structures."
Structures may be recursive and they need not be uniform. Thus it is
possible to have an up structure with three components: the first is a
number, the second is an up structure with two numerical components,
and the third is a down structure with two numerical components. Such
a structure has size (or length) 3, but it has five dimensions.
In Scmutils, the Scheme vectors are interpreted as up tuples, and
the down tuples are distinguished. The predicate "structure?" is true
of any down or up tuple, but the two can be distinguished by the
predicates "up?" and "down?".
(up?
(down?
) ==>
) ==>
(structure? ) = (or (down? ) (up? ))
In the following, stands for any structure-valued expression;
and will be used if necessary to make the distinction.
The generic type operation distinguishes the types:
(type )
==> *vector* or *down*
We reserve the right to change this implementation to distinguish
Scheme vectors from up tuples. Thus, we provide (null) conversions
between vectors and up tuples.
(vector->up )
==>
(vector->down ) ==>
(up->vector )
(down->vector )
==>
==>
Constructors are provided for these types, analogous to list and
vector.
(up . args)
(down . args)
==>
==> 12
The dimension of a structure is the number of entries, adding up the
numbers of entries from substructures. The dimension of any structure
can be determined by
(s:dimension
==> <+integer>
Processes that need to traverse a structure need to know the number of
components at the top level. This is the length of the structure,
(s:length )
==> <+integer>
The ith component (zero-based) can be accessed by
(s:ref i)
For example,
(s:ref (up 3 (up 5 6) (down 2 4)) 1)
(up 5 6)
As usual, the generic ref procedure can recursively access
substructure
(ref (up 3 (up 5 6) (down 2 4)) 1 0)
5
Given a structure we can make a new structure of the same type with
substituted for the th component of the given structure using
(s:with-substituted-coord )
We can construct an entirely new structure of length whose
components are the values of a procedure using s:generate:
(s:generate up/down )
The up/down argument may be either up or down. 13
The following generic arithmetic operations are defined for
structures.
(zero? ) ==>
is true if all of the components of the structure are zero.
(zero-like ) ==>
produces a new structure with the same shape as the given structure
but with all components being zero-like the corresponding component in
the given structure.
(negate )
(magnitude )
(abs )
(conjugate )
==>
==>
==>
==>
produce new structures which are the result of applying the generic
procedure elementwise to the given structure.
(= ... ) ==>
is true only when the corresponding components are =.
(+ ... ) ==>
(- ... ) ==>
These are componentwise addition and subtraction.
(* ) ==> or , a structure or a number
magically does what you want: If the structures are compatible for
contraction the product is the contraction (the sum of the products of
the corresponding components.) If the structures are not compatible
for contraction the product is the structure of the shape and length
of whose components are the products of with the
corresponding components of .
Structures are compatible for contraction if they are of the same
length, of opposite type, and if their corresponding elements are
compatible for contraction.
It is not obvious why this is what you want, but try it, you’ll like
it! 14
For example, the following are compatible for contraction:
(print-expression (* (up (up 2 3) (down 5 7 11))
(down (down 13 17) (up 19 23 29))))
652
Two up tuples are not compatible for contraction.
Their product is an outer product:
(print-expression (* (up 2 3) (up 5 7 11)))
(up (up 10 15) (up 14 21) (up 22 33))
(print-expression (* (up 5 7 11) (up 2 3)))
(up (up 10 14 22) (up 15 21 33))
This product is not generally associative or commutative. It is
commutative for structures that contract, and it is associative for
structures that represent linear transformations.
To yield additional flavor, the definition of square for structures is
inconsistent with the definition of product. It is possible to square
an up tuple or a down tuple. The result is the sum of the squares of
the components. This makes it convenient to write such things as
(/ (square p) (* 2 m)), but it is sometimes confusing.
Some structures, such as the ones that represent inertia tensors, must
be inverted. (The "m" above may be an inertia tensor!) Division is
arranged to make this work, when possible. The details are too hairy
to explain in this short document. We probably need to write a book
about this! 15
Matrices
There is an extensive set of operations for manipulating matrices.
Let , be matrix-valued expressions. The following operations
are provided
(matrix? )
(type )
(inexact? )
==>
==> *matrix*
==>
(m:num-rows )
==> ,
the number of rows in matrix M.
(m:num-cols )
==>