reduce-historical/r38/packages/numeric/numeric.txt

                               NUMERIC


                           Herbert Melenk

        Konrad--Zuse--Zentrum fuer Informationstechnik Berlin


                 E--mail:  Melenk@sc.zib--berlin.de


The  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.


1 Syntax


1.1 Intervals, Starting Points

Intervals  are  generally  coded  as  lower  bound  and  upper  bound
connected by the operator  `..', usually associated to a variable  in
an equation.  E.g.

     x= (2.5 .. 3.5)

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 PART  operator as the first and  second part respectively.   A
starting point is specified  by an equation with a numeric  righthand
side, e.g.

     x=3.0

If  for  multivariate   applications  several  coordinates  must   be
specified by intervals  or as a starting point, these  specifications


                                  1



2 MINIMA                                                            2


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 automatical style,  while the  flat form is  more
convenient for direct interactive input.


1.2 Accuracy Control

The keyword  parameters accuracy=a  and iterations=i,  where a and  i
must be positive  integer numbers, control the iterative  algorithms:
                                                              -a
the iteration is  continued until the local  error is below 10  ;  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.


1.3 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

    on trnumeric;


2 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:

NUMMIN (exp,var  [=val ][,var [=val  ]...]
               1      1      2     2
     [,accuracy=a][,iterations=i])
     or



3 ROOTS OF FUNCTIONS/ SOLUTIONS OF EQUATIONS                        3


NUMMIN (exp,{var  [=val ][,var [=val  ]...]}
                1      1      2     2
     [,accuracy=a][,iterations=i])

     where exp is a function expression,
     var ,var  ,... are the  variables in exp and  val ,val ,...  are
        1    2                                        1    2
     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.
Examples:

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



3 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



3 ROOTS OF FUNCTIONS/ SOLUTIONS OF EQUATIONS                        4


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:

NUM_SOLVE (exp ,var  [=val ][,accuracy=a][,iterations=i])
      or      1    1      1

NUM_SOLVE ({exp ,... ,exp },var  [=val ],... ,var [=val ]
	       1         n     1      1          1     n
         [,accuracy=a][,iterations=i])
     or

NUM_ SOLVE ({exp ,... ,exp },{var  [=val ],... ,var [=val ]}
	       1         n      1      1          1     n
         [,accuracy=a][,iterations=i])
     where exp ,... ,exp  are function expressions,
              1         n
     var ,... ,var  are the variables,
        1         n
     val ,... ,val  are optional start values.
        1         n
     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.

Example:
    num_solve({sin x=cos y, x + y = 1},{x=1,y=2});


    {X= - 1.8561957251,Y=2.856195584}

    jacobian;

    [COS(X)  SIN(Y)]
    [              ]



4 INTEGRALS                                                         5


    [  1       1   ]



4 Integrals

For the  numerical evaluation of univariate  integrals over a  finite
interval the following strategy is used:

  1. If  the function has  an antiderivative in  close form which  is
     bounded in the integration interval, this is used.
  2. 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.
  3. If none  of these methods is successful, an adaptive  multilevel
     quadrature algorithm is used.

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:

NUMINT (exp,var  =(l ..u  )[,var =(l  ..u )... ]
               1    1   1       2   2    2
       [,accuracy=a][,iterations=i])
     where exp is the function to be integrated,

     var ,var  ,... are the integration variables,
        1    2
     l ,l  ,... are the lower bounds,
      1  2
     u ,u  ,... are the upper bounds.
      1  2
     Result is the value of the integral.
Example:

    num_int(sin x,x=(0 .. pi));

    2.0000010334



5 ORDINARY DIFFERENTIAL EQUATIONS                                   6


5 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:
NUMODESOLVE (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.

Example:

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


6 BOUNDS OF A FUNCTION                                              7


  {0.8,2.2255208258},{1.0,2.7182511366}}


Remarks:

  -- If  in exp  the differential  is not  isolated  on the  lefthand
     side,  please ensure that the  dependent variable is  explicitly
     declared using a DEPEND statement, e.g.
         depend y,x;

     otherwise  the formal  derivative will  be computed  to zero  by
     REDUCE.

  -- 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.


6 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:

BOUNDS (exp,var  =(l ..u )[,var  =(l ..u  )...])
               1    1   1      2    2   2
BOUNDS (exp,{var  =(l ..u )[,var  =(l ..u  )...]})
                1    1   1      2    2   2
     where exp is the function to be investigated,



7 CHEBYSHEV CURVE FITTING                                           8


     var ,var  ,... are the variables of exp,
        1    2
     l ,l  ,... and u ,u ,...  specify the area (intervals).
      1  2           1  2
     BOUNDS  computes upper and  lower bounds for  the expression  in
     the given area.  An interval is returned.

Example:


    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




7 Chebyshev Curve Fitting

The  operator  family  Chebyshev-...  implements  approximation   and
                                                            (a,b)
evaluation  of functions  by the  Chebyshev method.    Let T      (x)

be the  Chebyshev polynomial of order  n transformed to the ninterval
(a,b).    Then a  function f(x)  can be approximated  in (a,b)  by  a

series

       N     (a,b)
f(x)?si=0 c T     (x)
           i i

The operator Chebyshev-fit computes this approximation and returns  a
list, which  has as first element the  sum expressed as a  polynomial



7 CHEBYSHEV CURVE FITTING                                           9


and as  second element  the sequence  of Chebyshev  coefficients c  .
                                                                  i
Chebyshevdf and  Chebyshevint 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 Chebysheveval 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.
CHEBYSHEVFIT (fcn,var=(lo..hi),n)

CHEBYSHEVEVAL (coeffs,var=(lo..hi),var=pt)
CHEBYSHEVDF (coeffs,var=(lo..hi))

CHEBYSHEVINT (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  CHEBYSHEVCOEFF, -DF
     or -INT.

Example:


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



8 GENERAL CURVE FITTING                                            10


8 General Curve Fitting


The operator NUMFIT 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:

NUMFIT (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.
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:


     % 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



8 GENERAL CURVE FITTING                                            11


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