This is maxima.info, produced by makeinfo version 4.7 from maxima.texi.

   This is a Texinfo Maxima Manual

   Copyright 1994,2001 William F. Schelter

START-INFO-DIR-ENTRY
* Maxima: (maxima).     A computer algebra system.
END-INFO-DIR-ENTRY


File: maxima.info,  Node: Functions and Variables for impdiff,  Prev: impdiff,  Up: impdiff

54.1 Functions and Variables for impdiff
========================================

 -- Function: implicit_derivative
          (<f>,<indvarlist>,<orderlist>,<depvar>)
     This subroutine computes implicit derivatives of multivariable
     functions.  <f> is an array function, the indexes are the
     derivative degree in the <indvarlist> order; <indvarlist> is the
     independent variable list; <orderlist> is the order desired; and
     <depvar> is the dependent variable.

     To use this function write first `load("impdiff")'.


File: maxima.info,  Node: implicit_plot,  Next: interpol,  Prev: impdiff,  Up: Top

55 implicit_plot
****************

* Menu:

* Functions and Variables for implicit_plot::


File: maxima.info,  Node: Functions and Variables for implicit_plot,  Prev: implicit_plot,  Up: implicit_plot

55.1 Functions and Variables for implicit_plot
==============================================

 -- Function: implicit_plot (<expr>, <x_range>, <y_range>)
 -- Function: implicit_plot ([<expr_1>, ..., <expr_n>], <x_range>,
          <y_range>)
     Displays a plot of one or more expressions in implicit form.
     <expr> is the expression to be plotted, <x_range> the range of the
     horizontal axis and <y_range> the range of vertical axis.
     `implicit_plot' only works with the gnuplot driver.
     `implicit_plot' respects global setting for the gnuplot driver set
     by the <set_plot_option> function. Options can also be passed to
     `implicit_plot' function as optional arguments.

     `implicit_plot' works by tracking sign changes on the area given
     by <x_range> and <y_range> and can fail for complicated
     expressions.

     `load(implicit_plot)' loads this function.

     Example:
          (%i1) implicit_plot (x^2 = y^3 - 3*y + 1, [x, -4, 4], [y, -4, 4],
           [gnuplot_preamble, "set zeroaxis"]);



File: maxima.info,  Node: interpol,  Next: lbfgs,  Prev: implicit_plot,  Up: Top

56 interpol
***********

* Menu:

* Introduction to interpol::
* Functions and Variables for interpol::


File: maxima.info,  Node: Introduction to interpol,  Next: Functions and Variables for interpol,  Prev: interpol,  Up: interpol

56.1 Introduction to interpol
=============================

Package `interpol' defines de Lagrangian, the linear and the cubic
splines methods for polynomial interpolation.

   For comments, bugs or suggestions, please contact me at <'mario AT
edu DOT xunta DOT es'>.


File: maxima.info,  Node: Functions and Variables for interpol,  Prev: Introduction to interpol,  Up: interpol

56.2 Functions and Variables for interpol
=========================================

 -- Function: lagrange (<points>)
 -- Function: lagrange (<points>, <option>)
     Computes the polynomial interpolation by the Lagrangian method.
     Argument <points> must be either:

        * a two column matrix, `p:matrix([2,4],[5,6],[9,3])',

        * a list of pairs, `p: [[2,4],[5,6],[9,3]]',

        * a list of numbers, `p: [4,6,3]', in which case the abscissas
          will be assigned automatically to 1, 2, 3, etc.

     In the first two cases the pairs are ordered with respect to the
     first coordinate before making computations.

     With the <option> argument it is possible to select the name for
     the independent variable, which is `'x' by default; to define
     another one, write something like `varname='z'.

     Note that when working with high degree polynomials, floating
     point evaluations are inestable.

     Examples:

          (%i1) load(interpol)$
          (%i2) p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$
          (%i3) lagrange(p);
                           4        3         2
                       73 x    701 x    8957 x    5288 x   186
          (%o3)        ----- - ------ + ------- - ------ + ---
                        420     210       420      105      5
          (%i4) f(x):=''%;
                               4        3         2
                           73 x    701 x    8957 x    5288 x   186
          (%o4)    f(x) := ----- - ------ + ------- - ------ + ---
                            420     210       420      105      5
          (%i5) /* Evaluate the polynomial at some points */
                map(f,[2.3,5/7,%pi]);
                                       919062
          (%o5)  [- 1.567534999999992, ------,
                                       84035
                                   4          3           2
                             73 %pi    701 %pi    8957 %pi    5288 %pi   186
                             ------- - -------- + --------- - -------- + ---]
                               420       210         420        105       5
          (%i6) %,numer;
          (%o6) [- 1.567534999999992, 10.9366573451538, 2.89319655125692]
          (%i7) load(draw)$  /* load draw package */
          (%i8) /* Plot the polynomial together with points */
                draw2d(
                  color      = red,
                  key        = "Lagrange polynomial",
                  explicit(f(x),x,0,10),
                  point_size = 3,
                  color      = blue,
                  key        = "Sample points",
                  points(p))$
          (%i9) /* Change variable name */
                lagrange(p, varname=w);
                           4        3         2
                       73 w    701 w    8957 w    5288 w   186
          (%o9)        ----- - ------ + ------- - ------ + ---
                        420     210       420      105      5


 -- Function: charfun2 (<x>, <a>, <b>)
     Returns `true' if number <x> belongs to the interval [a, b), and
     `false' otherwise.

 -- Function: linearinterpol (<points>)
 -- Function: linearinterpol (<points>, <option>)
     Computes the polynomial interpolation by the linear method.
     Argument <points> must be either:

        * a two column matrix, `p:matrix([2,4],[5,6],[9,3])',

        * a list of pairs, `p: [[2,4],[5,6],[9,3]]',

        * a list of numbers, `p: [4,6,3]', in which case the abscissas
          will be assigned automatically to 1, 2, 3, etc.

     In the first two cases the pairs are ordered with respect to the
     first coordinate before making computations.

     With the <option> argument it is possible to select the name for
     the independent variable, which is `'x' by default; to define
     another one, write something like `varname='z'.

     Examples:
          (%i1) load(interpol)$
          (%i2) p: matrix([7,2],[8,3],[1,5],[3,2],[6,7])$
          (%i3) linearinterpol(p);
                  13   3 x
          (%o3)  (-- - ---) charfun2(x, minf, 3)
                  2     2
           + (x - 5) charfun2(x, 7, inf) + (37 - 5 x) charfun2(x, 6, 7)
              5 x
           + (--- - 3) charfun2(x, 3, 6)
               3

          (%i4) f(x):=''%;
                          13   3 x
          (%o4)  f(x) := (-- - ---) charfun2(x, minf, 3)
                          2     2
           + (x - 5) charfun2(x, 7, inf) + (37 - 5 x) charfun2(x, 6, 7)
              5 x
           + (--- - 3) charfun2(x, 3, 6)
               3
          (%i5)  /* Evaluate the polynomial at some points */
                 map(f,[7.3,25/7,%pi]);
                                      62  5 %pi
          (%o5)                 [2.3, --, ----- - 3]
                                      21    3
          (%i6) %,numer;
          (%o6)  [2.3, 2.952380952380953, 2.235987755982989]
          (%i7) load(draw)$  /* load draw package */
          (%i8)  /* Plot the polynomial together with points */
                 draw2d(
                   color      = red,
                   key        = "Linear interpolator",
                   explicit(f(x),x,-5,20),
                   point_size = 3,
                   color      = blue,
                   key        = "Sample points",
                   points(args(p)))$
          (%i9)  /* Change variable name */
                 linearinterpol(p, varname='s);
                 13   3 s
          (%o9) (-- - ---) charfun2(s, minf, 3)
                 2     2
           + (s - 5) charfun2(s, 7, inf) + (37 - 5 s) charfun2(s, 6, 7)
              5 s
           + (--- - 3) charfun2(s, 3, 6)
               3


 -- Function: cspline (<points>)
 -- Function: cspline (<points>, <option1>, <option2>, ...)
     Computes the polynomial interpolation by the cubic splines method.
     Argument <points> must be either:

        * a two column matrix, `p:matrix([2,4],[5,6],[9,3])',

        * a list of pairs, `p: [[2,4],[5,6],[9,3]]',

        * a list of numbers, `p: [4,6,3]', in which case the abscissas
          will be assigned automatically to 1, 2, 3, etc.

     In the first two cases the pairs are ordered with respect to the
     first coordinate before making computations.

     There are three options to fit specific needs:
        * `'d1', default `'unknown', is the first derivative at x_1; if
          it is `'unknown', the second derivative at x_1 is made equal
          to 0 (natural cubic spline); if it is equal to a number, the
          second derivative is calculated based on this number.

        * `'dn', default `'unknown', is the first derivative at x_n; if
          it is `'unknown', the second derivative at x_n is made equal
          to 0 (natural cubic spline); if it is equal to a number, the
          second derivative is calculated based on this number.

        * `'varname', default `'x', is the name of the independent
          variable.

     Examples:
          (%i1) load(interpol)$
          (%i2) p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$
          (%i3) /* Unknown first derivatives at the extremes
                   is equivalent to natural cubic splines */
                cspline(p);
                        3         2
                  1159 x    1159 x    6091 x   8283
          (%o3)  (------- - ------- - ------ + ----) charfun2(x, minf, 3)
                   3288      1096      3288    1096
                      3         2
                2587 x    5174 x    494117 x   108928
           + (- ------- + ------- - -------- + ------) charfun2(x, 7, inf)
                 1644       137       1644      137
                    3          2
              4715 x    15209 x    579277 x   199575
           + (------- - -------- + -------- - ------) charfun2(x, 6, 7)
               1644       274        1644      274
                      3         2
                3287 x    2223 x    48275 x   9609
           + (- ------- + ------- - ------- + ----) charfun2(x, 3, 6)
                 4932       274      1644     274

          (%i4) f(x):=''%$
          (%i5) /* Some evaluations */
                map(f,[2.3,5/7,%pi]), numer;
          (%o5) [1.991460766423356, 5.823200187269903, 2.227405312429507]
          (%i6) load(draw)$  /* load draw package */
          (%i7) /* Plotting interpolating function */
                draw2d(
                  color      = red,
                  key        = "Cubic splines",
                  explicit(f(x),x,0,10),
                  point_size = 3,
                  color      = blue,
                  key        = "Sample points",
                  points(p))$
          (%i8) /* New call, but giving values at the derivatives */
                cspline(p,d1=0,dn=0);
                        3          2
                  1949 x    11437 x    17027 x   1247
          (%o8)  (------- - -------- + ------- + ----) charfun2(x, minf, 3)
                   2256       2256      2256     752
                      3          2
                1547 x    35581 x    68068 x   173546
           + (- ------- + -------- - ------- + ------) charfun2(x, 7, inf)
                  564       564        141      141
                   3          2
              607 x    35147 x    55706 x   38420
           + (------ - -------- + ------- - -----) charfun2(x, 6, 7)
               188       564        141      47
                      3         2
                3895 x    1807 x    5146 x   2148
           + (- ------- + ------- - ------ + ----) charfun2(x, 3, 6)
                 5076       188      141      47
          (%i8) /* Defining new interpolating function */
                g(x):=''%$
          (%i9) /* Plotting both functions together */
                draw2d(
                  color      = black,
                  key        = "Cubic splines (default)",
                  explicit(f(x),x,0,10),
                  color      = red,
                  key        = "Cubic splines (d1=0,dn=0)",
                  explicit(g(x),x,0,10),
                  point_size = 3,
                  color      = blue,
                  key        = "Sample points",
                  points(p))$



File: maxima.info,  Node: lbfgs,  Next: lindstedt,  Prev: interpol,  Up: Top

57 lbfgs
********

* Menu:

* Introduction to lbfgs::
* Functions and Variables for lbfgs::


File: maxima.info,  Node: Introduction to lbfgs,  Next: Functions and Variables for lbfgs,  Prev: Top,  Up: Top

57.1 Introduction to lbfgs
==========================

`lbfgs' is an implementation of the L-BFGS algorithm [1] to solve
unconstrained minimization problems via a limited-memory quasi-Newton
(BFGS) algorithm.  It is called a limited-memory method because a
low-rank approximation of the Hessian matrix inverse is stored instead
of the entire Hessian inverse.  The program was originally written in
Fortran [2] by Jorge Nocedal, incorporating some functions originally
written by Jorge J. More' and David J. Thuente, and translated into
Lisp automatically via the program `f2cl'.  The Maxima package `lbfgs'
comprises the translated code plus an interface function which manages
some details.

   References:

   [1] D. Liu and J. Nocedal. "On the limited memory BFGS method for
