eqSolvemt
==============================================
Purpose
----------------
Solves a system of nonlinear equations.
Format
----------------
.. function:: out = eqSolvemt(&fct, par, ..., c)
:param &fct: pointer to a procedure that computes the
function to be minimized. This procedure must have an instance of a :class:`PV` structure
containing the parameters as the first input. Extra data needed by this function may
be passed in after the :class:`PV` structure, but before the control structure.
This procedure must have one output argument, a column vector
containing the result of each equation.
:type &fct: scalar function pointer
:param par: The *par* instance is passed to the user-provided procedure pointed to by *&fct*. *par* is
constructed using the :func:`pvPack` functions.
:type par: an instance of a :class:`PV` structure
:param ....: Optional extra arguments.
These arguments are passed untouched to the user-provided objective function, by :func:`eqSolveMT`.
:param c: Optional. Normally an instance is initialized by calling
:func:`eqSolvemtControlCreate` and members of this instance
can be set to other values by the user. For an instance named
*c*, the members are:
.. list-table::
:widths: auto
* - *c.jacobianProc*
- pointer to a procedure which computes the analytical Jacobian. By default, :func:`eqSolvemt` will compute the Jacobian numerically.
* - *c.maxIters*
- scalar, the maximum number of iterations. Default = 100.
* - *c.stepTolerance*
- scalar, the step tolerance. Default = :math:`macheps^{2/3}`
* - *c.typicalF*
- Kx1 vector of the typical :code:`fct(x)` values at a point not near a root, used for scaling. This becomes important when the magnitudes of the components of :code:`fct(x)` are expected to be very different. By default, function values are not scaled.
* - *c.typicalX*
- Kx1 vector of the typical magnitude of *x*, used for scaling. This becomes important when the magnitudes of the components of *x* are expected to be very different. By default, variable values are not scaled.
* - *c.printIters*
- scalar, if nonzero, iteration information is printed. Default = 0.
* - *c.tolerance*
- scalar, the tolerance of the scalar function :math:`f = 0.5*||fct(X)||2` required to terminate the algorithm. That is, the condition that :math:`|f(x)| <= c.tolerance` must be met before that algorithm can terminate successfully. Default = 1e-5.
* - *c.altNames*
- Kx1 string array of alternate names to be used by the printed output. By default, the names :code:`X1,X2,X3...` will be used.
* - *c.title*
- string, printed as a title in output.
* - *c.output*
- scalar. If non-zero, final results are printed.
:type c: an instance of an :class:`eqSolvemtControl` structure.
:return out: instance of :class:`eqSolvemtOut` struct. For an instance named *out*, the members are:
.. list-table::
:widths: auto
* - *out.par*
- an instance of a :class:`PV` structure containing the parameter estimates.
* - *out.fct*
- scalar, function evaluated at *x*
* - *out.retcode*
- scalar, return code:
:-1: Jacobian is singular.
:1: Norm of the scaled function value is less than *c.tolerance*. *x* given is an approximate root of :code:`fct(x)` (unless c.tolerance is too large).
:2: The scaled distance between the last two steps is less than the step-tolerance (*c.stepTolerance*). *x* may be an approximate root of :code:`fct(x)`, but it is also possible that the algorithm is making very slow progress and is not near a root, or the step-tolerance is too large.
:3: The last global step failed to decrease :code:`norm2(fct(x))` sufficiently; either *x* is close to a root of :code:`fct(x)` and no more accuracy is possible, or an incorrectly coded analytic Jacobian is being used, or the secant approximation to the Jacobian is inaccurate, or the step-tolerance is too large.
:4: Iteration limit exceeded.
:5: Five consecutive steps of maximum step length have been taken; either :code:`norm2(fct(x))` asymptotes from above to a finite value in some direction or the maximum step length is too small.
:6: *x* seems to be an approximate local minimizer of :code:`norm2(fct(x))` that is not a root of :code:`fct(x)`. To find a root of :code:`fct(x)`, restart :func:`eqSolvemt` from a different region.
:rtype out: struct
Examples
----------------
Basic usage
+++++++++++
For example, consider a two-equation system given by:
.. math:: x_1^2 + x_2^2 - 5 = 0
.. math:: e^{x_1-1} + x_2^3 - 5 = 0
::
new;
// Declare 'par' to be an instance of a PV vector
struct PV par;
/*
** Create default PV struct and add a parameter
** named 'x1' with a starting value of 1
*/
par = pvPack(pvCreate(), 1, "x1");
/*
** Add a parameter named 'x2' to 'par'
** with a starting value of 1
*/
par = pvPack(par, 1, "x2");
/*
** Solve the system of equations
** and print the output to the screen
*/
call eqSolvemt(&fct, par);
// The definition of the function to be minimized
proc fct(struct PV p);
local x1, x2, z;
x1 = pvUnpack(p, "x1");
x2 = pvUnpack(p, "x2");
z = (x1^2 + x2^2 - 5) | (exp(x1 - 1) + x2^3 - 5);
retp(z);
endp;
After the code above, a short report will be printed to the **Command** window. Part of the output is displayed below:
::
--------------------------------------------------------------------
VARIABLE START ROOTS F(ROOTS)
--------------------------------------------------------------------
X1 1.00000 1.7146639625 0.0000000001
X2 1.00000 1.4352447511 0.0000000002
--------------------------------------------------------------------
Using control and output structures
+++++++++++++++++++++++++++++++++++
Again, consider a two-equation system given by:
.. math:: x_1^2 + x_2^2 - 5 = 0
.. math:: e^{x_1-1} + x_2^3 - 5 = 0
::
new;
// Declare control structure and fill with defaults
struct eqSolvemtControl c;
c = eqSolvemtControlCreate();
// Turn on printing of iteration information
c.printIters = 1;
// Assign variable names printed output
c.altNames = "alpha" $| "beta";
// Declare 'par' to be an instance of a PV vector
struct PV par;
/*
** Create default PV struct and add a parameter
** named 'x1' with a starting value of 1
*/
par = pvPack(pvCreate(),1, "x1");
/*
** Add a parameter named 'x2' to 'par'
** with a starting value of 1
*/
par = pvPack(par, 1, "x2");
// Declare output structure to hold results
struct eqSolvemtOut out;
// Solve the system of equations
out = eqSolvemt(&fct, par, c);
// The definition of the function to be minimized
proc fct(struct PV p);
local x1, x2, z;
x1 = pvUnpack(p, "x1");
x2 = pvUnpack(p, "x2");
z = (x1^2 + x2^2 - 5) | (exp(x1 - 1) + x2^3 - 5);
retp(z);
endp;
The code above will print out a report similar to the previous example. Notice that the variable names in the report are what we assigned to the *altNames* member of the control structure.
::
--------------------------------------------------------------------------------
VARIABLE START ROOTS F(ROOTS)
--------------------------------------------------------------------------------
alpha 1.00000 1.7146639625 0.0000000001
beta 1.00000 1.4352447511 0.0000000002
--------------------------------------------------------------------------------
The parameter values returned by :func:`eqSolveMT` are located in the par member of the :class:`eqsolveMTOut` struct. They can be accessed with :func:`pvGetParVector` or :func:`pvUnpack` like this:
::
// Return the values of 'x1' and 'x2' as a 2x1 vector
x_all = pvGetParVector(out.par);
// Return the value of 'x1'
x1 = pvUnpack(out.par, "x1");
// Return the value of 'x2'
x2 = pvUnpack(out.par, "x2");
Remarks
-------
The equation procedure should return a column vector containing the
result for each equation. For example, consider a two-equation system given by:
.. math:: x_1^2 + x_2^2 - 2 = 0
.. math:: e^{x_1-1} + x_2^3 - 2 = 0
::
proc (1) = f(struct PV p);
local x1, x2, eqns;
// Set x1
x1 = pvUnpack(p, "x1");
// Set x2
x2 = pvUnpack(p, "x2");
// Preallocate output vector
eqns = { 0, 0 };
// Equation 1
eqns[1] = x1^2 + x2^2 - 2;
// Equation 2
eqns[2] = exp(x1-1) + x2^3 - 2;
retp(eqns);
endp;
Note that the first equation in the system is contained in the first row of ``eqns`` and the second equation is contained in the second row of ``eqns``.
Source
------
eqsolvemt.src
.. seealso:: Functions :func:`eqSolvemtControlCreate`, :func:`eqSolvemtOutCreate`