This is Info file octave, produced by Makeinfo-1.64 from the input file
octave.tex.

START-INFO-DIR-ENTRY
* Octave: (octave).     Interactive language for numerical computations.
END-INFO-DIR-ENTRY

   Copyright (C) 1996 John W. Eaton.

   Permission is granted to make and distribute verbatim copies of this
manual provided the copyright notice and this permission notice are
preserved on all copies.

   Permission is granted to copy and distribute modified versions of
this manual under the conditions for verbatim copying, provided that
the entire resulting derived work is distributed under the terms of a
permission notice identical to this one.

   Permission is granted to copy and distribute translations of this
manual into another language, under the above conditions for modified
versions.


File: octave,  Node: Utility Functions,  Next: Complex Arithmetic,  Prev: Arithmetic,  Up: Arithmetic

Utility Functions
=================

   The following functions are available for working with complex
numbers.  Each expects a single argument.  They are called "mapping
functions" because when given a matrix argument, they apply the given
function to each element of the matrix.

 - Mapping Function:  ceil (X)
     Return the smallest integer not less than X.  If X is complex,
     return `ceil (real (X)) + ceil (imag (X)) * I'.

 - Mapping Function:  floor (X)
     Return the largest integer not greater than X.  If X is complex,
     return `floor (real (X)) + floor (imag (X)) * I'.

 - Mapping Function:  fix (X)
     Truncate X toward zero.  If X is complex, return `fix (real (X)) +
     fix (imag (X)) * I'.

 - Mapping Function:  round (X)
     Return the integer nearest to X.  If X is complex, return `round
     (real (X)) + round (imag (X)) * I'.

 - Mapping Function:  sign (X)
     Compute the "signum" function, which is defined as

                     -1, x < 0;
          sign (x) =  0, x = 0;
                      1, x > 0.

     For complex arguments, `sign' returns `x ./ abs (X)'.

 - Mapping Function:  exp (X)
     Compute the exponential of X.  To compute the matrix exponential,
     see *Note Linear Algebra::.

 - Mapping Function:  gcd (X, `...')
     Compute the greatest common divisor of the elements of X, or the
     list of all the arguments.  For example,

          gcd (a1, ..., ak)

     is the same as

          gcd ([a1, ..., ak])

     An optional second return value, V contains an integer vector such
     that

          g = v(1) * a(k) + ... + v(k) * a(k)

 - Mapping Function:  lcm (X, `...')
     Compute the least common multiple of the elements elements of X, or
     the list of all the arguments.  For example,

          lcm (a1, ..., ak)

     is the same as

          lcm ([a1, ..., ak]).

 - Mapping Function:  log (X)
     Compute the natural logarithm of X.  To compute the matrix
     logarithm, see *Note Linear Algebra::.

 - Mapping Function:  log2 (X)
     Compute the base-2 logarithm of X.

 - Mapping Function:  log10 (X)
     Compute the base-10 logarithm of X.

 - Mapping Function:  sqrt (X)
     Compute the square root of X.  To compute the matrix square root,
     see *Note Linear Algebra::.

 - Loadable Function:  max (X)
     For a vector argument, return the maximum value.  For a matrix
     argument, return the maximum value from each column, as a row
     vector.  Thus,

          max (max (X))

     returns the largest element of X.

     For complex arguments, the magnitude of the elements are used for
     comparison.

 - Loadable Function:  min (X)
     Like `max', but return the minimum value.

 - Mapping Function:  rem (X, Y)
     Return the remainder of `X / Y', computed using the expression

          x - y .* fix (x ./ y)

     An error message is printed if the dimensions of the arguments do
     not agree, or if either of the arguments is complex.


File: octave,  Node: Complex Arithmetic,  Next: Trigonometry,  Prev: Utility Functions,  Up: Arithmetic

Complex Arithmetic
==================

   The following functions are available for working with complex
numbers.  Each expects a single argument.  Given a matrix they work on
an element by element basis.

 - Mapping Function:  abs (Z)
     Compute the magnitude of Z.

 - Mapping Function:  arg (Z)
 - Mapping Function:  angle (Z)
     Compute the argument of Z.

 - Mapping Function:  conj (Z)
     Return the complex conjugate of Z.

 - Mapping Function:  imag (Z)
     Return the imaginary part of Z.

 - Mapping Function:  real (Z)
     Return the real part of Z.


File: octave,  Node: Trigonometry,  Next: Sums and Products,  Prev: Complex Arithmetic,  Up: Arithmetic

Trigonometry
============

   Octave provides the following trigonometric functions:

 - Mapping Function:  sin (Z)
 - Mapping Function:  cos (Z)
 - Mapping Function:  tan (Z)
 - Mapping Function:  sec (Z)
 - Mapping Function:  csc (Z)
 - Mapping Function:  cot (Z)
     The ordinary trigonometric functions.

 - Mapping Function:  asin (Z)
 - Mapping Function:  acos (Z)
 - Mapping Function:  atan (Z)
 - Mapping Function:  asec (Z)
 - Mapping Function:  acsc (Z)
 - Mapping Function:  acot (Z)
     The ordinary inverse trigonometric functions.

 - Mapping Function:  sinh (Z)
 - Mapping Function:  cosh (Z)
 - Mapping Function:  tanh (Z)
 - Mapping Function:  sech (Z)
 - Mapping Function:  csch (Z)
 - Mapping Function:  coth (Z)
     Hyperbolic trigonometric functions.

 - Mapping Function:  asinh (Z)
 - Mapping Function:  acosh (Z)
 - Mapping Function:  atanh (Z)
 - Mapping Function:  asech (Z)
 - Mapping Function:  acsch (Z)
 - Mapping Function:  acoth (Z)
     Inverse hyperbolic trigonometric functions.

   Each of these functions expect a single argument.  For matrix
arguments, they work on an element by element basis.  For example,

     sin ([1, 2; 3, 4])
          =>  0.84147   0.90930
              0.14112  -0.75680

 - Mapping Function:  atan2 (Y, X)
     Return the arctangent of Y/X.  The signs of the arguments are used
     to determine the quadrant of the result, which is in the range
     `pi' to -`pi'.


File: octave,  Node: Sums and Products,  Next: Special Functions,  Prev: Trigonometry,  Up: Arithmetic

Sums and Products
=================

 - Built-in Function:  sum (X)
     For a vector argument, return the sum of all the elements.  For a
     matrix argument, return the sum of the elements in each column, as
     a row vector.  The sum of an empty matrix is 0 if it has no
     columns, or a vector of zeros if it has no rows (*note Empty
     Matrices::.).

 - Built-in Function:  prod (X)
     For a vector argument, return the product of all the elements.
     For a matrix argument, return the product of the elements in each
     column, as a row vector.  The product of an empty matrix is 1 if
     it has no columns, or a vector of ones if it has no rows (*note
     Empty Matrices::.).

 - Built-in Function:  cumsum (X)
     Return the cumulative sum of each column of X.  For example,

          cumsum ([1, 2; 3, 4])
               =>  1  2
                   4  6

 - Built-in Function:  cumprod (X)
     Return the cumulative product of each column of X.  For example,

          cumprod ([1, 2; 3, 4])
               =>  1  2
                   3  8

 - Built-in Function:  sumsq (X)
     For a vector argument, return the sum of the squares of all the
     elements.  For a matrix argument, return the sum of the squares of
     the elements in each column, as a row vector.


File: octave,  Node: Special Functions,  Next: Mathematical Constants,  Prev: Sums and Products,  Up: Arithmetic

Special Functions
=================

 - Mapping Function:  beta (A, B)
     Returns the Beta function,

          beta (a, b) = gamma (a) * gamma (b) / gamma (a + b).

 - Mapping Function:  betai (A, B, X)
     Returns the incomplete Beta function,

                                              x
                                             /
          betai (a, b, x) = beta (a, b)^(-1) | t^(a-1) (1-t)^(b-1) dt.
                                             /
                                          t=0

     If x has more than one component, both A and B must be scalars.
     If X is a scalar, A and B must be of compatible dimensions.

 - Mapping Function:  erf (Z)
     Computes the error function,

                                   z
                                  /
          erf (z) = (2/sqrt (pi)) | e^(-t^2) dt
                                  /
                               t=0

 - Mapping Function:  erfc (Z)
     Computes the complementary error function, `1 - erf (Z)'.

 - Mapping Function:  gamma (Z)
     Computes the Gamma function,

                      infinity
                      /
          gamma (z) = | t^(z-1) exp (-t) dt.
                      /
                   t=0

 - Mapping Function:  gammai (A, X)
     Computes the incomplete gamma function,

                                        x
                              1        /
          gammai (a, x) = ---------    | exp (-t) t^(a-1) dt
                          gamma (a)    /
                                    t=0

     If A is scalar, then `gammai (A, X)' is returned for each element
     of X and vice versa.

     If neither A nor X is scalar, the sizes of A and X must agree, and
     GAMMAI is applied element-by-element.

 - Mapping Function:  lgamma (A, X)
     Returns the natural logarithm of the gamma function.


File: octave,  Node: Mathematical Constants,  Prev: Special Functions,  Up: Arithmetic

Mathematical Constants
======================

 - Built-in Variable: I
 - Built-in Variable: J
 - Built-in Variable: i
 - Built-in Variable: j
     A pure imaginary number, defined as   `sqrt (-1)'.  The `I' and
     `J' forms are true constants, and cannot be modified.  The `i' and
     `j' forms are like ordinary variables, and may be used for other
     purposes.  However, unlike other variables, they once again assume
     their special predefined values if they are cleared *Note
     Miscellaneous Utilities::.

 - Built-in Variable: Inf
 - Built-in Variable: inf
     Infinity.  This is the result of an operation like 1/0, or an
     operation that results in a floating point overflow.

 - Built-in Variable: NaN
 - Built-in Variable: nan
     Not a number.  This is the result of an operation like `0/0', or
     `Inf - Inf', or any operation with a NaN.

 - Built-in Variable: eps
     The machine precision.  More precisely, `eps' is the largest
     relative spacing between any two adjacent numbers in the machine's
     floating point system.  This number is obviously system-dependent.
     On machines that support 64 bit IEEE floating point arithmetic,
     `eps' is approximately  2.2204e-16.

 - Built-in Variable: pi
     The ratio of the circumference of a circle to its diameter.
     Internally, `pi' is computed as `4.0 * atan (1.0)'.

 - Built-in Variable: e
     The base of natural logarithms.  The constant  E  satisfies the
     equation  `log' (E) = 1.

 - Built-in Variable: realmax
     The largest floating point number that is representable.  The
     actual value is system-dependent.  On machines that support 64 bit
     IEEE floating point arithmetic, `realmax' is approximately
     1.7977e+308

 - Built-in Variable: realmin
     The smallest floating point number that is representable.  The
     actual value is system-dependent.  On machines that support 64 bit
     IEEE floating point arithmetic, `realmin' is approximately
     2.2251e-308


File: octave,  Node: Linear Algebra,  Next: Polynomial Manipulations,  Prev: Arithmetic,  Up: Top

Linear Algebra
**************

   This chapter documents the linear algebra functions of Octave.
Reference material for many of these options may be found in Golub and
Van Loan, `Matrix Computations, 2nd Ed.', Johns Hopkins, 1989, and in
`LAPACK Users' Guide', SIAM, 1992.

* Menu:

* Basic Matrix Functions::
* Matrix Factorizations::
* Functions of a Matrix::


File: octave,  Node: Basic Matrix Functions,  Next: Matrix Factorizations,  Prev: Linear Algebra,  Up: Linear Algebra

Basic Matrix Functions
======================

 - Loadable Function: AA = balance (A, OPT)
 - Loadable Function: [DD, AA] = balance (A, OPT)
 - Loadable Function: [CC, DD, AA, BB] = balance (A, B, OPT)
     `[dd, aa] = balance (a)' returns `aa = dd \ a * dd'.  `aa' is a
     matrix whose row/column norms are roughly equal in magnitude, and
     `dd' = `p * d', where `p' is a permutation matrix and `d' is a
     diagonal matrix of powers of two.  This allows the equilibration
     to be computed without roundoff.  Results of eigenvalue
     calculation are typically improved by balancing first.

     `[cc, dd, aa, bb] = balance (a, b)'  returns `aa' (`bb') `=
     cc*a*dd (cc*b*dd)'), where `aa' and `bb' have non-zero elements of
     approximately the same magnitude and `cc' and `dd'  are permuted
     diagonal matrices as in `dd' for the algebraic eigenvalue problem.

     The eigenvalue balancing option `opt' is selected as follows:

    `"N"', `"n"'
          No balancing; arguments copied, transformation(s) set to
          identity.

    `"P"', `"p"'
          Permute argument(s) to isolate eigenvalues where possible.

    `"S"', `"s"'
          Scale to improve accuracy of computed eigenvalues.

    `"B"', `"b"'
          Permute and scale, in that order. Rows/columns of a (and b)
          that are isolated by permutation are not scaled.  This is the
          default behavior.

     Algebraic eigenvalue balancing uses standard LAPACK routines.

     Generalized eigenvalue problem balancing uses Ward's algorithm
     (SIAM Journal on Scientific and Statistical Computing, 1981).

 - :  cond (A)
     Compute the (two-norm) condition number of a matrix. `cond (a)' is
     defined as `norm (a) * norm (inv (a))', and is computed via a
     singular value decomposition.

 - Loadable Function:  det (A)
     Compute the determinant of A using LINPACK.

 - Loadable Function: LAMBDA = eig (A)
 - Loadable Function: [V, LAMBDA] = eig (A)
     The eigenvalues (and eigenvectors) of a matrix are computed in a
     several step process which begins with a Hessenberg decomposition
     (see `hess'), followed by a Schur decomposition (see `schur'), from
     which the eigenvalues are apparent.  The eigenvectors, when
     desired, are computed by further manipulations of the Schur
     decomposition.

     See also: `hess', `schur'.

 - Loadable Function: G = givens (X, Y)
 - Loadable Function: [C, S] = givens (X, Y)
     `G = givens(x, y)' returns a 2 x 2 orthogonal matrix `G = [c s;
     -s' c]' such that `G [x; y] = [*; 0]'  (x, y scalars)

 - Loadable Function:  inv (A)
 - Loadable Function:  inverse (A)
     Compute the inverse of the square matrix A.

 - Function File:  norm (A, P)
     Compute the p-norm of the matrix A.  If the second argument is
     missing, `p = 2' is assumed.

     If A is a matrix:

    P = `1'
          1-norm, the largest column sum of A.

    P = `2'
          Largest singular value of A.

    P = `Inf'
          Infinity norm, the largest row sum of A.

    P = `"fro"'
          Frobenius norm of A, `sqrt (sum (diag (a' * a)))'.

     If A is a vector or a scalar:

    P = `Inf'
          `max (abs (a))'.

    P = `-Inf'
          `min (abs (a))'.

    other
          p-norm of A, `(sum (abs (a) .^ p)) ^ (1/p)'.

 - Function File:  null (A, TOL)
     Returns an orthonormal basis of the null space of A.

     The dimension of the null space is taken as the number of singular
     values of A not greater than TOL.  If the argument TOL is missing,
     it is computed as

          max (size (a)) * max (svd (a)) * eps

 - Function File:  orth (A, TOL)
     Returns an orthonormal basis of the range of A.

     The dimension of the range space is taken as the number of singular
     values of A greater than TOL.  If the argument TOL is missing, it
     is computed as

          max (size (a)) * max (svd (a)) * eps

 - Function File:  pinv (X, TOL)
     Returns the pseudoinverse of X.  Singular values less than TOL are
     ignored.

     If the second argument is omitted, it is assumed that

          tol = max (size (X)) * sigma_max (X) * eps,

     where `sigma_max (X)' is the maximal singular value of X.

 - Function File:  rank (A, TOL)
     Compute the rank of A, using the singular value decomposition.
     The rank is taken to be the number  of singular values of A that
     are greater than the specified tolerance TOL.  If the second
     argument is omitted, it is taken to be

          tol = max (size (a)) * sigma (1) * eps;

     where `eps' is machine precision and `sigma' is the largest
     singular value of `a'.

 - Function File:  trace (A)
     Compute the trace of A, `sum (diag (a))'.


File: octave,  Node: Matrix Factorizations,  Next: Functions of a Matrix,  Prev: Basic Matrix Functions,  Up: Linear Algebra

Matrix Factorizations
=====================

 - Loadable Function:  chol (A)
     Compute the Cholesky factor, R, of the symmetric positive definite
     matrix A, where

          r' * r = a.

 - Loadable Function: H = hess (A)
 - Loadable Function: [P, H] = hess (A)
     Compute the Hessenberg decomposition of the matrix A.

     The Hessenberg decomposition is usually used as the first step in
     an eigenvalue computation, but has other applications as well (see
     Golub, Nash, and Van Loan, IEEE Transactions on Automatic Control,
     1979.  The Hessenberg decomposition is `p * h * p' = a' where `p'
     is a square unitary matrix (`p' * p = I', using complex-conjugate
     transposition) and `h' is upper Hessenberg (`i >= j+1 => h (i, j)
     = 0').

 - Loadable Function: [L, U, P] = lu (A)
     Compute the LU decomposition of A, using subroutines from LAPACK.
     The result is returned in a permuted form, according to the
     optional return value P.  For example, given the matrix `a = [1,
     2; 3, 4]',

          [l, u, p] = lu (a)

     returns

          l =
          
            1.00000  0.00000
            0.33333  1.00000
          
          u =
          
            3.00000  4.00000
            0.00000  0.66667
          
          p =
          
            0  1
            1  0

 - Loadable Function: [Q, R] = qr (A)
     Compute the QR factorization of A, using standard LAPACK
     subroutines.  For example, given the matrix `a = [1, 2; 3, 4]',

          [q, r] = qr (a)

     returns

          q =
          
            -0.31623  -0.94868
            -0.94868   0.31623
          
          r =
          
            -3.16228  -4.42719
             0.00000  -0.63246

     The `qr' factorization has applications in the solution of least
     squares problems

          `min norm(A x - b)'

     for overdetermined systems of equations (i.e., `a'  is a tall,
     thin matrix).  The `qr' factorization is `q * r = a' where `q' is
     an orthogonal matrix and `r' is upper triangular.

     The permuted `qr' factorization `[q, r, pi] = qr (a)' forms the
     `qr' factorization such that the diagonal entries of `r' are
     decreasing in magnitude order.  For example, given the matrix `a =
     [1, 2; 3, 4]',

          [q, r, pi] = qr(a)

     returns

          q =
          
            -0.44721  -0.89443
            -0.89443   0.44721
          
          r =
          
            -4.47214  -3.13050
             0.00000   0.44721
          
          p =
          
             0  1
             1  0

     The permuted `qr' factorization `[q, r, pi] = qr (a)'
     factorization allows the construction of an orthogonal basis of
     `span (a)'.

 - Loadable Function: S schur (A)
 - Loadable Function: [U, S] = schur (A, OPT)
     The Schur decomposition is used to compute eigenvalues of a square
     matrix, and has applications in the solution of algebraic Riccati
     equations in control (see `are' and `dare').  `schur' always
     returns `s = u' * a * u' where `u'  is a unitary matrix (`u'* u'
     is identity) and `s' is upper triangular.  The eigenvalues of `a'
     (and `s') are the diagonal elements of `s' If the matrix `a' is
     real, then the real Schur decomposition is computed, in which the
     matrix `u' is orthogonal and `s' is block upper triangular with
     blocks of size at most `2 x 2' blocks along the diagonal.  The
     diagonal elements of `s' (or the eigenvalues of the `2 x 2'
     blocks, when appropriate) are the eigenvalues of `a' and `s'.

     The eigenvalues are optionally ordered along the diagonal
     according to the value of `opt'.  `opt = "a"' indicates that all
     eigenvalues with negative real parts should be moved to the leading
     block of `s' (used in `are'), `opt = "d"' indicates that all
     eigenvalues with magnitude less than one should be moved to the
     leading block of `s' (used in `dare'), and `opt = "u"', the
     default, indicates that no ordering of eigenvalues should occur.
     The leading `k' columns of `u' always span the `a'-invariant
     subspace corresponding to the `k' leading eigenvalues of `s'.

 - Loadable Function: S = svd (A)
 - Loadable Function: [U, S, V] = svd (A)
     Compute the singular value decomposition of A

          a = u * sigma * v'

     The function `svd' normally returns the vector of singular values.
     If asked for three return values, it computes U, S, and V.  For
     example,

          svd (hilb (3))

     returns

          ans =
          
            1.4083189
            0.1223271
            0.0026873

     and

          [u, s, v] = svd (hilb (3))

     returns

          u =
          
            -0.82704   0.54745   0.12766
            -0.45986  -0.52829  -0.71375
            -0.32330  -0.64901   0.68867
          
          s =
          
            1.40832  0.00000  0.00000
            0.00000  0.12233  0.00000
            0.00000  0.00000  0.00269
          
          v =
          
            -0.82704   0.54745   0.12766
            -0.45986  -0.52829  -0.71375
            -0.32330  -0.64901   0.68867

     If given a second argument, `svd' returns an economy-sized
     decomposition, eliminating the unnecessary rows or columns of U or
     V.


File: octave,  Node: Functions of a Matrix,  Prev: Matrix Factorizations,  Up: Linear Algebra

Functions of a Matrix
=====================

 - Loadable Function:  expm (A)
     Returns the exponential of a matrix, defined as the infinite
     Taylor series

          expm(a) = I + a + a^2/2! + a^3/3! + ...

     The Taylor series is *not* the way to compute the matrix
     exponential; see Moler and Van Loan, `Nineteen Dubious Ways to
     Compute the Exponential of a Matrix', SIAM Review, 1978.  This
     routine uses Ward's diagonal Pade' approximation method with three
     step preconditioning (SIAM Journal on Numerical Analysis, 1977).

     Diagonal Pade'  approximations are rational polynomials of matrices

               -1
          D (a)   N (a)

     whose Taylor series matches the first `2q+1' terms of the Taylor
     series above; direct evaluation of the Taylor series (with the
     same preconditioning steps) may be desirable in lieu of the Pade'
     approximation when `Dq(a)' is ill-conditioned.

 - Loadable Function:  logm (A)
     Compute the matrix logarithm of the square matrix A.  Note that
     this is currently implemented in terms of an eigenvalue expansion
     and needs to be improved to be more robust.

 - Loadable Function:  sqrtm (A)
     Compute the matrix square root of the square matrix A.  Note that
     this is currently implemented in terms of an eigenvalue expansion
     and needs to be improved to be more robust.

 - Function File:  kron (A, B)
     Form the kronecker product of two matrices, defined block by block
     as

          x = [a(i, j) b]

 - Function File:  qzhess (A, B)
     Compute the Hessenberg-triangular decomposition of the matrix
     pencil `(a, b)'.  This function returns `aa = q * a * z', `bb = q
     * b * z', `q', `z' orthogonal.  For example,

          [aa, bb, q, z] = qzhess (a, b)

     The Hessenberg-triangular decomposition is the first step in Moler
     and Stewart's QZ decomposition algorithm.  (The QZ decomposition
     will be included in a later release of Octave.)

     Algorithm taken from Golub and Van Loan, `Matrix Computations, 2nd
     edition'.

 - Loadable Function:  qzval (A, B)
     Compute generalized eigenvalues.

 - Loadable Function:  syl (A, B, C)
     Solve the Sylvester equation

          A X + X B + C = 0

     using standard LAPACK subroutines.


File: octave,  Node: Polynomial Manipulations,  Next: Nonlinear Equations,  Prev: Linear Algebra,  Up: Top

Polynomial Manipulations
************************

   In Octave, a polynomial is represented by its coefficients (arranged
in descending order).  For example, a vector  $c$ of length  N-th
order polynomial

     p(x) = C(1) x^N + ... + C(N) x + C(N+1).

 - Function File:  compan (C)
     Compute the companion matrix corresponding to polynomial
     coefficient vector C.

     The companion matrix is

               _                                                        _
              |  -c(2)/c(1)   -c(3)/c(1)  ...  -c(N)/c(1)  -c(N+1)/c(1)  |
              |       1            0      ...       0             0      |
              |       0            1      ...       0             0      |
          A = |       .            .   .            .             .      |
              |       .            .       .        .             .      |
              |       .            .           .    .             .      |
              |_      0            0      ...       1             0     _|

     The eigenvalues of the companion matrix are equal to the roots of
     the polynomial.

 - Function File:  conv (A, B)
     Convolve two vectors.

     `y = conv (a, b)' returns a vector of length equal to `length (a)
     + length (b) - 1'.  If A and B are polynomial coefficient vectors,
     `conv' returns the coefficients of the product polynomial.

 - Function File:  deconv (Y, A)
     Deconvolve two vectors.

     `[b, r] = deconv (y, a)' solves for B and R such that `y = conv
     (a, b) + r'.

     If Y and A are polynomial coefficient vectors, B will contain the
     coefficients of the polynomial quotient and R will be a remander
     polynomial of lowest order.

 - Function File:  poly (A)
     If A is a square N-by-N matrix, `poly (A)' is the row vector of
     the coefficients of `det (z * eye (N) - a)', the characteristic
     polynomial of A.  If X is a vector, `poly (X)' is a vector of
     coefficients of the polynomial whose roots are the elements of X.

 - Function File:  polyderiv (C)
     Returns the coefficients of the derivative of the polynomial whose
     coefficients are given by vector C.

 - Function File:  polyinteg (C)
     Returns the coefficients of the integral the polynomial whose
     coefficients are represented by the vector C.

     The constant of integration is set to zero.

 - Function File:  polyreduce (C)
     Reduces a polynomial coefficient vector to a minimum number of
     terms by stripping off any leading zeros.

 - Function File:  polyval (C, X)
     Evaluate a polynomial.

     `polyval (C, X)' will evaluate the polynomial at the specified
     value of X.

     If X is a vector or matrix, the polynomial is evaluated at each of
     the elements of X.

 - Function File:  polyvalm (C, X)
     Evaluate a polynomial in the matrix sense.

     `polyvalm (C, X)' will evaluate the polynomial in the matrix
     sense, i.e. matrix multiplication is used instead of element by
     element multiplication as is used in polyval.

     The argument X must be a square matrix.

 - Function File:  residue (B, A, TOL)
     If B and A are vectors of polynomial coefficients, then residue
     calculates the partial fraction expansion corresponding to the
     ratio of the two polynomials.

     The function `residue' returns R, P, K, and E, where the vector R
     contains the residue terms, P contains the pole values, K contains
     the coefficients of a direct polynomial term (if it exists) and E
     is a vector containing the powers of the denominators in the
     partial fraction terms.

     Assuming B and A represent polynomials  P (s) and Q(s)  we have:

           P(s)    M       r(m)         N
           ---- = SUM -------------  + SUM k(i)*s^(N-i)
           Q(s)   m=1 (s-p(m))^e(m)    i=1

     where M is the number of poles (the length of the R, P, and E
     vectors) and N is the length of the K vector.

     The argument TOL is optional, and if not specified, a default
     value of 0.001 is assumed.  The tolerance value is used to
     determine whether poles with small imaginary components are
     declared real.  It is also used to determine if two poles are
     distinct.  If the ratio of the imaginary part of a pole to the
     real part is less than TOL, the imaginary part is discarded.  If
     two poles are farther apart than TOL they are distinct.  For
     example,

     Example:

           b = [1, 1, 1];
           a = [1, -5, 8, -4];
          
           [r, p, k, e] = residue (b, a)

     returns

          r = [-2, 7, 3]
          
          p = [2, 2, 1]
          
          k = [](0x0)
          
          e = [1, 2, 1]

     which implies the following partial fraction expansion

                  s^2 + s + 1       -2        7        3
             ------------------- = ----- + ------- + -----
             s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)

 - Function File:  roots (V)
     For a vector V with N components, return the roots of the
     polynomial

          v(1) * z^(N-1) + ... + v(N-1) * z + v(N).


File: octave,  Node: Nonlinear Equations,  Next: Differential Equations,  Prev: Polynomial Manipulations,  Up: Top

Nonlinear Equations
*******************

   Octave can solve sets of nonlinear equations of the form

     F (x) = 0

using the function `fsolve', which is based on the MINPACK subroutine
`hybrd'.

 - Loadable Function: [X, INFO] = fsolve (FCN, X0)
     Given FCN, the name of a function of the form `f (X)' and an
     initial starting point X0, `fsolve' solves the set of equations
     such that `f(X) == 0'.

 - Loadable Function:  fsolve_options (OPT, VAL)
     When called with two arguments, this function allows you set
     options parameters for the function `fsolve'.  Given one argument,
     `fsolve_options' returns the value of the corresponding option.  If
     no arguments are supplied, the names of all the available options
     and their current values are displayed.

   Here is a complete example.  To solve the set of equations

     -2x^2 + 3xy   + 4 sin(y) = 6
      3x^2 - 2xy^2 + 3 cos(x) = -4

you first need to write a function to compute the value of the given
function.  For example:

     function y = f (x)
     
       y(1) = -2*x(1)^2 + 3*x(1)*x(2)   + 4*sin(x(2)) - 6;
       y(2) =  3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4;
     
     endfunction

   Then, call `fsolve' with a specified initial condition to find the
roots of the system of equations.  For example, given the function `f'
defined above,

     [x, info] = fsolve ("f", [1; 2])

results in the solution

     x =
     
       0.57983
       2.54621
     
     info = 1

   A value of `info = 1' indicates that the solution has converged.

   The function `perror' may be used to print English messages
corresponding to the numeric error codes.  For example,

     perror ("fsolve", 1)

prints

     solution converged to requested tolerance


File: octave,  Node: Differential Equations,  Next: Optimization,  Prev: Nonlinear Equations,  Up: Top

Differential Equations
**********************

   Octave has two built-in functions for solving differential equations.
Both are based on reliable ODE solvers written in Fortran.

* Menu:

* Ordinary Differential Equations::
* Differential-Algebraic Equations::


File: octave,  Node: Ordinary Differential Equations,  Next: Differential-Algebraic Equations,  Prev: Differential Equations,  Up: Differential Equations

Ordinary Differential Equations
===============================

   The function `lsode' can be used Solve ODEs of the form

     dx
     -- = f (x, t)
     dt

using Hindmarsh's ODE solver LSODE.

 - Lodable Function:  lsode (FCN, X0, T_OUT, T_CRIT)
     The first argument is the name of the function to call to compute
     the vector of right hand sides.  It must have the form

          XDOT = f (X, T)

     where XDOT and X are vectors and T is a scalar.

     The second argument specifies the initial condition, and the third
     specifies a vector of output times at which the solution is
     desired, including the time corresponding to the initial condition.

     The fourth argument is optional, and may be used to specify a set
     of times that the ODE solver should not integrate past.  It is
     useful for avoiding difficulties with singularities and points
     where there is a discontinuity in the derivative.

     Here is an example of solving a set of two differential equations
     using `lsode'.  The function

          function xdot = f (x, t)
          
            r = 0.25;
            k = 1.4;
            a = 1.5;
            b = 0.16;
            c = 0.9;
            d = 0.8;
          
            xdot(1) = r*x(1)*(1 - x(1)/k) - a*x(1)*x(2)/(1 + b*x(1));
            xdot(2) = c*a*x(1)*x(2)/(1 + b*x(1)) - d*x(2);
          
          endfunction

     is integrated with the command

          x = lsode ("f", [1; 2], (t = linspace (0, 50, 200)'));

     producing a set of 200 values stored in the variable X.  Note that
     this example takes advantage of the fact that an assignment
     produces a value to store the values of the output times in the
     variable T directly in the function call   The results can then be
     plotted using the command

          plot (t, x)

 - Lodable Function:  lsode_options (OPT, VAL)
     When called with two arguments, this function allows you set
     options parameters for the function `lsode'.  Given one argument,
     `lsode_options' returns the value of the corresponding option.  If
     no arguments are supplied, the names of all the available options
     and their current values are displayed.

   See Alan C. Hindmarsh, `ODEPACK, A Systematized Collection of ODE
Solvers', in Scientific Computing, R. S. Stepleman, editor, (1983) for
more information about this family of ODE solvers.


File: octave,  Node: Differential-Algebraic Equations,  Prev: Ordinary Differential Equations,  Up: Differential Equations

Differential-Algebraic Equations
================================

   The function `dassl' can be used Solve DAEs of the form

     0 = f (x-dot, x, t),    x(t=0) = x_0, x-dot(t=0) = x-dot_0

 - Loadable Function: [X, XDOT] = dassl (FCN, X_0, XDOT_0, T_OUT,
          T_CRIT)
     The first argument is the name of the function to call to compute
     the vector of residuals.  It must have the form

          RES = f (X, XDOT, T)

     where X, XDOT, and RES are vectors, and T is a scalar.

     The second and third arguments to `dassl' specify the initial
     condition of the states and their derivatives, and the fourth
     argument specifies a vector of output times at which the solution
     is desired, including the time corresponding to the initial
     condition.

     The set of initial states and derivatives are not strictly
     required to be consistent.  In practice, however, DASSL is not
     very good at determining a consistent set for you, so it is best
     if you ensure that the initial values result in the function
     evaluating to zero.

     The fifth argument is optional, and may be used to specify a set of
     times that the DAE solver should not integrate past.  It is useful
     for avoiding difficulties with singularities and points where
     there is a discontinuity in the derivative.

 - Loadable Function:  dassl_options (OPT, VAL)
     When called with two arguments, this function allows you set
     options parameters for the function `lsode'.  Given one argument,
     `dassl_options' returns the value of the corresponding option.  If
     no arguments are supplied, the names of all the available options
     and their current values are displayed.

   See K. E. Brenan, et al., `Numerical Solution of Initial-Value
Problems in Differential-Algebraic Equations', North-Holland (1989) for
more information about the implementation of DASSL.


File: octave,  Node: Optimization,  Next: Quadrature,  Prev: Differential Equations,  Up: Top

Optimization
************

* Menu:

* Quadratic Programming::
* Nonlinear Programming::
* Linear Least Squares::


File: octave,  Node: Quadratic Programming,  Next: Nonlinear Programming,  Prev: Optimization,  Up: Optimization

Quadratic Programming
=====================

 - Loadable Function: [X, OBJ, INFO, LAMBDA] = qpsol (X, H, C, LB, UB,
          LB, A, UB)
     Solve quadratic programs using Gill and Murray's QPSOL.  Because
     QPSOL is not freely redistributable, this function is only
     available if you have obtained your own copy of QPSOL.  *Note
     Installation::.

 - Loadable Function:  qpsol_options (OPT, VAL)
     When called with two arguments, this function allows you set
     options parameters for the function `qpsol'.  Given one argument,
     `qpsol_options' returns the value of the corresponding option.  If
     no arguments are supplied, the names of all the available options
     and their current values are displayed.


File: octave,  Node: Nonlinear Programming,  Next: Linear Least Squares,  Prev: Quadratic Programming,  Up: Optimization

Nonlinear Programming
=====================

 - Loadable Function: [X, OBJ, INFO, LAMBDA] = npsol (X, PHI, LB, UB,
          LB, A, UB, LB, G, UB)
     Solve nonlinear programs using Gill and Murray's NPSOL.  Because
     NPSOL is not freely redistributable, this function is only
     available if you have obtained your own copy of NPSOL.  *Note
     Installation::.

     The second argument is a string containing the name of the
     objective function to call.  The objective function must be of the
     form

          y = phi (x)

     where x is a vector and y is a scalar.

 - Loadable Function:  npsol_options (OPT, VAL)
     When called with two arguments, this function allows you set
     options parameters for the function `npsol'.  Given one argument,
     `npsol_options' returns the value of the corresponding option.  If
     no arguments are supplied, the names of all the available options
     and their current values are displayed.


File: octave,  Node: Linear Least Squares,  Prev: Nonlinear Programming,  Up: Optimization

Linear Least Squares
====================

 - Function File:  gls (Y, X, O)
     Generalized least squares (GLS) estimation for the multivariate
     model

          Y = X * B + E,  mean(E) = 0,  cov(vec(E)) = (s^2)*O

     with

          Y an T x p matrix
          X an T x k matrix
          B an k x p matrix
          E an T x p matrix
          O an Tp x Tp matrix

     Each row of Y and X is an observation and each column a variable.

     Returns BETA, v, and, R, where BETA is the GLS estimator for B, v
     is the GLS estimator for s^2, and R = Y - X*BETA is the matrix of
     GLS residuals.

 - Function File:  ols (Y, X)
     Ordinary Least Squares (OLS) estimation for the multivariate model

          Y = X*B + E,  mean (E) = 0,  cov (vec (E)) = kron (S, I)

     with

          Y an T x p matrix
          X an T x k matrix
          B an k x p matrix
          E an T x p matrix

     Each row of Y and X is an observation and each column a variable.

     Returns BETA, SIGMA, and R, where BETA is the OLS estimator for B,
     i.e.

          BETA = pinv(X)*Y,

     where pinv(X) denotes the pseudoinverse of X, SIGMA is the OLS
     estimator for the matrix S, i.e.

          SIGMA = (Y - X*BETA)'*(Y - X*BETA) / (T - rank(X))

     and R = Y - X*BETA is the matrix of OLS residuals.


File: octave,  Node: Quadrature,  Next: Control Theory,  Prev: Optimization,  Up: Top

Quadrature
**********

* Menu:

* Functions of one Variable::
* Orthogonal Collocation::


File: octave,  Node: Functions of one Variable,  Next: Orthogonal Collocation,  Prev: Quadrature,  Up: Quadrature

Functions of one Variable
=========================

 - Loadable Function: [V, IER, NFUN] = quad (F, A, B, TOL, SING)
     Integrate a nonlinear function of one variable using Quadpack.
     The first argument is the name of the  function to call to compute
     the value of the integrand.  It must have the form

          y = f (x)

     where Y and X are scalars.

     The second and third arguments are limits of integration.  Either
     or both may be infinite.

     The optional argument TOL is a vector that specifies the desired
     accuracy of the result.  The first element of the vector is the
     desired absolute tolerance, and the second element is the desired
     relative tolerance.  To choose a relative test only, set the
     absolute tolerance to zero.  To choose an absolute test only, set
     the relative tolerance to zero.

     The optional argument SING is a vector of values at which the
     integrand is known to be singular.

 - Loadable Function:  quad_options (OPT, VAL)
     When called with two arguments, this function allows you set
     options parameters for the function `quad'.  Given one argument,
     `quad_options' returns the value of the corresponding option.  If
     no arguments are supplied, the names of all the available options
     and their current values are displayed.


File: octave,  Node: Orthogonal Collocation,  Prev: Functions of one Variable,  Up: Quadrature

Orthogonal Collocation
======================

 - Loadable Function: [R, A, B, Q] = colloc (N, "left", "right")
     Compute derivative and integral weight matrices for orthogonal
     collocation using the subroutines given in J. Villadsen and M. L.
     Michelsen, `Solution of Differential Equation Models by Polynomial
     Approximation'.


File: octave,  Node: Control Theory,  Next: Signal Processing,  Prev: Quadrature,  Up: Top

Control Theory
**************

 - Function File:  abcddim (A, B, C, D)
     Check for compatibility of the dimensions of the matrices defining
     the linear system [A, B, C, D] corresponding to

          dx/dt = a x + b u
          y = c x + d u

     or a similar discrete-time system.

     If the matrices are compatibly dimensioned, then `abcddim' returns
     N = number of system states, M = number of system inputs, and P =
     number of system outputs.  Otherwise `abcddim' returns N = M = P =
     -1.

 - Function File:  are (A, B, C, OPT)
     Returns the solution, X, of the algebraic Riccati equation

          a' x + x a - x b x + c = 0

     for identically dimensioned square matrices A, B, C.  If B (C) is
     not square, then the function attempts to use `B*B'' (`C'*C')
     instead.

     Solution method: apply Laub's Schur method (IEEE Transactions on
     Automatic Control, 1979) to the appropriate Hamiltonian matrix.

     OPT is an option passed to the eigenvalue balancing routine.
     Default is `"B"'.

 - Function File:  c2d (A, B, T)
     Converts the continuous time system described by:

          dx/dt = a x + b u

     into a discrete time equivalent model

          x[k+1] = Ad x[k] + Bd u[k]

     via the matrix exponential assuming a zero-order hold on the input
     and sample time T.

 - Function File:  dare (A, B, C, R, OPT)
     Returns the solution, X of the discrete-time algebraic Riccati
     equation

          a' x a - x + a' x b (r + b' x b)^(-1) b' x a + c = 0

     for matrices with dimensions:

          A: N by N
          B: N by M
          C: N by N, symmetric positive semidefinite
          R: M by M, symmetric positive definite (invertible)

     If C is not square, then the function attempts to use `C'*C'
     instead.

     Solution method: Laub's Schur method (IEEE Transactions on
     Automatic Control, 1979) is applied to the appropriate symplectic
     matrix.

     See also: Ran and Rodman, `Stable Hermitian Solutions of Discrete
     Algebraic Riccati Equations', Mathematics of Control, Signals and
     Systems, Volume 5, Number 2 (1992).

     OPT is an option passed to the eigenvalue balancing routine.  The
     default is `"B"'.

 - Function File:  dgram (A, B)
     Returns the discrete controllability and observability gramian for
     the discrete time system described by

          x[k+1] = A x[k] + B u[k]
            y[k] = C x[k] + D u[k]

     `dgram (A, B)' returns the discrete controllability gramian and
     `dgram (A', C')' returns the observability gramian.

 - Function File:  dlqe (A, G, C, SIGW, SIGV [, Z])
     Linear quadratic estimator (Kalman filter) design for the discrete
     time system

          x[k+1] = A x[k] + B u[k] + G w[k]
            y[k] = C x[k] + D u[k] + w[k]

     where W, V are zero-mean gaussian noise processes with respective
     intensities `SIGW = cov (W, W)' and `SIGV = cov (V, V)'.

     If specified, Z is `cov (W, V)'.  Otherwise `cov (W, V) = 0'.

     The observer structure is

          z[k+1] = A z[k] + B u[k] + k(y[k] - C z[k] - D u[k])

     Returns:

     L is the observer gain, `(A - A L C)' is stable.

     M is the Ricatti equation solution.

     P is the estimate error covariance after the measurement update.

     E are the closed loop poles of `(A - A L C)'.

 - Function File:  dlqr (A, B, Q, R [, Z])
     Linear quadratic regulator design for the discrete time system

          x[k+1] = A x[k] + B u[k]

     to minimize the cost functional

          J = Sum [ x' Q x + u' R u ],              Z omitted

     or

          J = Sum [ x' Q x + u' R u +2 x' Z u ],    Z included

     Returns:

     K is the state feedback gain, `(A - B K)' is stable.

     P is the solution of algebraic Riccati equation.

     E are the closed loop poles of (A - B K).

 - Function File:  dlyap (A, B)
     Solve the discrete-time Lyapunov equation

          a x a' - x + b = 0

     for square matrices A, B.  If B is not square, then the function
     attempts to solve either

          a x a' - x + b b' = 0

     or

          a' x a - x + b' b = 0

     whichever is appropriate.

     Uses Schur decomposition method as in Kitagawa, International
     Journal of Control (1977); column-by-column solution method as
     suggested in Hammarling, IMA Journal of Numerical Analysis, (1982).

 - Function File:  is_controllable (A, B, TOL)
     If the pair (a, b) is controllable, then return value 1.
     Otherwise, returns a value of 0.

     TOL is a roundoff parameter, set to `2*eps' if omitted.

     Currently just constructs the controllability matrix and checks
     rank.  A better method is as follows (Boley and Golub, Systems and
     Control Letters, 1984):  Controllability is determined by applying
     Arnoldi iteration with complete re-orthogonalization to obtain an
     orthogonal basis of the Krylov subspace

                                n-1
           span ([b, a*b, ..., a^   *b]).

 - Function File:  is_observable (A, C, TOL)
     Returns 1 if the pair `(a, c)' is observable.  Otherwise, returns
     a value of 0.

 - Function File:  lqe (A, G, C, SIGW, SIGV, Z)
          [k, p, e] = lqe (a, g, c, sigw, sigv, z)

     Linear quadratic estimator (Kalman filter) design for the
     continuous time system

          dx
          -- = a x + b u
          dt
          
          y = c x + d u

     where W, V are zero-mean gaussian noise processes with respective
     intensities

          sigw = cov (w, w)
          sigv = cov (v, v)

     Z (if specified) is the cross-covariance `cov (W, V)'; the default
     value is `cov (W, V) = 0'.

     Observer structure is `dz/dt = A z + B u + k (y - C z - D u)'

     returns:

     K is observer gain:  `(A - K C)' is stable.

     P is solution of algebraic Riccati equation.

     E is the vector of closed loop poles of `(A - K C)'.

 - Function File:  lqr (A, B, Q, R, Z)
 - Function File: [K, P, E] = lqr (A, B, Q, R, Z)
     Linear quadratic regulator design for the continuous time system

          dx
          -- = A x + B u
          dt

     to minimize the cost functional

                infinity
                /
            J = |  x' Q x + u' R u
               /
              t=0

     Z omitted or

                infinity
                /
            J = |  x' Q x + u' R u + 2 x' Z u
               /
              t=0

     Z included

     Returns:

     K is state feedback gain:  `(A - B K)' is stable.

     P is the stabilizing solution of appropriate algebraic Riccati
     equation.

     E is the vector of the closed loop poles of `(A - B K)'.

 - Function File:  lyap (A, B, C)
     Solve the Lyapunov (or Sylvester) equation via the Bartels-Stewart
     algorithm (Communications of the ACM, 1972).

     If `(a, b, c)' are specified, then `lyap' returns the solution of
     the  Sylvester equation

          a x + x b + c = 0

     If only `(a, b)' are specified, then `lyap' returns the solution
     of the Lyapunov equation

          a' x + x a + b = 0

     If B is not square, then `lyap' returns the solution of either

          a' x + x a + b' b = 0

     or

          a x + x a' + b b' = 0

     whichever is appropriate.

     Solves by using the Bartels-Stewart algorithm (1972).

 - Function File:  tzero (A, B, C, D, BAL)
     Compute the transmission zeros of [A, B, C, D].

     BAL = balancing option (see balance); default is `"B"'.

     Needs to incorporate `mvzero' algorithm to isolate finite zeros;
     see Hodel, `Computation of System Zeros with Balancing', Linear
     Algebra and its Applications, July 1993.