large scale optimization". Mathematical Programming B 45:503-528 (1989)

   [2] `http://netlib.org/opt/lbfgs_um.shar'


File: maxima.info,  Node: Functions and Variables for lbfgs,  Prev: Introduction to lbfgs,  Up: Top

57.2 Functions and Variables for lbfgs
======================================

 -- Function: lbfgs (<FOM>, <X>, <X0>, <epsilon>, <iprint>)
     Finds an approximate solution of the unconstrained minimization of
     the figure of merit <FOM> over the list of variables <X>, starting
     from initial estimates <X0>, such that norm grad FOM < epsilon
     max(1, norm X).

     The algorithm applied is a limited-memory quasi-Newton (BFGS)
     algorithm [1].  It is called a limited-memory method because a
     low-rank approximation of the Hessian matrix inverse is stored
     instead of the entire Hessian inverse.  Each iteration of the
     algorithm is a line search, that is, a search along a ray in the
     variables <X>, with the search direction computed from the
     approximate Hessian inverse.  The FOM is always decreased by a
     successful line search.  Usually (but not always) the norm of the
     gradient of FOM also decreases.

     <iprint> controls progress messages printed by `lbfgs'.

    `iprint[1]'
          `<iprint>[1]' controls the frequency of progress messages.
         `iprint[1] < 0'
               No progress messages.

         `iprint[1] = 0'
               Messages at the first and last iterations.

         `iprint[1] > 0'
               Print a message every `<iprint>[1]' iterations.

    `iprint[2]'
          `<iprint>[2]' controls the verbosity of progress messages.
         `iprint[2] = 0'
               Print out iteration count, number of evaluations of
               <FOM>, value of <FOM>, norm of the gradient of <FOM>,
               and step length.

         `iprint[2] = 1'
               Same as `<iprint>[2] = 0', plus <X0> and the gradient of
               <FOM> evaluated at <X0>.

         `iprint[2] = 2'
               Same as `<iprint>[2] = 1', plus values of <X> at each
               iteration.

         `iprint[2] = 3'
               Same as `<iprint>[2] = 2', plus the gradient of <FOM> at
               each iteration.

     The columns printed by `lbfgs' are the following.

    `I'
          Number of iterations. It is incremented for each line search.

    `NFN'
          Number of evaluations of the figure of merit.

    `FUNC'
          Value of the figure of merit at the end of the most recent
          line search.

    `GNORM'
          Norm of the gradient of the figure of merit at the end of the
          most recent line search.

    `STEPLENGTH'
          An internal parameter of the search algorithm.

     Additional information concerning details of the algorithm are
     found in the comments of the original Fortran code [2].

     See also `lbfgs_nfeval_max' and `lbfgs_ncorrections'.

     References:

     [1] D. Liu and J. Nocedal. "On the limited memory BFGS method for
     large scale optimization". Mathematical Programming B 45:503-528
     (1989)

     [2] `http://netlib.org/opt/lbfgs_um.shar'

     Examples:

     The same FOM as computed by FGCOMPUTE in the program sdrive.f in
     the LBFGS package from Netlib.  Note that the variables in
     question are subscripted variables.  The FOM has an exact minimum
     equal to zero at u[k] = 1 for k = 1, ..., 8.

          (%i1) load (lbfgs);
          (%o1)   /usr/share/maxima/5.10.0cvs/share/lbfgs/lbfgs.mac
          (%i2) t1[j] := 1 - u[j];
          (%o2)                     t1  := 1 - u
                                      j         j
          (%i3) t2[j] := 10*(u[j + 1] - u[j]^2);
                                                    2
          (%o3)                t2  := 10 (u      - u )
                                 j         j + 1    j
          (%i4) n : 8;
          (%o4)                           8
          (%i5) FOM : sum (t1[2*j - 1]^2 + t2[2*j - 1]^2, j, 1, n/2);
                           2 2           2              2 2           2
          (%o5) 100 (u  - u )  + (1 - u )  + 100 (u  - u )  + (1 - u )
                      8    7           7           6    5           5
                               2 2           2              2 2           2
                  + 100 (u  - u )  + (1 - u )  + 100 (u  - u )  + (1 - u )
                          4    3           3           2    1           1
          (%i6) lbfgs (FOM, '[u[1],u[2],u[3],u[4],u[5],u[6],u[7],u[8]],
                 [-1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1], 1e-3, [1, 0]);
          *************************************************
            N=    8   NUMBER OF CORRECTIONS=25
                 INITIAL VALUES
           F=  9.680000000000000D+01   GNORM=  4.657353755084532D+02
          *************************************************

           I NFN   FUNC                    GNORM                   STEPLENGTH

           1   3   1.651479526340304D+01   4.324359291335977D+00   7.926153934390631D-04
           2   4   1.650209316638371D+01   3.575788161060007D+00   1.000000000000000D+00
           3   5   1.645461701312851D+01   6.230869903601577D+00   1.000000000000000D+00
           4   6   1.636867301275588D+01   1.177589920974980D+01   1.000000000000000D+00
           5   7   1.612153014409201D+01   2.292797147151288D+01   1.000000000000000D+00
           6   8   1.569118407390628D+01   3.687447158775571D+01   1.000000000000000D+00
           7   9   1.510361958398942D+01   4.501931728123680D+01   1.000000000000000D+00
           8  10   1.391077875774294D+01   4.526061463810632D+01   1.000000000000000D+00
           9  11   1.165625686278198D+01   2.748348965356917D+01   1.000000000000000D+00
          10  12   9.859422687859137D+00   2.111494974231644D+01   1.000000000000000D+00
          11  13   7.815442521732281D+00   6.110762325766556D+00   1.000000000000000D+00
          12  15   7.346380905773160D+00   2.165281166714631D+01   1.285316401779533D-01
          13  16   6.330460634066370D+00   1.401220851762050D+01   1.000000000000000D+00
          14  17   5.238763939851439D+00   1.702473787613255D+01   1.000000000000000D+00
          15  18   3.754016790406701D+00   7.981845727704576D+00   1.000000000000000D+00
          16  20   3.001238402309352D+00   3.925482944716691D+00   2.333129631296807D-01
          17  22   2.794390709718290D+00   8.243329982546473D+00   2.503577283782332D-01
          18  23   2.563783562918759D+00   1.035413426521790D+01   1.000000000000000D+00
          19  24   2.019429976377856D+00   1.065187312346769D+01   1.000000000000000D+00
          20  25   1.428003167670903D+00   2.475962450826961D+00   1.000000000000000D+00
          21  27   1.197874264861340D+00   8.441707983493810D+00   4.303451060808756D-01
          22  28   9.023848941942773D-01   1.113189216635162D+01   1.000000000000000D+00
          23  29   5.508226405863770D-01   2.380830600326308D+00   1.000000000000000D+00
          24  31   3.902893258815567D-01   5.625595816584421D+00   4.834988416524465D-01
          25  32   3.207542206990315D-01   1.149444645416472D+01   1.000000000000000D+00
          26  33   1.874468266362791D-01   3.632482152880997D+00   1.000000000000000D+00
          27  34   9.575763380706598D-02   4.816497446154354D+00   1.000000000000000D+00
          28  35   4.085145107543406D-02   2.087009350166495D+00   1.000000000000000D+00
          29  36   1.931106001379290D-02   3.886818608498966D+00   1.000000000000000D+00
          30  37   6.894000721499670D-03   3.198505796342214D+00   1.000000000000000D+00
          31  38   1.443296033051864D-03   1.590265471025043D+00   1.000000000000000D+00
          32  39   1.571766603154336D-04   3.098257063980634D-01   1.000000000000000D+00
          33  40   1.288011776581970D-05   1.207784183577257D-02   1.000000000000000D+00
          34  41   1.806140173752971D-06   4.587890233385193D-02   1.000000000000000D+00
          35  42   1.769004645459358D-07   1.790537375052208D-02   1.000000000000000D+00
          36  43   3.312164100763217D-10   6.782068426119681D-04   1.000000000000000D+00


           THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.
           IFLAG = 0
          (%o6) [u  = 1.000005339815974, u  = 1.000009942839805,
                  1                       2
          u  = 1.000005339815974, u  = 1.000009942839805,
           3                       4
          u  = 1.000005339815974, u  = 1.000009942839805,
           5                       6
          u  = 1.000005339815974, u  = 1.000009942839805]
           7                       8

     A regression problem.  The FOM is the mean square difference
     between the predicted value F(X[i]) and the observed value Y[i].
     The function F is a bounded monotone function (a so-called
     "sigmoidal" function).  In this example, `lbfgs' computes
     approximate values for the parameters of F and `plot2d' displays a
     comparison of F with the observed data.

          (%i1) load (lbfgs);
          (%o1)   /usr/share/maxima/5.10.0cvs/share/lbfgs/lbfgs.mac
          (%i2) FOM : '((1/length(X))*sum((F(X[i]) - Y[i])^2, i, 1,
                                                          length(X)));
                                         2
                         sum((F(X ) - Y ) , i, 1, length(X))
                                 i     i
          (%o2)          -----------------------------------
                                      length(X)
          (%i3) X : [1, 2, 3, 4, 5];
          (%o3)                    [1, 2, 3, 4, 5]
          (%i4) Y : [0, 0.5, 1, 1.25, 1.5];
          (%o4)                [0, 0.5, 1, 1.25, 1.5]
          (%i5) F(x) := A/(1 + exp(-B*(x - C)));
                                             A
          (%o5)            F(x) := ----------------------
                                   1 + exp((- B) (x - C))
          (%i6) ''FOM;
                          A               2            A                2
          (%o6) ((----------------- - 1.5)  + (----------------- - 1.25)
                    - B (5 - C)                  - B (4 - C)
                  %e            + 1            %e            + 1
                      A             2            A               2
           + (----------------- - 1)  + (----------------- - 0.5)
                - B (3 - C)                - B (2 - C)
              %e            + 1          %e            + 1
                       2
                      A
           + --------------------)/5
                - B (1 - C)     2
             (%e            + 1)
          (%i7) estimates : lbfgs (FOM, '[A, B, C], [1, 1, 1], 1e-4, [1, 0]);
          *************************************************
            N=    3   NUMBER OF CORRECTIONS=25
                 INITIAL VALUES
           F=  1.348738534246918D-01   GNORM=  2.000215531936760D-01
          *************************************************

          I  NFN  FUNC                    GNORM                   STEPLENGTH
          1    3  1.177820636622582D-01   9.893138394953992D-02   8.554435968992371D-01
          2    6  2.302653892214013D-02   1.180098521565904D-01   2.100000000000000D+01
          3    8  1.496348495303005D-02   9.611201567691633D-02   5.257340567840707D-01
          4    9  7.900460841091139D-03   1.325041647391314D-02   1.000000000000000D+00
          5   10  7.314495451266917D-03   1.510670810312237D-02   1.000000000000000D+00
          6   11  6.750147275936680D-03   1.914964958023047D-02   1.000000000000000D+00
          7   12  5.850716021108205D-03   1.028089194579363D-02   1.000000000000000D+00
          8   13  5.778664230657791D-03   3.676866074530332D-04   1.000000000000000D+00
          9   14  5.777818823650782D-03   3.010740179797255D-04   1.000000000000000D+00


           THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.
           IFLAG = 0
          (%o7) [A = 1.461933911464101, B = 1.601593973254802,
                                                     C = 2.528933072164854]
          (%i8) plot2d ([F(x), [discrete, X, Y]], [x, -1, 6]), ''estimates;
          (%o8)


 -- Variable: lbfgs_nfeval_max
     Default value: 100

     `lbfgs_nfeval_max' is the maximum number of evaluations of the
     figure of merit (FOM) in `lbfgs'.  When `lbfgs_nfeval_max' is
     reached, `lbfgs' returns the result of the last successful line
     search.


 -- Variable: lbfgs_ncorrections
     Default value: 25

     `lbfgs_ncorrections' is the number of corrections applied to the
     approximate inverse Hessian matrix which is maintained by `lbfgs'.



File: maxima.info,  Node: lindstedt,  Next: linearalgebra,  Prev: lbfgs,  Up: Top

58 lindstedt
************

* Menu:

* Functions and Variables for lindstedt::


File: maxima.info,  Node: Functions and Variables for lindstedt,  Prev: lindstedt,  Up: lindstedt

58.1 Functions and Variables for lindstedt
==========================================

 -- Function: Lindstedt (<eq>,<pvar>,<torder>,<ic>)
     This is a first pass at a Lindstedt code.  It can solve problems
     with initial conditions entered, which can be arbitrary constants,
     (just not <%k1> and <%k2>) where the initial conditions on the
     perturbation equations are z[i]=0, z'[i]=0 for i>0. <ic> is the
     list of initial conditions.

     Problems occur when initial conditions are not given, as the
     constants in the perturbation equations are the same as the zero
     order equation solution.  Also, problems occur when the initial
     conditions for the perturbation equations are not z[i]=0, z'[i]=0
     for i>0, such as the Van der Pol equation.

     Example:
          (%i1) load("makeOrders")$

          (%i2) load("lindstedt")$

          (%i3) Lindstedt('diff(x,t,2)+x-(e*x^3)/6,e,2,[1,0]);
                    2
                   e  (cos(5 T) - 24 cos(3 T) + 23 cos(T))
          (%o3) [[[---------------------------------------
                                    36864
             e (cos(3 T) - cos(T))
           - --------------------- + cos(T)],
                      192
                    2
                 7 e    e
          T = (- ---- - -- + 1) t]]
                 3072   16

     To use this function write first `load("makeOrders")' and
     `load("lindstedt")'.


File: maxima.info,  Node: linearalgebra,  Next: lsquares,  Prev: lindstedt,  Up: Top

59 linearalgebra
****************

* Menu:

* Introduction to linearalgebra::
* Functions and Variables for linearalgebra::


File: maxima.info,  Node: Introduction to linearalgebra,  Next: Functions and Variables for linearalgebra,  Prev: linearalgebra,  Up: linearalgebra

59.1 Introduction to linearalgebra
==================================

`linearalgebra' is a collection of functions for linear algebra.

   Example:

     (%i1) M : matrix ([1, 2], [1, 2]);
                                 [ 1  2 ]
     (%o1)                       [      ]
                                 [ 1  2 ]
     (%i2) nullspace (M);
                                    [  1  ]
                                    [     ]
     (%o2)                     span([   1 ])
                                    [ - - ]
                                    [   2 ]
     (%i3) columnspace (M);
                                     [ 1 ]
     (%o3)                      span([   ])
                                     [ 1 ]
     (%i4) ptriangularize (M - z*ident(2), z);
                              [ 1   2 - z   ]
     (%o4)                    [             ]
                              [           2 ]
                              [ 0  3 z - z  ]
     (%i5) M : matrix ([1, 2, 3], [4, 5, 6], [7, 8, 9]) - z*ident(3);
                          [ 1 - z    2      3   ]
                          [                     ]
     (%o5)                [   4    5 - z    6   ]
                          [                     ]
                          [   7      8    9 - z ]
     (%i6) MM : ptriangularize (M, z);
                   [ 4  5 - z            6            ]
                   [                                  ]
                   [                2                 ]
                   [     66        z    102 z   132   ]
                   [ 0   --      - -- + ----- + ---   ]
     (%o6)         [     49        7     49     49    ]
                   [                                  ]
                   [               3        2         ]
                   [           49 z    245 z    147 z ]
                   [ 0    0    ----- - ------ - ----- ]
                   [            264      88      44   ]
     (%i7) algebraic : true;
     (%o7)                         true
     (%i8) tellrat (MM [3, 3]);
                              3       2
     (%o8)                  [z  - 15 z  - 18 z]
     (%i9) MM : ratsimp (MM);
                    [ 4  5 - z           6           ]
                    [                                ]
                    [                2               ]
     (%o9)          [     66      7 z  - 102 z - 132 ]
                    [ 0   --    - ------------------ ]
                    [     49              49         ]
                    [                                ]
                    [ 0    0             0           ]
     (%i10) nullspace (MM);
                             [        1         ]
                             [                  ]
                             [   2              ]
                             [  z  - 14 z - 16  ]
                             [  --------------  ]
     (%o10)             span([        8         ])
                             [                  ]
                             [    2             ]
                             [   z  - 18 z - 12 ]
                             [ - -------------- ]
                             [         12       ]
     (%i11) M : matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12],
                        [13, 14, 15, 16]);
                            [ 1   2   3   4  ]
                            [                ]
                            [ 5   6   7   8  ]
     (%o11)                 [                ]
                            [ 9   10  11  12 ]
                            [                ]
                            [ 13  14  15  16 ]
     (%i12) columnspace (M);
                                [ 1  ]  [ 2  ]
                                [    ]  [    ]
                                [ 5  ]  [ 6  ]
     (%o12)                span([    ], [    ])
                                [ 9  ]  [ 10 ]
                                [    ]  [    ]
                                [ 13 ]  [ 14 ]
     (%i13) apply ('orthogonal_complement, args (nullspace (transpose (M))));
                                [ 0 ]  [  1  ]
                                [   ]  [     ]
                                [ 1 ]  [  0  ]
     (%o13)                span([   ], [     ])
                                [ 2 ]  [ - 1 ]
                                [   ]  [     ]
                                [ 3 ]  [ - 2 ]


File: maxima.info,  Node: Functions and Variables for linearalgebra,  Prev: Introduction to linearalgebra,  Up: linearalgebra

59.2 Functions and Variables for linearalgebra
==============================================

 -- Function: addmatrices (<f>, <M_1>, ..., <M_n>)
     Using the function <f> as the addition function, return the sum of
     the matrices <M_1>, ..., <M_n>. The function <f> must accept any
     number of arguments (a Maxima nary function).

     Examples:

          (%i1) m1 : matrix([1,2],[3,4])$
          (%i2) m2 : matrix([7,8],[9,10])$
          (%i3) addmatrices('max,m1,m2);
          (%o3) matrix([7,8],[9,10])
          (%i4) addmatrices('max,m1,m2,5*m1);
          (%o4) matrix([7,10],[15,20])


 -- Function: blockmatrixp (<M>)
     Return true if and only if <M> is a matrix and every entry of <M>
     is a matrix.


 -- Function: columnop (<M>, <i>, <j>, <theta>)
     If <M> is a matrix, return the matrix that results from doing the
     column operation `C_i <- C_i - <theta> * C_j'. If <M> doesn't have
     a row <i> or <j>, signal an error.


 -- Function: columnswap (<M>, <i>, <j>)
     If <M> is a matrix, swap columns <i> and <j>.  If <M> doesn't have
     a column <i> or <j>, signal an error.


 -- Function: columnspace (<M>)
     If <M> is a matrix, return `span (v_1, ..., v_n)', where the set
     `{v_1, ..., v_n}' is a basis for the column space of <M>.  The span
     of the empty set is `{0}'. Thus, when the column space has only
     one member, return `span ()'.


 -- Function: copy (<e>)
     Return a copy of the Maxima expression <e>. Although <e> can be any
     Maxima expression, the copy function is the most useful when <e>
     is either a list or a matrix; consider:
          (%i1) m : [1,[2,3]]$
          (%i2) mm : m$
          (%i3) mm[2][1] : x$
          (%i4) m;
          (%o4)                      [1,[x,3]]
          (%i5) mm;
          (%o5)                      [1,[x,3]]
     Let's try the same experiment, but this time let <mm> be a copy of
     <m>
          (%i6) m : [1,[2,3]]$
          (%i7) mm : copy(m)$
          (%i8) mm[2][1] : x$
          (%i9) m;
          (%o9)                     [1,[2,3]]
          (%i10) mm;
          (%o10)                    [1,[x,3]]
     This time, the assignment to <mm> does not change the value of <m>.


 -- Function: cholesky (<M>)
 -- Function: cholesky (<M>, <field>)
     Return the Cholesky factorization of the matrix selfadjoint (or
     hermitian) matrix <M>. The second argument defaults to
     'generalring.' For a description of the possible values for
     <field>, see `lu_factor'.


 -- Function: ctranspose (<M>)
     Return the complex conjugate transpose of the matrix <M>. The
     function `ctranspose' uses `matrix_element_transpose' to transpose
     each matrix element.


 -- Function: diag_matrix (<d_1>, <d_2>,...,<d_n>)
     Return a diagonal matrix with diagonal entries <d_1>,
     <d_2>,...,<d_n>.  When the diagonal entries are matrices, the zero
     entries of the returned matrix are zero matrices of the
     appropriate size; for example:
          (%i1) diag_matrix(diag_matrix(1,2),diag_matrix(3,4));

                                      [ [ 1  0 ]  [ 0  0 ] ]
                                      [ [      ]  [      ] ]
                                      [ [ 0  2 ]  [ 0  0 ] ]
          (%o1)                       [                    ]
                                      [ [ 0  0 ]  [ 3  0 ] ]
                                      [ [      ]  [      ] ]
                                      [ [ 0  0 ]  [ 0  4 ] ]
          (%i2) diag_matrix(p,q);

                                             [ p  0 ]
          (%o2)                              [      ]
                                             [ 0  q ]

 -- Function: dotproduct (<u>, <v>)
     Return the dotproduct of vectors <u> and <v>.  This is the same as
     `conjugate (transpose (<u>)) . <v>'.  The arguments <u> and <v>
     must be column vectors.


 -- Function: eigens_by_jacobi (<A>)
 -- Function: eigens_by_jacobi (<A>, <field_type>)
     Computes the eigenvalues and eigenvectors of <A> by the method of
     Jacobi rotations.  <A> must be a symmetric matrix (but it need not
     be positive definite nor positive semidefinite).  <field_type>
     indicates the computational field, either `floatfield' or
     `bigfloatfield'.  If <field_type> is not specified, it defaults to
     `floatfield'.

     The elements of <A> must be numbers or expressions which evaluate
     to numbers via `float' or `bfloat' (depending on <field_type>).

     Examples:

          (%i1) S: matrix([1/sqrt(2), 1/sqrt(2)],[-1/sqrt(2), 1/sqrt(2)]);
                               [     1         1    ]
                               [  -------   ------- ]
                               [  sqrt(2)   sqrt(2) ]
          (%o1)                [                    ]
                               [      1        1    ]
                               [ - -------  ------- ]
                               [   sqrt(2)  sqrt(2) ]
          (%i2) L : matrix ([sqrt(3), 0], [0, sqrt(5)]);
                                [ sqrt(3)     0    ]
          (%o2)                 [                  ]
                                [    0     sqrt(5) ]
          (%i3) M : S . L . transpose (S);
                      [ sqrt(5)   sqrt(3)  sqrt(5)   sqrt(3) ]
                      [ ------- + -------  ------- - ------- ]
                      [    2         2        2         2    ]
          (%o3)       [                                      ]
                      [ sqrt(5)   sqrt(3)  sqrt(5)   sqrt(3) ]
                      [ ------- - -------  ------- + ------- ]
                      [    2         2        2         2    ]
          (%i4) eigens_by_jacobi (M);
          The largest percent change was 0.1454972243679
          The largest percent change was 0.0
          number of sweeps: 2
          number of rotations: 1
          (%o4) [[1.732050807568877, 2.23606797749979],
                                  [  0.70710678118655   0.70710678118655 ]
                                  [                                      ]]
                                  [ - 0.70710678118655  0.70710678118655 ]
          (%i5) float ([[sqrt(3), sqrt(5)], S]);
          (%o5) [[1.732050807568877, 2.23606797749979],
                                  [  0.70710678118655   0.70710678118655 ]
                                  [                                      ]]
                                  [ - 0.70710678118655  0.70710678118655 ]
          (%i6) eigens_by_jacobi (M, bigfloatfield);
          The largest percent change was 1.454972243679028b-1
          The largest percent change was 0.0b0
          number of sweeps: 2
          number of rotations: 1
          (%o6) [[1.732050807568877b0, 2.23606797749979b0],
                          [  7.071067811865475b-1   7.071067811865475b-1 ]
                          [                                              ]]
                          [ - 7.071067811865475b-1  7.071067811865475b-1 ]


 -- Function: get_lu_factors (<x>)
     When `<x> = lu_factor (<A>)', then `get_lu_factors' returns a list
     of the form `[P, L, U]', where <P> is a permutation matrix, <L> is
     lower triangular with ones on the diagonal, and <U> is upper
     triangular, and `<A> = <P> <L> <U>'.


 -- Function: hankel (<col>)
 -- Function: hankel (<col>, <row>)
     Return a Hankel matrix <H>. The first column of <H> is <col>;
     except for the first entry, the last row of <H> is <row>. The
     default for <row> is the zero vector with the same length as <col>.


 -- Function: hessian (<f>,<vars>)
     Return the hessian matrix of <f> with respect to the variables in
     the list <vars>.  The <i>,<j> entry of the hessian matrix is
     <diff(f vars[i],1,vars[j],1)>.


 -- Function: hilbert_matrix (<n>)
     Return the <n> by <n> Hilbert matrix. When <n> isn't a positive
     integer, signal an error.


 -- Function: identfor (<M>)
 -- Function: identfor (<M>, <fld>)
     Return an identity matrix that has the same shape as the matrix
     <M>.  The diagonal entries of the identity matrix are the
     multiplicative identity of the field <fld>; the default for <fld>
     is <generalring>.

     The first argument <M> should be a square matrix or a non-matrix.
     When <M> is a matrix, each entry of <M> can be a square matrix -
     thus <M> can be a blocked Maxima matrix. The matrix can be blocked
     to any (finite) depth.

     See also `zerofor'


 -- Function: invert_by_lu (<M>, <(rng generalring)>)
     Invert a matrix <M> by using the LU factorization.  The LU
     factorization is done using the ring <rng>.


 -- Function: kronecker_product (<A>, <B>)
     Return the Kronecker product of the matrices <A> and <B>.


 -- Function: listp (<e>, <p>)
 -- Function: listp (<e>)
     Given an optional argument <p>, return `true' if <e> is a Maxima
     list and <p> evaluates to `true' for every list element.  When
     `listp' is not given the optional argument, return `true' if <e> is
     a Maxima list.  In all other cases, return `false'.


 -- Function: locate_matrix_entry (<M>, <r_1>, <c_1>, <r_2>, <c_2>,
          <f>, <rel>)
     The first argument must be a matrix; the arguments <r_1> through
     <c_2> determine a sub-matrix of <M> that consists of rows <r_1>
     through <r_2> and columns <c_1> through <c_2>.

     Find a entry in the sub-matrix <M> that satisfies some property.
     Three cases:

     (1) `<rel> = 'bool' and <f> a predicate:

     Scan the sub-matrix from left to right then top to bottom, and
     return the index of the first entry that satisfies the predicate
     <f>. If no matrix entry satisfies <f>, return `false'.

     (2) `<rel> = 'max' and <f> real-valued:

     Scan the sub-matrix looking for an entry that maximizes <f>.
     Return the index of a maximizing entry.

     (3) `<rel> = 'min' and <f> real-valued:

     Scan the sub-matrix looking for an entry that minimizes <f>.
     Return the index of a minimizing entry.


 -- Function: lu_backsub (<M>, <b>)
     When `<M> = lu_factor (<A>, <field>)', then `lu_backsub (<M>,
     <b>)' solves the linear system `<A> <x> = <b>'.


 -- Function: lu_factor (<M>, <field>)
     Return a list of the form `[<LU>, <perm>, <fld>]', or `[<LU>,
     <perm>, <fld>, <lower-cnd> <upper-cnd>]', where

     (1) The matrix <LU> contains the factorization of <M> in a packed
     form. Packed       form means three things: First, the rows of
     <LU> are permuted according to the       list <perm>.  If, for
     example, <perm> is the list `[3,2,1]', the actual first row
     of the <LU> factorization is the third row of the matrix <LU>.
     Second,       the lower triangular factor of m is the lower
     triangular part of <LU> with the       diagonal entries replaced
     by all ones. Third, the upper triangular factor of       <M> is
     the upper triangular part of <LU>.

     (2) When the field is either `floatfield' or `complexfield',
     the numbers <lower-cnd> and <upper-cnd> are lower and upper bounds
     for the       infinity norm condition number of <M>.  For all
     fields, the condition number       might not be estimated; for
     such fields, `lu_factor' returns a two item list.        Both the
     lower and upper bounds can differ from their true values by
     arbitrarily large factors. (See also `mat_cond'.)

     The argument <M> must be a square matrix.

     The optional argument <fld> must be a symbol that determines a
     ring or field. The pre-defined   fields and rings are:

     (a) `generalring' - the ring of Maxima expressions,     (b)
     `floatfield' -  the field of floating point numbers of the type
     double,     (c) `complexfield' -  the field of complex floating
     point numbers of the         type double,     (d) `crering'  - the
     ring of Maxima CRE expressions,     (e) `rationalfield' - the
     field of rational numbers,     (f) `runningerror' - track the all
     floating point rounding errors,     (g) `noncommutingring' - the
     ring of Maxima expressions where multiplication is the
     non-commutative dot operator.

     When the field is `floatfield', `complexfield', or `runningerror',
     the algorithm uses partial pivoting; for all other fields, rows
     are switched only when needed to avoid a zero pivot.

     Floating point addition arithmetic isn't associative, so the
     meaning of 'field' differs from the mathematical definition.

     A member of the field `runningerror' is a two member Maxima list
     of the form `[x,n]',where <x> is a floating point number and `n'
     is an integer. The relative difference between the 'true' value of
     `x' and `x' is approximately bounded by the machine epsilon times
     `n'. The running error bound drops some terms that of the order
     the square of the machine epsilon.

     There is no user-interface for defining a new field. A user that is
     familiar with Common Lisp should be able to define a new field.
     To do this, a user must define functions for the arithmetic
     operations and functions for converting from the field
     representation to Maxima and back. Additionally, for ordered
     fields (where partial pivoting will be used), a user must define
     functions for the magnitude and for comparing field members.
     After that all that remains is to define a Common Lisp structure
     `mring'.  The file `mring' has many examples.

     To compute the factorization, the first task is to convert each
     matrix entry to a member of the indicated field. When conversion
     isn't possible, the factorization halts with an error message.
     Members of the field needn't be Maxima expressions.  Members of the
     `complexfield', for example, are Common Lisp complex numbers. Thus
     after computing the factorization, the matrix entries must be
     converted to Maxima expressions.

     See also  `get_lu_factors'.

     Examples:

          (%i1) w[i,j] := random (1.0) + %i * random (1.0);
          (%o1)          w     := random(1.) + %i random(1.)
                          i, j
          (%i2) showtime : true$
          Evaluation took 0.00 seconds (0.00 elapsed)
          (%i3) M : genmatrix (w, 100, 100)$
          Evaluation took 7.40 seconds (8.23 elapsed)
          (%i4) lu_factor (M, complexfield)$
          Evaluation took 28.71 seconds (35.00 elapsed)
          (%i5) lu_factor (M, generalring)$
          Evaluation took 109.24 seconds (152.10 elapsed)
          (%i6) showtime : false$

          (%i7) M : matrix ([1 - z, 3], [3, 8 - z]);
                                  [ 1 - z    3   ]
          (%o7)                   [              ]
                                  [   3    8 - z ]
          (%i8) lu_factor (M, generalring);
                    [ 1 - z         3        ]
                    [                        ]
          (%o8)    [[   3            9       ], [1, 2], generalring]
                    [ -----  - z - ----- + 8 ]
                    [ 1 - z        1 - z     ]
          (%i9) get_lu_factors (%);
                            [   1    0 ]  [ 1 - z         3        ]
                  [ 1  0 ]  [          ]  [                        ]
          (%o9)  [[      ], [   3      ], [                9       ]]
                  [ 0  1 ]  [ -----  1 ]  [   0    - z - ----- + 8 ]
                            [ 1 - z    ]  [              1 - z     ]
          (%i10) %[1] . %[2] . %[3];
                                  [ 1 - z    3   ]
          (%o10)                  [              ]
                                  [   3    8 - z ]


 -- Function: mat_cond (<M>, 1)
 -- Function: mat_cond (<M>, inf)
     Return the <p>-norm matrix condition number of the matrix <m>. The
     allowed values for <p> are 1 and <inf>.  This function uses the LU
     factorization to invert the matrix <m>. Thus the running time for
     `mat_cond' is proportional to the cube of the matrix size;
     `lu_factor' determines lower and upper bounds for the infinity
     norm condition number in time proportional to the square of the
     matrix size.


 -- Function: mat_norm (<M>, 1)
 -- Function: mat_norm (<M>, inf)
 -- Function: mat_norm (<M>, frobenius)
     Return the matrix <p>-norm of the matrix <M>.  The allowed values
     for <p> are 1, `inf', and `frobenius' (the Frobenius matrix norm).
     The matrix <M> should be an unblocked matrix.


 -- Function: matrixp (<e>, <p>)
 -- Function: matrixp (<e>)
     Given an optional argument <p>, return `true' if <e> is a matrix
     and <p> evaluates to `true' for every matrix element.  When
     `matrixp' is not given an optional argument, return `true' if `e'
     is a matrix.  In all other cases, return `false'.

     See also `blockmatrixp'


 -- Function: matrix_size (<M>)
     Return a two member list that gives the number of rows and
     columns, respectively of the matrix <M>.


 -- Function: mat_fullunblocker (<M>)
     If <M> is a block matrix, unblock the matrix to all levels. If <M>
     is a matrix, return <M>; otherwise, signal an error.


 -- Function: mat_trace (<M>)
     Return the trace of the matrix <M>. If <M> isn't a matrix, return a
     noun form. When <M> is a block matrix, `mat_trace(M)' returns the
     same value as does `mat_trace(mat_unblocker(m))'.


 -- Function: mat_unblocker (<M>)
     If <M> is a block matrix, unblock <M> one level. If <M> is a
     matrix, `mat_unblocker (M)' returns <M>; otherwise, signal an
     error.

     Thus if each entry of <M> is matrix, `mat_unblocker (M)' returns an
     unblocked matrix, but if each entry of <M> is a block matrix,
     `mat_unblocker (M)' returns a block matrix with one less level of
     blocking.

     If you use block matrices, most likely you'll want to set
     `matrix_element_mult' to `"."' and `matrix_element_transpose' to
     `'transpose'. See also `mat_fullunblocker'.

     Example:

          (%i1) A : matrix ([1, 2], [3, 4]);
                                      [ 1  2 ]
          (%o1)                       [      ]
                                      [ 3  4 ]
          (%i2) B : matrix ([7, 8], [9, 10]);
                                      [ 7  8  ]
          (%o2)                       [       ]
                                      [ 9  10 ]
          (%i3) matrix ([A, B]);
                               [ [ 1  2 ]  [ 7  8  ] ]
          (%o3)                [ [      ]  [       ] ]
                               [ [ 3  4 ]  [ 9  10 ] ]
          (%i4) mat_unblocker (%);
                                   [ 1  2  7  8  ]
          (%o4)                    [             ]
                                   [ 3  4  9  10 ]


 -- Function: nonnegintegerp (<n>)
     Return `true' if and only if `<n> >= 0' and <n> is an integer.


 -- Function: nullspace (<M>)
     If <M> is a matrix, return `span (v_1, ..., v_n)', where the set
     `{v_1, ..., v_n}' is a basis for the nullspace of <M>.  The span
     of the empty set is `{0}'.  Thus, when the nullspace has only one
     member, return `span ()'.


 -- Function: nullity (<M>)
     If <M> is a matrix, return the dimension of the nullspace of <M>.


 -- Function: orthogonal_complement (<v_1>, ..., <v_n>)
     Return `span (u_1, ..., u_m)', where the set `{u_1, ..., u_m}' is a
     basis for the orthogonal complement of the set `(v_1, ..., v_n)'.

     Each vector <v_1> through <v_n> must be a column vector.


 -- Function: polynomialp (<p>, <L>, <coeffp>, <exponp>)
 -- Function: polynomialp (<p>, <L>, <coeffp>)
 -- Function: polynomialp (<p>, <L>)
     Return `true' if <p> is a polynomial in the variables in the list
     <L>, The predicate <coeffp> must evaluate to `true' for each
     coefficient, and the predicate <exponp> must evaluate to `true'
     for all exponents of the variables in <L>. If you want to use a
     non-default value for <exponp>, you must supply <coeffp> with a
     value even if you want to use the default for <coeffp>.

     `polynomialp (<p>, <L>, <coeffp>)' is equivalent to `polynomialp
     (<p>, <L>, <coeffp>, 'nonnegintegerp)'.

     `polynomialp (<p>, <L>)' is equivalent to `polynomialp (<p>, L<,>
     'constantp, 'nonnegintegerp)'.

     The polynomial needn't be expanded:

          (%i1) polynomialp ((x + 1)*(x + 2), [x]);
          (%o1)                         true
          (%i2) polynomialp ((x + 1)*(x + 2)^a, [x]);
          (%o2)                         false

     An example using non-default values for coeffp and exponp:

          (%i1) polynomialp ((x + 1)*(x + 2)^(3/2), [x], numberp, numberp);
          (%o1)                         true
          (%i2) polynomialp ((x^(1/2) + 1)*(x + 2)^(3/2), [x], numberp,
                                                                  numberp);
          (%o2)                         true

     Polynomials with two variables:

          (%i1) polynomialp (x^2 + 5*x*y + y^2, [x]);
          (%o1)                         false
          (%i2) polynomialp (x^2 + 5*x*y + y^2, [x, y]);
          (%o2)                         true


 -- Function: polytocompanion (<p>, <x>)
     If <p> is a polynomial in <x>, return the companion matrix of <p>.
     For a monic polynomial <p> of degree <n>, we have `<p> = (-1)^<n>
     charpoly (polytocompanion (<p>, <x>))'.

     When <p> isn't a polynomial in <x>, signal an error.


 -- Function: ptriangularize (<M>, <v>)
     If <M> is a matrix with each entry a polynomial in <v>, return a
     matrix <M2> such that

     (1) <M2> is upper triangular,

     (2) `<M2> = <E_n> ... <E_1> <M>', where <E_1> through <E_n> are
     elementary matrices whose entries are polynomials in <v>,

     (3) `|det (<M>)| = |det (<M2>)|',

     Note: This function doesn't check that every entry is a polynomial
     in <v>.


 -- Function: rowop (<M>, <i>, <j>, <theta>)
     If <M> is a matrix, return the matrix that results from doing the
     row operation `R_i <- R_i - theta * R_j'. If <M> doesn't have a row
     <i> or <j>, signal an error.


 -- Function: rank (<M>)
     Return the rank of that matrix <M>. The rank is the dimension of
     the column space. Example:
          (%i1) rank(matrix([1,2],[2,4]));
          (%o1)                                  1
          (%i2) rank(matrix([1,b],[c,d]));
          Proviso:  {d - b c # 0}
          (%o2)                                  2


 -- Function: rowswap (<M>, <i>, <j>)
     If <M> is a matrix, swap rows <i> and <j>. If <M> doesn't have a
     row <i> or <j>, signal an error.


 -- Function: toeplitz (<col>)
 -- Function: toeplitz (<col>, <row>)
     Return a Toeplitz matrix <T>. The first first column of <T> is
     <col>; except for the first entry, the first row of <T> is <row>.
     The default for <row> is complex conjugate of <col>. Example:
          (%i1)  toeplitz([1,2,3],[x,y,z]);

                                            [ 1  y  z ]
                                            [         ]
          (%o1)                             [ 2  1  y ]
                                            [         ]
                                            [ 3  2  1 ]
          (%i2)  toeplitz([1,1+%i]);

                                        [   1     1 - %I ]
          (%o2)                         [                ]
                                        [ %I + 1    1    ]


 -- Function: vandermonde_matrix ([<x_1>, ..., <x_n>])
     Return a <n> by <n> matrix whose <i>-th row is `[1, <x_i>,
     <x_i>^2, ... <x_i>^(<n>-1)]'.


 -- Function: zerofor (<M>)
 -- Function: zerofor (<M>, <fld>)
     Return a zero  matrix that has the same shape as the matrix <M>.
     Every entry of the zero matrix is the additive identity of the
     field <fld>; the default for <fld> is <generalring>.

     The first argument <M> should be a square matrix or a non-matrix.
     When <M> is a matrix, each entry of <M> can be a square matrix -
     thus <M> can be a blocked Maxima matrix. The matrix can be blocked
     to any (finite) depth.

     See also `identfor'


 -- Function: zeromatrixp (<M>)
     If <M> is not a block matrix, return `true' if `is (equal (<e>,
     0))' is true for each element <e> of the matrix <M>.  If <M> is a
     block matrix, return `true' if `zeromatrixp' evaluates to `true'
     for each element of <e>.



File: maxima.info,  Node: lsquares,  Next: makeOrders,  Prev: linearalgebra,  Up: Top

60 lsquares
***********

* Menu:

* Functions and Variables for lsquares::


File: maxima.info,  Node: Functions and Variables for lsquares,  Prev: lsquares,  Up: lsquares

60.1 Functions and Variables for lsquares
=========================================

 -- Function: lsquares_estimates (<D>, <x>, <e>, <a>)
 -- Function: lsquares_estimates (<D>, <x>, <e>, <a>, initial = <L>,
          tol = <t>)
     Estimate parameters <a> to best fit the equation <e> in the
     variables <x> and <a> to the data <D>, as determined by the method
     of least squares.  `lsquares_estimates' first seeks an exact
     solution, and if that fails, then seeks an approximate solution.

     The return value is a list of lists of equations of the form `[a =
     ..., b = ..., c = ...]'.  Each element of the list is a distinct,
     equivalent minimum of the mean square error.

     The data <D> must be a matrix.  Each row is one datum (which may
     be called a `record' or `case' in some contexts), and each column
     contains the values of one variable across all data.  The list of
     variables <x> gives a name for each column of <D>, even the
     columns which do not enter the analysis.  The list of parameters
     <a> gives the names of the parameters for which estimates are
     sought.  The equation <e> is an expression or equation in the
     variables <x> and <a>; if <e> is not an equation, it is treated
     the same as `<e> = 0'.

     Additional arguments to `lsquares_estimates' are specified as
     equations and passed on verbatim to the function `lbfgs' which is
     called to find estimates by a numerical method when an exact
     result is not found.

     If some exact solution can be found (via `solve'), the data <D>
     may contain non-numeric values.  However, if no exact solution is
     found, each element of <D> must have a numeric value.  This
     includes numeric constants such as `%pi' and `%e' as well as
     literal numbers (integers, rationals, ordinary floats, and
     bigfloats).  Numerical calculations are carried out with ordinary
     floating-point arithmetic, so all other kinds of numbers are
     converted to ordinary floats for calculations.

     `load(lsquares)' loads this function.

     See also `lsquares_estimates_exact',
     `lsquares_estimates_approximate', `lsquares_mse',
     `lsquares_residuals', and `lsquares_residual_mse'.

     Examples:

     A problem for which an exact solution is found.

          (%i1) load (lsquares)$
          (%i2) M : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]);
                                     [ 1  1  1 ]
                                     [         ]
                                     [ 3       ]
                                     [ -  1  2 ]
                                     [ 2       ]
                                     [         ]
          (%o2)                      [ 9       ]
                                     [ -  2  1 ]
                                     [ 4       ]
                                     [         ]
                                     [ 3  2  2 ]
                                     [         ]
                                     [ 2  2  1 ]
          (%i3) lsquares_estimates (M, [z, x, y], (z + D)^2 = A*x + B*y + C, [A, B, C, D]);
                            59        27      10921        107
          (%o3)     [[A = - --, B = - --, C = -----, D = - ---]]
                            16        16      1024         32

     A problem for which no exact solution is found, so
     `lsquares_estimates' resorts to numerical approximation.

          (%i1) load (lsquares)$
          (%i2) M : matrix ([1, 1], [2, 7/4], [3, 11/4], [4, 13/4]);
                                      [ 1  1  ]
                                      [       ]
                                      [    7  ]
                                      [ 2  -  ]
                                      [    4  ]
                                      [       ]
          (%o2)                       [    11 ]
                                      [ 3  -- ]
                                      [    4  ]
                                      [       ]
                                      [    13 ]
                                      [ 4  -- ]
                                      [    4  ]
          (%i3) lsquares_estimates (M, [x, y], y = a*x^b + c, [a, b, c], initial = [3, 3, 3], iprint = [-1, 0]);
          (%o3) [[a = 1.387365874920637, b = .7110956639593767,
                                                  c = - .4142705622439105]]


 -- Function: lsquares_estimates_exact (<MSE>, <a>)
     Estimate parameters <a> to minimize the mean square error <MSE>,
     by constructing a system of equations and attempting to solve them
     symbolically via `solve'.  The mean square error is an expression
     in the parameters <a>, such as that returned by `lsquares_mse'.

     The return value is a list of lists of equations of the form `[a =
     ..., b = ..., c = ...]'.  The return value may contain zero, one,
     or two or more elements.  If two or more elements are returned,
     each represents a distinct, equivalent minimum of the mean square
     error.

     See also `lsquares_estimates', `lsquares_estimates_approximate',
     `lsquares_mse', `lsquares_residuals', and `lsquares_residual_mse'.

     Example:

          (%i1) load (lsquares)$
          (%i2) M : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]);
                                     [ 1  1  1 ]
                                     [         ]
                                     [ 3       ]
                                     [ -  1  2 ]
                                     [ 2       ]
                                     [         ]
          (%o2)                      [ 9       ]
                                     [ -  2  1 ]
                                     [ 4       ]
                                     [         ]
                                     [ 3  2  2 ]
                                     [         ]
                                     [ 2  2  1 ]
          (%i3) mse : lsquares_mse (M, [z, x, y], (z + D)^2 = A*x + B*y + C);
                     5
                    ====
                    \                 2                         2
                     >    ((D + M    )  - C - M     B - M     A)
                    /            i, 1          i, 3      i, 2
                    ====
                    i = 1
          (%o3)     ---------------------------------------------
                                          5
          (%i4) lsquares_estimates_exact (mse, [A, B, C, D]);
                            59        27      10921        107
          (%o4)     [[A = - --, B = - --, C = -----, D = - ---]]
                            16        16      1024         32


 -- Function: lsquares_estimates_approximate (<MSE>, <a>, initial =
          <L>, tol = <t>)
     Estimate parameters <a> to minimize the mean square error <MSE>,
     via the numerical minimization function `lbfgs'.  The mean square
     error is an expression in the parameters <a>, such as that
     returned by `lsquares_mse'.

     The solution returned by `lsquares_estimates_approximate' is a
     local (perhaps global) minimum of the mean square error.  For
     consistency with `lsquares_estimates_exact', the return value is a
     nested list which contains one element, namely a list of equations
     of the form `[a = ..., b = ..., c = ...]'.

     Additional arguments to `lsquares_estimates_approximate' are
     specified as equations and passed on verbatim to the function
     `lbfgs'.

     <MSE> must evaluate to a number when the parameters are assigned
     numeric values.  This requires that the data from which <MSE> was
     constructed comprise only numeric constants such as `%pi' and `%e'
     and literal numbers (integers, rationals, ordinary floats, and
     bigfloats).  Numerical calculations are carried out with ordinary
     floating-point arithmetic, so all other kinds of numbers are
     converted to ordinary floats for calculations.

     `load(lsquares)' loads this function.

     See also `lsquares_estimates', `lsquares_estimates_exact',
     `lsquares_mse', `lsquares_residuals', and `lsquares_residual_mse'.

     Example:

          (%i1) load (lsquares)$
          (%i2) M : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]);
                                     [ 1  1  1 ]
                                     [         ]
                                     [ 3       ]
                                     [ -  1  2 ]
                                     [ 2       ]
                                     [         ]
          (%o2)                      [ 9       ]
                                     [ -  2  1 ]
                                     [ 4       ]
                                     [         ]
                                     [ 3  2  2 ]
                                     [         ]
                                     [ 2  2  1 ]
          (%i3) mse : lsquares_mse (M, [z, x, y], (z + D)^2 = A*x + B*y + C);
                     5
                    ====
                    \                 2                         2
                     >    ((D + M    )  - C - M     B - M     A)
                    /            i, 1          i, 3      i, 2
                    ====
                    i = 1
          (%o3)     ---------------------------------------------
                                          5
          (%i4) lsquares_estimates_approximate (mse, [A, B, C, D], iprint = [-1, 0]);
          (%o4) [[A = - 3.67850494740174, B = - 1.683070351177813,
                           C = 10.63469950148635, D = - 3.340357993175206]]


 -- Function: lsquares_mse (<D>, <x>, <e>)
     Returns the mean square error (MSE), a summation expression, for
     the equation <e> in the variables <x>, with data <D>.

     The MSE is defined as:

                              n
                             ====
                             \                        2
                              >    (lhs(e ) - rhs(e ))
                             /           i         i
                             ====
                             i = 1
                             --------------------------
                                         n

     where <n> is the number of data and `<e>[i]' is the equation <e>
     evaluated with the variables in <x> assigned values from the
     `i'-th datum, `<D>[i]'.

     `load(lsquares)' loads this function.

     Example:

          (%i1) load (lsquares)$
          (%i2) M : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]);
                                     [ 1  1  1 ]
                                     [         ]
                                     [ 3       ]
                                     [ -  1  2 ]
                                     [ 2       ]
                                     [         ]
          (%o2)                      [ 9       ]
                                     [ -  2  1 ]
                                     [ 4       ]
                                     [         ]
                                     [ 3  2  2 ]
                                     [         ]
                                     [ 2  2  1 ]
          (%i3) mse : lsquares_mse (M, [z, x, y], (z + D)^2 = A*x + B*y + C);
                     5
                    ====
                    \                 2                         2
                     >    ((D + M    )  - C - M     B - M     A)
                    /            i, 1          i, 3      i, 2
                    ====
                    i = 1
          (%o3)     ---------------------------------------------
                                          5
          (%i4) diff (mse, D);
                   5
                  ====
                  \                             2
                4  >    (D + M    ) ((D + M    )  - C - M     B - M     A)
                  /           i, 1         i, 1          i, 3      i, 2
                  ====
                  i = 1
          (%o4) ----------------------------------------------------------
                                            5
          (%i5) ''mse, nouns;
                         2                 2         9 2               2
          (%o5) (((D + 3)  - C - 2 B - 2 A)  + ((D + -)  - C - B - 2 A)
                                                     4
                     2               2         3 2               2
           + ((D + 2)  - C - B - 2 A)  + ((D + -)  - C - 2 B - A)
                                               2
                     2             2
           + ((D + 1)  - C - B - A) )/5


 -- Function: lsquares_residuals (<D>, <x>, <e>, <a>)
     Returns the residuals for the equation <e> with specified
     parameters <a> and data <D>.

     <D> is a matrix, <x> is a list of variables, <e> is an equation or
     general expression; if not an equation, <e> is treated as if it
     were `<e> = 0'.  <a> is a list of equations which specify values
     for any free parameters in <e> aside from <x>.

     The residuals are defined as:

                                  lhs(e ) - rhs(e )
                                       i         i

     where `<e>[i]' is the equation <e> evaluated with the variables in
     <x> assigned values from the `i'-th datum, `<D>[i]', and assigning
     any remaining free variables from <a>.

     `load(lsquares)' loads this function.

     Example:

          (%i1) load (lsquares)$
          (%i2) M : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]);
                                     [ 1  1  1 ]
                                     [         ]
                                     [ 3       ]
                                     [ -  1  2 ]
                                     [ 2       ]
                                     [         ]
          (%o2)                      [ 9       ]
                                     [ -  2  1 ]
                                     [ 4       ]
                                     [         ]
                                     [ 3  2  2 ]
                                     [         ]
                                     [ 2  2  1 ]
          (%i3) a : lsquares_estimates (M, [z, x, y], (z + D)^2 = A*x + B*y + C, [A, B, C, D]);
                            59        27      10921        107
          (%o3)     [[A = - --, B = - --, C = -----, D = - ---]]
                            16        16      1024         32
          (%i4) lsquares_residuals (M, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a));
                               13    13    13  13  13
          (%o4)               [--, - --, - --, --, --]
                               64    64    32  64  64


 -- Function: lsquares_residual_mse (<D>, <x>, <e>, <a>)
     Returns the residual mean square error (MSE) for the equation <e>
     with specified parameters <a> and data <D>.

     The residual MSE is defined as:

                              n
                             ====
                             \                        2
                              >    (lhs(e ) - rhs(e ))
                             /           i         i
                             ====
                             i = 1
                             --------------------------
                                         n

     where `<e>[i]' is the equation <e> evaluated with the variables in
     <x> assigned values from the `i'-th datum, `<D>[i]', and assigning
     any remaining free variables from <a>.

     `load(lsquares)' loads this function.

     Example:

          (%i1) load (lsquares)$
          (%i2) M : matrix ([1, 1, 1], [3/2, 1, 2], [9/4, 2, 1], [3, 2, 2], [2, 2, 1]);
                                     [ 1  1  1 ]
                                     [         ]
                                     [ 3       ]
                                     [ -  1  2 ]
                                     [ 2       ]
                                     [         ]
          (%o2)                      [ 9       ]
                                     [ -  2  1 ]
                                     [ 4       ]
                                     [         ]
                                     [ 3  2  2 ]
                                     [         ]
                                     [ 2  2  1 ]
          (%i3) a : lsquares_estimates (M, [z, x, y], (z + D)^2 = A*x + B*y + C, [A, B, C, D]);
                            59        27      10921        107
          (%o3)     [[A = - --, B = - --, C = -----, D = - ---]]
                            16        16      1024         32
          (%i4) lsquares_residual_mse (M, [z, x, y], (z + D)^2 = A*x + B*y + C, first (a));
                                        169
          (%o4)                         ----
                                        2560


 -- Function: plsquares (<Mat>,<VarList>,<depvars>)
 -- Function: plsquares (<Mat>,<VarList>,<depvars>,<maxexpon>)
 -- Function: plsquares
          (<Mat>,<VarList>,<depvars>,<maxexpon>,<maxdegree>)
     Multivariable polynomial adjustment of a data table by the "least
     squares" method. <Mat> is a matrix containing the data, <VarList>
     is a list of variable names (one for each Mat column, but use "-"
     instead of varnames to ignore Mat columns), <depvars> is the name
     of a dependent variable or a list with one or more names of
     dependent variables (which names should be in <VarList>),
     <maxexpon> is the optional maximum exponent for each independent
     variable (1 by default), and <maxdegree> is the optional maximum
     polynomial degree (<maxexpon> by default); note that the sum of
     exponents of each term must be equal or smaller than <maxdegree>,
     and if `maxdgree = 0' then no limit is applied.

     If <depvars> is the name of a dependent variable (not in a list),
     `plsquares' returns the adjusted polynomial. If <depvars> is a
     list of one or more dependent variables, `plsquares' returns a
     list with the adjusted polynomial(s). The Coefficients of
     Determination  are displayed in order to inform about the goodness
     of fit, which ranges from 0 (no correlation) to 1 (exact
     correlation). These values are also stored in the global variable
     <DETCOEF> (a list if <depvars> is a list).

     A simple example of multivariable linear adjustment:
          (%i1) load("plsquares")$

          (%i2) plsquares(matrix([1,2,0],[3,5,4],[4,7,9],[5,8,10]),
                          [x,y,z],z);
               Determination Coefficient for z = .9897039897039897
                                 11 y - 9 x - 14
          (%o2)              z = ---------------
                                        3

     The same example without degree restrictions:
          (%i3) plsquares(matrix([1,2,0],[3,5,4],[4,7,9],[5,8,10]),
                          [x,y,z],z,1,0);
               Determination Coefficient for z = 1.0
                              x y + 23 y - 29 x - 19
          (%o3)           z = ----------------------
                                        6

     How many diagonals does a N-sides polygon have? What polynomial
     degree should be used?
          (%i4) plsquares(matrix([3,0],[4,2],[5,5],[6,9],[7,14],[8,20]),
                          [N,diagonals],diagonals,5);
               Determination Coefficient for diagonals = 1.0
                                          2
                                         N  - 3 N
          (%o4)              diagonals = --------
                                            2
          (%i5) ev(%, N=9);   /* Testing for a 9 sides polygon */
          (%o5)                 diagonals = 27

     How many ways do we have to put two queens without they are
     threatened into a n x n chessboard?
          (%i6) plsquares(matrix([0,0],[1,0],[2,0],[3,8],[4,44]),
                          [n,positions],[positions],4);
               Determination Coefficient for [positions] = [1.0]
                                   4       3      2
                                3 n  - 10 n  + 9 n  - 2 n
          (%o6)    [positions = -------------------------]
                                            6
          (%i7) ev(%[1], n=8); /* Testing for a (8 x 8) chessboard */
          (%o7)                positions = 1288

     An example with six dependent variables:
          (%i8) mtrx:matrix([0,0,0,0,0,1,1,1],[0,1,0,1,1,1,0,0],
                            [1,0,0,1,1,1,0,0],[1,1,1,1,0,0,0,1])$
          (%i8) plsquares(mtrx,[a,b,_And,_Or,_Xor,_Nand,_Nor,_Nxor],
                               [_And,_Or,_Xor,_Nand,_Nor,_Nxor],1,0);
                Determination Coefficient for
          [_And, _Or, _Xor, _Nand, _Nor, _Nxor] =
          [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
          (%o2) [_And = a b, _Or = - a b + b + a,
          _Xor = - 2 a b + b + a, _Nand = 1 - a b,
          _Nor = a b - b - a + 1, _Nxor = 2 a b - b - a + 1]

     To use this function write first `load("lsquares")'.


File: maxima.info,  Node: makeOrders,  Next: mnewton,  Prev: lsquares,  Up: Top

61 makeOrders
*************

* Menu:

* Functions and Variables for makeOrders::


File: maxima.info,  Node: Functions and Variables for makeOrders,  Prev: makeOrders,  Up: makeOrders

61.1 Functions and Variables for makeOrders
===========================================

 -- Function: makeOrders (<indvarlist>,<orderlist>)
     Returns a list of all powers for a polynomial up to and including
     the arguments.

          (%i1) load("makeOrders")$

          (%i2) makeOrders([a,b],[2,3]);
          (%o2) [[0, 0], [0, 1], [0, 2], [0, 3], [1, 0], [1, 1],
                      [1, 2], [1, 3], [2, 0], [2, 1], [2, 2], [2, 3]]
          (%i3) expand((1+a+a^2)*(1+b+b^2+b^3));
                 2  3      3    3    2  2      2    2    2
          (%o3) a  b  + a b  + b  + a  b  + a b  + b  + a  b + a b
                                                            2
                                                     + b + a  + a + 1
     where `[0, 1]' is associated with the term b and `[2, 3]' with a^2
     b^3.

     To use this function write first `load("makeOrders")'.


File: maxima.info,  Node: mnewton,  Next: numericalio,  Prev: makeOrders,  Up: Top

62 mnewton
**********

* Menu:

* Functions and Variables for mnewton::


File: maxima.info,  Node: Functions and Variables for mnewton,  Prev: mnewton,  Up: mnewton

62.1 Functions and Variables for mnewton
========================================

 -- Option variable: newtonepsilon
     Default value: `10.0^(-fpprec/2)'

     Precision to determine when the `mnewton' function has converged
     towards the solution.

     See also `mnewton'.

 -- Option variable: newtonmaxiter
     Default value: `50'

     Maximum number of iterations to stop the `mnewton' function if it
     does not converge or if it converges too slowly.

     See also `mnewton'.

 -- Function: mnewton (<FuncList>,<VarList>,<GuessList>)
     Multiple nonlinear functions solution using the Newton method.
     <FuncList> is the list of functions to solve, <VarList> is the
     list of variable names, and <GuessList> is the list of initial
     approximations.

     The solution is returned in the same format that `solve()' returns.
     If the solution isn't found, `[]' is returned.

     This function is controlled by global variables `newtonepsilon'
     and `newtonmaxiter'.

          (%i1) load("mnewton")$

          (%i2) mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1],
                        [x1, x2], [5, 5]);
          (%o2) [[x1 = 3.756834008012769, x2 = 2.779849592817897]]
          (%i3) mnewton([2*a^a-5],[a],[1]);
          (%o3)             [[a = 1.70927556786144]]
          (%i4) mnewton([2*3^u-v/u-5, u+2^v-4], [u, v], [2, 2]);
          (%o4) [[u = 1.066618389595407, v = 1.552564766841786]]

     To use this function write first `load("mnewton")'. See also
     `newtonepsilon' and `newtonmaxiter'.


