reduce-historical/r36/doc/NUMERIC.TEX

\documentstyle[11pt,reduce]{article}
\date{}
\title{NUMERIC}
\author{Herbert Melenk \\
Konrad--Zuse--Zentrum f\"ur Informationstechnik Berlin \\
Heilbronner Strasse 10 \\
D--10711 Berlin -- Wilmersdorf \\
Federal Republic of Germany \\[0.05in]
E--mail: melenk@sc.zib--berlin.de}
\begin{document}
\maketitle

\index{NUMERIC package}
The {\small NUMERIC} package implements some numerical (approximative)
algorithms for \REDUCE, based on the \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 list form is more appropriate when the
parameters are built from other REDUCE calculations in an
automatic 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.

NUM\_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.

NUM\_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 a side effect in 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}

For the numerical evaluation of univariate integrals over a finite
interval the following strategy is used:
\begin{enumerate}
\item If the function has an antiderivative in close form
    which is bounded in the integration interval, this
    is used.
\item Otherwise a Chebyshev approximation is computed,
    starting with order 20, eventually up to order 80.
    If that is recognized as sufficiently convergent
    it is used for computing the integral by directly
    integrating the coefficient sequence.
\item If none of these methods is successful, an
    adaptive multilevel quadrature algorithm is used.
\end{enumerate}
For multivariate integrals only the adaptive quadrature is used.
This 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.25 .. 0.75

\end{verbatim}

\section{Chebyshev Curve Fitting}

The operator family $Chebyshev\_\ldots$ implements approximation
and evaluation of functions by the Chebyshev method.
Let $T_n^{(a,b)}(x)$ be the Chebyshev polynomial of order $n$
transformed to the interval $(a,b)$. Then a function $f(x)$ can be
approximated in $(a,b)$ by a series

$f(x) \approx \sum_{i=0}^N c_i T_i^{(a,b)}(x)$

The operator $Chebyshev\_fit$ computes this approximation and
returns a list, which has as first element the sum expressed
as a polynomial and as second element the sequence
of Chebyshev coefficients ${c_i}$.
$Chebyshev\_df$ and $Chebyshev\_int$ transform a Chebyshev
coefficient list into the coefficients of the corresponding
derivative or integral respectively. For evaluating a Chebyshev
approximation at a given point in the basic interval the
operator $Chebyshev\_eval$ can be used. Note that
$Chebyshev\_eval$ is based on a recurrence relation which is
in general more stable than a direct evaluation of the
complete polynomial.

\begin{description}
\item[CHEBYSHEV\_FIT] $(fcn,var=(lo .. hi),n)$

\item[CHEBYSHEV\_EVAL] $(coeffs,var=(lo .. hi),var=pt)$

\item[CHEBYSHEV\_DF] $(coeffs,var=(lo .. hi))$

\item[CHEBYSHEV\_INT] $(coeffs,var=(lo .. hi))$

where $fcn$ is an algebraic expression (the function to be
fitted), $var$ is the variable of $fcn$, $lo$ and $hi$ are
numerical real values which describe an interval ($lo < hi$),
$n$ is the approximation order,an integer $>0$, set to 20 if missing,
$pt$ is a numerical value in the interval and $coeffs$ is
a series of Chebyshev coefficients, computed by one of
$CHEBYSHEV\_COEFF$, $\_DF$ or $\_INT$.
\end{description}

Example:

\begin{verbatim}

on rounded;

w:=chebyshev_fit(sin x/x,x=(1 .. 3),5);

               3           2
w := {0.03824*x  - 0.2398*x  + 0.06514*x + 0.9778,

      {0.8991,-0.4066,-0.005198,0.009464,-0.00009511}}

chebyshev_eval(second w, x=(1 .. 3), x=2.1);

0.4111

\end{verbatim}

\section{General 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
     Chebyshev_base_T(x,n)    Chebyshev polynomials first kind
     Chebyshev_base_U(x,n)    Chebyshev polynomials second kind

\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}