Chapter 5
ROOTS and OPTIMA
Summary
The functions in this library segment are used to
find the roots of nonlinear equations and systems
of such equations and to locate local extrema of
functions. The areas covered are:
o Roots of Equations
o Optimization
-------------------------------------------------------------------------------
Notes on Contents
Functions in this section of the library are used to extract the roots of
nonlinear equations and to find local optima.
o Roots of Equations:
solnl -------- solve a system of nonlinear equations.
solnlx ------- solve a nonlinear system using an initial
Jacobian.
secrt -------- solve a nonlinear equation using the secant
method.
plrt --------- find the roots (real and complex) of a polynomial.
polyc -------- evaluate a polynomial for complex arguments.
o Optimization:
optmiz ------- perform an unconstrained search for the
optimal parameter vector.
optsh -------- perform a golden section search for optima in
one dimension.
-------------------------------------------------------------------------------
General Technical Comments
The general root finding functions employ a matrix update algorithm
developed by Broyden. The second form offers an enlarged domain of
convergence at the expense of requiring the Jacobian matrix of the system at
the initial point of the search. The secant method is a simple and standard
root extraction technique, provided for convenience. Roots of a polynomial
with real coefficients are computed by plrt.
The nonlinear optimal search algorithm uses the BFGS matrix update
method. The robustness of this procedure has been confirmed in numerous tests
of optimization techniques. The golden section method is useful and efficient
for bounded search in one-dimension.
Functions for finding the roots of a nonlinear system of equations
and for nonlinear optimization are based on efficient quasi-Newton matrix
update algorithms. This approach has been shown to combine high execution
efficiency (ie. a small number of function evaluations) with a large
convergence domain. Simplicity of use was also an important criterion in the
selection of the algorithms. The functions implemented do not require user
supplied procedures for computation of derivatives. Fortunately, this does
not entail a performance sacrifice, since the quasi-Newton procedures exhibit
excellent convergence properties when operating with numerical derivatives.
The algorithms employed for finding roots and optima are iterative.
They require an initial estimate of the solution as input. While the
implemented algorithms have robust convergence, refinement of these initial
estimates is often the most effective way to improve algorithm performance.
A few hours spent estimating appropriate scales for variables and refining
initial point estimates may save days of agonizing search for a "golden
algorithm". Nonlinear analysis is a field so rich that vast bodies of code
are no substitute for a thoughtful understanding of the problem at hand.
-------------------------------------------------------------------------------
FUNCTION SYNOPSES
-------------------------------------------------------------------------------
Roots of Equations:
-------------------------------------------------------------------------------
solnl
Solve a system of nonlinear equations.
int solnl(double *x,double *f,double (*fvec[])(),int n,double test)
x = pointer to array containing initial solution estimate,
converted to solution vector at exit
f = pointer to array containing final values of the system functions
{ f[k]=(*fvec[k])(x) }
fvec = array of pointers to the system functions
{ (*fvec[k])(x)=0 for k=0,1,2,--,n-1 }
n = dimension of system
test = convergence threshold (f~*f normal exit
0 -> convergence failure
-------------------------------------------------------------------
solnx
Variant of solnl, which permits input of an initial system
Jacobian. (Has an enlarged convergence domain.)
int solnx(double *x,double *f,double (*fvec)(),double *jm,int n, \
double *test)
x = pointer to array containing initial solution estimate,
converted to solution vector at exit
f = pointer to array containing final values of the system
functions { f[k]=(*fvec[k])(x) }
fvec = array of pointers to the system functions
{ (*fvec[k])(x)=0 for k=0,1,2,--,n-1 }
jm = pointer to array of dimension n*n containing an initial system
Jacobian ( row order: JM[i,j] = df[i]/dx[j] )
n = dimension of system
test = convergence threshold (f~*f normal exit
0 -> convergence failure
--------------------------------------------------------------------
plrt
Find the real and complex roots of a polynomial with real
coefficients.
typedef struct complex {double re,im;} Cpx;
int plrt(double *cof,int n,Cpx *root,double ra,double rb)
cof = pointer to array containing polynomial coefficients, with
the leading coefficient cof[n] != 0
root = pointer to array of complex structures, containing the
polynomial roots at exit (dimension n)
n = degree of polynomial
ra,rb = initial root estimates, with
( rb>0 -> real = ra, imag = rb and
rb<0 -> real roots ra+rb, ra-rb )
return: 0 -> normal exit
m>0 -> convergence failure for roots k0 -> success (iteration count)
m=0 -> convergence failure
---------------------------------------------------------------------
optsh
Conduct a golden section line search for a function minima.
double optsch(double (*func)(),double a,double b,double test)
func = pointer to function to be minimized
a,b = bounds on function argument (a<=x<=b)
test = convergence threshold (length of section
containing minima