Octave Higher Math Functions
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)
-| solution converged to requested tolerance
Quadrature
**********
Functions of One Variable
=========================
- Loadable Function: [V, IER, NFUN, ERR] = 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.
The result of the integration is returned in V and IER contains an
integer error code (0 indicates a successful integration). The
value of NFUN indicates how many function evaluations were
required, and ERR contains an estimate of the error in the
solution.
- 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.
Here is an example of using `quad' to integrate the function
F(X) = X * sin (1/X) * sqrt (abs (1 - X))
from X = 0 to X = 3.
This is a fairly difficult integration (plot the function over the
range of integration to see why).
The first step is to define the function:
function y = f (x)
y = x .* sin (1 ./ x) .* sqrt (abs (1 - x));
endfunction
Note the use of the `dot' forms of the operators. This is not
necessary for the call to `quad', but it makes it much easier to
generate a set of points for plotting (because it makes it possible to
call the function with a vector argument to produce a vector result).
Then we simply call quad:
[v, ier, nfun, err] = quad ("f", 0, 3)
=> 1.9819
=> 1
=> 5061
=> 1.1522e-07
Although `quad' returns a nonzero value for IER, the result is
reasonably accurate (to see why, examine what happens to the result if
you move the lower bound to 0.1, then 0.01, then 0.001, etc.).
Orthogonal Collocation
======================
- Loadable Function: [R, AMAT, BMAT, 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'.
Here is an example of using `colloc' to generate weight matrices for
solving the second order differential equation U' - ALPHA * U" = 0 with
the boundary conditions U(0) = 0 and U(1) = 1.
First, we can generate the weight matrices for N points (including
the endpoints of the interval), and incorporate the boundary conditions
in the right hand side (for a specific value of ALPHA).
n = 7;
alpha = 0.1;
[r, a, b] = colloc (n-2, "left", "right");
at = a(2:n-1,2:n-1);
bt = b(2:n-1,2:n-1);
rhs = alpha * b(2:n-1,n) - a(2:n-1,n);
Then the solution at the roots R is
u = [ 0; (at - alpha * bt) \ rhs; 1]
=> [ 0.00; 0.004; 0.01 0.00; 0.12; 0.62; 1.00 ]
Differential Equations
**********************
Octave has two built-in functions for solving differential equations.
Both are based on reliable ODE solvers written in Fortran.
Ordinary Differential Equations
===============================
The function `lsode' can be used to solve ODEs of the form
dx
-- = f (x, t)
dt
using Hindmarsh's ODE solver LSODE.
- Loadable Function: lsode (FCN, X0, T, T_CRIT)
Return a matrix of X as a function of T, given the initial state
of the system X0. Each row in the result matrix corresponds to
one of the elements in the vector T. The first element of T
corresponds to the initial state X0, so that the first row of the
output is X0.
The first argument, FCN, is a string that names the function to
call to compute the vector of right hand sides for the set of
equations. It must have the form
XDOT = f (X, T)
where XDOT and X are vectors and T is a scalar.
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 three differential equations
using `lsode'. Given the function
function xdot = f (x, t)
xdot = zeros (3,1);
xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) \
- 8.375e-06*x(1)^2);
xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
xdot(3) = 0.161*(x(1) - x(3));
endfunction
and the initial condition `x0 = [ 4; 1.1; 4 ]', the set of equations
can be integrated using the command
t = linspace (0, 500, 1000);
y = lsode ("f", x0, t);
If you try this, you will see that the value of the result changes
dramatically between T = 0 and 5, and again around T = 305. A more
efficient set of output points might be
t = [0, logspace (-1, log10(303), 150), \
logspace (log10(304), log10(500), 150)];
- Loadable 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 the inner workings of `lsode'.
Differential-Algebraic Equations
================================
The function `dassl' can be used to solve DAEs of the form
0 = f (x-dot, x, t), x(t=0) = x_0, x-dot(t=0) = x-dot_0
using Petzold's DAE solver DASSL.
- Loadable Function: [X, XDOT] = dassl (FCN, X0, XDOT0, T, T_CRIT)
Return a matrix of states and their first derivatives with respect
to T. Each row in the result matrices correspond to one of the
elements in the vector T. The first element of T corresponds to
the initial state X0 and derivative XDOT0, so that the first row
of the output X is X0 and the first row of the output XDOT is
XDOT0.
The first argument, FCN, is a string that names the function to
call to compute the vector of residuals for the set of equations.
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.
Optimization
************
Quadratic Programming
=====================
Nonlinear Programming
=====================
Linear Least Squares
====================
- Function File: [BETA, V, R] = gls (Y, X, O)
Generalized least squares estimation for the multivariate model y
= x b + e with mean (e) = 0 and cov (vec (e)) = (s^2) o, where y
is a t by p matrix, x is a t by k matrix, b is a k by p matrix, e
is a t by p matrix, and o is a t p by t p matrix.
Each row of Y and X is an observation and each column a variable.
The return values BETA, V, and R are defined as follows.
BETA
The GLS estimator for b.
V
The GLS estimator for s^2.
R
The matrix of GLS residuals, r = y - x beta.
- Function File: [BETA, SIGMA, R] = ols (Y, X)
Ordinary least squares estimation for the multivariate model y = x
b + e with mean (e) = 0 and cov (vec (e)) = kron (s, I). where y
is a t by p matrix, x is a t by k matrix, b is a k by p matrix, and
e is a t by p matrix.
Each row of Y and X is an observation and each column a variable.
The return values BETA, SIGMA, and R are defined as follows.
BETA
The OLS estimator for B, `BETA = pinv (X) * Y', where `pinv
(X)' denotes the pseudoinverse of X.
SIGMA
The OLS estimator for the matrix S,
SIGMA = (Y-X*BETA)'
* (Y-X*BETA)
/ (T-rank(X))
R
The matrix of OLS residuals, `R = Y - X * BETA'.