File: maxima.info,  Node: numericalio,  Next: opsubst,  Prev: mnewton,  Up: Top

63 numericalio
**************

* Menu:

* Introduction to numericalio::
* Functions and Variables for numericalio::


File: maxima.info,  Node: Introduction to numericalio,  Next: Functions and Variables for numericalio,  Prev: numericalio,  Up: numericalio

63.1 Introduction to numericalio
================================

`numericalio' is a collection of functions to read and write data files
and streams.

   It is assumed that each item to read or write is an atom: an
integer, float, bigfloat, string, or symbol, and not a rational or
complex number or any other kind of nonatomic expression.  The
`numericalio' functions may attempt to do something sensible faced with
nonatomic expressions, but the results are not specified here and
subject to change.

   Atoms in both input and output files have the same format as in
Maxima batch files or the interactive console.  In particular, strings
are enclosed in double quotes, backslash `\' prevents any special
interpretation of the next character, and the question mark `?' is
recognized at the beginning of a symbol to mean a Lisp symbol (as
opposed to a Maxima symbol).  No continuation character (to join broken
lines) is recognized.

   <separator_flag> tells which character separates elements.  It is an
optional argument for all read and write functions.

63.1.1 Separator flag values for input
--------------------------------------

For input, these values of <separator_flag> are recognized: `comma' for
comma separated values, `pipe' for values separated by the vertical bar
character `|', `semicolon' for values separated by semicolon `;', and
`space' for values separated by space or tab characters.  If the file
name ends in `.csv' and <separator_flag> is not specified, `comma' is
assumed.  If the file name ends in something other than `.csv' and
`separator_flag' is not specified, `space' is assumed.

   In input, multiple successive space and tab characters count as a
single separator.  However, multiple comma, pipe, or semicolon
characters are significant.  Successive comma, pipe, or semicolon
characters (with or without intervening spaces or tabs) are considered
to have `false' between the separators.  For example, `1234,,Foo' is
treated the same as `1234,false,Foo'.

63.1.2 Separator flag values for output
---------------------------------------

For output, `tab', for values separated by the tab character, is
recognized as a value of <separator_flag>, as well as `comma', `pipe',
`semicolon', and `space'.

   In output, `false' atoms are written as such; a list `[1234, false,
Foo]' is written `1234,false,Foo', and there is no attempt to collapse
the output to `1234,,Foo'.


File: maxima.info,  Node: Functions and Variables for numericalio,  Prev: Introduction to numericalio,  Up: numericalio

63.2 Functions and Variables for numericalio
============================================

 -- Function: read_matrix (<file_name>)
 -- Function: read_matrix (<file_name>, <separator_flag>)
 -- Function: read_matrix (<S>)
 -- Function: read_matrix (<S>, <separator_flag>)
     Reads the file <file_name> or file stream <S> and returns its
     entire content as a matrix.  If <separator_flag> is not specified,
     the file is assumed space-delimited.

     `read_matrix' infers the size of the matrix from the input data.
     Each line of the file becomes one row of the matrix.  If some
     lines have different lengths, `read_matrix' complains.


 -- Function: read_lisp_array (<file_name>, <A>)
 -- Function: read_lisp_array (<file_name>, <A>, <separator_flag>)
 -- Function: read_lisp_array (<S>, <A>)
 -- Function: read_lisp_array (<S>, <A>, <separator_flag>)
     Reads the file <file_name> or file stream <S> and returns its
     entire content as a Lisp array.  If <separator_flag> is not
     specified, the file is assumed space-delimited.

     `read_lisp_array' requires that the array be declared by
     `make_array' before calling the read function. (This obviates the
     need to infer the array dimensions, which could be a problem for
     arrays with multiple dimensions.)

     `read_lisp_array' does not check to see that the input file
     conforms in some way to the array dimensions; the input is read as
     a flat list, then the array is filled using `fillarray'.


 -- Function: read_maxima_array (<file_name>, <A>)
 -- Function: read_maxima_array (<file_name>, <A>, <separator_flag>)
 -- Function: read_maxima_array (<S>, <A>)
 -- Function: read_maxima_array (<S>, <A>, <separator_flag>)
     Reads the file <file_name> or file stream <S> and returns its
     entire content as a Maxima array.  If <separator_flag> is not
     specified, the file is assumed space-delimited.

     `read_maxima_array' requires that the array be declared by `array'
     before calling the read function. (This obviates the need to infer
     the array dimensions, which could be a hassle for arrays with
     multiple dimensions.)

     `read_maxima_array' does not check to see that the input file
     conforms in some way to the array dimensions; the input is read as
     a flat list, then the array is filled using `fillarray'.


 -- Function: read_hashed_array (<file_name>, <A>)
 -- Function: read_hashed_array (<file_name>, <A>, <separator_flag>)
 -- Function: read_hashed_array (<S>, <A>)
 -- Function: read_hashed_array (<S>, <A>, <separator_flag>)
     Reads the file <file_name> or file stream <S> and returns its
     entire content as a hashed array.  If <separator_flag> is not
     specified, the file is assumed space-delimited.

     `read_hashed_array' treats the first item on a line as a hash key,
     and associates the remainder of the line (as a list) with the key.
     For example, the line `567 12 17 32 55' is equivalent to `A[567]:
     [12, 17, 32, 55]$'.  Lines need not have the same numbers of
     elements.


 -- Function: read_nested_list (<file_name>)
 -- Function: read_nested_list (<file_name>, <separator_flag>)
 -- Function: read_nested_list (<S>)
 -- Function: read_nested_list (<S>, <separator_flag>)
     Reads the file <file_name> or file stream <S> and returns its
     entire content as a nested list.  If <separator_flag> is not
     specified, the file is assumed space-delimited.

     `read_nested_list' returns a list which has a sublist for each
     line of input. Lines need not have the same numbers of elements.
     Empty lines are not ignored: an empty line yields an empty sublist.


 -- Function: read_list (<file_name>)
 -- Function: read_list (<file_name>, <separator_flag>)
 -- Function: read_list (<S>)
 -- Function: read_list (<S>, <separator_flag>)
     Reads the file <file_name> or file stream <S> and returns its
     entire content as a flat list.  If <separator_flag> is not
     specified, the file is assumed space-delimited.


 -- Function: write_data (<X>, <file_name>)
 -- Function: write_data (<object>, <file_name>, <separator_flag>)
 -- Function: write_data (<X>, <S>)
 -- Function: write_data (<object>, <S>, <separator_flag>)
     `write_data' writes the object <X> to the file <file_name> or file
     stream <S>.

     `write_data' writes matrices in row-major form, with one line per
     row.

     `write_data' writes Lisp and Maxima declared arrays in row-major
     form, with a new line at the end of every slab.
     Higher-dimensional slabs are separated by additional new lines.

     `write_data' writes hashed arrays with a key followed by the
     associated list on each line.

     `write_data' writes a nested list with each sublist on one line.

     `write_data' writes a flat list all on one line.

     When `write_data' writes to a file, the global variable
     `file_output_append' governs whether the output file is appended
     or truncated.

     No special action concerning appending or truncating is taken when
     `write_data' writes to a file stream.  The file stream is not
     closed by `write_data'.



File: maxima.info,  Node: opsubst,  Next: orthopoly,  Prev: numericalio,  Up: Top

64 opsubst
**********

* Menu:

* Functions and Variables for opsubst::


File: maxima.info,  Node: Functions and Variables for opsubst,  Prev: opsubst,  Up: opsubst

64.1 Functions and Variables for opsubst
========================================

 -- Function: opsubst (<f>,<g>,<e>)
 -- Function: opsubst (<g>=<f>,<e>)
 -- Function: opsubst ([<g1>=<f1>,<g2>=<f2>,..., <gn>=<fn>],<e>)
     The function `opsubst' is similar to the function `subst', except
     that `opsubst' only makes substitutions for the operators in an
     expression. In general, When <f> is an operator in the expression
     <e>, substitute <g> for <f> in the expression <e>.

     To determine the operator, `opsubst' sets `inflag' to true. This
     means `opsubst' substitutes for the internal, not the displayed,
     operator in the expression.

     Examples:
          (%i1) load (opsubst)$

          (%i2) opsubst(f,g,g(g(x)));
          (%o2)                     f(f(x))
          (%i3) opsubst(f,g,g(g));
          (%o3)                       f(g)
          (%i4) opsubst(f,g[x],g[x](z));
          (%o4)                       f(z)
          (%i5) opsubst(g[x],f, f(z));
          (%o5)                      g (z)
                                      x
          (%i6) opsubst(tan, sin, sin(sin));
          (%o6)                     tan(sin)
          (%i7) opsubst([f=g,g=h],f(x));
          (%o7)                       h(x)

     Internally, Maxima does not use the unary negation, division, or
     the subtraction operators; thus:
          (%i8) opsubst("+","-",a-b);
          (%o8)                     a - b
          (%i9) opsubst("f","-",-a);
          (%o9)                      - a
          (%i10) opsubst("^^","//",a/b);
                                       a
          (%o10)                       -
                                       b

     The internal representation of -a*b is *(-1,a,b); thus
          (%i11) opsubst("[","*", -a*b);
          (%o11)                  [- 1, a, b]

     When either operator isn't a Maxima symbol, generally some other
     function will signal an error:
          (%i12) opsubst(a+b,f, f(x));

          Improper name or value in functional position:
          b + a
           -- an error.  Quitting.  To debug this try debugmode(true);

     However, subscripted operators are allowed:
          (%i13) opsubst(g[5],f, f(x));
          (%o13)                     g (x)
                                      5

     To use this function write first `load("opsubst")'.


File: maxima.info,  Node: orthopoly,  Next: plotdf,  Prev: opsubst,  Up: Top

65 orthopoly
************

* Menu:

* Introduction to orthogonal polynomials::
* Functions and Variables for orthogonal polynomials::


File: maxima.info,  Node: Introduction to orthogonal polynomials,  Next: Functions and Variables for orthogonal polynomials,  Prev: orthopoly,  Up: orthopoly

65.1 Introduction to orthogonal polynomials
===========================================

`orthopoly' is a package for symbolic and numerical evaluation of
several kinds of orthogonal polynomials, including Chebyshev, Laguerre,
Hermite, Jacobi, Legendre, and ultraspherical (Gegenbauer) polynomials.
Additionally, `orthopoly' includes support for the spherical Bessel,
spherical Hankel, and spherical harmonic functions.

   For the most part, `orthopoly' follows the conventions of Abramowitz
and Stegun Handbook of Mathematical Functions, Chapter 22 (10th
printing, December 1972); additionally, we use Gradshteyn and Ryzhik,
Table of Integrals, Series, and Products (1980 corrected and enlarged
edition), and Eugen Merzbacher Quantum Mechanics (2nd edition, 1970).

   Barton Willis of the University of Nebraska at Kearney (UNK) wrote
the `orthopoly' package and its documentation. The package is released
under the GNU General Public License (GPL).

65.1.1 Getting Started with orthopoly
-------------------------------------

`load (orthopoly)' loads the `orthopoly' package.

   To find the third-order Legendre polynomial,

     (%i1) legendre_p (3, x);
                           3             2
                  5 (1 - x)    15 (1 - x)
     (%o1)      - ---------- + ----------- - 6 (1 - x) + 1
                      2             2

   To express this as a sum of powers of <x>, apply <ratsimp> or <rat>
to the result.

     (%i2) [ratsimp (%), rat (%)];
                             3           3
                          5 x  - 3 x  5 x  - 3 x
     (%o2)/R/            [----------, ----------]
                              2           2

   Alternatively, make the second argument to `legendre_p' (its "main"
variable) a canonical rational expression (CRE).

     (%i1) legendre_p (3, rat (x));
                                   3
                                5 x  - 3 x
     (%o1)/R/                   ----------
                                    2

   For floating point evaluation, `orthopoly' uses a running error
analysis to estimate an upper bound for the error. For example,

     (%i1) jacobi_p (150, 2, 3, 0.2);
     (%o1) interval(- 0.062017037936715, 1.533267919277521E-11)

   Intervals have the form `interval (<c>, <r>)', where <c> is the
center and <r> is the radius of the interval. Since Maxima does not
support arithmetic on intervals, in some situations, such as graphics,
you want to suppress the error and output only the center of the
interval. To do this, set the option variable
`orthopoly_returns_intervals' to `false'.

     (%i1) orthopoly_returns_intervals : false;
     (%o1)                         false
     (%i2) jacobi_p (150, 2, 3, 0.2);
     (%o2)                  - 0.062017037936715

   Refer to the section *note Floating point Evaluation:: for more
information.

   Most functions in `orthopoly' have a `gradef' property; thus

     (%i1) diff (hermite (n, x), x);
     (%o1)                     2 n H     (x)
                                    n - 1
     (%i2) diff (gen_laguerre (n, a, x), x);
                   (a)               (a)
                n L   (x) - (n + a) L     (x) unit_step(n)
                   n                 n - 1
     (%o2)      ------------------------------------------
                                    x

   The unit step function in the second example prevents an error that
