\documentstyle[11pt,reduce]{article}
\date{}
\title{A Vector Algebra and Calculus Package for REDUCE}
\author{David Harper \\
Astronomy Unit \\
Queen Mary and Westfield College \\
University of London \\
Mile End Road \\
London E1 4NS \\
England \\[0.05in]
Electronic mail: {\it adh@star.qmw.ac.uk}}
\begin{document}
\maketitle
\index{AVECTOR package}
\section{Introduction}
This package \footnote{Reference: Computer Physics Communications,
{\bf 54}, 295-305 (1989)}
is written in RLISP (the LISP meta-language) and is
intended for use with REDUCE 3.4. \index{vector algebra} It provides
REDUCE with the ability to perform vector algebra using the same
notation as scalar algebra. The basic algebraic operations are
supported, as are differentiation and integration of vectors with
respect to scalar variables, cross product and dot product, component
manipulation and application of scalar functions ({\em e.g.} cosine)
to a vector to yield a vector result.
A set of vector calculus operators are provided for use with any
orthogonal curvilinear coordinate system. These operators are
gradient, divergence, curl and del-squared (Laplacian). The Laplacian
operator can take scalar or vector arguments.
Several important coordinate systems are pre-defined and can be
invoked by name. It is also possible to create new coordinate systems
by specifying the names of the coordinates and the values of the scale
factors.
\section{Vector declaration and initialisation}
Any name may be declared to be a vector, provided that it has
not previously been declared as a matrix or an array. To
declare a list of names to be vectors use the VEC command:
\index{VEC command}
\begin{verbatim}
VEC A,B,C;
\end{verbatim}
declares the variables {\tt A}, {\tt B} and {\tt C} to be vectors.
If they have already been assigned (scalar) values, these will be lost.
When a vector is declared using the {\tt VEC} command, it does not
have an initial value.
If a vector value is assigned to a scalar variable, then that
variable will automatically be declared as a vector and the
user will be notified that this has happened.
\index{AVEC function}
A vector may be initialised using the {\tt AVEC} function which
takes three scalar arguments and returns a vector made up
from those scalars. For example
\begin{verbatim}
A := AVEC(A1, A2, A3);
\end{verbatim}
sets the components of the vector {\tt A} to {\tt A1}, {\tt A2} and {\tt A3}.
\section{Vector algebra}
(In the examples which follow, {\tt V}, {\tt V1}, {\tt V2} {\em etc}
are assumed to be vectors while {\tt S}, {\tt S1}, {\tt S2} etc are scalars.)
\index{+ ! vector} \index{- ! vector} \index{* ! vector} \index{/ ! vector}
The scalar algebra operators +,-,* and / may be used with
vector operands according to the rules of vector algebra.
Thus multiplication and division of a vector by a scalar
are both allowed, but it is an error to multiply or
divide one vector by another.
\begin{tabular}{l l}
{\tt V := V1 + V2 - V3;} & Addition and subtraction \\
{\tt V := S1*3*V1;} & Scalar multiplication \\
{\tt V := V1/S;} & Scalar division \\
{\tt V := -V1;} & Negation \\
\end{tabular}
\index{DOT ! vector} \index{Dot product} \index{CROSS ! vector}
\index{cross product}
\noindent Vector multiplication is carried out using the infix
operators {\tt DOT} and {\tt CROSS}. These are defined to have
higher precedence than scalar multiplication and
division.
\begin{tabular}{l l}
{\tt V := V1 CROSS V2;} & Cross product \\
{\tt S := V1 DOT V2;} & Dot product \\
{\tt V := V1 CROSS V2 + V3;} & \\
{\tt V := (V1 CROSS V2) + V3;} & \\
\end{tabular}
The last two expressions are equivalent due to the precedence of
the {\tt CROSS} operator.
\index{VMOD operator}
The modulus of a vector may be calculated using the {\tt VMOD} operator.
\begin{verbatim}
S := VMOD V;
\end{verbatim}
A unit vector may be generated from any vector using the {\tt VMOD}
operator.
\begin{verbatim}
V1 := V/(VMOD V);
\end{verbatim}
Components may be extracted from any vector using index notation
in the same way as an array.
\begin{tabular}{l l}
{\tt V := AVEC(AX, AY, AZ);} & \\
{\tt V(0);} & yields AX \\
{\tt V(1);} & yields AY \\
{\tt V(2);} & yields AZ \\
\end{tabular}
It is also possible to set values of individual components. Following
from above:
\begin{verbatim}
V(1) := B;
\end{verbatim}
The vector {\tt V} now has components {\tt AX}, {\tt B}, {\tt AZ}.
\index{vector ! differentiation} \index{vector ! integration}
\index{differentiation ! vector} \index{differentiation ! vector}
Vectors may be used as arguments in the differentiation and
integration routines in place of the dependent expression.
\begin{tabular}{l l}
{\tt V := AVEC(X**2, SIN(X), Y);} & \\
{\tt DF(V,X);} & yields (2*X, COS(X), 0) \\
{\tt INT(V,X);} & yields (X**3/3, -COS(X), Y*X) \\
\end{tabular}
Vectors may be given as arguments to monomial functions such as {\tt
SIN}, {\tt LOG} and {\tt TAN}. The result is a vector obtained by
applying the function component-wise to the argument vector.
\begin{tabular}{l l}
{\tt V := AVEC(A1, A2, A3);} & \\
{\tt SIN(V);} & yields (SIN(A1), SIN(A2), SIN(A3)) \\
\end{tabular}
\section{ Vector calculus}
\index{DIV ! operator} \index{divergence ! vector field}
\index{GRAD ! operator} \index{gradient ! vector field}
\index{CURL ! operator} \index{curl ! vector field}
\index{DELSQ ! operator} \index{Laplacian ! vector field}
The vector calculus operators div, grad and curl are recognised.
The Laplacian operator is also available and may be applied to
scalar and vector arguments.
\begin{tabular}{l l}
{\tt V := GRAD S;} & Gradient of a scalar field \\
{\tt S := DIV V;} & Divergence of a vector field \\
{\tt V := CURL V1;} & Curl of a vector field \\
{\tt S := DELSQ S1;} & Laplacian of a scalar field \\
{\tt V := DELSQ V1;} & Laplacian of a vector field \\
\end{tabular}
These operators may be used in any orthogonal curvilinear coordinate
system. The user may alter the names of the coordinates and the values
of the scale factors. Initially the coordinates are {\tt X}, {\tt Y}
and {\tt Z} and the scale factors are all unity.
\index{COORDS vector} \index{HFACTORS scale factors}
There are two special vectors : {\tt COORDS} contains the names
of the coordinates in the current system and {\tt HFACTORS}
contains the values of the scale factors.
\index{COORDINATES operator}
The coordinate names may be changed using the {\tt COORDINATES}
operator.
\begin{verbatim}
COORDINATES R,THETA,PHI;
\end{verbatim}
This command changes the coordinate names to {\tt R}, {\tt THETA} and
{\tt PHI}.
\index{SCALEFACTORS operator}
The scale factors may be altered using the {\tt SCALEFACTORS} operator.
\begin{verbatim}
SCALEFACTORS(1,R,R*SIN(THETA));
\end{verbatim}
This command changes the scale factors to {\tt 1}, {\tt R} and {\tt R
SIN(THETA)}.
Note that the arguments of {\tt SCALEFACTORS} must be enclosed in
parentheses. This is not necessary with {\tt COORDINATES}.
When vector differential operators are applied to an expression,
the current set of coordinates are used as the independent
variables and the scale factors are employed in the calculation.
(See, for example, Batchelor G.K. 'An Introduction to Fluid
Mechanics', Appendix 2.)
\index{"!*CSYSTEMS global (AVECTOR)}
Several coordinate systems are pre-defined and may be invoked by
name. To see a list of valid names enter
\begin{verbatim}
SYMBOLIC !*CSYSTEMS;
\end{verbatim}
and REDUCE will respond with something like
\begin{verbatim}
(CARTESIAN SPHERICAL CYLINDRICAL)
\end{verbatim}
\index{GETCSYSTEM command}
To choose a coordinate system by name, use the command {\tt GETCSYSTEM}.
To choose the Cartesian coordinate system :
\begin{verbatim}
GETCSYSTEM 'CARTESIAN;
\end{verbatim}
\index{PUTCSYSTEM command}
Note the quote which prefixes the name of the coordinate system. This
is required because {\tt GETCSYSTEM} (and its complement {\tt
PUTCSYSTEM}) is a {\tt SYMBOLIC} procedure which requires a literal
argument.
REDUCE responds by typing a list of the coordinate names in that
coordinate system. The example above would produce the response
\begin{verbatim}
(X Y Z)
\end{verbatim}
whilst
\begin{verbatim}
GETCSYSTEM 'SPHERICAL;
\end{verbatim}
would produce
\begin{verbatim}
(R THETA PHI)
\end{verbatim}
Note that any attempt to invoke a coordinate system is subject to the
same restrictions as the implied calls to {\tt COORDINATES} and {\tt
SCALEFACTORS}. In particular, {\tt GETCSYSTEM} fails if any of the
coordinate names has been assigned a value and the previous coordinate
system remains in effect.
A user-defined coordinate system can be assigned a name using the
command {\tt PUTCSYSTEM}. It may then be re-invoked at a later stage using
{\tt GETCSYSTEM}.
\example\index{AVECTOR package ! example}
We define a general coordinate system with coordinate names {\tt
X},{\tt Y},{\tt Z} and scale factors {\tt H1},{\tt H2},{\tt H3} :
\begin{verbatim}
COORDINATES X,Y,Z;
SCALEFACTORS(H1,H2,H3);
PUTCSYSTEM 'GENERAL;
\end{verbatim}
This system may later be invoked by entering
\begin{verbatim}
GETCSYSTEM 'GENERAL;
\end{verbatim}
\section{Volume and Line Integration}
Several functions are provided to perform volume and line integrals.
These operate in any orthogonal curvilinear coordinate system and
make use of the scale factors described in the previous section.
Definite integrals of scalar and vector expressions may be calculated
using the {\tt DEFINT} function.
\example\index{AVECTOR package ! example}
\index{DEFINT function} \index{integration ! definite (simple)}
\index{definite integration (simple)}
\noindent To calculate the definite integral of $\sin(x)^2$ between 0 and
2$\pi$ we enter
\begin{verbatim}
DEFINT(SIN(X)**2,X,0,2*PI);
\end{verbatim}
This function is a simple extension of the {\tt INT} function taking
two extra arguments, the lower and upper bounds of integration
respectively.
\index{VOLINTEGRAL function} \index{integration ! volume}
Definite volume integrals may be calculated using the {\tt
VOLINTEGRAL} function whose syntax is as follows :
\noindent {\tt VOLINTEGRAL}({\tt integrand}, vector {\tt lower-bound},
vector {\tt upper-bound});
\example\index{AVECTOR package ! example}
\noindent In spherical polar coordinates we may calculate the volume of a
sphere by integrating unity over the range $r$=0 to {\tt RR}, $\theta$=0 to
{\tt PI}, $\phi$=0 to 2*$\pi$ as follows :
\begin{tabular}{l l}
{\tt VLB := AVEC(0,0,0);} & Lower bound \\
{\tt VUB := AVEC(RR,PI,2*PI);} & Upper bound in $r, \theta, \phi$
respectively \\
{\tt VOLINTORDER := (0,1,2);} & The order of integration \\
{\tt VOLINTEGRAL(1,VLB,VUB);} & \\
\end{tabular}
\index{VOLINTORDER vector}
Note the use of the special vector {\tt VOLINTORDER} which controls
the order in which the integrations are carried out. This vector
should be set to contain the number 0, 1 and 2 in the required order.
The first component of {\tt VOLINTORDER} contains the index of the
first integration variable, the second component is the index of the
second integration variable and the third component is the index of
the third integration variable.
\example\index{AVECTOR package ! example}
Suppose we wish to calculate the volume of a right circular cone. This
is equivalent to integrating unity over a conical region with the
bounds:
\begin{tabular}{l l}
z = 0 to H & (H = the height of the cone) \\
r = 0 to pZ & (p = ratio of base diameter to height) \\
phi = 0 to 2*PI & \\
\end{tabular}
We evaluate the volume by integrating a series of infinitesimally thin
circular disks of constant z-value. The integration is thus performed
in the order : d($\phi$) from 0 to $2\pi$, dr from 0 to p*Z, dz from 0 to H.
The order of the indices is thus 2, 0, 1.
\begin{verbatim}
VOLINTORDER := AVEC(2,0,1);
VLB := AVEC(0,0,0);
VUB := AVEC(P*Z,H,2*PI);
VOLINTEGRAL(1,VLB,VUB);
\end{verbatim}
(At this stage, we replace {\tt P*H} by {\tt RR}, the base radius of the cone,
to obtain the result in its more familiar form.)
\index{LINEINT function} \index{DEFLINEINT function}
\index{integration ! line} \index{line integrals}
Line integrals may be calculated using the {\tt LINEINT} and {\tt DEFLINEINT}
functions. Their general syntax is
\noindent {\tt LINEINT}({\tt vector-function}, {\tt vector-curve},
{\tt variable});
\noindent{\tt DEFLINENINT}({\tt vector-function}, {\tt vector-curve},
{\tt variable}, {\tt lower-bound}, {\tt upper-bound});
\noindent where
\begin{description}
\item[{\tt vector-function}] is any vector-valued expression;
\item[{\tt vector-curve}] is a vector expression which describes the path of
integration in terms of the independent variable;
\item[{\tt variable}] is the independent variable;
\item[{\tt lower-bound}]
\item[{\tt upper-bound}] are the bounds of integration in terms of the
independent variable.
\end{description}
\example\index{AVECTOR package ! example}
In spherical polar coordinates, we may integrate round a line of
constant theta (`latitude') to find the length of such a line. The
vector function is thus the tangent to the `line of latitude', (0,0,1)
and the path is {\tt (0,LAT,PHI)} where {\tt PHI} is the independent
variable. We show how to obtain the definite integral {\em i.e.} from
$\phi=0$ to $2 \pi$ :
\begin{verbatim}
DEFLINEINT(AVEC(0,0,1),AVEC(0,LAT,PHI),PHI,0,2*PI);
\end{verbatim}
\section{Defining new functions and procedures}
Most of the procedures in this package are defined in symbolic mode
and are invoked by the REDUCE expression-evaluator when a vector
expression is encountered. It is not generally possible to define
procedures which accept or return vector values in algebraic mode.
This is a consequence of the way in which the REDUCE interpreter
operates and it affects other non-scalar data types as well : arrays
cannot be passed as algebraic procedure arguments, for example.
\section{Acknowledgements}
This package was written whilst the author was the U.K. Computer
Algebra Support Officer at the University of Liverpool Computer Laboratory.
\end{document}