\documentstyle[11pt,reduce]{article}
\date{}
\title{NUMERIC}
\author{Herbert Melenk \\
Konrad--Zuse--Zentrum f\"ur Informationstechnik Berlin \\
E--mail: Melenk@sc.zib--berlin.de}
\begin{document}
\maketitle
\index{NUMERIC package}
The {\small NUMERIC} package implements some numerical (approximative)
algorithms for {\small REDUCE}, based on the {\small REDUCE} rounded mode
arithmetic. These algorithms are implemented for standard cases.
They should not be called for ill-conditioned problems;
please use standard mathematical libraries for these.
\section{Syntax}
\subsection{Intervals, Starting Points}
Intervals are generally coded as lower bound and
upper bound connected by the operator \verb+`..'+, usually
associated to a variable in an
equation. E.g.
\begin{verbatim}
x= 2.5 .. 3.5
\end{verbatim}
means that the variable x is taken in the range from 2.5 up to
3.5. Note, that the bounds can be algebraic
expressions, which, however, must evaluate to numeric results.
In cases where an interval is returned as the result, the lower
and upper bounds can be extracted by the \verb+PART+ operator
as the first and second part respectively.
A starting point is specified by an equation with a numeric
righthand side, e.g.
\begin{verbatim}
x=3.0
\end{verbatim}
If for multivariate applications several coordinates must be
specified by intervals or as a starting point, these
specifications can be collected in one parameter (which is then
a list) or they can be given as separate parameters
alternatively. The listified form is more appropriate when the
parameters are built from other REDUCE calculations in an
automatical style, while the flat form is more convenient
for direct interactive input.
\subsection{Accuracy Control}
The keyword parameters $accuracy=a$ and $iterations=i$, where
$a$ and $i$ must be positive integer numbers, control the
iterative algorithms: the iteration is continued until
the local error is below $10^{-a}$; if that is impossible
within $i$ steps, the iteration is terminated with an
error message. The values reached so far are then returned
as the result.
\subsection{tracing}
Normally the algorithms produce only a minimum of printed
output during their operation. In cases of an unsuccessful
or unexpected long operation a trace of the iteration can be
printed by setting
\begin{verbatim}
on trnumeric;
\end{verbatim}
\section{Minima}
The Fletcher Reeves version of the $steepest\ descent$
algorithms is used to find the minimum of a
function of one or more variables. The
function must have continuous partial derivatives with respect to all
variables. The starting point of the search can be
specified; if not, random values are taken instead.
The steepest descent algorithms in general find only local
minima.
Syntax:
\begin{description}
\item[NUM\_MIN] $(exp, var_1[=val_1] [,var_2[=val_2] \ldots]$
$ [,accuracy=a][,iterations=i]) $
or
\item[NUM\_MIN] $(exp, \{ var_1[=val_1] [,var_2[=val_2] \ldots] \}$
$ [,accuracy=a][,iterations=i]) $
where $exp$ is a function expression,
$var_1, var_2, \ldots$ are the variables in $exp$ and
$val_1,val_2, \ldots$ are the (optional) start values.
MIN tries to find the next local minimum along the descending
path starting at the given point. The result is a list
with the minimum function value as first element followed by a list
of equations, where the variables are equated to the coordinates
of the result point.
\end{description}
Examples:
\begin{verbatim}
num_min(sin(x)+x/5, x);
{4.9489585606,{X=29.643767785}}
num_min(sin(x)+x/5, x=0);
{ - 1.3342267466,{X= - 1.7721582671}}
% Rosenbrock function (well known as hard to minimize).
fktn := 100*(x1**2-x2)**2 + (1-x1)**2;
num_min(fktn, x1=-1.2, x2=1, iterations=200);
{0.00000021870228295,{X1=0.99953284494,X2=0.99906807238}}
\end{verbatim}
\section{Roots of Functions/ Solutions of Equations}
An adaptively damped Newton iteration is used to find
an approximative zero of a function, a function vector or the solution
of an equation or an equation system. Equations are
internally converted to a difference of lhs and rhs such
that the Newton method (=zero detection) can be applied. The expressions
must have continuous derivatives for all variables.
A starting point for the iteration can be given. If not given,
random values are taken instead. If the number of
forms is not equal to the number of variables, the
Newton method cannot be applied. Then the minimum
of the sum of absolute squares is located instead.
With ON COMPLEX solutions with imaginary parts can be
found, if either the expression(s) or the starting point
contain a nonzero imaginary part.
Syntax:
\begin{description}
\item[NUM\_SOLVE] $(exp_1, var_1[=val_1][,accuracy=a][,iterations=i])$
or
\item[NUM\_SOLVE] $(\{exp_1,\ldots,exp_n\},
var_1[=val_1],\ldots,var_1[=val_n]$
\item[\ \ \ \ \ \ \ \ ]$[,accuracy=a][,iterations=i])$
or
\item[NUM\_SOLVE] $(\{exp_1,\ldots,exp_n\},
\{var_1[=val_1],\ldots,var_1[=val_n]\}$
\item[\ \ \ \ \ \ \ \ ]$[,accuracy=a][,iterations=i])$
where $exp_1, \ldots,exp_n$ are function expressions,
$var_1, \ldots, var_n$ are the variables,
$val_1, \ldots, val_n$ are optional start values.
SOLVE tries to find a zero/solution of the expression(s).
Result is a list of equations, where the variables are
equated to the coordinates of the result point.
The Jacobian matrix is stored as side effect the shared
variable JACOBIAN.
\end{description}
Example:
\begin{verbatim}
num_solve({sin x=cos y, x + y = 1},{x=1,y=2});
{X= - 1.8561957251,Y=2.856195584}
jacobian;
[COS(X) SIN(Y)]
[ ]
[ 1 1 ]
\end{verbatim}
\section{Integrals}
An adaptive multilevel integration algorithnm is used
to compute univariate integrals over an interval or
multivariate integrals over a rectangular domain.
The algorithm tolerates isolated singularities.
The value $iterations$ here limits the number of
local interval intersection levels.
$Accuracy$ is a measure for the relative total discretization
error (comparison of order 1 and order 2 approximations).
Syntax:
\begin{description}
\item[NUM\_INT] $(exp,var_1=l_1 .. u_1 [,var_2=l_2 .. u_2 \ldots]$
\item[\ \ \ \ \ \ ]$[,accuracy=a][,iterations=i])$
where $exp$ is the function to be integrated,
$var_1, var_2 , \ldots$ are the integration variables,
$l_1, l_2 , \ldots$ are the lower bounds,
$u_1, u_2 , \ldots$ are the upper bounds.
Result is the value of the integral.
\end{description}
Example:
\begin{verbatim}
num_int(sin x,x=0 .. pi);
2.0000010334
\end{verbatim}
\section{Ordinary Differential Equations}
A Runge-Kutta method of order 3 finds an approximate graph for
the solution of a ordinary differential equation
real initial value problem.
Syntax:
\begin{description}
\item[NUM\_ODESOLVE]($exp$,$depvar=dv$,$indepvar$=$from .. to$
$ [,accuracy=a][,iterations=i]) $
where
$exp$ is the differential expression/equation,
$depvar$ is an identifier representing the dependent variable
(function to be found),
$indepvar$ is an identifier representing the independent variable,
$exp$ is an equation (or an expression implicitly set to zero) which
contains the first derivative of $depvar$ wrt $indepvar$,
$from$ is the starting point of integration,
$to$ is the endpoint of integration (allowed to be below $from$),
$dv$ is the initial value of $depvar$ in the point $indepvar=from$.
The ODE $exp$ is converted into an explicit form, which then is
used for a Runge Kutta iteration over the given range. The
number of steps is controlled by the value of $i$
(default: 20).
If the steps are too coarse to reach the desired
accuracy in the neighborhood of the starting point, the number is
increased automatically.
Result is a list of pairs, each representing a point of the
approximate solution of the ODE problem.
\end{description}
Example:
\begin{verbatim}
num_odesolve(df(y,x)=y,y=1,x=0 .. 1, iterations=5);
{{0.0,1.0},{0.2,1.2214},{0.4,1.49181796},{0.6,1.8221064563},
{0.8,2.2255208258},{1.0,2.7182511366}}
\end{verbatim}
Remarks:
\begin{enumerate}
\item[--] If in $exp$ the differential is not isolated on the lefthand side,
please ensure that the dependent variable is explicitly declared
using a \verb+DEPEND+ statement, e.g.
\begin{verbatim}
depend y,x;
\end{verbatim}
otherwise the formal derivative will be computed to zero by REDUCE.
\item[--] The REDUCE package SOLVE is used to convert the form into
an explicit ODE. If that process fails or has no unique result,
the evaluation is stopped with an error message.
\end{enumerate}
\section{Bounds of a Function}
Upper and lower bounds of a real valued function over an
interval or a rectangular multivariate domain are computed
by the operator BOUNDS. The algorithmic basis is the computation
with inequalities: starting from the interval(s) of the
variables, the bounds are propagated in the expression
using the rules for inequality computation. Some knowledge
about the behavior of special functions like ABS, SIN, COS, EXP, LOG,
fractional exponentials etc. is integrated and can be evaluated
if the operator BOUNDS is called with rounded mode on
(otherwise only algebraic evaluation rules are available).
If BOUNDS finds a singularity within an interval, the evaluation
is stopped with an error message indicating the problem part
of the expression.
Syntax:
\begin{description}
\item[BOUNDS]$(exp,var_1=l_1 .. u_1 [,var_2=l_2 .. u_2 \ldots])$
\item[{\it BOUNDS}]$(exp,\{var_1=l_1 .. u_1 [,var_2=l_2 .. u_2 \ldots]\})$
where $exp$ is the function to be investigated,
$var_1, var_2 , \ldots$ are the variables of exp,
$l_1, l_2 , \ldots$ and $u_1, u_2 , \ldots$ specify the area (intervals).
$BOUNDS$ computes upper and lower bounds for the expression in the
given area. An interval is returned.
\end{description}
Example:
\begin{verbatim}
bounds(sin x,x=1 .. 2);
{-1,1}
on rounded;
bounds(sin x,x=1 .. 2);
0.84147098481 .. 1
bounds(x**2+x,x=-0.5 .. 0.5);
- 0.5 .. 0.75
bounds((x-10)**2+x-10,x=9.5 .. 10.5);
- 0.5 .. 0.75
\end{verbatim}
\section{Curve Fitting}
The operator $NUM\_FIT$ finds for a set of
points the linear combination of a given set of
functions (function basis) which approximates the
points best under the objective of the least squares
criterion (minimum of the sum of the squares of the deviation).
The solution is found as zero of the
gradient vector of the sum of squared errors.
Syntax:
\begin{description}
\item[NUM\_FIT] $(vals,basis,var=pts)$
where $vals$ is a list of numeric values,
$var$ is a variable used for the approximation,
$pts$ is a list of coordinate values which correspond to $var$,
$basis$ is a set of functions varying in $var$ which is used
for the approximation.
\end{description}
The result is a list containing as first element the
function which approximates the given values, and as
second element a list of coefficients which were used
to build this function from the basis.
Example:
\begin{verbatim}
% approximate a set of factorials by a polynomial
pts:=for i:=1 step 1 until 5 collect i$
vals:=for i:=1 step 1 until 5 collect
for j:=1:i product j$
num_fit(vals,{1,x,x**2},x=pts);
2
{14.571428571*X - 61.428571429*X + 54.6,{54.6,
- 61.428571429,14.571428571}}
num_fit(vals,{1,x,x**2,x**3,x**4},x=pts);
4 3
{2.2083333234*X - 20.249999879*X
2
+ 67.791666154*X - 93.749999133*X
+ 44.999999525,
{44.999999525, - 93.749999133,67.791666154,
- 20.249999879,2.2083333234}}
\end{verbatim}
\section{Function Bases}
The following procedures compute sets of functions
e.g. to be used for approximation.
All procedures have
two parameters, the expression to be used as $variable$
(an identifier in most cases) and the
order of the desired system.
The functions are not scaled to a specific interval, but
the $variable$ can be accompanied by a scale factor
and/or a translation
in order to map the generic interval of orthogonality to another
(e.g. $(x- 1/2 ) * 2 pi$).
The result is a function list with ascending order, such that
the first element is the function of order zero and (for
the polynomial systems) the function of order $n$ is the $n+1$-th
element.
\begin{verbatim}
monomial_base(x,n) {1,x,...,x**n}
trigonometric_base(x,n) {1,sin x,cos x,sin(2x),cos(2x)...}
Bernstein_base(x,n) Bernstein polynomials
Legendre_base(x,n) Legendre polynomials
Laguerre_base(x,n) Laguerre polynomials
Hermite_base(x,n) Hermite polynomials
Chebyshew_base(x,n) Chebyshew polynomials
\end{verbatim}
Example:
\begin{verbatim}
Bernstein_base(x,5);
5 4 3 2
{ - X + 5*X - 10*X + 10*X - 5*X + 1,
4 3 2
5*X*(X - 4*X + 6*X - 4*X + 1),
2 3 2
10*X *( - X + 3*X - 3*X + 1),
3 2
10*X *(X - 2*X + 1),
4
5*X *( - X + 1),
5
X }
\end{verbatim}
\end{document}