would otherwise arise by evaluating with <n> equal to 0.

     (%i3) ev (%, n = 0);
     (%o3)                           0

   The `gradef' property only applies to the "main" variable;
derivatives with respect other arguments usually result in an error
message; for example

     (%i1) diff (hermite (n, x), x);
     (%o1)                     2 n H     (x)
                                    n - 1
     (%i2) diff (hermite (n, x), n);

     Maxima doesn't know the derivative of hermite with respect the first
     argument
      -- an error.  Quitting.  To debug this try debugmode(true);

   Generally, functions in `orthopoly' map over lists and matrices. For
the mapping to fully evaluate, the option variables `doallmxops' and
`listarith' must both be `true' (the defaults).  To illustrate the
mapping over matrices, consider

     (%i1) hermite (2, x);
                                          2
     (%o1)                    - 2 (1 - 2 x )
     (%i2) m : matrix ([0, x], [y, 0]);
                                 [ 0  x ]
     (%o2)                       [      ]
                                 [ y  0 ]
     (%i3) hermite (2, m);
                    [                             2  ]
                    [      - 2        - 2 (1 - 2 x ) ]
     (%o3)          [                                ]
                    [             2                  ]
                    [ - 2 (1 - 2 y )       - 2       ]

   In the second example, the `i, j' element of the value is `hermite
(2, m[i,j])'; this is not the same as computing `-2 + 4 m . m', as seen
in the next example.

     (%i4) -2 * matrix ([1, 0], [0, 1]) + 4 * m . m;
                         [ 4 x y - 2      0     ]
     (%o4)               [                      ]
                         [     0      4 x y - 2 ]

   If you evaluate a function at a point outside its domain, generally
`orthopoly' returns the function unevaluated. For example,

     (%i1) legendre_p (2/3, x);
     (%o1)                        P   (x)
                                   2/3

   `orthopoly' supports translation into TeX; it also does
two-dimensional output on a terminal.

     (%i1) spherical_harmonic (l, m, theta, phi);
                               m
     (%o1)                    Y (theta, phi)
                               l
     (%i2) tex (%);
     $$Y_{l}^{m}\left(\vartheta,\varphi\right)$$
     (%o2)                         false
     (%i3) jacobi_p (n, a, a - b, x/2);
                               (a, a - b) x
     (%o3)                    P          (-)
                               n          2
     (%i4) tex (%);
     $$P_{n}^{\left(a,a-b\right)}\left({{x}\over{2}}\right)$$
     (%o4)                         false

65.1.2 Limitations
------------------

When an expression involves several orthogonal polynomials with
symbolic orders, it's possible that the expression actually vanishes,
yet Maxima is unable to simplify it to zero. If you divide by such a
quantity, you'll be in trouble. For example, the following expression
vanishes for integers <n> greater than 1, yet Maxima is unable to
simplify it to zero.

     (%i1) (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x)
           + (1 - n) * legendre_p (n - 2, x);
     (%o1)  (2 n - 1) P     (x) x - n P (x) + (1 - n) P     (x)
                       n - 1           n               n - 2

   For a specific <n>, we can reduce the expression to zero.

     (%i2) ev (% ,n = 10, ratsimp);
     (%o2)                           0

   Generally, the polynomial form of an orthogonal polynomial is
ill-suited for floating point evaluation. Here's an example.

     (%i1) p : jacobi_p (100, 2, 3, x)$

     (%i2) subst (0.2, x, p);
     (%o2)                3.4442767023833592E+35
     (%i3) jacobi_p (100, 2, 3, 0.2);
     (%o3)  interval(0.18413609135169, 6.8990300925815987E-12)
     (%i4) float(jacobi_p (100, 2, 3, 2/10));
     (%o4)                   0.18413609135169

   The true value is about 0.184; this calculation suffers from extreme
subtractive cancellation error. Expanding the polynomial and then
evaluating, gives a better result.
     (%i5) p : expand(p)$
     (%i6) subst (0.2, x, p);
     (%o6) 0.18413609766122982

   This isn't a general rule; expanding the polynomial does not always
result in an expression that is better suited for numerical evaluation.
By far, the best way to do numerical evaluation is to make one or more
of the function arguments floating point numbers. By doing that,
specialized floating point algorithms are used for evaluation.

   Maxima's `float' function is somewhat indiscriminate; if you apply
`float' to an expression involving an orthogonal polynomial with a
symbolic degree or order parameter, these parameters may be converted
into floats; after that, the expression will not evaluate fully.
Consider

     (%i1) assoc_legendre_p (n, 1, x);
                                    1
     (%o1)                         P (x)
                                    n
     (%i2) float (%);
                                   1.0
     (%o2)                        P   (x)
                                   n
     (%i3) ev (%, n=2, x=0.9);
                                  1.0
     (%o3)                       P   (0.9)
                                  2

   The expression in (%o3) will not evaluate to a float; `orthopoly'
doesn't recognize floating point values where it requires an integer.
Similarly, numerical evaluation of the `pochhammer' function for orders
that exceed `pochhammer_max_index' can be troublesome; consider

     (%i1) x :  pochhammer (1, 10), pochhammer_max_index : 5;
     (%o1)                         (1)
                                      10

   Applying `float' doesn't evaluate <x> to a float

     (%i2) float (x);
     (%o2)                       (1.0)
                                      10.0

   To evaluate <x> to a float, you'll need to bind
`pochhammer_max_index' to 11 or greater and apply `float' to <x>.

     (%i3) float (x), pochhammer_max_index : 11;
     (%o3)                       3628800.0

   The default value of `pochhammer_max_index' is 100; change its value
after loading `orthopoly'.

   Finally, be aware that reference books vary on the definitions of the
orthogonal polynomials; we've generally used the conventions of
conventions of Abramowitz and Stegun.

   Before you suspect a bug in orthopoly, check some special cases to
determine if your definitions match those used by `orthopoly'.
Definitions often differ by a normalization; occasionally, authors use
"shifted" versions of the functions that makes the family orthogonal on
an interval other than (-1, 1). To define, for example, a Legendre
polynomial that is orthogonal on (0, 1), define

     (%i1) shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$

     (%i2) shifted_legendre_p (2, rat (x));
                                 2
     (%o2)/R/                 6 x  - 6 x + 1
     (%i3) legendre_p (2, rat (x));
                                    2
                                 3 x  - 1
     (%o3)/R/                    --------
                                    2

65.1.3 Floating point Evaluation
--------------------------------

Most functions in `orthopoly' use a running error analysis to estimate
the error in floating point evaluation; the exceptions are the
spherical Bessel functions and the associated Legendre polynomials of
the second kind. For numerical evaluation, the spherical Bessel
functions call SLATEC functions. No specialized method is used for
numerical evaluation of the associated Legendre polynomials of the
second kind.

   The running error analysis ignores errors that are second or higher
order in the machine epsilon (also known as unit roundoff). It also
ignores a few other errors. It's possible (although unlikely) that the
actual error exceeds the estimate.

   Intervals have the form `interval (<c>, <r>)', where <c> is the
center of the interval and <r> is its radius. The center of an interval
can be a complex number, and the radius is always a positive real
number.

   Here is an example.

     (%i1) fpprec : 50$

     (%i2) y0 : jacobi_p (100, 2, 3, 0.2);
     (%o2) interval(0.1841360913516871, 6.8990300925815987E-12)
     (%i3) y1 : bfloat (jacobi_p (100, 2, 3, 1/5));
     (%o3) 1.8413609135168563091370224958913493690868904463668b-1

   Let's test that the actual error is smaller than the error estimate

     (%i4) is (abs (part (y0, 1) - y1) < part (y0, 2));
     (%o4)                         true

   Indeed, for this example the error estimate is an upper bound for the
true error.

   Maxima does not support arithmetic on intervals.

     (%i1) legendre_p (7, 0.1) + legendre_p (8, 0.1);
     (%o1) interval(0.18032072148437508, 3.1477135311021797E-15)
             + interval(- 0.19949294375000004, 3.3769353084291579E-15)

   A user could define arithmetic operators that do interval math. To
define interval addition, we can define

     (%i1) infix ("@+")$

     (%i2) "@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2)
           + part (y, 2))$

     (%i3) legendre_p (7, 0.1) @+ legendre_p (8, 0.1);
     (%o3) interval(- 0.019172222265624955, 6.5246488395313372E-15)

   The special floating point routines get called when the arguments
are complex.  For example,

     (%i1) legendre_p (10, 2 + 3.0*%i);
     (%o1) interval(- 3.876378825E+7 %i - 6.0787748E+7,
                                                1.2089173052721777E-6)

   Let's compare this to the true value.

     (%i1) float (expand (legendre_p (10, 2 + 3*%i)));
     (%o1)          - 3.876378825E+7 %i - 6.0787748E+7

   Additionally, when the arguments are big floats, the special
floating point routines get called; however, the big floats are
converted into double floats and the final result is a double.

     (%i1) ultraspherical (150, 0.5b0, 0.9b0);
     (%o1) interval(- 0.043009481257265, 3.3750051301228864E-14)

65.1.4 Graphics and `orthopoly'
-------------------------------

To plot expressions that involve the orthogonal polynomials, you must
do two things:
  1. Set the option variable `orthopoly_returns_intervals' to `false',

  2. Quote any calls to `orthopoly' functions.
        If function calls aren't quoted, Maxima evaluates them to
polynomials before plotting; consequently, the specialized floating
point code doesn't get called.  Here is an example of how to plot an
expression that involves a Legendre polynomial.

     (%i1) plot2d ('(legendre_p (5, x)), [x, 0, 1]),
                             orthopoly_returns_intervals : false;
     (%o1)

   The entire expression `legendre_p (5, x)' is quoted; this is
different than just quoting the function name using `'legendre_p (5,
<x>)'.

65.1.5 Miscellaneous Functions
------------------------------

The `orthopoly' package defines the Pochhammer symbol and a unit step
function. `orthopoly' uses the Kronecker delta function and the unit
step function in `gradef' statements.

   To convert Pochhammer symbols into quotients of gamma functions, use
`makegamma'.

     (%i1) makegamma (pochhammer (x, n));
                               gamma(x + n)
     (%o1)                     ------------
                                 gamma(x)
     (%i2) makegamma (pochhammer (1/2, 1/2));
                                     1
     (%o2)                       ---------
                                 sqrt(%pi)

   Derivatives of the Pochhammer symbol are given in terms of the `psi'
function.

     (%i1) diff (pochhammer (x, n), x);
     (%o1)             (x)  (psi (x + n) - psi (x))
                          n     0             0
     (%i2) diff (pochhammer (x, n), n);
     (%o2)                   (x)  psi (x + n)
                                n    0

   You need to be careful with the expression in (%o1); the difference
of the `psi' functions has polynomials when `<x> = -1, -2, .., -<n>'.
These polynomials cancel with factors in `pochhammer (<x>, <n>)' making
the derivative a degree `<n> - 1' polynomial when <n> is a positive
integer.

   The Pochhammer symbol is defined for negative orders through its
representation as a quotient of gamma functions. Consider

     (%i1) q : makegamma (pochhammer (x, n));
                               gamma(x + n)
     (%o1)                     ------------
                                 gamma(x)
     (%i2) sublis ([x=11/3, n= -6], q);
                                    729
     (%o2)                        - ----
                                    2240

   Alternatively, we can get this result directly.

     (%i1) pochhammer (11/3, -6);
                                    729
     (%o1)                        - ----
                                    2240

   The unit step function is left-continuous; thus

     (%i1) [unit_step (-1/10), unit_step (0), unit_step (1/10)];
     (%o1)                       [0, 0, 1]

   If you need a unit step function that is neither left or right
continuous at zero, define your own using `signum'; for example,

     (%i1) xunit_step (x) := (1 + signum (x))/2$

     (%i2) [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)];
                                     1
     (%o2)                       [0, -, 1]
                                     2

   Do not redefine `unit_step' itself; some code in `orthopoly'
requires that the unit step function be left-continuous.

65.1.6 Algorithms
-----------------

Generally, `orthopoly' does symbolic evaluation by using a hypergeometic
representation of the orthogonal polynomials. The hypergeometic
functions are evaluated using the (undocumented) functions `hypergeo11'
and `hypergeo21'. The exceptions are the half-integer Bessel functions
and the associated Legendre function of the second kind. The
half-integer Bessel functions are evaluated using an explicit
representation, and the associated Legendre function of the second kind
is evaluated using recursion.

   For floating point evaluation, we again convert most functions into
a hypergeometic form; we evaluate the hypergeometic functions using
forward recursion. Again, the exceptions are the half-integer Bessel
functions and the associated Legendre function of the second kind.
Numerically, the half-integer Bessel functions are evaluated using the
SLATEC code.


File: maxima.info,  Node: Functions and Variables for orthogonal polynomials,  Prev: Introduction to orthogonal polynomials,  Up: orthopoly

65.2 Functions and Variables for orthogonal polynomials
=======================================================

 -- Function: assoc_legendre_p (<n>, <m>, <x>)
     The associated Legendre function of the first kind.

     Reference: Abramowitz and Stegun, equations 22.5.37, page 779,
     8.6.6 (second equation), page 334, and 8.2.5, page 333.

 -- Function: assoc_legendre_q (<n>, <m>, <x>)
     The associated Legendre function of the second kind.

     Reference: Abramowitz and Stegun, equation 8.5.3 and 8.1.8.

 -- Function: chebyshev_t (<n>, <x>)
     The Chebyshev function of the first kind.

     Reference: Abramowitz and Stegun, equation 22.5.47, page 779.

 -- Function: chebyshev_u (<n>, <x>)
     The Chebyshev function of the second kind.

     Reference: Abramowitz and Stegun, equation 22.5.48, page 779.

 -- Function: gen_laguerre (<n>, <a>, <x>)
     The generalized Laguerre polynomial.

     Reference: Abramowitz and Stegun, equation 22.5.54, page 780.

 -- Function: hermite (<n>, <x>)
     The Hermite polynomial.

     Reference: Abramowitz and Stegun, equation 22.5.55, page 780.

 -- Function: intervalp (<e>)
     Return `true' if the input is an interval and return false if it
     isn't.

 -- Function: jacobi_p (<n>, <a>, <b>, <x>)
     The Jacobi polynomial.

     The Jacobi polynomials are actually defined for all <a> and <b>;
     however, the Jacobi polynomial weight `(1 - <x>)^<a> (1 +
     <x>)^<b>' isn't integrable for `<a> <= -1' or `<b> <= -1'.

     Reference: Abramowitz and Stegun, equation 22.5.42, page 779.

 -- Function: laguerre (<n>, <x>)
     The Laguerre polynomial.

     Reference: Abramowitz and Stegun, equations 22.5.16 and 22.5.54,
     page 780.

 -- Function: legendre_p (<n>, <x>)
     The Legendre polynomial of the first kind.

     Reference: Abramowitz and Stegun, equations 22.5.50 and 22.5.51,
     page 779.

 -- Function: legendre_q (<n>, <x>)
     The Legendre polynomial of the first kind.

     Reference: Abramowitz and Stegun, equations 8.5.3 and 8.1.8.

 -- Function: orthopoly_recur (<f>, <args>)
     Returns a recursion relation for the orthogonal function family
     <f> with arguments <args>. The recursion is with respect to the
     polynomial degree.

          (%i1) orthopoly_recur (legendre_p, [n, x]);
                          (2 n - 1) P     (x) x + (1 - n) P     (x)
                                     n - 1                 n - 2
          (%o1)   P (x) = -----------------------------------------
                   n                          n

     The second argument to `orthopoly_recur' must be a list with the
     correct number of arguments for the function <f>; if it isn't,
     Maxima signals an error.

          (%i1) orthopoly_recur (jacobi_p, [n, x]);

          Function jacobi_p needs 4 arguments, instead it received 2
           -- an error.  Quitting.  To debug this try debugmode(true);

     Additionally, when <f> isn't the name of one of the families of
     orthogonal polynomials, an error is signalled.

          (%i1) orthopoly_recur (foo, [n, x]);

          A recursion relation for foo isn't known to Maxima
           -- an error.  Quitting.  To debug this try debugmode(true);

 -- Variable: orthopoly_returns_intervals
     Default value: `true'

     When `orthopoly_returns_intervals' is `true', floating point
     results are returned in the form `interval (<c>, <r>)', where <c>
     is the center of an interval and <r> is its radius. The center can
     be a complex number; in that case, the interval is a disk in the
     complex plane.

 -- Function: orthopoly_weight (<f>, <args>)
     Returns a three element list; the first element is the formula of
     the weight for the orthogonal polynomial family <f> with arguments
     given by the list <args>; the second and third elements give the
     lower and upper endpoints of the interval of orthogonality. For
     example,

          (%i1) w : orthopoly_weight (hermite, [n, x]);
                                      2
                                   - x
          (%o1)                 [%e    , - inf, inf]
          (%i2) integrate(w[1]*hermite(3, x)*hermite(2, x), x, w[2], w[3]);
          (%o2)                           0

     The main variable of <f> must be a symbol; if it isn't, Maxima
     signals an error.


 -- Function: pochhammer (<n>, <x>)
     The Pochhammer symbol. For nonnegative integers <n> with `<n> <=
     pochhammer_max_index', the expression `pochhammer (<x>, <n>)'
     evaluates to the product `<x> (<x> + 1) (<x> + 2) ... (<x> + n -
     1)' when `<n> > 0' and to 1 when `<n> = 0'. For negative <n>,
     `pochhammer (<x>, <n>)' is defined as `(-1)^<n> / pochhammer (1 -
     <x>, -<n>)'.  Thus

          (%i1) pochhammer (x, 3);
          (%o1)                   x (x + 1) (x + 2)
          (%i2) pochhammer (x, -3);
                                           1
          (%o2)               - -----------------------
                                (1 - x) (2 - x) (3 - x)

     To convert a Pochhammer symbol into a quotient of gamma functions,
     (see Abramowitz and Stegun, equation 6.1.22) use `makegamma'; for
     example

          (%i1) makegamma (pochhammer (x, n));
                                    gamma(x + n)
          (%o1)                     ------------
                                      gamma(x)

     When <n> exceeds `pochhammer_max_index' or when <n> is symbolic,
     `pochhammer' returns a noun form.

          (%i1) pochhammer (x, n);
          (%o1)                         (x)
                                           n

 -- Variable: pochhammer_max_index
     Default value: 100

     `pochhammer (<n>, <x>)' expands to a product if and only if `<n>
     <= pochhammer_max_index'.

     Examples:

          (%i1) pochhammer (x, 3), pochhammer_max_index : 3;
          (%o1)                   x (x + 1) (x + 2)
          (%i2) pochhammer (x, 4), pochhammer_max_index : 3;
          (%o2)                         (x)
                                           4

     Reference: Abramowitz and Stegun, equation 6.1.16, page 256.

 -- Function: spherical_bessel_j (<n>, <x>)
     The spherical Bessel function of the first kind.

     Reference: Abramowitz and Stegun, equations 10.1.8, page 437 and
     10.1.15, page 439.

 -- Function: spherical_bessel_y (<n>, <x>)
     The spherical Bessel function of the second kind.

     Reference: Abramowitz and Stegun, equations 10.1.9, page 437 and
     10.1.15, page 439.

 -- Function: spherical_hankel1 (<n>, <x>)
     The spherical Hankel function of the first kind.

     Reference: Abramowitz and Stegun, equation 10.1.36, page 439.

 -- Function: spherical_hankel2 (<n>, <x>)
     The spherical Hankel function of the second kind.

     Reference: Abramowitz and Stegun, equation 10.1.17, page 439.

 -- Function: spherical_harmonic (<n>, <m>, <x>, <y>)
     The spherical harmonic function.

     Reference: Merzbacher 9.64.

 -- Function: unit_step (<x>)
     The left-continuous unit step function; thus `unit_step (<x>)'
     vanishes for `<x> <= 0' and equals 1 for `<x> > 0'.

     If you want a unit step function that takes on the value 1/2 at
     zero, use `(1 + signum (<x>))/2'.

 -- Function: ultraspherical (<n>, <a>, <x>)
     The ultraspherical polynomial (also known as the Gegenbauer
     polynomial).

     Reference: Abramowitz and Stegun, equation 22.5.46, page 779.


File: maxima.info,  Node: plotdf,  Next: romberg,  Prev: orthopoly,  Up: Top

66 plotdf
*********

* Menu:

* Introduction to plotdf::
* Functions and Variables for plotdf::


File: maxima.info,  Node: Introduction to plotdf,  Next: Functions and Variables for plotdf,  Prev: plotdf,  Up: plotdf

66.1 Introduction to plotdf
===========================

The function `plotdf' creates a plot of the direction field of a
first-order Ordinary Differential Equation (ODE) or a system of two
autonomous first-order ODE's.

   Since this is an additional package, in order to use it you must
first load it with `load("plotdf")'. Plotdf requires Openmath, which is
provided by the package Xmaxima (Xmaxima is not only used as a graphical
console for Maxima but also to plot graphs in the Openmath format).

   To plot the direction field of a single ODE, the ODE must be written
in the form:
            dy
            -- = F(x,y)
            dx

   and the function <F> should be given as the argument for `plotdf'.
If the independent and dependent variables are not <x>, and <y>, as in
the equation above, then those two variables should be named explicitly
in a list given as an argument to the plotdf command (see the examples).

   To plot the direction field of a set of two autonomous ODE's, they
must be written in the form
            dx             dy
            -- = G(x,y)    -- = F(x,y)
            dt             dt

   and the argument for `plotdf' should be a list with the two
functions <G> and <F>, in that order; namely, the first expression in
the list will be taken to be the time derivative of the variable
represented on the horizontal axis, and the second expression will be
the time derivative of the variable represented on the vertical axis.
Those two variables do not have to be <x> and <y>, but if they are not,
then the second argument given to plotdf must be another list naming
the two variables, first the one on the horizontal axis and then the
one on the vertical axis.

   If only one ODE is given, `plotdf' will implicitly admit `x=t', and
`G(x,y)=1', transforming the non-autonomous equation into a system of
two autonomous equations.


File: maxima.info,  Node: Functions and Variables for plotdf,  Prev: Introduction to plotdf,  Up: plotdf

66.2 Functions and Variables for plotdf
=======================================

 -- Function: plotdf (<dydx>, ...options...)
 -- Function: plotdf (<dvdu>, `['<u>,<v>`]', ...options...)
 -- Function: plotdf (`['<dxdt>,<dydt>`]', ...options...)
 -- Function: plotdf (`['<dudt>,<dvdt>`]', `['<u>,<v>`]', ...options...)
     Displays a direction field in two dimensions <x> and <y>.

     <dydx>, <dxdt> and <dydt> are expressions that depend on <x> and
     <y>. <dvdu>, <dudt> and <dvdt> are expressions that depend on <u>
     and <v>. In addition to those two variables, the expressions can
     also depend on a set of parameters, with numerical values given
     with the `parameters' option (the option syntax is given below),
     or with a range of allowed values specified by a <sliders> option.

     Several other options can be given within the command, or selected
     in the menu. Integral curves can be obtained by clicking on the
     plot, or with the option `trajectory_at'. The direction of the
     integration can be controlled with the `direction' option, which
     can have values of _forward_, _backward_ or _both_. The number of
     integration steps is given by `nsteps' and the time interval
     between them is set up with the `tstep' option. The Adams Moulton
     method is used for the integration; it is also possible to switch
     to an adaptive Runge-Kutta 4th order method.

     Plot window menu:

     The menu in the plot window has the following options: _Zoom_, will
     change the behavior of the mouse so that it will allow you to zoom
     in on a region of the plot by clicking with the left button. Each
     click near a point magnifies the plot, keeping the center at the
     point where you clicked. Holding the <Shift> key while clicking,
     zooms out to the previous magnification. To resume computing
     trajectories when you click on a point, select _Integrate_ from
     the menu.

     The option _Config_ in the menu can be used to change the ODE(s) in
     use and various other settings. After configuration changes are
     made, the menu option _Replot_ should be selected, to activate the
     new settings.  If a pair of coordinates are entered in the field
     _Trajectory at_ in the _Config_ dialog menu, and the <enter> key
     is pressed, a new integral curve will be shown, in addition to the
     ones already shown. When _Replot_ is selected, only the last
     integral curve entered will be shown.

     Holding the right mouse button down while the cursor is moved, can
     be used to drag the plot sideways or up and down. Additional
     parameters such as the number of steps, the initial value of <t>
     and the x and y centers and radii, may be set in the Config menu.

     A copy of the plot can be saved as a postscript file, using the
     menu option _Save_.

     Plot options:

     The `plotdf' command may include several commands, each command is
     a list of two or more items.  The first item is the name of the
     option, and the remainder comprises the value or values assigned
     to the option.

     The options which are recognized by `plotdf' are the following:

        * "tstep" defines the length of the increments on the
          independent variable <t>, used to compute an integral curve.
          If only one expression <dydx> is given to `plotdf', the <x>
          variable will be directly proportional to <t>.  The default
          value is 0.1.

        * "nsteps" defines the number of steps of length `tstep' that
          will be used for the independent variable, to compute an
          integral curve.  The default value is 100.

        * "direction" defines the direction of the independent variable
          that will be followed to compute an integral curve. Possible
          values are `forward', to make the independent variable
          increase `nsteps' times, with increments `tstep', `backward',
          to make the independent variable decrease, or `both' that
          will lead to an integral curve that extends `nsteps' forward,
          and `nsteps' backward. The keywords `right' and `left' can be
          used as synonyms for `forward' and `backward'.  The default
          value is `both'.

        * "tinitial" defines the initial value of variable <t> used to
          compute integral curves. Since the differential equations are
          autonomous, that setting will only appear in the plot of the
          curves as functions of <t>.  The default value is 0.

        * "versus_t" is used to create a second plot window, with a
          plot of an integral curve, as two functions <x>, <y>, of the
          independent variable <t>. If `versus_t' is given any value
          different from 0, the second plot window will be displayed.
          The second plot window includes another menu, similar to the
          menu of the main plot window.  The default value is 0.

        * "trajectory_at" defines the coordinates <xinitial> and
          <yinitial> for the starting point of an integral curve.  The
          option is empty by default.

        * "parameters" defines a list of parameters, and their
          numerical values, used in the definition of the differential
          equations. The name and values of the parameters must be
          given in a string with a comma-separated sequence of pairs
          `name=value'.

        * "sliders" defines a list of parameters that will be changed
          interactively using slider buttons, and the range of
          variation of those parameters. The names and ranges of the
          parameters must be given in a string with a comma-separated
          sequence of elements `name=min:max'

        * "xfun" defines a string with semi-colon-separated sequence of
          functions of <x> to be displayed, on top of the direction
          field.  Those functions will be parsed by Tcl and not by
          Maxima.

        * "x" should be followed by two numbers, which will set up the
          minimum and maximum values shown on the horizontal axis. If
          the variable on the horizontal axis is not <x>, then this
          option should have the name of the variable on the horizontal
          axis.  The default horizontal range is from -10 to 10.

        * "y" should be followed by two numbers, which will set up the
          minimum and maximum values shown on the vertical axis. If the
          variable on the vertical axis is not <y>, then this option
          should have the name of the variable on the vertical axis.
          The default vertical range is from -10 to 10.

     Examples:

        * To show the direction field of the differential equation y' =
          exp(-x) + y and the solution that goes through (2, -0.1):
               (%i1) load("plotdf")$

               (%i2) plotdf(exp(-x)+y,[trajectory_at,2,-0.1])$

        * To obtain the direction field for the equation diff(y,x) = x
          - y^2 and the solution with initial condition y(-1) = 3, we
          can use the command:
               (%i3) plotdf(x-y^2,[xfun,"sqrt(x);-sqrt(x)"],
                         [trajectory_at,-1,3], [direction,forward],
                         [y,-5,5], [x,-4,16])$
          The graph also shows the function y = sqrt(x).

        * The following example shows the direction field of a harmonic
          oscillator, defined by the two equations dz/dt = v and dv/dt
          = -k*z/m, and the integral curve through (z,v) = (6,0), with
          a slider that will allow you to change the value of m
          interactively (k is fixed at 2):
               (%i4) plotdf([v,-k*z/m], [z,v], [parameters,"m=2,k=2"],
                           [sliders,"m=1:5"], [trajectory_at,6,0])$

        * To plot the direction field of the Duffing equation,
          m*x"+c*x'+k*x+b*x^3 = 0, we introduce the variable y=x' and
          use:
               (%i5) plotdf([y,-(k*x + c*y + b*x^3)/m],
                             [parameters,"k=-1,m=1.0,c=0,b=1"],
                             [sliders,"k=-2:2,m=-1:1"],[tstep,0.1])$

        * The direction field for a damped pendulum, including the
          solution for the given initial conditions, with a slider that
          can be used to change the value of the mass m, and with a
          plot of the two state variables as a function of time:

               (%i6) plotdf([w,-g*sin(a)/l - b*w/m/l], [a,w],
                        [parameters,"g=9.8,l=0.5,m=0.3,b=0.05"],
                        [trajectory_at,1.05,-9],[tstep,0.01],
                        [a,-10,2], [w,-14,14], [direction,forward],
                        [nsteps,300], [sliders,"m=0.1:1"], [versus_t,1])$




File: maxima.info,  Node: romberg,  Next: simplex,  Prev: plotdf,  Up: Top

67 romberg
**********

* Menu:

* Functions and Variables for romberg::


File: maxima.info,  Node: Functions and Variables for romberg,  Prev: Top,  Up: Top

67.1 Functions and Variables for romberg
========================================

 -- Function: romberg (<expr>, <x>, <a>, <b>)
 -- Function: romberg (<F>, <a>, <b>)
     Computes a numerical integration by Romberg's method.

     `romberg(<expr>, <x>, <a>, <b>)' returns an estimate of the
     integral `integrate(<expr>, <x>, <a>, <b>)'.  <expr> must be an
     expression which evaluates to a floating point value when <x> is
     bound to a floating point value.

     `romberg(<F>, <a>, <b>)' returns an estimate of the integral
     `integrate(<F>(x), x, <a>, <b>)' where `x' represents the unnamed,
     sole argument of <F>; the actual argument is not named `x'.  <F>
     must be a Maxima or Lisp function which returns a floating point
     value when the argument is a floating point value.  <F> may name a
     translated or compiled Maxima function.

     The accuracy of `romberg' is governed by the global variables
     `rombergabs' and `rombergtol'.  `romberg' terminates successfully
     when the absolute difference between successive approximations is
     less than `rombergabs', or the relative difference in successive
     approximations is less than `rombergtol'.  Thus when `rombergabs'
     is 0.0 (the default) only the relative error test has any effect
     on `romberg'.

     `romberg' halves the stepsize at most `rombergit' times before it
     gives up; the maximum number of function evaluations is therefore
     `2^rombergit'.  If the error criterion established by `rombergabs'
     and `rombergtol' is not satisfied, `romberg' prints an error
     message.  `romberg' always makes at least `rombergmin' iterations;
     this is a heuristic intended to prevent spurious termination when
     the integrand is oscillatory.

     `romberg' repeatedly evaluates the integrand after binding the
     variable of integration to a specific value (and not before).
     This evaluation policy makes it possible to nest calls to
     `romberg', to compute multidimensional integrals.  However, the
     error calculations do not take the errors of nested integrations
     into account, so errors may be underestimated.  Also, methods
     devised especially for multidimensional problems may yield the
     same accuracy with fewer function evaluations.

     `load(romberg)' loads this function.

     See also `QUADPACK', a collection of numerical integration
     functions.

     Examples:

     A 1-dimensional integration.

          (%i1) load (romberg);
          (%o1)    /usr/share/maxima/5.11.0/share/numeric/romberg.lisp
          (%i2) f(x) := 1/((x - 1)^2 + 1/100) + 1/((x - 2)^2 + 1/1000)
                        + 1/((x - 3)^2 + 1/200);
                              1                 1                1
          (%o2) f(x) := -------------- + --------------- + --------------
                               2    1           2    1            2    1
                        (x - 1)  + ---   (x - 2)  + ----   (x - 3)  + ---
                                   100              1000              200
          (%i3) rombergtol : 1e-6;
          (%o3)                 9.9999999999999995E-7
          (%i4) rombergit : 15;
          (%o4)                          15
          (%i5) estimate : romberg (f(x), x, -5, 5);
          (%o5)                   173.6730736617464
          (%i6) exact : integrate (f(x), x, -5, 5);
          (%o6) 10 sqrt(10) atan(70 sqrt(10))
           + 10 sqrt(10) atan(30 sqrt(10)) + 10 sqrt(2) atan(80 sqrt(2))
           + 10 sqrt(2) atan(20 sqrt(2)) + 10 atan(60) + 10 atan(40)
          (%i7) abs (estimate - exact) / exact, numer;
          (%o7)                7.5527060865060088E-11

     A 2-dimensional integration, implemented by nested calls to
     `romberg'.

          (%i1) load (romberg);
          (%o1)    /usr/share/maxima/5.11.0/share/numeric/romberg.lisp
          (%i2) g(x, y) := x*y / (x + y);
                                              x y
          (%o2)                   g(x, y) := -----
                                             x + y
          (%i3) rombergtol : 1e-6;
          (%o3)                 9.9999999999999995E-7
          (%i4) estimate : romberg (romberg (g(x, y), y, 0, x/2), x, 1, 3);
          (%o4)                   0.81930239628356
          (%i5) assume (x > 0);
          (%o5)                        [x > 0]
          (%i6) integrate (integrate (g(x, y), y, 0, x/2), x, 1, 3);
                                                    3
                                              2 log(-) - 1
                              9                     2        9
          (%o6)       - 9 log(-) + 9 log(3) + ------------ + -
                              2                    6         2
          (%i7) exact : radcan (%);
                              26 log(3) - 26 log(2) - 13
          (%o7)             - --------------------------
                                          3
          (%i8) abs (estimate - exact) / exact, numer;
          (%o8)                1.3711979871851024E-10


 -- Option variable: rombergabs
     Default value: 0.0

     The accuracy of `romberg' is governed by the global variables
     `rombergabs' and `rombergtol'.  `romberg' terminates successfully
     when the absolute difference between successive approximations is
     less than `rombergabs', or the relative difference in successive
     approximations is less than `rombergtol'.  Thus when `rombergabs'
     is 0.0 (the default) only the relative error test has any effect
     on `romberg'.

     See also `rombergit' and `rombergmin'.


 -- Option variable: rombergit
     Default value: 11

     `romberg' halves the stepsize at most `rombergit' times before it
     gives up; the maximum number of function evaluations is therefore
     `2^rombergit'.  `romberg' always makes at least `rombergmin'
     iterations; this is a heuristic intended to prevent spurious
     termination when the integrand is oscillatory.

     See also `rombergabs' and `rombergtol'.


 -- Option variable: rombergmin
     Default value: 0

     `romberg' always makes at least `rombergmin' iterations; this is a
     heuristic intended to prevent spurious termination when the
     integrand is oscillatory.

     See also `rombergit', `rombergabs', and `rombergtol'.


 -- Option variable: rombergtol
     Default value: 1e-4

     The accuracy of `romberg' is governed by the global variables
     `rombergabs' and `rombergtol'.  `romberg' terminates successfully
     when the absolute difference between successive approximations is
     less than `rombergabs', or the relative difference in successive
     approximations is less than `rombergtol'.  Thus when `rombergabs'
     is 0.0 (the default) only the relative error test has any effect
     on `romberg'.

     See also `rombergit' and `rombergmin'.



File: maxima.info,  Node: simplex,  Next: simplification,  Prev: romberg,  Up: Top

68 simplex
**********

* Menu:

* Introduction to simplex::
* Functions and Variables for simplex::


File: maxima.info,  Node: Introduction to simplex,  Next: Functions and Variables for simplex,  Prev: simplex,  Up: simplex

68.1 Introduction to simplex
============================

`simplex' is a package for linear optimization using the simplex
algorithm.

   Example:

     (%i1) load("simplex")$
     (%i2) minimize_lp(x+y, [3*x+2*y>2, x+4*y>3]);
                       9        7       1
     (%o2)            [--, [y = --, x = -]]
                       10       10      5


File: maxima.info,  Node: Functions and Variables for simplex,  Prev: Introduction to simplex,  Up: simplex

68.2 Functions and Variables for simplex
========================================

 -- Option variable: epsilon_lp
     Default value: `10^-8'

     Epsilon used for numerical computations in `linear_program'.

     See also: `linear_program'.


 -- Function: linear_program (<A>, <b>, <c>)
     `linear_program' is an implementation of the simplex algorithm.
     `linear_program(A, b, c)' computes a vector <x> for which `c.x' is
     minimum possible among vectors for which `A.x = b' and `x >= 0'.
     Argument <A> is a matrix and arguments <b> and <c> are lists.

     `linear_program' returns a list which contains the minimizing
     vector <x> and the minimum value `c.x'. If the problem is not
     bounded, it returns "Problem not bounded!" and if the problem is
     not feasible, it returns "Problem not feasible!".

     To use this function first load the `simplex' package with
     `load(simplex);'.

     Example:

          (%i2) A: matrix([1,1,-1,0], [2,-3,0,-1], [4,-5,0,0])$
          (%i3) b: [1,1,6]$
          (%i4) c: [1,-2,0,0]$
          (%i5) linear_program(A, b, c);
                             13     19        3
          (%o5)            [[--, 4, --, 0], - -]
                             2      2         2

     See also: `minimize_lp', `scale_lp', and `epsilon_lp'.


 -- Function: maximize_lp (<obj>, <cond>, [<pos>])
     Maximizes linear objective function <obj> subject to some linear
     constraints <cond>. See `minimize_lp' for detailed description of
     arguments and return value.

     See also: `minimize_lp'.


 -- Function: minimize_lp (<obj>, <cond>, [<pos>])
     Minimizes a linear objective function <obj> subject to some linear
     constraints <cond>. <cond> a list of linear equations or
     inequalities. In strict inequalities `>' is replaced by `>=' and
     `<' by `<='. The optional argument <pos> is a list of decision
     variables which are assumed to be positive.

     If the minimum exists, `minimize_lp' returns a list which contains
     the minimum value of the objective function and a list of decision
     variable values for which the minimum is attained. If the problem
     is not bounded, `minimize_lp' returns "Problem not bounded!" and
     if the problem is not feasible, it returns "Ploblem not feasible!".

     The decision variables are not assumed to be nonegative by
     default. If all decision variables are nonegative, set
     `nonegative_lp' to `true'.  If only some of decision variables are
     positive, list them in the optional argument <pos> (note that this
     is more efficient than adding constraints).

     `minimize_lp' uses the simplex algorithm which is implemented in
     maxima `linear_program' function.

     To use this function first load the `simplex' package with
     `load(simplex);'.

     Examples:

          (%i1) minimize_lp(x+y, [3*x+y=0, x+2*y>2]);
                                4       6        2
          (%o1)                [-, [y = -, x = - -]]
                                5       5        5
          (%i2) minimize_lp(x+y, [3*x+y>0, x+2*y>2]), nonegative_lp=true;
          (%o2)                [1, [y = 1, x = 0]]
          (%i3) minimize_lp(x+y, [3*x+y=0, x+2*y>2]), nonegative_lp=true;
          (%o3)                Problem not feasible!
          (%i4) minimize_lp(x+y, [3*x+y>0]);
          (%o4)                Problem not bounded!

     See also: `maximize_lp', `nonegative_lp', `epsilon_lp'.


 -- Option variable: nonegative_lp
     Default value: `false'

     If `nonegative_lp' is true all decision variables to `minimize_lp'
     and `maximize_lp' are assumed to be positive.

     See also: `minimize_lp'.



File: maxima.info,  Node: simplification,  Next: solve_rec,  Prev: simplex,  Up: Top

69 simplification
*****************

* Menu:

* Introduction to simplification::
* Functions and Variables for simplification::


File: maxima.info,  Node: Introduction to simplification,  Next: Functions and Variables for simplification,  Prev: simplification,  Up: simplification

69.1 Introduction to simplification
===================================

The directory `maxima/share/simplification' contains several scripts
which implement simplification rules and functions, and also some
functions not related to simplification.


File: maxima.info,  Node: Functions and Variables for simplification,  Prev: Introduction to simplification,  Up: simplification

69.2 Functions and Variables for simplification
===============================================

69.2.1 Package absimp
---------------------

The `absimp' package contains pattern-matching rules that extend the
built-in simplification rules for the `abs' and `signum' functions.
`absimp' respects relations established with the built-in `assume'
function and by declarations such as  `modedeclare (m, even, n, odd)'
for even or odd integers.

   `absimp' defines `unitramp' and `unitstep' functions in terms of
`abs' and `signum'.

   `load (absimp)' loads this package.  `demo (absimp)' shows a
demonstration of this package.

   Examples:

     (%i1) load (absimp)$
     (%i2) (abs (x))^2;
                                            2
     (%o2)                                 x
     (%i3) diff (abs (x), x);
                                           x
     (%o3)                               ------
                                         abs(x)
     (%i4) cosh (abs (x));
     (%o4)                               cosh(x)

69.2.2 Package facexp
---------------------

The `facexp' package contains several related  functions that provide
the user with the ability to structure expressions by controlled
expansion.   This capability  is especially  useful when  the expression
contains variables that have physical meaning, because it is  often true
that the most economical form  of such an expression can be  obtained by
fully expanding the expression with respect to those variables, and then
factoring their coefficients.  While it is  true that this  procedure is
not difficult to carry out using standard Maxima  functions, additional
fine-tuning may also  be desirable, and  these finishing touches  can be
more  difficult to  apply.

   The  function `facsum'  and its  related forms provide a convenient
means for controlling the structure  of expressions in this way.
Another function, `collectterms', can be used to add  two or more
expressions that have already been simplified to this form, without
resimplifying the whole expression again.  This function may be useful
when the expressions are very large.

   `load (facexp)' loads this package.  `demo (facexp)' shows a
demonstration of this package.

 -- Function: facsum (<expr>, <arg_1>, ..., <arg_n>)
     Returns  a form  of <expr>  which depends  on the arguments
     <arg_1>, ..., <arg_n>.  The arguments can be any form suitable for
     `ratvars', or they can be lists  of such  forms.  If  the
     arguments  are not  lists, then  the form returned is  fully
     expanded with respect  to the arguments,  and the coefficients of
     the arguments are factored.  These  coefficients are free of the
     arguments, except perhaps in a non-rational sense.

     If any of the arguments are  lists, then all such lists are
     combined into  a  single  list,   and  instead  of  calling
     `factor'   on  the coefficients  of  the  arguments,  `facsum'
     calls  itself   on  these coefficients, using  this newly
     constructed  single list as  the new argument list  for this
     recursive  call.  This  process can  be  repeated to arbitrary
     depth by nesting the desired elements in lists.

     It is possible that one may wish to `facsum' with respect  to more
     complicated subexpressions,  such as  `log (x + y)'.  Such
     arguments are also  permissible.   With  no  variable
     specification,  for example `facsum (<expr>)', the  result
     returned  is the same  as that  returned by `ratsimp (<expr>)'.

     Occasionally the user may wish to obtain any of the  above forms
     for expressions which are specified only by their leading
     operators.  For example, one may wish  to `facsum' with respect to
     all  `log''s.  In this situation, one may  include among the
     arguments either  the specific `log''s which are to be treated in
     this way, or  alternatively, either the expression  `operator
     (log)' or `'operator (log)'.   If one  wished to `facsum' the
     expression <expr> with respect to the operators <op_1>, ...,
     <op_n>, one   would  evaluate  `facsum (<expr>, operator (<op_1>,
     ..., <op_n>))'.  The `operator' form may also appear inside list
     arguments.

     In  addition,  the  setting  of  the  switches   `facsum_combine'
     and `nextlayerfactor' may affect the result of `facsum'.

 -- Global variable: nextlayerfactor
     Default value: `false'

     When `nextlayerfactor' is `true', recursive calls  of `facsum' are
     applied  to  the  factors  of  the  factored  form   of  the
     coefficients of the arguments.

     When  `false', `facsum' is applied to each coefficient as a whole
     whenever recusive calls to  `facsum' occur.

     Inclusion   of   the  atom `nextlayerfactor' in  the argument
     list of `facsum'  has the  effect of `nextlayerfactor: true', but
     for the next level of the expression only.  Since
     `nextlayerfactor' is  always bound to  either `true' or  `false',
     it must be presented single-quoted whenever it appears in the
     argument list of `facsum'.

 -- Global variable: facsum_combine
     Default value: `true'

     `facsum_combine' controls the form  of the final result  returned
     by `facsum'  when  its  argument  is  a  quotient  of
     polynomials.   If `facsum_combine' is `false'  then the form will
     be returned as  a fully expanded  sum  as described  above,  but
     if  `true',  then  the expression returned is a ratio of
     polynomials, with each polynomial in the form described above.

     The `true' setting of this switch is useful when one wants to
     `facsum' both  the numerator and  denominator of  a rational
     expression,  but  does not  want  the denominator  to  be
     multiplied through the terms of the numerator.

 -- Function: factorfacsum (<expr>, <arg_1>, ... <arg_n>)
     Returns a  form of <expr>  which is obtained by calling  `facsum'
     on the factors  of <expr> with <arg_1>, ... <arg_n> as arguments.
     If any of the factors of <expr> is raised to a  power, both the
     factor and the exponent will be processed in this way.

 -- Function: collectterms (<arg_1>, ..., <arg_n>)
     If several  expressions  have been simplified  with  `facsum',
     `factorfacsum',  `factenexpand',  `facexpten' or
     `factorfacexpten',  and  they are  to  be added  together,  it
     may be desirable   to  combine   them  using   the   function
     `collecterms'.  `collecterms' can take as arguments  all of the
     arguments that  can be given  to these  other associated
     functions with  the  exception of `nextlayerfactor', which has no
     effect on `collectterms'.  The advantage of `collectterms'  is
     that it  returns a form  similar to  `facsum', but since it is
     adding forms that have already been processed by `facsum', it
     does  not  need  to  repeat  that  effort.   This  capability is
     especially useful when the expressions to be summed are very large.

69.2.3 Package functs
---------------------

 -- Function: rempart (<expr>, <n>)
     Removes part <n> from the expression <expr>.

     If <n> is a list of the form `[<l>, <m>]' then parts <l> thru <m>
     are removed.

     To use this function write first `load(functs)'.

 -- Function: wronskian ([<f_1>, ..., <f_n>], <x>)
     Returns the Wronskian matrix of the functions <f_1>, ..., <f_n> in
     the variable <x>.

     <f_1>, ..., <f_n> may be the names of user-defined functions, or
     expressions in the variable <x>.

     The determinant of the Wronskian matrix is the Wronskian
     determinant of the set of functions.  The functions are linearly
     dependent if this determinant is zero.

     To use this function write first `load(functs)'.

 -- Function: tracematrix (<M>)
     Returns the trace (sum of the diagonal elements) of matrix <M>.

     To use this function write first `load(functs)'.

 -- Function: rational (`z')
     Multiplies numerator and denominator of <z> by the complex
     conjugate of denominator, thus rationalizing the denominator.
     Returns canonical rational expression (CRE) form if given one,
     else returns general form.

     To use this function write first `load(functs)'.

 -- Function: logand (`x',`y')
     Returns logical (bit-wise) "and" of arguments x and y.

     To use this function write first `load(functs)'.

 -- Function: logor (`x',`y')
     Returns logical (bit-wise) "or" of arguments x and y.

     To use this function write first `load(functs)'.

 -- Function: logxor (`x',`y')
     Returns logical (bit-wise) exclusive-or of arguments x and y.

     To use this function write first `load(functs)'.

 -- Function: nonzeroandfreeof (<x>, <expr>)
     Returns `true' if <expr> is nonzero and `freeof (<x>, <expr>)'
     returns `true'.  Returns `false' otherwise.

     To use this function write first `load(functs)'.

 -- Function: linear (<expr>, <x>)
     When <expr> is an expression linear in variable <x>, `linear'
     returns `<a>*<x> + <b>' where <a> is nonzero, and <a> and <b> are
     free of <x>.  Otherwise, `linear' returns <expr>.

     To use this function write first `load(functs)'.

 -- Function: gcdivide (<p>, <q>)
     When `takegcd' is `true', `gcdivide' divides the polynomials <p>
     and <q> by their greatest common divisor and returns the ratio of
     the results.

     When `takegcd' is `false', `gcdivide' returns the ratio `<p>/<q>'.

     To use this function write first `load(functs)'.

 -- Function: arithmetic (<a>, <d>, <n>)
     Returns the <n>-th term of the arithmetic series `<a>, <a> + <d>,
     <a> + 2*<d>, ..., <a> + (<n> - 1)*<d>'.

     To use this function write first `load(functs)'.

 -- Function: geometric (<a>, <r>, <n>)
     Returns the <n>-th term of the geometric series `<a>, <a>*<r>,
     <a>*<r>^2, ..., <a>*<r>^(<n> - 1)'.

     To use this function write first `load(functs)'.

 -- Function: harmonic (<a>, <b>, <c>, <n>)
     Returns the <n>-th term of the harmonic series `<a>/<b>, <a>/(<b>
     + <c>), <a>/(<b> + 2*<c>), ..., <a>/(<b> + (<n> - 1)*<c>)'.

     To use this function write first `load(functs)'.

 -- Function: arithsum (<a>, <d>, <n>)
     Returns the sum of the arithmetic series from 1 to <n>.

     To use this function write first `load(functs)'.

 -- Function: geosum (<a>, <r>, <n>)
     Returns the sum of the geometric series from 1 to <n>.  If <n> is
     infinity (`inf') then a sum is finite only if the absolute value
     of <r> is less than 1.

     To use this function write first `load(functs)'.

 -- Function: gaussprob (<x>)
     Returns the Gaussian probability function `%e^(-<x>^2/2) /
     sqrt(2*%pi)'.

     To use this function write first `load(functs)'.

 -- Function: gd (<x>)
     Returns the Gudermannian function `2 * atan(%e^<x> - %pi/2)'.

     To use this function write first `load(functs)'.

 -- Function: agd (<x>)
     Returns the inverse Gudermannian function `log (tan (%pi/4 +
     x/2)))'.

     To use this function write first `load(functs)'.

 -- Function: vers (<x>)
     Returns the versed sine `1 - cos (x)'.

     To use this function write first `load(functs)'.

 -- Function: covers (<x>)
     Returns the coversed sine `1 - sin (<x>)'.

     To use this function write first `load(functs)'.

 -- Function: exsec (<x>)
     Returns the exsecant `sec (<x>) - 1'.

     To use this function write first `load(functs)'.

 -- Function: hav (<x>)
     Returns the haversine `(1 - cos(x))/2'.

     To use this function write first `load(functs)'.

 -- Function: combination (<n>, <r>)
     Returns the number of combinations of <n> objects taken <r> at a
     time.

     To use this function write first `load(functs)'.

 -- Function: permutation (<n>, <r>)
     Returns the number of permutations of <r> objects selected from a
     set of <n> objects.

     To use this function write first `load(functs)'.

69.2.4 Package ineq
-------------------

The `ineq' package contains simplification rules for inequalities.

   Example session:

     (%i1) load(ineq)$
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     Warning: Putting rules on '+' or '*' is inefficient, and may not work.
     (%i2) a>=4;  /* a sample inequality */
     (%o2)                               a >= 4
     (%i3) (b>c)+%; /* add a second, strict inequality */
     (%o3)                            b + a > c + 4
     (%i4) 7*(x<y); /* multiply by a positive number */
     (%o4)                              7 x < 7 y
     (%i5) -2*(x>=3*z); /* multiply by a negative number */
     (%o5)                           - 2 x <= - 6 z
     (%i6) (1+a^2)*(1/(1+a^2)<=1); /* Maxima knows that 1+a^2 > 0 */
                                             2
     (%o6)                             1 <= a  + 1
     (%i7) assume(x>0)$ x*(2<3); /* assuming x>0 */
     (%o7)                              2 x < 3 x
     (%i8) a>=b; /* another inequality */
     (%o8)                               a >= b
     (%i9) 3+%; /* add something */
     (%o9)                           a + 3 >= b + 3
     (%i10) %-3; /* subtract it out */
     (%o10)                              a >= b
     (%i11) a>=c-b; /* yet another inequality */
     (%o11)                            a >= c - b
     (%i12) b+%; /* add b to both sides */
     (%o12)                            b + a >= c
     (%i13) %-c; /* subtract c from both sides */
     (%o13)                         - c + b + a >= 0
     (%i14) -%;  /* multiply by -1 */
     (%o14)                          c - b - a <= 0
     (%i15) (z-1)^2>-2*z; /* determining truth of assertion */
                                           2
     (%o15)                         (z - 1)  > - 2 z
     (%i16) expand(%)+2*z; /* expand this and add 2*z to both sides */
                                        2
     (%o16)                            z  + 1 > 0
     (%i17) %,pred;
     (%o17)                               true

   Be careful about using parentheses around the inequalities: when the
user types in `(A > B) + (C = 5)' the result is `A + C > B + 5', but `A
> B + C = 5' is a syntax error, and `(A > B + C) = 5' is something else
entirely.

   Do `disprule (all)' to see a complete listing of the rule
definitions.

   The user will be queried if Maxima is unable to decide the sign of a
quantity multiplying an inequality.

   The most common mis-feature is illustrated by:

     (%i1) eq: a > b;
     (%o1)                              a > b
     (%i2) 2*eq;
     (%o2)                            2 (a > b)
     (%i3) % - eq;
     (%o3)                              a > b

   Another problem is 0 times an inequality; the default to have this
turn into 0 has been left alone. However, if you type
`X*<some_inequality>' and Maxima asks about the sign of `X' and you
respond `zero' (or `z'), the program returns `X*<some_inequality>' and
not use the information that `X' is 0. You should do `ev (%, x: 0)' in
such a case, as the database will only be used for comparison purposes
in decisions, and not for the purpose of evaluating `X'.

   The user may note a slower response when this package is loaded, as
the simplifier is forced to examine more rules than without the
package, so you might wish to remove the rules after making use of
them. Do `kill (rules)' to eliminate all of the rules (including any
that you might have defined); or you may be more selective by killing
only some of them; or use `remrule' on a specific rule.

   Note that if you load this package after defining your own rules you
will clobber your rules that have the same name. The rules in this
package are: `*rule1', ..., `*rule8', `+rule1', ..., `+rule18', and you
must enclose the rulename in quotes to refer to it, as in `remrule
("+", "+rule1")' to specifically remove the first rule on `"+"' or
`disprule ("*rule2")' to display the definition of the second
multiplicative rule.

69.2.5 Package rducon
---------------------

 -- Function: reduce_consts (<expr>)
     Replaces constant subexpressions of <expr> with constructed
     constant atoms, saving the definition of all these constructed
     constants in the list of equations `const_eqns', and returning the
     modified <expr>.  Those parts of <expr> are constant which return
     `true' when operated on by the function `constantp'.  Hence,
     before invoking `reduce_consts', one should do

          declare ([<objects to be given the constant property>], constant)$

     to set up a database of the constant quantities occurring in your
     expressions.

     If you are planning to generate Fortran output after these symbolic
     calculations, one of the first code sections should be the
     calculation of all constants.  To generate this code segment, do

          map ('fortran, const_eqns)$

     Variables besides `const_eqns' which affect `reduce_consts' are:

     `const_prefix' (default value: `xx') is the string of characters
     used to prefix all symbols generated by `reduce_consts' to
     represent constant subexpressions.

     `const_counter' (default value: 1) is the integer index used to
     generate unique symbols to represent each constant subexpression
     found by `reduce_consts'.

     `load (rducon)' loads this function.  `demo (rducon)' shows a
     demonstration of this function.

69.2.6 Package scifac
---------------------

 -- Function: gcfac (<expr>)
     `gcfac' is a factoring function that attempts to apply the same
     heuristics which scientists apply in trying to make expressions
     simpler.  `gcfac' is limited to monomial-type factoring.  For a
     sum, `gcfac' does the following:

       1. Factors over the integers.

       2. Factors out the largest powers of terms occurring as
          coefficients, regardless of the complexity of the terms.

       3. Uses (1) and (2) in factoring adjacent pairs of terms.

       4. Repeatedly and recursively applies these techniques until the
          expression no longer changes.

     Item (3) does not necessarily do an optimal job of pairwise
     factoring because of the combinatorially-difficult nature of
     finding which of all possible rearrangements of the pairs yields
     the most compact pair-factored result.

     `load (scifac)' loads this function.  `demo (scifac)' shows a
     demonstration of this function.

69.2.7 Package sqdnst
---------------------

 -- Function: sqrtdenest (<expr>)
     Denests `sqrt' of simple, numerical, binomial surds, where
     possible.  E.g.

          (%i1) load (sqdnst)$
          (%i2) sqrt(sqrt(3)/2+1)/sqrt(11*sqrt(2)-12);
                                              sqrt(3)
                                         sqrt(------- + 1)
                                                 2
          (%o2)                        ---------------------
                                       sqrt(11 sqrt(2) - 12)
          (%i3) sqrtdenest(%);
                                            sqrt(3)   1
                                            ------- + -
                                               2      2
          (%o3)                            -------------
                                              1/4    3/4
                                           3 2    - 2

     Sometimes it helps to apply `sqrtdenest' more than once, on such as
     `(19601-13860 sqrt(2))^(7/4)'.

     `load (sqdnst)' loads this function.


File: maxima.info,  Node: solve_rec,  Next: stats,  Prev: simplification,  Up: Top

70 solve_rec
************

* Menu:

* Introduction to solve_rec::
* Functions and Variables for solve_rec::


File: maxima.info,  Node: Introduction to solve_rec,  Next: Functions and Variables for solve_rec,  Prev: solve_rec,  Up: solve_rec

70.1 Introduction to solve_rec
==============================

`solve_rec' is a package for solving linear recurrences with polynomial
coefficients.

   A demo is available with `demo(solve_rec);'.

   Example:

     (%i1) load("solve_rec")$
     (%i2) solve_rec((n+4)*s[n+2] + s[n+1] - (n+1)*s[n], s[n]);
                                         n
                      %k  (2 n + 3) (- 1)          %k
                        1                            2
     (%o2)       s  = -------------------- + ---------------
                  n     (n + 1) (n + 2)      (n + 1) (n + 2)


File: maxima.info,  Node: Functions and Variables for solve_rec,  Prev: Introduction to solve_rec,  Up: solve_rec

70.2 Functions and Variables for solve_rec
==========================================

 -- Function: reduce_order (<rec>, <sol>, <var>)
     Reduces the order of linear recurrence <rec> when a particular
     solution <sol> is known. The reduced reccurence can be used to get
     other solutions.

     Example:

          (%i3) rec: x[n+2] = x[n+1] + x[n]/n;
                                                x
                                                 n
          (%o3)               x      = x      + --
                               n + 2    n + 1   n
          (%i4) solve_rec(rec, x[n]);
          WARNING: found some hypergeometrical solutions!
          (%o4)                    x  = %k  n
                                    n     1
          (%i5) reduce_order(rec, n, x[n]);
          (%t5)                    x  = n %z
                                    n       n

                                     n - 1
                                     ====
                                     \
          (%t6)                %z  =  >     %u
                                 n   /        %j
                                     ====
                                     %j = 0

          (%o6)             (- n - 2) %u     - %u
                                        n + 1     n
          (%i6) solve_rec((n+2)*%u[n+1] + %u[n], %u[n]);
                                               n
                                      %k  (- 1)
                                        1
          (%o6)                 %u  = ----------
                                  n    (n + 1)!

          So the general solution is

                       n - 1
                       ====        n
                       \      (- 1)
                 %k  n  >    -------- + %k  n
                   2   /     (n + 1)!     1
                       ====
                       n = 0


 -- Option variable: simplify_products
     Default value: `true'

     If `simplify_products' is `true', `solve_rec' will try to simplify
     products in result.

     See also: `solve_rec'.


 -- Function: simplify_sum (<expr>)
     Tries to simplify all sums appearing in <expr> to a closed form.

     `simplify_sum' uses Gosper and Zeilberger algorithms to simplify
     sums.

     To use this function first load the `simplify_sum' package with
     `load(simplify_sum)'.

     Example:

          (%i1) load("simplify_sum")$
          (%i2) sum(binom(n+k,k)/2^k, k, 0, n)
                                       + sum(binom(2*n, 2*k), k, 0, n);
                   n                            n
                  ====                         ====
                  \      binomial(n + k, k)    \
          (%o2)    >     ------------------ +   >    binomial(2 n, 2 k)
                  /               k            /
                  ====           2             ====
                  k = 0                        k = 0
          (%i3) simplify_sum(%);
                                         n
                                        4     n
          (%o3)                         -- + 2
                                        2


 -- Function: solve_rec (<eqn>, <var>, [<init>])
     Solves for hypergeometrical solutions to linear recurrence <eqn>
     with polynomials coefficient in variable <var>. Optional arguments
     <init> are initial conditions.

     `solve_rec' can solve linear recurrences with constant
     coefficients, finds hypergeometrical solutions to homogeneous
     linear recurrences with polynomial coefficients, rational
     solutions to linear recurrences with polynomial coefficients and
     can solve Ricatti type recurrences.

     Note that the running time of the algorithm used to find
     hypergeometrical solutions is exponential in the degree of the
     leading and trailing coefficient.

     To use this function first load the `solve_rec' package with
     `load(solve_rec);'.

     Example of linear recurrence with constant coefficients:

          (%i2) solve_rec(a[n]=a[n-1]+a[n-2]+n/2^n, a[n]);
                                  n          n
                     (sqrt(5) - 1)  %k  (- 1)
                                      1           n
          (%o2) a  = ------------------------- - ----
                 n               n                  n
                                2                5 2
                                                          n
                                             (sqrt(5) + 1)  %k
                                                              2    2
                                           + ------------------ - ----
                                                      n              n
                                                     2            5 2

     Example of linear recurrence with polynomial coefficients:

          (%i7) 2*x*(x+1)*y[x] - (x^2+3*x-2)*y[x+1] + (x-1)*y[x+2];
                                   2
          (%o7) (x - 1) y      - (x  + 3 x - 2) y      + 2 x (x + 1) y
                         x + 2                   x + 1                x
          (%i8) solve_rec(%, y[x], y[1]=1, y[3]=3);
                                        x
                                     3 2    x!
          (%o9)                 y  = ---- - --
                                 x    4     2

     Example of Ricatti type recurrence:

          (%i2) x*y[x+1]*y[x] - y[x+1]/(x+2) + y[x]/(x-1) = 0;
                                      y         y
                                       x + 1     x
          (%o2)         x y  y      - ------ + ----- = 0
                           x  x + 1   x + 2    x - 1
          (%i3) solve_rec(%, y[x], y[3]=5)$
          (%i4) ratsimp(minfactorial(factcomb(%)));
                                             3
                                         30 x  - 30 x
          (%o4) y  = - -------------------------------------------------
                 x        6      5       4       3       2
                       5 x  - 3 x  - 25 x  + 15 x  + 20 x  - 12 x - 1584

     See also: `solve_rec_rat', `simplify_products', and
     `product_use_gamma'.


 -- Function: solve_rec_rat (<eqn>, <var>, [<init>])
     Solves for rational solutions to linear recurrences. See solve_rec
     for description of arguments.

     To use this function first load the `solve_rec' package with
     `load(solve_rec);'.

     Example:

          (%i1) (x+4)*a[x+3] + (x+3)*a[x+2] - x*a[x+1] + (x^2-1)*a[x];
          (%o1)  (x + 4) a      + (x + 3) a      - x a
                          x + 3            x + 2      x + 1
                                                             2
                                                         + (x  - 1) a
                                                                      x
          (%i2) solve_rec_rat(% = (x+2)/(x+1), a[x]);
                                 1
          (%o2)      a  = ---------------
                      x   (x - 1) (x + 1)

     See also: `solve_rec'.


 -- Option variable: product_use_gamma
     Default value: `true'

     When simplifying products, `solve_rec' introduces gamma function
     into the expression if `product_use_gamma' is `true'.

     See also: `simplify_products', `solve_rec'.


 -- Function: summand_to_rec (<summand>, <k>, <n>)
 -- Function: summand_to_rec (<summand>, [<k>, <lo>, <hi>], <n>)
     Returns the recurrence sattisfied by the sum

               hi
              ====
              \
               >     summand
              /
              ====
            k = lo

     where summand is hypergeometrical in <k> and <n>. If <lo> and <hi>
     are omited, they are assumed to be `lo = -inf' and `hi = inf'.

     To use this function first load the `simplify_sum' package with
     `load(simplify_sum)'.

     Example:

          (%i1) load("simplify_sum")$
          (%i2) summand: binom(n,k);
          (%o2)                           binomial(n, k)
          (%i3) summand_to_rec(summand,k,n);
          (%o3)                      2 sm  - sm      = 0
                                         n     n + 1
          (%i7) summand: binom(n, k)/(k+1);
                                          binomial(n, k)
          (%o7)                           --------------
                                              k + 1
          (%i8) summand_to_rec(summand, [k, 0, n], n);
          (%o8)               2 (n + 1) sm  - (n + 2) sm      = - 1
                                          n             n + 1



File: maxima.info,  Node: stats,  Next: stirling,  Prev: solve_rec,  Up: Top

71 stats
********

* Menu:

* Introduction to stats::
* Functions and Variables for inference_result::
* Functions and Variables for stats::
* Functions and Variables for special distributions::


File: maxima.info,  Node: Introduction to stats,  Next: Functions and Variables for inference_result,  Prev: Top,  Up: Top

71.1 Introduction to stats
==========================

Package `stats' contains a set of classical statistical inference and
hypothesis testing procedures.

   All these functions return an `inference_result' Maxima object which
contains the necessary results for population inferences and decision
making.

   Global variable `stats_numer' controls whether results are given in
floating point or symbolic and rational format; its default value is
`true' and results are returned in floating point format.

   Package `descriptive' contains some utilities to manipulate data
structures (lists and matrices); for example, to extract subsamples. It
also contains some examples on how to use package `numericalio' to read
data from plain text files. See `descriptive' and `numericalio' for
more details.

   Package `stats' loads packages `descriptive', `distrib' and
`inference_result'.

   For comments, bugs or suggestions, please contact the author at

   <'mario AT edu DOT xunta DOT es'>.


File: maxima.info,  Node: Functions and Variables for inference_result,  Next: Functions and Variables for stats,  Prev: Introduction to stats,  Up: Top

71.2 Functions and Variables for inference_result
=================================================

 -- Function: inference_result (<title>, <values>, <numbers>)
     Constructs an `inference_result' object of the type returned by the
     stats functions. Argument <title> is a Maxima string with the name
     of the procedure; <values> is a list with elements of the form
     `symbol = value' and <numbers> is a list with positive integer
     numbers ranging from one to `length(<values>)', indicating which
     values will be shown by default.

     Example:

     This is a simple example showing results concerning a rectangle.
     The title of this object is the string `"Rectangle"', it stores
     five results, named `'base', `'height', `'diagonal', `'area', and
     `'perimeter', but only the first, second, fifth, and fourth will
     be displayed. The `'diagonal' is stored in this object, but it is
     not displayed; to access its value, make use of function
     `take_inference'.

          (%i1) load(inference_result)$
          (%i2) b: 3$ h: 2$
          (%i3) inference_result("Rectangle",
                                  ['base=b,
                                   'height=h,
                                   'diagonal=sqrt(b^2+h^2),
                                   'area=b*h,
                                   'perimeter=2*(b+h)],
                                  [1,2,5,4] );
                                  |   Rectangle
                                  |
                                  |    base = 3
                                  |
          (%o3)                   |   height = 2
                                  |
                                  | perimeter = 10
                                  |
                                  |    area = 6
          (%i4) take_inference('diagonal,%);
          (%o4)                        sqrt(13)

     See also `take_inference'.

 -- Function: inferencep (<obj>)
     Returns `true' or `false', depending on whether <obj> is an
     `inference_result' object or not.


 -- Function: items_inference (<obj>)
     Returns a list with the names of the items stored in <obj>, which
     must be an `inference_result' object.

     Example:

     The `inference_result' object stores two values, named `'pi' and
     `'e', but only the second is displayed. The `items_inference'
     function returns the names of all items, no matter they are
     displayed or not.

          (%i1) load(inference_result)$
          (%i2) inference_result("Hi", ['pi=%pi,'e=%e],[2]);
                                      |   Hi
          (%o2)                       |
                                      | e = %e
          (%i3) items_inference(%);
          (%o3)                        [pi, e]

 -- Function: take_inference (<n>, <obj>)
 -- Function: take_inference (<name>, <obj>)
 -- Function: take_inference (<list>, <obj>)
     Returns the <n>-th value stored in <obj> if <n> is a positive
     integer, or the item named <name> if this is the name of an item.
     If the first argument is a list of numbers and/or symbols,
     function `take_inference' returns a list with the corresponding
     results.

     Example:

     Given an `inference_result' object, function `take_inference' is
     called in order to extract some information stored in it.

          (%i1) load(inference_result)$
          (%i2) b: 3$ h: 2$
          (%i3) sol: inference_result("Rectangle",
                                      ['base=b,
                                       'height=h,
                                       'diagonal=sqrt(b^2+h^2),
                                       'area=b*h,
                                       'perimeter=2*(b+h)],
                                      [1,2,5,4] );
                                  |   Rectangle
                                  |
                                  |    base = 3
                                  |
          (%o3)                   |   height = 2
                                  |
                                  | perimeter = 10
                                  |
                                  |    area = 6
          (%i4) take_inference('base,sol);
          (%o4)                           3
          (%i5) take_inference(5,sol);
          (%o5)                          10
          (%i6) take_inference([1,'diagonal],sol);
          (%o6)                     [3, sqrt(13)]
          (%i7) take_inference(items_inference(sol),sol);
          (%o7)                [3, 2, sqrt(13), 6, 10]

     See also `inference_result' and `take_inference'.


File: maxima.info,  Node: Functions and Variables for stats,  Next: Functions and Variables for special distributions,  Prev: Functions and Variables for inference_result,  Up: Top

71.3 Functions and Variables for stats
======================================

 -- Option variable: stats_numer
     Default value: `true'

     If `stats_numer' is `true', inference statistical functions return
     their results in floating point numbers. If it is `false', results
     are given in symbolic and rational format.


 -- Function: test_mean (<x>)
 -- Function: test_mean (<x>, <option_1>, <option_2>, ...)
     This is the mean <t>-test. Argument <x> is a list or a column
     matrix containing a one dimensional sample. It also performs an
     asymptotic test based on the Central Limit Theorem if option
     `'asymptotic' is `true'.

     Options:

        * `'mean', default `0', is the mean value to be checked.

        * `'alternative', default `'twosided', is the alternative
          hypothesis; valid values are: `'twosided', `'greater' and
          `'less'.

        * `'dev', default `'unknown', this is the value of the standard
          deviation when it is known; valid values are: `'unknown' or a
          positive expression.

        * `'conflevel', default `95/100', confidence level for the
          confidence interval; it must be an expression which takes a
          value in (0,1).

        * `'asymptotic', default `false', indicates whether it performs
          an exact <t>-test or an asymptotic one based on the Central
          Limit Theorem; valid values are `true' and `false'.


     The output of function `test_mean' is an `inference_result' Maxima
     object showing the following results:

       1. `'mean_estimate': the sample mean.

       2. `'conf_level': confidence level selected by the user.

       3. `'conf_interval': confidence interval for the population mean.

       4. `'method': inference procedure.

       5. `'hypotheses': null and alternative hypotheses to be tested.

       6. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

       7. `'distribution': distribution of the sample statistic,
          together with its parameter(s).

       8. `'p_value': p-value of the test.


     Examples:

     Performs an exact <t>-test with unknown variance. The null
     hypothesis is H_0: mean=50 against the one sided alternative H_1:
     mean<50; according to the results, the p-value is too great, there
     are no evidence for rejecting H_0.

          (%i1) load("stats")$
          (%i2) data: [78,64,35,45,45,75,43,74,42,42]$
          (%i3) test_mean(data,'conflevel=0.9,'alternative='less,'mean=50);
                    |                 MEAN TEST
                    |
                    |            mean_estimate = 54.3
                    |
                    |              conf_level = 0.9
                    |
                    | conf_interval = [minf, 61.51314273502712]
                    |
          (%o3)     |  method = Exact t-test. Unknown variance.
                    |
                    | hypotheses = H0: mean = 50 , H1: mean < 50
                    |
                    |       statistic = .8244705235071678
                    |
                    |       distribution = [student_t, 9]
                    |
                    |        p_value = .7845100411786889

     This time Maxima performs an asymptotic test, based on the Central
     Limit Theorem.  The null hypothesis is H_0: equal(mean, 50)
     against the two sided alternative H_1: not equal(mean, 50);
     according to the results, the p-value is very small, H_0 should be
     rejected in favor of the alternative H_1. Note that, as indicated
     by the `Method' component, this procedure should be applied to
     large samples.

          (%i1) load("stats")$
          (%i2) test_mean([36,118,52,87,35,256,56,178,57,57,89,34,25,98,35,
                        98,41,45,198,54,79,63,35,45,44,75,42,75,45,45,
                        45,51,123,54,151],
                        'asymptotic=true,'mean=50);
                    |                       MEAN TEST
                    |
                    |           mean_estimate = 74.88571428571429
                    |
                    |                   conf_level = 0.95
                    |
                    | conf_interval = [57.72848600856194, 92.04294256286663]
                    |
          (%o2)     |    method = Large sample z-test. Unknown variance.
                    |
                    |       hypotheses = H0: mean = 50 , H1: mean # 50
                    |
                    |             statistic = 2.842831192874313
                    |
                    |             distribution = [normal, 0, 1]
                    |
                    |             p_value = .004471474652002261


 -- Function: test_means_difference (<x1>, <x2>)
 -- Function: test_means_difference (<x1>, <x2>, <option_1>,
          <option_2>, ...)
     This is the difference of means <t>-test for two samples.
     Arguments <x1> and <x2> are lists or column matrices containing
     two independent samples. In case of different unknown variances
     (see options `'dev1', `'dev2' and `'varequal' bellow), the degrees
     of freedom are computed by means of the Welch approximation.  It
     also performs an asymptotic test based on the Central Limit
     Theorem if option `'asymptotic' is set to `true'.

     Options:

        *        * `'alternative', default `'twosided', is the alternative
          hypothesis; valid values are: `'twosided', `'greater' and
          `'less'.

        * `'dev1', default `'unknown', this is the value of the
          standard deviation of the <x1> sample when it is known; valid
          values are: `'unknown' or a positive expression.

        * `'dev2', default `'unknown', this is the value of the
          standard deviation of the <x2> sample when it is known; valid
          values are: `'unknown' or a positive expression.

        * `'varequal', default `false', whether variances should be
          considered to be equal or not; this option takes effect only
          when `'dev1' and/or `'dev2' are  `'unknown'.

        * `'conflevel', default `95/100', confidence level for the
          confidence interval; it must be an expression which takes a
          value in (0,1).

        * `'asymptotic', default `false', indicates whether it performs
          an exact <t>-test or an asymptotic one based on the Central
          Limit Theorem; valid values are `true' and `false'.


     The output of function `test_means_difference' is an
     `inference_result' Maxima object showing the following results:

       1. `'diff_estimate': the difference of means estimate.

       2. `'conf_level': confidence level selected by the user.

       3. `'conf_interval': confidence interval for the difference of
          means.

       4. `'method': inference procedure.

       5. `'hypotheses': null and alternative hypotheses to be tested.

       6. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

       7. `'distribution': distribution of the sample statistic,
          together with its parameter(s).

       8. `'p_value': p-value of the test.


     Examples:

     The equality of means is tested with two small samples <x> and <y>,
     against the alternative H_1: m_1>m_2, being m_1 and m_2 the
     populations means; variances are unknown and supposed to be
     different.

          (%i1) load("stats")$
          (%i2) x: [20.4,62.5,61.3,44.2,11.1,23.7]$
          (%i3) y: [1.2,6.9,38.7,20.4,17.2]$
          (%i4) test_means_difference(x,y,'alternative='greater);
                      |              DIFFERENCE OF MEANS TEST
                      |
                      |         diff_estimate = 20.31999999999999
                      |
                      |                 conf_level = 0.95
                      |
                      |    conf_interval = [- .04597417812882298, inf]
                      |
          (%o4)       |        method = Exact t-test. Welch approx.
                      |
                      | hypotheses = H0: mean1 = mean2 , H1: mean1 > mean2
                      |
                      |           statistic = 1.838004300728477
                      |
                      |    distribution = [student_t, 8.62758740184604]
                      |
                      |            p_value = .05032746527991905

     The same test as before, but now variances are supposed to be
     equal.

          (%i1) load("stats")$
          (%i2) x: [20.4,62.5,61.3,44.2,11.1,23.7]$
          (%i3) y: matrix([1.2],[6.9],[38.7],[20.4],[17.2])$
          (%i4) test_means_difference(x,y,'alternative='greater,
                                                           'varequal=true);
                      |              DIFFERENCE OF MEANS TEST
                      |
                      |         diff_estimate = 20.31999999999999
                      |
                      |                 conf_level = 0.95
                      |
                      |     conf_interval = [- .7722627696897568, inf]
                      |
          (%o4)       |   method = Exact t-test. Unknown equal variances
                      |
                      | hypotheses = H0: mean1 = mean2 , H1: mean1 > mean2
                      |
                      |           statistic = 1.765996124515009
                      |
                      |           distribution = [student_t, 9]
                      |
                      |            p_value = .05560320992529344


 -- Function: test_variance (<x>)
 -- Function: test_variance (<x>, <option_1>, <option_2>, ...)
     This is the variance <chi^2>-test. Argument <x> is a list or a
     column matrix containing a one dimensional sample taken from a
     normal population.

     Options:

        * `'mean', default `'unknown', is the population's mean, when
          it is known.

        * `'alternative', default `'twosided', is the alternative
          hypothesis; valid values are: `'twosided', `'greater' and
          `'less'.

        * `'variance', default `1', this is the variance value
          (positive) to be checked.

        * `'conflevel', default `95/100', confidence level for the
          confidence interval; it must be an expression which takes a
          value in (0,1).


     The output of function `test_variance' is an `inference_result'
     Maxima object showing the following results:

       1. `'var_estimate': the sample variance.

       2. `'conf_level': confidence level selected by the user.

       3. `'conf_interval': confidence interval for the population
          variance.

       4. `'method': inference procedure.

       5. `'hypotheses': null and alternative hypotheses to be tested.

       6. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

       7. `'distribution': distribution of the sample statistic,
          together with its parameter.

       8. `'p_value': p-value of the test.


     Examples:

     It is tested whether the variance of a population with unknown mean
     is equal to or greater than 200.

          (%i1) load("stats")$
          (%i2) x: [203,229,215,220,223,233,208,228,209]$
          (%i3) test_variance(x,'alternative='greater,'variance=200);
                       |                  VARIANCE TEST
                       |
                       |              var_estimate = 110.75
                       |
                       |                conf_level = 0.95
                       |
                       |     conf_interval = [57.13433376937479, inf]
                       |
          (%o3)        | method = Variance Chi-square test. Unknown mean.
                       |
                       |    hypotheses = H0: var = 200 , H1: var > 200
                       |
                       |                 statistic = 4.43
                       |
                       |             distribution = [chi2, 8]
                       |
                       |           p_value = .8163948512777689


 -- Function: test_variance_ratio (<x1>, <x2>)
 -- Function: test_variance_ratio (<x1>, <x2>, <option_1>, <option_2>,
          ...)
     This is the variance ratio <F>-test for two normal populations.
     Arguments <x1> and <x2> are lists or column matrices containing
     two independent samples.

     Options:

        * `'alternative', default `'twosided', is the alternative
          hypothesis; valid values are: `'twosided', `'greater' and
          `'less'.

        * `'mean1', default `'unknown', when it is known, this is the
          mean of the population from which <x1> was taken.

        * `'mean2', default `'unknown', when it is known, this is the
          mean of the population from which <x2> was taken.

        * `'conflevel', default `95/100', confidence level for the
          confidence interval of the ratio; it must be an expression
          which takes a value in (0,1).


     The output of function `test_variance_ratio' is an
     `inference_result' Maxima object showing the following results:

       1. `'ratio_estimate': the sample variance ratio.

       2. `'conf_level': confidence level selected by the user.

       3. `'conf_interval': confidence interval for the variance ratio.

       4. `'method': inference procedure.

       5. `'hypotheses': null and alternative hypotheses to be tested.

       6. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

       7. `'distribution': distribution of the sample statistic,
          together with its parameters.

       8. `'p_value': p-value of the test.


     Examples:

     The equality of the variances of two normal populations is checked
     against the alternative that the first is greater than the second.

          (%i1) load("stats")$
          (%i2) x: [20.4,62.5,61.3,44.2,11.1,23.7]$
          (%i3) y: [1.2,6.9,38.7,20.4,17.2]$
          (%i4) test_variance_ratio(x,y,'alternative='greater);
                        |              VARIANCE RATIO TEST
                        |
                        |       ratio_estimate = 2.316933391522034
                        |
                        |               conf_level = 0.95
                        |
                        |    conf_interval = [.3703504689507268, inf]
                        |
          (%o4)         | method = Variance ratio F-test. Unknown means.
                        |
                        | hypotheses = H0: var1 = var2 , H1: var1 > var2
                        |
                        |         statistic = 2.316933391522034
                        |
                        |            distribution = [f, 5, 4]
                        |
                        |          p_value = .2179269692254457


 -- Function: test_sign (<x>)
 -- Function: test_sign (<x>, <option_1>, <option_2>, ...)
     This is the non parametric sign test for the median of a
     continuous population.  Argument <x> is a list or a column matrix
     containing a one dimensional sample.

     Options:

        * `'alternative', default `'twosided', is the alternative
          hypothesis; valid values are: `'twosided', `'greater' and
          `'less'.

        * `'median', default `0', is the median value to be checked.


     The output of function `test_sign' is an `inference_result' Maxima
     object showing the following results:

       1. `'med_estimate': the sample median.

       2. `'method': inference procedure.

       3. `'hypotheses': null and alternative hypotheses to be tested.

       4. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

       5. `'distribution': distribution of the sample statistic,
          together with its parameter(s).

       6. `'p_value': p-value of the test.


     Examples:

     Checks whether the population from which the sample was taken has
     median 6, against the alternative H_1: median > 6.

          (%i1) load("stats")$
          (%i2) x: [2,0.1,7,1.8,4,2.3,5.6,7.4,5.1,6.1,6]$
          (%i3) test_sign(x,'median=6,'alternative='greater);
                         |                  SIGN TEST
                         |
                         |              med_estimate = 5.1
                         |
                         |      method = Non parametric sign test.
                         |
          (%o3)          | hypotheses = H0: median = 6 , H1: median > 6
                         |
                         |                statistic = 7
                         |
                         |      distribution = [binomial, 10, 0.5]
                         |
                         |         p_value = .05468749999999989


 -- Function: test_signed_rank (<x>)
 -- Function: test_signed_rank (<x>, <option_1>, <option_2>, ...)
     This is the Wilcoxon signed rank test to make inferences about the
     median of a continuous population. Argument <x> is a list or a
     column matrix containing a one dimensional sample. Performs normal
     approximation if the sample size is greater than 20, or if there
     are zeroes or ties.

     See also `pdf_rank_test' and `cdf_rank_test'.

     Options:

        * `'median', default `0', is the median value to be checked.

        * `'alternative', default `'twosided', is the alternative
          hypothesis; valid values are: `'twosided', `'greater' and
          `'less'.


     The output of function `test_signed_rank' is an `inference_result'
     Maxima object with the following results:

       1. `'med_estimate': the sample median.

       2. `'method': inference procedure.

       3. `'hypotheses': null and alternative hypotheses to be tested.

       4. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

       5. `'distribution': distribution of the sample statistic,
          together with its parameter(s).

       6. `'p_value': p-value of the test.


     Examples:

     Checks the null hypothesis H_0: median = 15 against the
     alternative H_1: median > 15. This is an exact test, since there
     are no ties.

          (%i1) load("stats")$
          (%i2) x: [17.1,15.9,13.7,13.4,15.5,17.6]$
          (%i3) test_signed_rank(x,median=15,alternative=greater);
                           |             SIGNED RANK TEST
                           |
                           |           med_estimate = 15.7
                           |
                           |           method = Exact test
                           |
          (%o3)            | hypotheses = H0: med = 15 , H1: med > 15
                           |
                           |              statistic = 14
                           |
                           |     distribution = [signed_rank, 6]
                           |
                           |            p_value = 0.28125

     Checks the null hypothesis H_0: equal(median, 2.5) against the
     alternative H_1: not equal(median, 2.5). This is an approximated
     test, since there are ties.

          (%i1) load("stats")$
          (%i2) y:[1.9,2.3,2.6,1.9,1.6,3.3,4.2,4,2.4,2.9,1.5,3,2.9,4.2,3.1]$
          (%i3) test_signed_rank(y,median=2.5);
                       |                 SIGNED RANK TEST
                       |
                       |                med_estimate = 2.9
                       |
                       |          method = Asymptotic test. Ties
                       |
          (%o3)        |    hypotheses = H0: med = 2.5 , H1: med # 2.5
                       |
                       |                 statistic = 76.5
                       |
                       | distribution = [normal, 60.5, 17.58195097251724]
                       |
                       |           p_value = .3628097734643669


 -- Function: test_rank_sum (<x1>, <x2>)
 -- Function: test_rank_sum (<x1>, <x2>, <option_1>)
     This is the Wilcoxon-Mann-Whitney test for comparing the medians
     of two continuous populations. The first two arguments <x1> and
     <x2> are lists or column matrices with the data of two independent
     samples. Performs normal approximation if any of the sample sizes
     is greater than 10, or if there are ties.

     Option:

        * `'alternative', default `'twosided', is the alternative
          hypothesis; valid values are: `'twosided', `'greater' and
          `'less'.


     The output of function `test_rank_sum' is an `inference_result'
     Maxima object with the following results:

       1. `'method': inference procedure.

       2. `'hypotheses': null and alternative hypotheses to be tested.

       3. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

       4. `'distribution': distribution of the sample statistic,
          together with its parameters.

       5. `'p_value': p-value of the test.


     Examples:

     Checks whether populations have similar medians. Samples sizes are
     small and an exact test is made.

          (%i1) load("stats")$
          (%i2) x:[12,15,17,38,42,10,23,35,28]$
          (%i3) y:[21,18,25,14,52,65,40,43]$
          (%i4) test_rank_sum(x,y);
                        |                 RANK SUM TEST
                        |
                        |              method = Exact test
                        |
                        | hypotheses = H0: med1 = med2 , H1: med1 # med2
          (%o4)         |
                        |                 statistic = 22
                        |
                        |        distribution = [rank_sum, 9, 8]
                        |
                        |          p_value = .1995886466474702

     Now, with greater samples and ties, the procedure makes normal
     approximation. The alternative hypothesis is H_1: median1 <
     median2.

          (%i1) load("stats")$
          (%i2) x: [39,42,35,13,10,23,15,20,17,27]$
          (%i3) y: [20,52,66,19,41,32,44,25,14,39,43,35,19,56,27,15]$
          (%i4) test_rank_sum(x,y,'alternative='less);
                       |                  RANK SUM TEST
                       |
                       |          method = Asymptotic test. Ties
                       |
                       |  hypotheses = H0: med1 = med2 , H1: med1 < med2
          (%o4)        |
                       |                 statistic = 48.5
                       |
                       | distribution = [normal, 79.5, 18.95419580097078]
                       |
                       |           p_value = .05096985666598441


 -- Function: test_normality (<x>)
     Shapiro-Wilk test for normality. Argument <x> is a list of
     numbers, and sample size must be greater than 2 and less or equal
     than 5000, otherwise, function `test_normality' signals an error
     message.

     Reference:

     [1] Algorithm AS R94, Applied Statistics (1995), vol.44, no.4,
     547-551

     The output of function `test_normality' is an `inference_result'
     Maxima object with the following results:

       1. `'statistic': value of the <W> statistic.

       2. `'p_value': p-value under normal assumption.


     Examples:

     Checks for the normality of a population, based on a sample of
     size 9.

          (%i1) load("stats")$
          (%i2) x:[12,15,17,38,42,10,23,35,28]$
          (%i3) test_normality(x);
                                 |      SHAPIRO - WILK TEST
                                 |
          (%o3)                  | statistic = .9251055695162436
                                 |
                                 |  p_value = .4361763918860381


 -- Function: simple_linear_regression (<x>)
 -- Function: simple_linear_regression (<x> <option_1>)
     Simple linear regression, y_i=a+b x_i+e_i, where e_i are N(0,sigma)
     independent random variables. Argument <x> must be a two column
     matrix or a list of pairs.

     Options:

        * `'conflevel', default `95/100', confidence level for the
          confidence interval; it must be an expression which takes a
          value in (0,1).

        * `'regressor', default `'x', name of the independent variable.


     The output of function `simple_linear_regression' is an
     `inference_result' Maxima object with the following results:

       1. `'model': the fitted equation. Useful to make new
          predictions. See examples bellow.

       2. `'means': bivariate mean.

       3. `'variances': variances of both variables.

       4. `'correlation': correlation coefficient.

       5. `'adc': adjusted determination coefficient.

       6. `'a_estimation': estimation of parameter <a>.

       7. `'a_conf_int': confidence interval of parameter <a>.

       8. `'b_estimation': estimation of parameter <b>.

       9. `'b_conf_int': confidence interval of parameter <b>.

      10. `'hypotheses': null and alternative hypotheses about
          parameter <b>.

      11. `'statistic': value of the sample statistic used for testing
          the null hypothesis.

      12. `'distribution': distribution of the sample statistic,
          together with its parameter.

      13. `'p_value': p-value of the test about <b>.

      14. `'v_estimation': unbiased variance estimation, or residual
          variance.

      15. `'v_conf_int': variance confidence interval.

      16. `'cond_mean_conf_int': confidence interval for the
          conditioned mean. See examples bellow.

      17. `'new_pred_conf_int': confidence interval for a new
          prediction. See examples bellow.

      18. `'residuals': list of pairs (prediction, residual), ordered
          with respect to predictions.  This is useful for goodness of
          fit analysis. See examples bellow.


     Only items 1, 4, 14, 9, 10, 11, 12, and 13 above, in this order,
     are shown by default. The rest remain hidden until the user makes
     use of functions `items_inference' and `take_inference'.

     Example:

     Fitting a linear model to a bivariate sample. Input `%i4' plots
     the sample together with the regression line; input `%i5' computes
     `y' given `x=113'; the means and the confidence interval for a new
     prediction when `x=113' are also calculated.

          (%i1) load("stats")$
          (%i2) s:[[125,140.7], [130,155.1], [135,160.3], [140,167.2],
                                                          [145,169.8]]$
          (%i3) z:simple_linear_regression(s,conflevel=0.99);
                     |               SIMPLE LINEAR REGRESSION
                     |
                     |   model = 1.405999999999985 x - 31.18999999999804
                     |
                     |           correlation = .9611685255255155
                     |
                     |           v_estimation = 13.57966666666665
                     |
          (%o3)      | b_conf_int = [.04469633662525263, 2.767303663374718]
                     |
                     |          hypotheses = H0: b = 0 ,H1: b # 0
                     |
                     |            statistic = 6.032686683658114
                     |
                     |            distribution = [student_t, 3]
                     |
                     |             p_value = 0.0038059549413203
          (%i4) plot2d([[discrete, s], take_inference(model,z)],
                  [x,120,150],
                  [gnuplot_curve_styles, ["with points","with lines"]] )$
          (%i5) take_inference(model,z), x=133;
          (%o5)                         155.808
          (%i6) take_inference(means,z);
          (%o6)                     [135.0, 158.62]
          (%i7) take_inference(new_pred_conf_int,z), x=133;
          (%o7)              [132.0728595995113, 179.5431404004887]



File: maxima.info,  Node: Functions and Variables for special distributions,  Prev: Functions and Variables for stats,  Up: Top

71.4 Functions and Variables for special distributions
======================================================

 -- Function: pdf_signed_rank (<x>, <n>)
     Probability density function of the exact distribution of the
     signed rank statistic. Argument <x> is a real number and <n> a
     positive integer.

     See also `test_signed_rank'.

 -- Function: cdf_signed_rank (<x>, <n>)
     Cumulative density function of the exact distribution of the
     signed rank statistic. Argument <x> is a real number and <n> a
     positive integer.

     See also `test_signed_rank'.

 -- Function: pdf_rank_sum (<x>, <n>, <m>)
     Probability density function of the exact distribution of the rank
     sum statistic. Argument <x> is a real number and <n> and <m> are
     both positive integers.

     See also `test_rank_sum'.

 -- Function: cdf_rank_sum (<x>, <n>, <m>)
     Cumulative density function of the exact distribution of the rank
     sum statistic. Argument <x> is a real number and <n> and <m> are
     both positive integers.

     See also `test_rank_sum'.


File: maxima.info,  Node: stirling,  Next: stringproc,  Prev: stats,  Up: Top

72 stirling
***********

* Menu:

* Functions and Variables for stirling::


File: maxima.info,  Node: Functions and Variables for stirling,  Prev: stirling,  Up: stirling

72.1 Functions and Variables for stirling
=========================================

 -- Function: stirling (<z>,<n>)
     Replace `gamma(x)' with the O(1/x^(2n-1)) Stirling formula. when
     <n> isn't a nonnegative integer, signal an error.

     Reference: Abramowitz & Stegun, " Handbook of mathematical
     functions", 6.1.40.

     Examples:
          (%i1) load (stirling)$

          (%i2) stirling(gamma(%alpha+x)/gamma(x),1);
                 1/2 - x             x + %alpha - 1/2
          (%o2) x        (x + %alpha)
                                             1           1
                                      --------------- - ---- - %alpha
                                      12 (x + %alpha)   12 x
                                    %e
          (%i3) taylor(%,x,inf,1);
                              %alpha       2    %alpha
                    %alpha   x       %alpha  - x       %alpha
          (%o3)/T/ x       + -------------------------------- + . . .
                                           2 x
          (%i4) map('factor,%);
                                                 %alpha - 1
                   %alpha   (%alpha - 1) %alpha x
          (%o4)   x       + -------------------------------
                                            2

     The function `stirling' knows the difference between the variable
     'gamma' and the function gamma:

          (%i5) stirling(gamma + gamma(x),0);
                                              x - 1/2   - x
          (%o5)    gamma + sqrt(2) sqrt(%pi) x        %e
          (%i6) stirling(gamma(y) + gamma(x),0);
                                   y - 1/2   - y
          (%o6) sqrt(2) sqrt(%pi) y        %e
                                                        x - 1/2   - x
                                   + sqrt(2) sqrt(%pi) x        %e

     To use this function write first `load("stirling")'.


File: maxima.info,  Node: stringproc,  Next: unit,  Prev: stirling,  Up: Top

73 stringproc
*************

* Menu:

* Introduction to string processing::
* Functions and Variables for input and output::
* Functions and Variables for characters::
* Functions and Variables for strings::


File: maxima.info,  Node: Introduction to string processing,  Next: Functions and Variables for input and output,  Prev: stringproc,  Up: stringproc

73.1 Introduction to string processing
======================================

`stringproc.lisp' enlarges Maximas capabilities of working with strings
and adds some useful functions for file in/output.

   For questions and bugs please mail to van.nek at arcor.de .

   In Maxima a string is easily constructed by typing "text".  Note
that Maxima-strings are no Lisp-strings and vice versa.  `stringp'
tests for Maxima strings, and `lstringp' tests for Lisp strings.  If
for some reasons you have a value, that is a Lisp-string, maybe when
using Maxima-function `sconcat', you can convert via `sunlisp'.

     (%i1) m: "text";
     (%o1)                         text
     (%i2) [stringp(m),lstringp(m)];
     (%o2)                     [true, false]
     (%i3) l: sconcat("text");
     (%o3)                         text
     (%i4) [stringp(l),lstringp(l)];
     (%o4)                     [false, true]
     (%i5) stringp( sunlisp(l) );
     (%o5)                         true

   All functions in `stringproc.lisp', that return strings, return
Maxima-strings.

   Characters are introduced as Maxima-strings of length 1.  Of course,
these are no Lisp-characters.  Tests can be done with `charp'
(respectively `lcharp' and conversion from Lisp to Maxima with
`cunlisp').

     (%i1) c: "e";
     (%o1)                           e
     (%i2) [charp(c),lcharp(c)];
     (%o2)                     [true, false]
     (%i3) supcase(c);
     (%o3)                           E
     (%i4) charp(%);
     (%o4)                         true

   Again, all functions in `stringproc.lisp', that return characters,
return Maxima-characters.  Due to the fact, that the introduced
characters are strings of length 1, you can use a lot of string
functions also for characters.  As seen, `supcase' is one example.

   It is important to know, that the first character in a Maxima-string
is at position 1.  This is designed due to the fact that the first
element in a Maxima-list is at position 1 too.  See definitions of
`charat' and `charlist' for examples.

   In applications string-functions are often used when working with
files.  You will find some useful stream- and print-functions in
`stringproc.lisp'.  The following example shows some of the here
introduced functions at work.

   Example:

   `openw' returns an output stream to a file, `printf' then allows
formatted writing to this file. See `printf' for details.

     (%i1) s: openw("E:/file.txt");
     (%o1)                    #<output stream E:/file.txt>
     (%i2) for n:0 thru 10 do printf( s, "~d ", fib(n) );
     (%o2)                                done
     (%i3) printf( s, "~%~d ~f ~a ~a ~f ~e ~a~%",
                   42,1.234,sqrt(2),%pi,1.0e-2,1.0e-2,1.0b-2 );
     (%o3)                                false
     (%i4) close(s);
     (%o4)                                true

   After closing the stream you can open it again, this time with input
direction.  `readline' returns the entire line as one string. The
`stringproc' package now offers a lot of functions for manipulating
strings. Tokenizing can be done by `split' or `tokens'.

     (%i5) s: openr("E:/file.txt");
     (%o5)                     #<input stream E:/file.txt>
     (%i6) readline(s);
     (%o6)                     0 1 1 2 3 5 8 13 21 34 55
     (%i7) line: readline(s);
     (%o7)               42 1.234 sqrt(2) %pi 0.01 1.0E-2 1.0b-2
     (%i8) list: tokens(line);
     (%o8)           [42, 1.234, sqrt(2), %pi, 0.01, 1.0E-2, 1.0b-2]
     (%i9) map( parsetoken, list );
     (%o9)           [42, 1.234, false, false, 0.01, 0.01, false]

   `parsetoken' only parses integer and float numbers. Parsing symbols
or bigfloats needs `parse_string', which will be autoloaded from
`eval_string.lisp'.

     (%i10) map( parse_string, list );
     (%o10)           [42, 1.234, sqrt(2), %pi, 0.01, 0.01, 1.0b-2]
     (%i11) float(%);
     (%o11) [42.0, 1.234, 1.414213562373095, 3.141592653589793, 0.01,
                                                          0.01, 0.01]
     (%i12) readline(s);
     (%o12)                               false
     (%i13) close(s)$

   `readline' returns `false' when the end of file occurs.


File: maxima.info,  Node: Functions and Variables for input and output,  Next: Functions and Variables for characters,  Prev: Introduction to string processing,  Up: stringproc

73.2 Functions and Variables for input and output
=================================================

Example:

     (%i1) s: openw("E:/file.txt");
     (%o1)                     #<output stream E:/file.txt>
     (%i2) control:
     "~2tAn atom: ~20t~a~%~2tand a list: ~20t~{~r ~}~%~2t\
     and an integer: ~20t~d~%"$
     (%i3) printf( s,control, 'true,[1,2,3],42 )$
     (%o3)                                false
     (%i4) close(s);
     (%o4)                                true
     (%i5) s: openr("E:/file.txt");
     (%o5)                     #<input stream E:/file.txt>
     (%i6) while stringp( tmp:readline(s) ) do print(tmp)$
       An atom:          true
       and a list:       one two three
       and an integer:   42
     (%i7) close(s)$

 -- Function: close (<stream>)
     Closes <stream> and returns `true' if <stream> had been open.


 -- Function: flength (<stream>)
     Returns the number of elements in <stream>.


 -- Function: fposition (<stream>)
 -- Function: fposition (<stream>, <pos>)
     Returns the current position in <stream>, if <pos> is not used.
     If <pos> is used, `fposition' sets the position in <stream>.
     <pos> has to be a positive number, the first element in <stream>
     is in position 1.


 -- Function: freshline ()
 -- Function: freshline (<stream>)
     Writes a new line (to <stream>), if the position is not at the
     beginning of a line.  See also `newline'.

 -- Function: newline ()
 -- Function: newline (<stream>)
     Writes a new line (to <stream>).  See `sprint' for an example of
     using `newline()'.  Note that there are some cases, where
     `newline()' does not work as expected.


 -- Function: opena (<file>)
     Returns an output stream to <file>.  If an existing file is
     opened, `opena' appends elements at the end of file.


 -- Function: openr (<file>)
     Returns an input stream to <file>.  If <file> does not exist, it
     will be created.


 -- Function: openw (<file>)
     Returns an output stream to <file>.  If <file> does not exist, it
     will be created.  If an existing file is opened, `openw'
     destructively modifies <file>.


 -- Function: printf (<dest>, <string>)
 -- Function: printf (<dest>, <string>, <expr_1>, ..., <expr_n>)
     Makes the Common Lisp function FORMAT available in Maxima.  (From
     gcl.info: "format produces formatted output by outputting the
     characters of control-string string and observing that a tilde
     introduces a directive.  The character after the tilde, possibly
     preceded by prefix parameters and modifiers, specifies what kind
     of formatting is desired.  Most directives use one or more
     elements of args to create their output.")

     The following description and the examples may give an idea of
     using `printf'.  See a Lisp reference for more information.

             ~%       new line
             ~&       fresh line
             ~t       tab
             ~$       monetary
             ~d       decimal integer
             ~b       binary integer
             ~o       octal integer
             ~x       hexadecimal integer
             ~br      base-b integer
             ~r       spell an integer
             ~p       plural
             ~f       floating point
             ~e       scientific notation
             ~g       ~f or ~e, depending upon magnitude
             ~a       as printed by Maxima function print
             ~s       strings enclosed in "double quotes"
             ~~       ~
             ~<       justification, ~> terminates
             ~(       case conversion, ~) terminates
             ~[       selection, ~] terminates
             ~{       iteration, ~} terminates

     Please note that there is no format specifier for bigfloats.
     However bigfloats can simply be printed by using the `~a'
     directive.  `~s' prints strings enclosed in "double quotes", you
     can avoid this by using `~a'.  Note that the selection directive
     `~[' is zero-indexed.  Also note that there are some directives,
     which do not work in Maxima.  For example, `~:[' fails.

          (%i1) printf( false, "~a ~a ~4f ~a ~@r",
                        "String",sym,bound,sqrt(12),144), bound = 1.234;
          (%o1)                 String sym 1.23 2*sqrt(3) CXLIV
          (%i2) printf( false,"~{~a ~}",["one",2,"THREE"] );
          (%o2)                          one 2 THREE
          (%i3) printf(true,"~{~{~9,1f ~}~%~}",mat ),
                    mat = args(matrix([1.1,2,3.33],[4,5,6],[7,8.88,9]))$
                1.1       2.0       3.3
                4.0       5.0       6.0
                7.0       8.9       9.0
          (%i4) control: "~:(~r~) bird~p ~[is~;are~] singing."$
          (%i5) printf( false,control, n,n,if n=1 then 0 else 1 ), n=2;
          (%o5)                    Two birds are singing.

     If <dest> is a stream or `true', then `printf' returns `false'.
     Otherwise, `printf' returns a string containing the output.


 -- Function: readline (<stream>)
     Returns a string containing the characters from the current
     position in <stream> up to the end of the line or <false> if the
     end of the file is encountered.


 -- Function: sprint (<expr_1>, ..., <expr_n>)
     Evaluates and displays its arguments one after the other `on a
     line' starting at the leftmost position.  The numbers are printed
     with the '-' right next to the number, and it disregards line
     length. `newline()', which will be autoloaded from
     `stringproc.lisp' might be useful, if you whish to place
     intermediate line breaking.

          (%i1) for n:0 thru 20 do sprint( fib(n) )$
          0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765
          (%i2) for n:0 thru 22 do (
                   sprint(fib(n)), if mod(n,10)=9 then newline() )$
          0 1 1 2 3 5 8 13 21 34
          55 89 144 233 377 610 987 1597 2584 4181
          6765 10946 17711



File: maxima.info,  Node: Functions and Variables for characters,  Next: Functions and Variables for strings,  Prev: Functions and Variables for input and output,  Up: stringproc

73.3 Functions and Variables for characters
===========================================

 -- Function: alphacharp (<char>)
     Returns `true' if <char> is an alphabetic character.


 -- Function: alphanumericp (<char>)
     Returns `true' if <char> is an alphabetic character or a digit.


 -- Function: ascii (<int>)
     Returns the character corresponding to the ASCII number <int>.  (
     -1 < int < 256 )

          (%i1) for n from 0 thru 255 do (
             tmp: ascii(n), if alphacharp(tmp) then sprint(tmp),
                if n=96 then newline() )$
          A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
          a b c d e f g h i j k l m n o p q r s t u v w x y z


 -- Function: cequal (<char_1>, <char_2>)
     Returns `true' if <char_1> and <char_2> are the same.


 -- Function: cequalignore (<char_1>, <char_2>)
     Like `cequal' but ignores case.


 -- Function: cgreaterp (<char_1>, <char_2>)
     Returns `true' if the ASCII number of <char_1> is greater than the
     number of <char_2>.


 -- Function: cgreaterpignore (<char_1>, <char_2>)
     Like `cgreaterp' but ignores case.


 -- Function: charp (<obj>)
     Returns `true' if <obj> is a Maxima-character.  See introduction
     for example.


 -- Function: cint (<char>)
     Returns the ASCII number of <char>.


 -- Function: clessp (<char_1>, <char_2>)
     Returns `true' if the ASCII number of <char_1> is less than the
     number of <char_2>.


 -- Function: clesspignore (<char_1>, <char_2>)
     Like `clessp' but ignores case.


 -- Function: constituent (<char>)
     Returns `true' if <char> is a graphic character and not the space
     character.  A graphic character is a character one can see, plus
     the space character.  (`constituent' is defined by Paul Graham,
     ANSI Common Lisp, 1996, page 67.)

          (%i1) for n from 0 thru 255 do (
          tmp: ascii(n), if constituent(tmp) then sprint(tmp) )$
          ! " #  %  ' ( ) * + , - . / 0 1 2 3 4 5 6 7 8 9 : ; < = > ? @ A B
          C D E F G H I J K L M N O P Q R S T U V W X Y Z [ \ ] ^ _ ` a b c
          d e f g h i j k l m n o p q r s t u v w x y z { | } ~


 -- Function: cunlisp (<lisp_char>)
     Converts a Lisp-character into a Maxima-character.  (You won't
     need it.)


 -- Function: digitcharp (<char>)
     Returns `true' if <char> is a digit.


 -- Function: lcharp (<obj>)
     Returns `true' if <obj> is a Lisp-character.  (You won't need it.)


 -- Function: lowercasep (<char>)
     Returns `true' if <char> is a lowercase character.


 -- Variable: newline
     The newline character.


 -- Variable: space
     The space character.


 -- Variable: tab
     The tab character.


 -- Function: uppercasep (<char>)
     Returns `true' if <char> is an uppercase character.



File: maxima.info,  Node: Functions and Variables for strings,  Prev: Functions and Variables for characters,  Up: stringproc

73.4 Functions and Variables for strings
========================================

 -- Function: sunlisp (<lisp_string>)
     Converts a Lisp-string into a Maxima-string.  (In general you
     won't need it.)


 -- Function: lstringp (<obj>)
     Returns `true' if <obj> is a Lisp-string.  (In general you won't
     need it.)


 -- Function: stringp (<obj>)
     Returns `true' if <obj> is a Maxima-string.  See introduction for
     example.


 -- Function: charat (<string>, <n>)
     Returns the <n>-th character of <string>.  The first character in
     <string> is returned with <n> = 1.

          (%i1) charat("Lisp",1);
          (%o1)                           L


 -- Function: charlist (<string>)
     Returns the list of all characters in <string>.

          (%i1) charlist("Lisp");
          (%o1)                     [L, i, s, p]
          (%i2) %[1];
          (%o2)                           L


 -- Function: parsetoken (<string>)
     `parsetoken' converts the first token in <string> to the
     corresponding number or returns `false' if the number cannot be
     determined .  The delimiter set for tokenizing is `{space, comma,
     semicolon, tab, newline}'.

          (%i1) 2*parsetoken("1.234 5.678");
          (%o1)                         2.468

     For parsing you can also use function parse_string.  See
     description in file 'share\contrib\eval_string.lisp'.


 -- Function: sconc (<expr_1>, ..., <expr_n>)
     Evaluates its arguments and concatenates them into a string.
     `sconc' is like `sconcat' but returns a Maxima string.

          (%i1) sconc("xx[",3,"]:",expand((x+y)^3));
          (%o1)             xx[3]:y^3+3*x*y^2+3*x^2*y+x^3
          (%i2) stringp(%);
          (%o2)                         true


 -- Function: scopy (<string>)
     Returns a copy of <string> as a new string.


 -- Function: sdowncase (<string>)
 -- Function: sdowncase (<string>, <start>)
 -- Function: sdowncase (<string>, <start>, <end>)
     Like `supcase', but uppercase characters are converted to
     lowercase.


 -- Function: sequal (<string_1>, <string_2>)
     Returns `true' if <string_1> and <string_2> are the same length
     and contain the same characters.


 -- Function: sequalignore (<string_1>, <string_2>)
     Like `sequal' but ignores case.


 -- Function: sexplode (<string>)
     `sexplode' is an alias for function `charlist'.


 -- Function: simplode (<list>)
 -- Function: simplode (<list>, <delim>)
     `simplode' takes a list of expressions and concatenates them into
     a string.  If no delimiter <delim> is used, `simplode' is like
     `sconc' and uses no delimiter.  <delim> can be any string.

          (%i1) simplode(["xx[",3,"]:",expand((x+y)^3)]);
          (%o1)             xx[3]:y^3+3*x*y^2+3*x^2*y+x^3
          (%i2) simplode( sexplode("stars")," * " );
          (%o2)                   s * t * a * r * s
          (%i3) simplode( ["One","more","coffee."]," " );
          (%o3)                   One more coffee.


 -- Function: sinsert (<seq>, <string>, <pos>)
     Returns a string that is a concatenation of `substring (<string>,
     1, <pos> - 1)', the string <seq> and `substring (<string>, <pos>)'.
     Note that the first character in <string> is in position 1.

          (%i1) s: "A submarine."$
          (%i2) sconc( substring(s,1,3),"yellow ",substring(s,3) );
          (%o2)                  A yellow submarine.
          (%i3) sinsert("hollow ",s,3);
          (%o3)                  A hollow submarine.


 -- Function: sinvertcase (<string>)
 -- Function: sinvertcase (<string>, <start>)
 -- Function: sinvertcase (<string>, <start>, <end>)
     Returns <string> except that each character from position <start>
     to <end> is inverted.  If <end> is not given, all characters from
     start to the <end> of <string> are replaced.

          (%i1) sinvertcase("sInvertCase");
          (%o1)                      SiNVERTcASE


 -- Function: slength (<string>)
     Returns the number of characters in <string>.


 -- Function: smake (<num>, <char>)
     Returns a new string with a number of <num> characters <char>.

          (%i1) smake(3,"w");
          (%o1)                          www


 -- Function: smismatch (<string_1>, <string_2>)
 -- Function: smismatch (<string_1>, <string_2>, <test>)
     Returns the position of the first character of <string_1> at which
     <string_1> and <string_2> differ or `false'.  Default test
     function for matching is `sequal'.  If `smismatch' should ignore
     case, use `sequalignore' as test.

          (%i1) smismatch("seven","seventh");
          (%o1)                           6


 -- Function: split (<string>)
 -- Function: split (<string>, <delim>)
 -- Function: split (<string>, <delim>, <multiple>)
     Returns the list of all tokens in <string>.  Each token is an
     unparsed string.  `split' uses <delim> as delimiter.  If <delim>
     is not given, the space character is the default delimiter.
     <multiple> is a boolean variable with `true' by default.  Multiple
     delimiters are read as one.  This is useful if tabs are saved as
     multiple space characters.  If <multiple> is set to `false', each
     delimiter is noted.

          (%i1) split("1.2   2.3   3.4   4.5");
          (%o1)                 [1.2, 2.3, 3.4, 4.5]
          (%i2) split("first;;third;fourth",";",false);
          (%o2)               [first, , third, fourth]


 -- Function: sposition (<char>, <string>)
     Returns the position of the first character in <string> which
     matches <char>.  The first character in <string> is in position 1.
     For matching characters ignoring case see `ssearch'.


 -- Function: sremove (<seq>, <string>)
 -- Function: sremove (<seq>, <string>, <test>)
 -- Function: sremove (<seq>, <string>, <test>, <start>)
 -- Function: sremove (<seq>, <string>, <test>, <start>, <end>)
     Returns a string like <string> but without all substrings matching
     <seq>.  Default test function for matching is `sequal'.  If
     `sremove' should ignore case while searching for <seq>, use
     `sequalignore' as test.  Use <start> and <end> to limit searching.
     Note that the first character in <string> is in position 1.

          (%i1) sremove("n't","I don't like coffee.");
          (%o1)                   I do like coffee.
          (%i2) sremove ("DO ",%,'sequalignore);
          (%o2)                    I like coffee.


 -- Function: sremovefirst (<seq>, <string>)
 -- Function: sremovefirst (<seq>, <string>, <test>)
 -- Function: sremovefirst (<seq>, <string>, <test>, <start>)
 -- Function: sremovefirst (<seq>, <string>, <test>, <start>, <end>)
     Like `sremove' except that only the first substring that matches
     `seq' is removed.


 -- Function: sreverse (<string>)
     Returns a string with all the characters of <string> in reverse
     order.


 -- Function: ssearch (<seq>, <string>)
 -- Function: ssearch (<seq>, <string>, <test>)
 -- Function: ssearch (<seq>, <string>, <test>, <start>)
 -- Function: ssearch (<seq>, <string>, <test>, <start>, <end>)
     Returns the position of the first substring of <string> that
     matches the string <seq>.  Default test function for matching is
     `sequal'.  If `ssearch' should ignore case, use `sequalignore' as
     test.  Use <start> and <end> to limit searching.  Note that the
     first character in <string> is in position 1.

          (%i1) ssearch("~s","~{~S ~}~%",'sequalignore);
          (%o1)                                  4


 -- Function: ssort (<string>)
 -- Function: ssort (<string>, <test>)
     Returns a string that contains all characters from <string> in an
     order such there are no two successive characters <c> and <d> such
     that `test (<c>, <d>)' is `false' and `test (<d>, <c>)' is `true'.
     Default test function for sorting is <clessp>.  The set of test
     functions is `{clessp, clesspignore, cgreaterp, cgreaterpignore,
     cequal, cequalignore}'.

          (%i1) ssort("I don't like Mondays.");
          (%o1)                    '.IMaddeiklnnoosty
          (%i2) ssort("I don't like Mondays.",'cgreaterpignore);
          (%o2)                 ytsoonnMlkIiedda.'


 -- Function: ssubst (<new>, <old>, <string>)
 -- Function: ssubst (<new>, <old>, <string>, <test>)
 -- Function: ssubst (<new>, <old>, <string>, <test>, <start>)
 -- Function: ssubst (<new>, <old>, <string>, <test>, <start>, <end>)
     Returns a string like <string> except that all substrings matching
     <old> are replaced by <new>.  <old> and <new> need not to be of
     the same length.  Default test function for matching is `sequal'.
     If `ssubst' should ignore case while searching for old, use
     `sequalignore' as test.  Use <start> and <end> to limit searching.
     Note that the first character in <string> is in position 1.

          (%i1) ssubst("like","hate","I hate Thai food. I hate green tea.");
          (%o1)          I like Thai food. I like green tea.
          (%i2) ssubst("Indian","thai",%,'sequalignore,8,12);
          (%o2)         I like Indian food. I like green tea.


 -- Function: ssubstfirst (<new>, <old>, <string>)
 -- Function: ssubstfirst (<new>, <old>, <string>, <test>)
 -- Function: ssubstfirst (<new>, <old>, <string>, <test>, <start>)
 -- Function: ssubstfirst (<new>, <old>, <string>, <test>, <start>,
          <end>)
     Like `subst' except that only the first substring that matches
     <old> is replaced.


 -- Function: strim (<seq>,<string>)
     Returns a string like <string>, but with all characters that
     appear in <seq> removed from both ends.

          (%i1) "/* comment */"$
          (%i2) strim(" /*",%);
          (%o2)                        comment
          (%i3) slength(%);
          (%o3)                           7


 -- Function: striml (<seq>, <string>)
     Like `strim' except that only the left end of <string> is trimmed.


 -- Function: strimr (<seq>, <string>)
     Like `strim' except that only the right end of string is trimmed.


 -- Function: substring (<string>, <start>)
 -- Function: substring (<string>, <start>, <end>)
     Returns the substring of <string> beginning at position <start>
     and ending at position <end>.  The character at position <end> is
     not included.  If <end> is not given, the substring contains the
     rest of the string.  Note that the first character in <string> is
     in position 1.

          (%i1) substring("substring",4);
          (%o1)                        string
          (%i2) substring(%,4,6);
          (%o2)                          in


 -- Function: supcase (<string>)
 -- Function: supcase (<string>, <start>)
 -- Function: supcase (<string>, <start>, <end>)
     Returns <string> except that lowercase characters from position
     <start> to <end> are replaced by the corresponding uppercase ones.
     If <end> is not given, all lowercase characters from <start> to
     the end of <string> are replaced.

          (%i1) supcase("english",1,2);
          (%o1)                        English


 -- Function: tokens (<string>)
 -- Function: tokens (<string>, <test>)
     Returns a list of tokens, which have been extracted from <string>.
     The tokens are substrings whose characters satisfy a certain test
     function.  If test is not given, <constituent> is used as the
     default test.  `{constituent, alphacharp, digitcharp, lowercasep,
     uppercasep, charp, characterp, alphanumericp}' is the set of test
     functions.  (The Lisp-version of `tokens' is written by Paul
     Graham. ANSI Common Lisp, 1996, page 67.)

          (%i1) tokens("24 October 2005");
          (%o1)                  [24, October, 2005]
          (%i2) tokens("05-10-24",'digitcharp);
          (%o2)                     [05, 10, 24]
          (%i3) map(parsetoken,%);
          (%o3)                      [5, 10, 24]



File: maxima.info,  Node: unit,  Next: zeilberger,  Prev: stringproc,  Up: Top

74 unit
*******

* Menu:

* Introduction to Units::
* Functions and Variables for Units::


File: maxima.info,  Node: Introduction to Units,  Next: Functions and Variables for Units,  Prev: unit,  Up: unit

74.1 Introduction to Units
==========================

The _unit_ package enables the user to convert between arbitrary units
and work with dimensions in equations. The functioning of this package
is radically different from the original Maxima units package - whereas
the original was a basic list of definitions, this package uses
rulesets to allow the user to chose, on a per dimension basis, what
unit final answers should be rendered in.  It will separate units
instead of intermixing them in the display, allowing the user to
readily identify the units associated with a particular answer.  It
will allow a user to simplify an expression to its fundamental Base
Units, as well as providing fine control over simplifying to derived
units.  Dimensional analysis is possible, and a variety of tools are
available to manage conversion and simplification options. In addition
to customizable automatic conversion, _units_ also provides a
traditional manual conversion option.

   Note - when unit conversions are inexact Maxima will make
approximations resulting in fractions. This is a consequence of the
techniques used to simplify units.  The messages warning of this type
of substitution are disabled by default in the case of units (normally
they are on) since this situation occurs frequently and the warnings
clutter the output.  (The existing state of ratprint is restored after
unit conversions, so user changes to that setting will be preserved
otherwise.)  If  the user needs this information for units, they can set
_unitverbose:on_ to reactivate the printing of warnings from the unit
conversion process.

   _unit_ is included in Maxima in the share/contrib/unit directory. It
obeys normal Maxima package loading conventions:

     (%i1) load("unit")$
     *******************************************************************
     *                       Units version 0.50                        *
     *          Definitions based on the NIST Reference on             *
     *              Constants, Units, and Uncertainty                  *
     *       Conversion factors from various sources including         *
     *                   NIST and the GNU units package                *
     *******************************************************************

     Redefining necessary functions...
     WARNING: DEFUN/DEFMACRO: redefining function TOPLEVEL-MACSYMA-EVAL ...
     WARNING: DEFUN/DEFMACRO: redefining function MSETCHK ...
     WARNING: DEFUN/DEFMACRO: redefining function KILL1 ...
     WARNING: DEFUN/DEFMACRO: redefining function NFORMAT ...
     Initializing unit arrays...
     Done.

   The WARNING messages are expected and not a cause for concern - they
indicate the _unit_ package is redefining functions already defined in
Maxima proper.  This is necessary in order to properly handle units.
The user should be aware that if other changes have been made to these
functions by other packages those changes will be overwritten by this
loading process.

   The _unit.mac_ file also loads a lisp file _unit-functions.lisp_
which contains the lisp functions needed for the package.

   Clifford Yapp is the primary author.  He has received valuable
assistance from Barton Willis of the University of Nebraska at Kearney
(UNK), Robert Dodier, and other intrepid folk of the Maxima mailing
list.

   There are probably lots of bugs.  Let me know.  `float' and `numer'
don't do what is expected.

   TODO : dimension functionality, handling of temperature, showabbr
and friends.  Show examples with addition of quantities containing
units.


File: maxima.info,  Node: Functions and Variables for Units,  Prev: Introduction to Units,  Up: unit

74.2 Functions and Variables for Units
======================================

 -- Function: setunits (<list>)
     By default, the _unit_ package does not use any derived
     dimensions, but will convert all units to the seven fundamental
     dimensions using MKS units.
          (%i2) N;
                                               kg m
          (%o2)                                ----
                                                 2
                                                s
          (%i3) dyn;
                                             1      kg m
          (%o3)                           (------) (----)
                                           100000     2
                                                     s
          (%i4) g;
                                              1
          (%o4)                             (----) (kg)
                                             1000
          (%i5) centigram*inch/minutes^2;
                                            127        kg m
          (%o5)                       (-------------) (----)
                                       1800000000000     2
                                                        s

     In some cases this is the desired behavior.  If the user wishes to
     use other units, this is achieved with the `setunits' command:
          (%i6) setunits([centigram,inch,minute]);
          (%o6)                                done
          (%i7) N;
                                      1800000000000   %in cg
          (%o7)                      (-------------) (------)
                                           127            2
                                                      %min
          (%i8) dyn;
                                         18000000   %in cg
          (%o8)                         (--------) (------)
                                           127          2
                                                    %min
          (%i9) g;
          (%o9)                             (100) (cg)
          (%i10) centigram*inch/minutes^2;
                                              %in cg
          (%o10)                              ------
                                                  2
                                              %min

     The setting of units is quite flexible.  For example, if we want to
     get back to kilograms, meters, and seconds as defaults for those
     dimensions we can do:
          (%i11) setunits([kg,m,s]);
          (%o11)                               done
          (%i12) centigram*inch/minutes^2;
                                            127        kg m
          (%o12)                      (-------------) (----)
                                       1800000000000     2
                                                        s

     Derived units are also handled by this command:
          (%i17) setunits(N);
          (%o17)                               done
          (%i18) N;
          (%o18)                                 N
          (%i19) dyn;
                                              1
          (%o19)                           (------) (N)
                                            100000
          (%i20) kg*m/s^2;
          (%o20)                                 N
          (%i21) centigram*inch/minutes^2;
                                              127
          (%o21)                        (-------------) (N)
                                         1800000000000

     Notice that the _unit_ package recognized the non MKS combination
     of mass, length, and inverse time squared as a force, and
     converted it to Newtons.  This is how Maxima works in general.
     If, for example, we prefer dyne to Newtons, we simply do the
     following:
          (%i22) setunits(dyn);
          (%o22)                               done
          (%i23) kg*m/s^2;
          (%o23)                          (100000) (dyn)
          (%i24) centigram*inch/minutes^2;
                                            127
          (%o24)                         (--------) (dyn)
                                          18000000

     To discontinue simplifying to any force, we use the uforget
     command:
          (%i26) uforget(dyn);
          (%o26)                               false
          (%i27) kg*m/s^2;
                                               kg m
          (%o27)                               ----
                                                 2
                                                s
          (%i28) centigram*inch/minutes^2;
                                            127        kg m
          (%o28)                      (-------------) (----)
                                       1800000000000     2
                                                        s
     This would have worked equally well with `uforget(N)' or
     `uforget(%force)'.

     See also `uforget'. To use this function write first
     `load("unit")'.

 -- Function: uforget (<list>)
     By default, the _unit_ package converts all units to the seven
     fundamental dimensions using MKS units. This behavior can be
     changed with the `setunits' command. After that, the user can
     restore the default behavior for a particular dimension by means
     of the `uforget' command:
          (%i13) setunits([centigram,inch,minute]);
          (%o13)                               done
          (%i14) centigram*inch/minutes^2;
                                              %in cg
          (%o14)                              ------
                                                  2
                                              %min
          (%i15) uforget([cg,%in,%min]);
          (%o15)                      [false, false, false]
          (%i16) centigram*inch/minutes^2;
                                            127        kg m
          (%o16)                      (-------------) (----)
                                       1800000000000     2
                                                        s

     `uforget' operates on dimensions, not units, so any unit of a
     particular dimension will work.  The dimension itself is also a
     legal argument.

     See also `setunits'. To use this function write first
     `load("unit")'.

 -- Function: convert (<expr>, <list>)
     When resetting the global environment is overkill, there is the
     `convert' command, which allows one time conversions.  It can
     accept either a single argument or a list of units to use in
     conversion.  When a convert operation is done, the normal global
     evaluation system is bypassed, in order to avoid the desired
     result being converted again.  As a consequence, for inexact
     calculations "rat" warnings will be visible if the global
     environment controlling this behavior (`ratprint') is true.  This
     is also useful for spot-checking the accuracy of a global
     conversion.  Another feature is `convert' will allow a user to do
     Base Dimension conversions even if the global environment is set to
     simplify to a Derived Dimension.

          (%i2) kg*m/s^2;
                                               kg m
          (%o2)                                ----
                                                 2
                                                s
          (%i3) convert(kg*m/s^2,[g,km,s]);
                                               g km
          (%o3)                                ----
                                                 2
                                                s
          (%i4) convert(kg*m/s^2,[g,inch,minute]);

          `rat' replaced 39.37007874015748 by 5000//127 = 39.37007874015748
                                        18000000000   %in g
          (%o4)                        (-----------) (-----)
                                            127           2
                                                      %min
          (%i5) convert(kg*m/s^2,[N]);
          (%o5)                                  N
          (%i6) convert(kg*m^2/s^2,[N]);
          (%o6)                                 m N
          (%i7) setunits([N,J]);
          (%o7)                                done
          (%i8) convert(kg*m^2/s^2,[N]);
          (%o8)                                 m N
          (%i9) convert(kg*m^2/s^2,[N,inch]);

          `rat' replaced 39.37007874015748 by 5000//127 = 39.37007874015748
                                           5000
          (%o9)                           (----) (%in N)
                                           127
          (%i10) convert(kg*m^2/s^2,[J]);
          (%o10)                                 J
          (%i11) kg*m^2/s^2;
          (%o11)                                 J
          (%i12) setunits([g,inch,s]);
          (%o12)                               done
          (%i13) kg*m/s^2;
          (%o13)                                 N
          (%i14) uforget(N);
          (%o14)                               false
          (%i15) kg*m/s^2;
                                          5000000   %in g
          (%o15)                         (-------) (-----)
                                            127       2
                                                     s
          (%i16) convert(kg*m/s^2,[g,inch,s]);

          `rat' replaced 39.37007874015748 by 5000//127 = 39.37007874015748
                                          5000000   %in g
          (%o16)                         (-------) (-----)
                                            127       2
                                                     s

     See also `setunits' and `uforget'. To use this function write
     first `load("unit")'.

 -- Optional variable: usersetunits
     Default value: none

     If a user wishes to have a default unit behavior other than that
     described, they can make use of _maxima-init.mac_ and the
     _usersetunits_ variable.  The _unit_ package will check on startup
     to see if this variable has been assigned a list.  If it has, it
     will use setunits on that list and take the units from that list
     to be defaults.  `uforget' will revert to the behavior defined by
     usersetunits over its own defaults.  For example, if we have a
     _maxima-init.mac_ file containing:
          usersetunits : [N,J];
     we would see the following behavior:
          (%i1) load("unit")$
          *******************************************************************
          *                       Units version 0.50                        *
          *          Definitions based on the NIST Reference on             *
          *              Constants, Units, and Uncertainty                  *
          *       Conversion factors from various sources including         *
          *                   NIST and the GNU units package                *
          *******************************************************************

          Redefining necessary functions...
          WARNING: DEFUN/DEFMACRO: redefining function
           TOPLEVEL-MACSYMA-EVAL ...
          WARNING: DEFUN/DEFMACRO: redefining function MSETCHK ...
          WARNING: DEFUN/DEFMACRO: redefining function KILL1 ...
          WARNING: DEFUN/DEFMACRO: redefining function NFORMAT ...
          Initializing unit arrays...
          Done.
          User defaults found...
          User defaults initialized.
          (%i2) kg*m/s^2;
          (%o2)                                  N
          (%i3) kg*m^2/s^2;
          (%o3)                                  J
          (%i4) kg*m^3/s^2;
          (%o4)                                 J m
          (%i5) kg*m*km/s^2;
          (%o5)                             (1000) (J)
          (%i6) setunits([dyn,eV]);
          (%o6)                                done
          (%i7) kg*m/s^2;
          (%o7)                           (100000) (dyn)
          (%i8) kg*m^2/s^2;
          (%o8)                     (6241509596477042688) (eV)
          (%i9) kg*m^3/s^2;
          (%o9)                    (6241509596477042688) (eV m)
          (%i10) kg*m*km/s^2;
          (%o10)                   (6241509596477042688000) (eV)
          (%i11) uforget([dyn,eV]);
          (%o11)                           [false, false]
          (%i12) kg*m/s^2;
          (%o12)                                 N
          (%i13) kg*m^2/s^2;
          (%o13)                                 J
          (%i14) kg*m^3/s^2;
          (%o14)                                J m
          (%i15) kg*m*km/s^2;
          (%o15)                            (1000) (J)
     Without `usersetunits', the initial inputs would have been
     converted to MKS, and uforget would have resulted in a return to
     MKS rules.  Instead, the user preferences are respected in both
     cases.  Notice these can still be overridden if desired.  To
     completely eliminate this simplification - i.e.  to have the user
     defaults reset to factory defaults - the `dontusedimension'
     command can be used.  `uforget' can restore user settings again,
     but only if `usedimension' frees it for use.  Alternately,
     `kill(usersetunits)' will completely remove all knowledge of the
     user defaults from the session.  Here are some examples of how
     these various options work.
          (%i2) kg*m/s^2;
          (%o2)                                  N
          (%i3) kg*m^2/s^2;
          (%o3)                                  J
          (%i4) setunits([dyn,eV]);
          (%o4)                                done
          (%i5) kg*m/s^2;
          (%o5)                           (100000) (dyn)
          (%i6) kg*m^2/s^2;
          (%o6)                     (6241509596477042688) (eV)
          (%i7) uforget([dyn,eV]);
          (%o7)                          [false, false]
          (%i8) kg*m/s^2;
          (%o8)                                  N
          (%i9) kg*m^2/s^2;
          (%o9)                                  J
          (%i10) dontusedimension(N);
          (%o10)                             [%force]
          (%i11) dontusedimension(J);
          (%o11)                         [%energy, %force]
          (%i12) kg*m/s^2;
                                               kg m
          (%o12)                               ----
                                                 2
                                                s
          (%i13) kg*m^2/s^2;
                                                   2
                                               kg m
          (%o13)                               -----
                                                 2
                                                s
          (%i14) setunits([dyn,eV]);
          (%o14)                               done
          (%i15) kg*m/s^2;
                                               kg m
          (%o15)                               ----
                                                 2
                                                s
          (%i16) kg*m^2/s^2;
                                                   2
                                               kg m
          (%o16)                               -----
                                                 2
                                                s
          (%i17) uforget([dyn,eV]);
          (%o17)                         [false, false]
          (%i18) kg*m/s^2;
                                               kg m
          (%o18)                               ----
                                                 2
                                                s
          (%i19) kg*m^2/s^2;
                                                   2
                                               kg m
          (%o19)                               -----
                                                 2
                                                s
          (%i20) usedimension(N);
          Done.  To have Maxima simplify to this dimension, use
          setunits([unit]) to select a unit.
          (%o20)                               true
          (%i21) usedimension(J);
          Done.  To have Maxima simplify to this dimension, use
          setunits([unit]) to select a unit.
          (%o21)                               true
          (%i22) kg*m/s^2;
                                               kg m
          (%o22)                               ----
                                                 2
                                                s
          (%i23) kg*m^2/s^2;
                                                   2
                                               kg m
          (%o23)                               -----
                                                 2
                                                s
          (%i24) setunits([dyn,eV]);
          (%o24)                               done
          (%i25) kg*m/s^2;
          (%o25)                          (100000) (dyn)
          (%i26) kg*m^2/s^2;
          (%o26)                    (6241509596477042688) (eV)
          (%i27) uforget([dyn,eV]);
          (%o27)                           [false, false]
          (%i28) kg*m/s^2;
          (%o28)                                 N
          (%i29) kg*m^2/s^2;
          (%o29)                                 J
          (%i30) kill(usersetunits);
          (%o30)                               done
          (%i31) uforget([dyn,eV]);
          (%o31)                          [false, false]
          (%i32) kg*m/s^2;
                                               kg m
          (%o32)                               ----
                                                 2
                                                s
          (%i33) kg*m^2/s^2;
                                                   2
                                               kg m
          (%o33)                               -----
                                                 2
                                                s
     Unfortunately this wide variety of options is a little confusing
     at first, but once the user grows used to them they should find
     they have very full control over their working environment.


 -- Function: metricexpandall (<x>)
     Rebuilds global unit lists automatically creating all desired
     metric units.  <x> is a numerical argument which is used to
     specify how many metric prefixes the user wishes defined.  The
     arguments are as follows, with each higher number defining all
     lower numbers' units:
                     0 - none. Only base units
                     1 - kilo, centi, milli
          (default)  2 - giga, mega, kilo, hecto, deka, deci, centi, milli,
                         micro, nano
                     3 - peta, tera, giga, mega, kilo, hecto, deka, deci,
                         centi, milli, micro, nano, pico, femto
                     4 - all
     Normally, Maxima will not define the full expansion since this
     results in a very large number of units, but `metricexpandall' can
     be used to rebuild the list in a more or less complete fashion.
     The relevant variable in the _unit.mac_ file is <%unitexpand>.

 -- Variable: %unitexpand
     Default value: `2'

     This is the value supplied to `metricexpandall' during the initial
     loading of _unit_.



File: maxima.info,  Node: zeilberger,  Next: Function and Variable Index,  Prev: unit,  Up: Top

75 zeilberger
*************

* Menu:

* Introduction to zeilberger::
* Functions and Variables for zeilberger::


File: maxima.info,  Node: Introduction to zeilberger,  Next: Functions and Variables for zeilberger,  Prev: zeilberger,  Up: zeilberger

75.1 Introduction to zeilberger
===============================

`zeilberger' is a implementation of Zeilberger's algorithm for definite
hypergeometric summation, and also Gosper's algorithm for indefinite
hypergeometric summation.

   `zeilberger' makes use of the "filtering" optimization method
developed by Axel Riese.

   `zeilberger' was developed by Fabrizio Caruso.

   `load (zeilberger)' loads this package.

75.1.0.1 The indefinite summation problem
.........................................

`zeilberger' implements Gosper's algorithm for indefinite
hypergeometric summation.  Given a hypergeometric term F_k in k we want
to find its hypergeometric anti-difference, that is, a hypergeometric
term f_k such that F_k = f_(k+1) - f_k.

75.1.0.2 The definite summation problem
.......................................

`zeilberger' implements Zeilberger's algorithm for definite
hypergeometric summation.  Given a proper hypergeometric term (in n and
k) F_(n,k) and a positive integer d we want to find a d-th order linear
recurrence with polynomial coefficients (in n) for F_(n,k) and a
rational function R in n and k such that

   a_0 F_(n,k) + ... + a_d F_(n+d),k = Delta_K(R(n,k) F_(n,k))

   where Delta_k is the k-forward difference operator, i.e.,
Delta_k(t_k) := t_(k+1) - t_k.

75.1.1 Verbosity levels
-----------------------

There are also verbose versions of the commands which are called by
adding one of the following prefixes:

`Summary'
     Just a summary at the end is shown

`Verbose'
     Some information in the intermidiate steps

`VeryVerbose'
     More information

`Extra'
     Even more information including information on the linear system
     in Zeilberger's algorithm

   For example: `GosperVerbose', `parGosperVeryVerbose',
`ZeilbergerExtra', `AntiDifferenceSummary'.


File: maxima.info,  Node: Functions and Variables for zeilberger,  Prev: Introduction to zeilberger,  Up: zeilberger

75.2 Functions and Variables for zeilberger
===========================================

 -- Function: AntiDifference (<F_k>, <k>)
     Returns the hypergeometric anti-difference of <F_k>, if it exists.
     Otherwise `AntiDifference' returns `no_hyp_antidifference'.

 -- Function: Gosper (<F_k>, <k>)
     Returns the rational certificate <R(k)> for <F_k>, that is, a
     rational function such that

     F_k = R(k+1) F_(k+1) - R(k) F_k

     if it exists.  Otherwise, `Gosper' returns `no_hyp_sol'.

 -- Function: GosperSum (<F_k>, <k>, <a>, <b>)
     Returns the summmation of <F_k> from <k> = <a> to <k> = <b> if
     <F_k> has a hypergeometric anti-difference.  Otherwise,
     `GosperSum' returns `nongosper_summable'.

     Examples:

          (%i1) load (zeilberger);
          (%o1)  /usr/share/maxima/share/contrib/Zeilberger/zeilberger.mac
          (%i2) GosperSum ((-1)^k*k / (4*k^2 - 1), k, 1, n);

          Dependent equations eliminated:  (1)
                                     3       n + 1
                                (n + -) (- 1)
                                     2               1
          (%o2)               - ------------------ - -
                                            2        4
                                2 (4 (n + 1)  - 1)
          (%i3) GosperSum (1 / (4*k^2 - 1), k, 1, n);
                                          3
                                    - n - -
                                          2       1
          (%o3)                  -------------- + -
                                          2       2
                                 4 (n + 1)  - 1
          (%i4) GosperSum (x^k, k, 1, n);
                                    n + 1
                                   x          x
          (%o4)                    ------ - -----
                                   x - 1    x - 1
          (%i5) GosperSum ((-1)^k*a! / (k!*(a - k)!), k, 1, n);
                                          n + 1
                          a! (n + 1) (- 1)              a!
          (%o5)       - ------------------------- - ----------
                        a (- n + a - 1)! (n + 1)!   a (a - 1)!
          (%i6) GosperSum (k*k!, k, 1, n);

          Dependent equations eliminated:  (1)
          (%o6)                     (n + 1)! - 1
          (%i7) GosperSum ((k + 1)*k! / (k + 1)!, k, 1, n);
                            (n + 1) (n + 2) (n + 1)!
          (%o7)             ------------------------ - 1
                                    (n + 2)!
          (%i8) GosperSum (1 / ((a - k)!*k!), k, 1, n);
          (%o8)                  nonGosper_summable

 -- Function: parGosper (<F_{n,k}>, <k>, <n>, <d>)
     Attempts to find a <d>-th order recurrence for <F_{n,k}>.

     The algorithm yields a sequence [s_1, s_2, ..., s_m] of solutions.
     Each solution has the form

     [R(n, k), [a_0, a_1, ..., a_d]]

     `parGosper' returns `[]' if it fails to find a recurrence.

 -- Function: Zeilberger (<F_{n,k}>, <k>, <n>)
     Attempts to compute the indefinite hypergeometric summation of
     <F_{n,k}>.

     `Zeilberger' first invokes `Gosper', and if that fails to find a
     solution, then invokes `parGosper' with order 1, 2, 3, ..., up to
     `MAX_ORD'.  If Zeilberger finds a solution before reaching
     `MAX_ORD', it stops and returns the solution.

     The algorithms yields a sequence [s_1, s_2, ..., s_m] of solutions.
     Each solution has the form

     [R(n,k), [a_0, a_1, ..., a_d]]

     `Zeilberger' returns `[]' if it fails to find a solution.

     `Zeilberger' invokes `Gosper' only if `gosper_in_zeilberger' is
     `true'.

75.3 General global variables
=============================

 -- Global variable: MAX_ORD
     Default value: 5

     `MAX_ORD' is the maximum recurrence order attempted by
     `Zeilberger'.

 -- Global variable: simplified_output
     Default value: `false'

     When `simplified_output' is `true', functions in the `zeilberger'
     package attempt further simplification of the solution.

 -- Global variable: linear_solver
     Default value: `linsolve'

     `linear_solver' names the solver which is used to solve the system
     of equations in Zeilberger's algorithm.

 -- Global variable: warnings
     Default value: `true'

     When `warnings' is `true', functions in the `zeilberger' package
     print warning messages during execution.

 -- Global variable: gosper_in_zeilberger
     Default value: `true'

     When `gosper_in_zeilberger' is `true', the `Zeilberger' function
     calls `Gosper' before calling `parGosper'.  Otherwise,
     `Zeilberger' goes immediately to `parGosper'.

 -- Global variable: trivial_solutions
     Default value: `true'

     When `trivial_solutions' is `true', `Zeilberger' returns solutions
     which have certificate equal to zero, or all coefficients equal to
     zero.

75.4 Variables related to the modular test
==========================================

 -- Global variable: mod_test
     Default value: `false'

     When `mod_test' is `true', `parGosper' executes a modular test for
     discarding systems with no solutions.

 -- Global variable: modular_linear_solver
     Default value: `linsolve'

     `modular_linear_solver' names the linear solver used by the
     modular test in `parGosper'.

 -- Global variable: ev_point
     Default value: `big_primes[10]'

     `ev_point' is the value at which the variable <n> is evaluated
     when executing the modular test in `parGosper'.

 -- Global variable: mod_big_prime
     Default value: `big_primes[1]'

     `mod_big_prime' is the modulus used by the modular test in
     `parGosper'.

 -- Global variable: mod_threshold
     Default value: 4

     `mod_threshold' is the greatest order for which the modular test
     in `parGosper' is attempted.